Universality of the Turbulent Magnetic Field in Hypermassive Neutron Stars Produced by Binary Mergers

The GW170817 event (Abbott et al. 2017b, 2017c) was arguably one of the most important astrophysical findings of the last decade for several reasons. First, it demonstrated that binary neutron star (BNS) mergers can produce strong gravitational wave (GW) signals and power bright electromagnetic emissions on a broad range of the spectrum (Goldstein et al. 2017; Savchenko et al. 2017; Abbott et al. 2017a, 2017d; Metzger 2017; D'avanzo et al. 2018; Fong et al. 2019; Dobie et al. 2018; Mooley et al. 2018). Second, the comparison of theoretical models with the observations allowed one to narrow some of the physical properties of neutron stars (NSs), such as their radius, maximum mass, tidal deformability, and equation of state (EoS; see, e.g., Margalit & Metzger 2017; Shibata et al. 2017; Abbott et al. 2018). In addition, it also served to measure the GW speed with unprecedented precision, setting strong constraints in many alternative theories of gravity (see, e.g., Baker et al. 2017; Gong et al. 2018).

Although the dynamics of BNS mergers is fairly well understood (e.g., Ciolfi 2020a), there are still many open questions regarding the details of the physical processes taking place during and after the merger. Here we focus on one of them, namely the amplification and large-scale reorganization of the magnetic field, which is thought to be necessary in order to launch a magnetically dominated jet associated with the short gamma-ray burst (SGRB; see, e.g., McKinney & Blandford 2009; Liska et al. 2020). Although the BNS merger and post-merger evolution of the remnant has been extensively studied through general-relativistic magnetohydrodynamics (GRMHD) simulations (Palenzuela et al. 2013b; Kiuchi et al. 2014; Neilsen et al. 2014; Kiuchi et al. 2015; Giacomazzo et al. 2015; Palenzuela et al. 2015; Ruiz et al. 2016; Kiuchi et al. 2018; Ciolfi et al. 2019; Ciolfi 2020b; Ruiz et al. 2020; Mösta et al. 2020), the problem has not yet been fully resolved. The impossibility of capturing the small (but dynamically important) scales induced by the MHD instabilities at play, possibly combined with other instabilities, obscures a definite answer as to the topology and intensity of the magnetic field in the remnant, when the turbulent amplification approaches saturation. Moreover, these strong magnetic fields are thought to enhance the angular momentum transport (see, e.g., Ciolfi 2020a and references therein), a key factor in the fate of the remnant. In the presence of large-scale magnetic fields, the formation of a jet can be favored, and the amount mass ejecta could also change, compared to a nonmagnetized remnant. Moreover, the presence of a jet or magnetically driven winds also depends on the field topology.

Different GRMHD simulations have attempted to resolve all the relevant scales of the problem. The highest-resolution simulations so far (Kiuchi et al. 2018) showed that the average magnetic field of the remnant can be amplified from 10¹³ to 10¹⁶ G during the first milliseconds after the merger. Unfortunately, even the very fine spatial grid size used there (i.e., Δ ∼ 12.5 m) is still much larger than the estimated wavelength of the fastest-growing unstable modes of the Kelvin–Helmholtz instability (KHI): insufficient to well capture the small, highly dynamical scales and the magnetic field amplification at this stage. As a result, the saturation level of the magnetic field intensity did not converge to any clear value.

At present (and arguably for the foreseeable future) it is not possible to perform direct numerical simulations of this scenario, since the range of the relevant scales (from hundreds/thousands of kilometers of the domain to the submeter shear layer thickness) makes the computational cost unfeasible, even employing the most efficient numerical and parallelization methods currently available. In order to overcome this limitation, different strategies have been implemented to reproduce the under-resolved amplified small-scale magnetic field. A commonly used strategy is to impose a purely poloidal magnetic field, with unrealistically large strengths ∼10¹⁴⁻¹⁶ G, either before or after the merger (e.g., Ruiz et al. 2016; Kiuchi et al. 2018; Ciolfi et al. 2019; Ciolfi 2020b; Mösta et al. 2020; Ruiz et al. 2020). This choice is hardly comparable to the real effects of the dynamo operating at small scales, for which a purely large-scale ordered field is not an expected outcome.

Other alternatives involve the use of large-eddy simulations (LESs), which consist of including extra terms in the discretized version of the evolution equations to account for the unresolved subgrid scale (SGS) dynamics (see, e.g., Yang 2015). The main idea of a LES is to reproduce the imprints of the SGS dynamics on the large-scale (numerically resolved) fields, thus providing a magnetic field growth with a realistic topology and spectrum. Following this line, we have recently extended and implemented the so-called gradient SGS model used in nonrelativistic fluid dynamics (Leonard 1975; Müller & Carati 2002) to the nonrelativistic (Viganò et al. 2019), special (Carrasco et al. 2020), and general-relativistic (Viganò et al. 2020) MHD with excellent results in capturing the small-scale effects of turbulent flow. We have performed LES of BNS coalescence (Aguilera-Miret et al. 2020), and found that the average magnetic field in the remnant is amplified with much less computational resources than the higher-resolution simulations leading to comparable results. In an accompanying paper (Palenzuela et al. 2021), we show that high-resolution LESs provide an amplification of the average magnetic field from 10¹¹ G to 10¹⁶ G. More importantly, for the first time the magnetic field strength and its energy spectral distribution converge to the same saturated level.

The work shown in this Letter, which uses the same techniques as in Palenzuela et al. (2021), focuses on the role of the initial topology and intensity of the magnetic field in the post-merger remnant. Most (if not all) simulations of magnetized BNS mergers up to date start with mainly dipolar magnetic fields, often confined in each star, for simplicity. However, magnetic topologies are expected to be much richer, with relevant contributions from small scales and from the magnetospheric currents. This holds throughout a NS's life: at birth after the core-collapse supernova (Reboul-Salze et al. 2021), for middle (Myr) ages, as shown by magnetar observations (Tiengo et al. 2013; Borghese et al. 2015) and simulations (Gourgouliatos et al. 2016), and for late (Gyr) ages similar to what merging NSs should have (as shown by NICER studies of old millisecond pulsars; Riley et al. 2019). The key question we want to address here is the following: how does the choice of the pre-merger configuration affect the final magnetic field of the remnant? The results presented in this Letter indicate that the memory of the initial magnetic field configuration is lost during the amplification phase induced by the small-scale dynamo. For all the initial topologies considered, the bulk of the remnant (i.e., the regions with ρ ≥ 10¹³ g cm⁻³) is endowed with a very similar isotropic, turbulence-like configuration with an average magnetic field of approximately 10¹⁶ G.

The paper is organized as follows. The evolution equations and the numerical setup are described briefly in Section 2. The results of the simulations are presented and analyzed in Section 3. Conclusions are drawn in Section 4.

Both the Einstein and the full set of filtered GRMHD evolution equations, including all the gradient SGS terms, can be found in Viganò et al. (2020) and Aguilera-Miret et al. (2020). Following those works, we include only the SGS term (i.e., ${\tau }_{M}^{{ki}}$) appearing in the induction equation, which in a flat spacetime can be written as

$$\begin{eqnarray*}&&\begin{array}{l}{\partial }_{t}{B}^{i}+{\partial }_{k}[{B}^{i}{v}^{k}-{B}^{k}{v}^{i}-{\tau }_{M}^{{ki}}]=0,{\tau }_{M}^{{ki}}=-\displaystyle \frac{{\rm{\Delta }}{x}^{2}}{24}{{ \mathcal C }}_{{ \mathcal M }}\,{H}_{M}^{{ki}},\end{array}\end{eqnarray*}$$

where the explicit expressions for the tensor ${{H}}_{{M}}^{{ki}}$ in terms of field gradients can be found in Aguilera-Miret et al. (2020). As we found in that study, the value ${{ \mathcal C }}_{{ \mathcal M }}=8$ reproduces the magnetic field amplification more accurately for our numerical schemes, and we will set the same value for all our simulations in this work. Note that less-dissipative numerical schemes, or higher-resolution simulations, would require smaller values of ${{ \mathcal C }}_{{ \mathcal M }}$ closer to the theoretical expectation ${{ \mathcal C }}_{{ \mathcal M }}\simeq 1$.

As in our previous works, we use the MHDuet code, generated by the Simflowny platform (Arbona et al. 2013, 2018; Palenzuela et al. 2021) and based on the SAMRAI infrastructure (Hornung & Kohn 2002; Gunney & Anderson 2016). MHDuet employs at least fourth-order-accurate operators to discretize both the spatial and time derivatives, with a high-resolution shock-capturing method for the fluid based on the Lax–Friedrichs flux-splitting formula (Shu 1998) and the fifth-order reconstruction method (MP5; Suresh & Huynh 1997). A complete assessment of the implemented numerical methods can be found in Palenzuela et al. (2018). We employ the same hybrid EoS as in Palenzuela et al. (2021) for the evolution, with a cold contribution given by a tabulated polytrope fit to the APR4 zero-temperature EoS (Read et al. 2009), and thermal effects modeled by the ideal gas EoS with adiabatic index Γ_th = 1.8.

The conversion from the conserved fields to the primitive ones is performed by using the robust procedure given in Kastaun et al. (2020). To minimize further failures in the very tenuous regions outside the star, we impose a minimum density of 6.2 × 10⁵ g cm⁻³, with the regions having such values referred to hereafter as atmosphere. Moreover, we apply the SGS terms only in regions where the density is higher than 6.2 × 10¹³ g cm⁻³ in order to avoid spurious effects near the stellar surface. Since the remnant's maximum density is above 10¹⁵ g cm⁻³, the SGS model is applied only in the most dense regions of the star.

The initial data is created with the Lorene package (2010), using the same tabulated polytropic EoS described above. We consider an equal-mass BNS in a quasi-circular orbit with an irrotational configuration. The total mass of the system is M = 2.7 M_⊙ and the initial separation is 45 km, corresponding to an initial angular frequency of 1775 rad s⁻¹.

The binary is solved in a cubic domain of side $\left[-1204,1204\right]$ km. The inspiral is fully covered by seven fixed-mesh refinement levels, each being a cube doubling the resolution of the previous one, and one adaptive-mesh refinement level, achieving a maximum resolution of 60 m in a domain covering at least the bulk of the remnant.

For each star, we consider a commonly used initial axially symmetric magnetic field, confined in the region where the fluid pressure P is larger than a value P_cut, set to 100 times the atmospheric pressure. The azimuthal (ϕ) component of the vector potential has a radial dependence A_ϕ ∝ r²(P − P_cut), where r is the distance from the center of the star. We have considered four initial configurations that differ in the intensity and the colatitude (θ) dependence of the magnetic field, as follows (see also Table 1):

• 1.  
  Aligned dipole-like (Dip): a very ordered (large-scale) poloidal field ${A}_{\phi }={A}_{0}{r}^{2}{\sin }^{2}\theta (P-{P}_{\mathrm{cut}})$, similar to a dipole (which would go $\propto \sin \theta$), with the magnetic moment aligned to the orbital axis, and a normalization value A₀ such that the maximum intensity (at the center) is 10¹² G, orders of magnitude lower than the large initial fields of other simulations (e.g., Kiuchi et al. 2015; Ruiz et al. 2016; Kiuchi et al. 2018; Ciolfi et al. 2019; Ciolfi 2020b; Ruiz et al. 2020).
• 2.  
  Highly magnetized (BHigh): the same as the above Dip model, except that the intensity is 1000 times larger, reaching a maximum of 10¹⁵ G.
• 3.  
  Misaligned dipole (Misal): the same as the Dip model, except that the magnetic moment is orthogonal to the orbital axis.
• 4.  
  Multipole (Mult): a more complex topology containing high multipolar structures, with ${A}_{\phi }\propto {r}^{2}{\sin }^{6}\theta \left(1+\cos \theta \right)(P-{P}_{\mathrm{cut}})$, with the same maximum intensity as Dip.

Our initial binary system evolves for five orbits before merging. The merger produces a differentially rotating remnant that relaxes to a hypermassive NS in a few milliseconds. Before the merger occurs (hereafter, t = 0), we set a magnetic field topology on each star corresponding to one of the cases described before. Therefore, we have then considered four different simulations, as summarized in Table 1.

Figure 1 displays some snapshots, in the orbital plane z = 0, of the Dip, BHigh, Misal, and Mult simulations (from top to bottom) at t = {2, 5, 10, 20} ms (from left to right) after the merger. The orange scale represents the intensity of the magnetic field, while the two thin black lines are mass density contours corresponding to ρ = 10¹⁴ g cm⁻³ (inner line) and 10¹³ g cm⁻³ (outer line). The shape of the remnant varies at the initial times, but clearly restructures itself at later ones, approaching an almost axisymmetric structure. There we can also see the KHI appearing at the merging layers, amplifying local values of the magnetic field (MF) up to a maximum of 10¹⁷ G. Thus, the MF changes from fully turbulent to partially ordered at ∼20 ms after the merger, where we can see the systematic formation of azimuthal/spiraling filaments. This is due to the winding that from now on rule the rising of the magnetic field (see Palenzuela et al. 2021 for an in-depth discussion about the mechanisms contributing to the amplification).

Zoom In Zoom Out Reset image size

In Figure 2 we represent the evolution of the volume-integrated thermal energy (top, solid line), kinetic rotational (top, dashed line), and magnetic energy (bottom) for the four models. We can see the energies of all cases rise soon after merger, at the expense of the large gravitational energy available. For all cases, the thermal energy still rises monotonically after the merger while the rotational one is losing energy, all becoming remnants objects with higher temperature but rotating more slowly. We obtain comparable values between all cases at 30 ms after the merger for the rotational energy. For the thermal one, differences are around a factor of ∼1.2 between BHigh and Mult at the same time.

Zoom In Zoom Out Reset image size

After about 5 ms, the magnetic energy growth saturates, with a maximum factor difference of ∼3 between Dip and Misal. The KHI that appears during the merger between the two stars, possibly combined with the Rayleigh–Taylor instability near the surface, is responsible for the amplification of the magnetic energy, which for all cases (except by the BHigh) increases 10 orders of magnitude, from 10⁴⁰ to 10⁵⁰ erg. For the BHigh case, the initially large value implies a smaller amplification by four orders of magnitude from 10⁴⁶ erg to a similar value of 10⁵⁰ erg. Overall, differences in energies lie within a factor ∼3, much less than both the orders-of-magnitude differences in the initial magnetic energy and across the evolution.

Figure 3 shows the intensity of the poloidal (solid lines) and toroidal (dashed lines) components of the magnetic field, as a function of time, averaged in the bulk of the remnant. There we can see the amplification for both components in the first 5 ms after the merger. What is important to remark upon here is that (i) during the exponential growth, both components are similar due to the isotropic character of turbulence, and (ii) after the rising, the toroidal component is the dominant one, as its shape is practically the same as the magnetic energy plot from Figure 2. The toroidal component, in the saturation phase, maintains its value as merely constant for all cases, close to 10¹⁶ G. The poloidal component, on the other hand, is slowly decreasing, probably due to energy transfer to the toroidal component, showing values around 10¹⁵ G. For the toroidal component of the magnetic field we can see almost the same behavior for all cases where, again, the misaligned and the multipolar ones are below the others by a factor ∼3. The differences are less drastic when focusing on the poloidal component of the magnetic field, where at the beginning the same cases (i.e., misaligned and multipolar) do not rise as much as the others but at 30 ms after the merger the differences are only about a factor ∼2.

Zoom In Zoom Out Reset image size

Besides the volume-integrated quantities, we can analyze the evolving spectral energy distribution ${ \mathcal E }(k)$ (for details on how it is calculated, see the Appendix of Viganò et al. 2020). This will allow us to see whether the different scales of the problem behave similarly or not between the cases we here consider. The spectral distribution of the kinetic (solid lines) and magnetic (dotted lines) energies, as a function of the wavenumber k, is displayed in Figure 4. The four plots correspond to, from left to right, t = {5, 10, 20, 30} ms after the merger. As a reference, Kolmogorov (k^−5/3, thin solid line) and Kazantsev (k^3/2, thin dotted line) slopes are also included in all plots. Large dots indicate the spectra-weighted average of the wavenumber, $\bar{k}=\tfrac{{\int }_{k}[k{ \mathcal E }(k)]}{{\int }_{k}{ \mathcal E }(k)}$, which gives the typical size $\bar{\lambda }\,=2\pi /\bar{k}$ of either the fluid or the magnetic (${\bar{\lambda }}_{B}$) structures.

Zoom In Zoom Out Reset image size

For all times represented, the kinetic energy spectra behave in the same way for all cases (having a Kolmogorov slope in the inertial range), regardless of the scale considered. For the magnetic energy spectra, they initially follow the Kazantsev slope up to the numerical dissipative scale (intrinsically set by the discretization scheme). This is a property of the kinematic phase of the dynamo, until the dynamo approaches saturation at small scales (large k). At t = 10 ms all cases have roughly the same amount of magnetic energy spectra. At t = 20 ms after the merger, small differences begin to appear, although they may be in part due to stochastic variations. At t = 30 ms after the merger, such differences are less evident. Moreover, the amplification has saturated and magnetic spectra appear to be compatible with a Kolmogorov slope at intermediate scales.

We find that ${\bar{\lambda }}_{B}\sim 800$ m soon after the merger, increasing to almost 2 km at t = 30 ms. This confirms that larger coherent magnetic field structures are being formed in the remnant. Clearly, there is no significant difference (less than 7% in ${\bar{\lambda }}_{B}$) between the simulations with different initial topologies considered here.

In Figure 5 we further analyze the magnetic energy spectra, identifying the contributions coming from the toroidal (dashed lines) and poloidal (solid) components. At 5 ms after the merger, both components are similar for all simulations. As time passes, the poloidal component of all cases decreases two orders of magnitude, while the toroidal one increases by about one order of magnitude. However, for all times, both components are still comparable among the different simulations. We note that, comparing both components in different models, the differences are up to a factor 3 only, much less than the relative differences and the overall changes in time. Finally, at t = 30 ms, the slope of the toroidal component (the dominant contribution to the magnetic energy) approaches a Kolmogorov slope at intermediate scales. An interesting difference at these times is that, starting with an (unrealistically) high magnetic field, there is an excess of large-scale magnetic fields (low k), an effect probably due to the winding acting on the already quite organized field.

Zoom In Zoom Out Reset image size

Finally, we note that the resemblance of integrated magnetic energy and the spectra for the different models point to a universality of the magnetic field not only in the bulk, but in the entire domain.

We have studied the influence of different initial magnetic configurations on the evolution of BNS mergers, using high-resolution LESs, employed also in an accompanying paper (Palenzuela et al. 2021), to which we refer for further details on the methods and amplification mechanisms at work. In particular, we have considered initial magnetic fields confined within the stars, varying the intensity, the magnetic moment misalignment with the orbital axis, and the poloidal topology. Looking at the evolution of integrated energy and spectral distribution, we proved that the differences lie within a factor 3 at most, which could be even smaller in more accurate simulations. This, then, ensures that the initial topology of the stars is not relevant at all because the small-scale turbulence induced in the remnant will erase any memory of realistic magnetic fields of B ≤ 10¹² G within only a few milliseconds after the merger.

In this work we have explored only some of the infinite possible magnetic configurations. More choices could be explored, in particular the presence of a toroidal field, nonaxisymmetric topology, a magnetic field extended outside the star, and more complex, small-scale-dominated configurations. However, the results shown here already suggest that the expected dependence on the initial topology is basically negligible, compared to other much more uncertain issues. Among the latter, we mention the importance of the numerical capability to resolve the small-scale magnetic amplification and the physics involved in the post-merger phase (neutrino transport, temperature-dependent equation of state, etc.).

This universality of the magnetic field outcome after the merger sets serious doubts as to how could we infer the initial magnetic field of the stars in a BNS merger through multimessenger astronomy. The only foreseeable possibilities are through the presence of a precursor electromagnetic signal before the merger, or other kinds of outflows that may appear during the first milliseconds after the merger, which would have information mainly on the initial topology and intensity of the magnetic field (Palenzuela et al. 2013b, 2013a; Ponce et al. 2014).

The final message is that the commonly used simplification on the topology of the magnetic field is acceptable in BNS mergers, as long as the magnetic field is not too large and enough turbulence is developed to erase the seed and produce the correct spectra distribution. In those cases, the system will tend to be practically the same remnant regardless of its initial configuration. However, if one wants to focus on the realistic generation of a large-scale field in the post-merger, the rule-of-thumb is basically the following: be sure that the large-scale initial magnetic field is much smaller than the amplified average field that the numerical scheme is able to reproduce.

We thank Riccardo Ciolfi, Wolfgang Kastaun, and Jay Vijay Kalinani for providing us the EoS and for the useful discussions. This work was supported by European Union FEDER funds, the Ministry of Science, Innovation and Universities, and the Spanish Agencia Estatal de Investigación grant No. PID2019-110301GB-I00. D.V. is funded by the European Research Council (ERC) Starting grant IMAGINE (grant agreement No. [948582]) under the European Union's Horizon 2020 research and innovation program. D.V.'s work was also partially supported by the program Unidad de Excelencia María de Maeztu CEX2020-001058-M. The authors thankfully acknowledge the computer resources at MareNostrum and the technical support provided by the Barcelona Supercomputing Center (BSC) through grant project Long-LESBNS by the 22nd PRACE regular call (Proposal 2019215177, PIs: C.P. and D.V.).