Universality of the turbulent magnetic field in hypermassive neutron stars produced by binary mergers

The detection of a binary neutron star merger in 2017 through both gravitational waves and electromagnetic emission opened a new era of multimessenger astronomy. The understanding of the magnetic field amplification triggered by the Kelvin-Helmholtz instability during the merger is still a numerically unresolved problem because of the relevant small scales involved. One of the uncertainties comes from the simplifications usually assumed in the initial magnetic topology of merging neutron stars. We perform high-resolution, convergent large-eddy simulations of binary neutron star mergers, following the newly formed remnant for up to $30$ milliseconds. Here we specifically focus on the comparison between simulations with different initial magnetic configurations, going beyond the widespread-used aligned dipole confined within each star. The results obtained show that the initial topology is quickly forgotten, in a timescale of few miliseconds after the merger. Moreover, at the end of the simulations, the average intensity ($B\sim 10^{16}$ G) and the spectral distribution of magnetic energy over spatial scales barely depend on the initial configuration. This is expected due to the small-scale efficient dynamo involved, and thus it holds as long as: (i) the initial large-scale magnetic field is not unrealistically high (as often imposed in mergers studies); (ii) the turbulent instability is numerically (at least partially) resolved, so that the amplified magnetic energy is distributed across a wide range of scales and becomes orders-of-magnitude larger than the initial one.


I. INTRODUCTION
The GW170817 event [1,2] 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 [3][4][5][6][7][8][9][10][11]. Second, the comparison of theoretical models with the observations allowed 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., [12][13][14]). 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., [15,16]).
Although the dynamics of BNS mergers is fairly well understood (e.g., [17]), there are still many open questions regarding the details of the physical processes taking place during and after the merger. Here we focus in 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 to the short gamma-ray burst (SGRB) (see e.g. [18,19]). Although the BNS merger and post-merger evolution of the remnant has been extensively studied through general-relativistic magnetohydrodynamics (GRMHD) simulations [20][21][22][23][24][25][26][27][28][29][30][31], the problem is not fully resolved yet. 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 on what is the topol-ogy 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. [17] and references therein), a key factor in the fate of the remnant. In presence of large-scale magnetic fields, the formation of the jet can be favored, and the amount mass ejecta could also change, compared to a non-magnetized 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 highestresolution simulations so far [27] showed that the average magnetic field of the remnant can be amplified from 10 13 G to 10 16 G during the first miliseconds after the merger. Unfortunately, even the very fine spatial grid size used there (i.e. ∆ ∼ 12.5m) 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 was not converging to any clear value.
Nowadays (and arguably in 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 km of the domain to the sub-meter 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 one is to impose a purely poloidal magnetic field, with unrealistically large strengths ∼ 10 14−16 G, either before or after the merger(e.g., [26][27][28][29][30][31]). This choice is hardly comparable to the real effects of the dynamo operating at small scales, for which purely large-scale ordered field is not an expected outcome.
Other alternatives involve the use of large-eddy simulations (LES), which consist on including extra terms in the discretized version of the evolution equations to account for the unresolved sub-grid scale (SGS) dynamics (see e.g. [32]). The main idea of the LES is to reproduce the imprints of the sub-grid 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 non-relativistic fluid dynamics [33,34] to the non-relativistic [35], special [36] and general relativistic [37] MHD with excellent results on capturing the small scale effects of turbulent flow. We have performed LES of BNS coalescence [38], 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 [39], we show that high-resolution LES provide an amplification of the average magnetic field from 10 11 G to 10 16 G. More importantly, for the first time the magnetic field strength and its energy spectral distribution are converging to the same saturated level.
The work shown in this letter, which uses the same techniques as in [39], focuses on the role of the initial topology and intensity of the magnetic field in the postmerger 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 neutron star's life: at birth after the core-collapse supernova [40], for middle (Myr) ages (as shown by magnetars' observations [41,42] and simulations [43]) and for late (Gyr) ages similar to what merging NSs should have (as shown by NICER studies of old millisecond pulsars [44]). The key question we want to address here is the following: how 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 13 g cm −3 ) is endowed with a very similar isotropic, turbulent-like configuration with an average magnetic field of approximately 10 16 G.
The paper is organized as follows: the evolution equations and the numerical setup are described briefly in §II. The results of the simulations are presented and analysed in §III. Conclusions are drawn in §IV.

II. INITIAL MODELS
Both the Einstein and the full set of filtered GRMHD evolution equations, including all the gradient SGS terms, can be found in [37,38]. Following those works, we include only the SGS term (i.e., τ ki M ) appearing on the induction equation, that in a flat spacetime can be written as where the explicit expressions for the tensor H ki M in terms of field gradients can be found in [38]. As we found in that study, the value C 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 C M closer to the theoretical expectation C M 1.
As in our previous works, we use the code MHDuet, generated by the platform Simflowny [45][46][47] and based on the SAMRAI infrastructure [48,49]. MHDuet employs at least fourth-order-accurate operators to discretize both the spatial and time derivatives, with highresolution shock-capturing method for the fluid based on the Lax-Friedrich flux splitting formula [50] and the fifthorder reconstruction method MP5 [51]. A complete assessment of the implemented numerical methods can be found in [52]. We employ the same hybrid EoS as in [39] for the evolution, with a cold contribution given by a tabulated polytrope fit to the APR4 zero-temperature EoS [53], 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 [54]. To minimize further failures in the very tenuous regions outside the star, we impose a minimum density of 6.2 × 10 5 g cm −3 , with the regions having such values referred hereafter as atmosphere. Moreover, we apply the SGS terms only in regions where the density is higher than 6.2 × 10 13 g cm −3 in order to avoid spurious effects near the stellar surface. Since the remnant's maximum density is above 10 15 g cm −3 , the SGS model is applied only in the most dense regions of the star.
The initial data is created with the Lorene package [55], using the same tabulated polytropic EoS described above. We consider an equal-mass BNS in quasicircular orbit with an irrotational configuration. The total mass of the system is M = 2.7 M and the initial separation is 45km, corresponding to an initial angular frequency of 1775 rad s −1 .
The binary is solved in a cubic domain of side [−1204, 1204] km. The inspiral is fully covered by 7 Fixed Mesh Refinement (FMR) levels, each being a cube doubling the resolution of the previous one, and one Adaptive Mesh Refinement (AMR) 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 2 (P − P cut ), where r is the distance from the center of the star. We have considered four initial configurations that differ among themselves in the intensity and the co-latitude (θ) dependence of the magnetic field, as follows (see also Table I): • Aligned Dipole-like (Dip): A very ordered (large scale) poloidal field A φ = A 0 r 2 sin 2 θ(P − P cut ), similar to a dipole (which would go ∝ sin θ), with the magnetic moment aligned to the orbital axis, and a normalization value A 0 such that the maximum intensity (at the centre) is 10 12 G, orders of magnitude lower than the large initial fields of other simulations (e.g., [23,[26][27][28][29][30]).
• Highly magnetized (BHigh): The same as the above Dip model, except that the intensity is 1000 times larger, reaching a maximum of 10 15 G.
• Misaligned dipole (Misal): The same as the Dip model, except that the magnetic moment is orthogonal to the orbital axis.
• Multipole (Mult): A more complex topology containing high multi-polar structures, with A φ ∝ r 2 sin 6 θ (1 + cos θ) (P − P cut ), with the same maximum intensity as Dip. We indicate the initial values of the maximum intensity of the magnetic field, the angular misalignment between the orbital and magnetic field axes, and the initial topology of the magnetic field. Cells remarked in grey correspond to the differences with respect to the reference model Dip.

III. RESULTS
Our initial binary system evolves for 5 orbits before merging. The merger produces a differentially rotating remnant that relaxes to an hypermassive-neutron star (HMNS) 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 I. Fig. 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 black thin lines are mass density contours corresponding to ρ = 10 14 g cm −3 (inner line) and 10 13 g cm −3 (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 MF up to a maximum of 10 17 G. Thus, 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 [39] for an in-depth discussion about the mechanisms contributing to the amplification). In Fig. 2 we represent the evolution of the volumeintegrated 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, becoming all 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.
After about 5 milliseconds, 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 the responsible of the amplification of the magnetic energy, which for all cases (except by the BHigh) increases 10 orders of magnitude, from 10 40 erg to 10 50 erg. For the BHigh case, the initially large value implies a smaller amplification by 4 orders of magnitude from 10 46 erg to a similar value of 10 50 erg. Overall, differences in energies lie within a factor ∼ 3, much less than both the orders-ofmagnitude differences in the initial magnetic energy and across the evolution. Fig. 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 due to the KHI in the firsts 5 milliseconds after the merger. What is important to remark 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 Fig. 2. The toroidal component, in the saturation phase, maintains its value merely constant for all cases, close to 10 16 G. The poloidal component, on the other hand, is slowly decreasing, probably due to energy transfer to the toroidal component, getting values around 10 15 G. For the toroidal component of the magnetic field we can see almost the same behaviour for all cases where, again, the misaligned and the multipolar ones are below the others by a factor ∼ 3. Differences are less notorious 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.
Besides the volume-integrated quantities, we can analyze the evolving spectral energy distribution E(k) (for details on how it is calculated, see appendix of [37]). This will allow us to see whether the different scales of the problem behave similar or not between the cases we are considering. The spectral distribution of the kinetic (solid lines) and magnetic (dotted lines) energies, as a function of the wavenumber k, is displayed in Fig. 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,k = k [kE(k)] k E(k) , which gives the typical sizeλ = 2π/k of either the fluid or the magnetic (λ B ) structures. 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 we are considering. 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 found thatλ B ∼ 800 m soon after the merger, increasing to almost 2 km at t = 30 ms. This confirm that larger coherent magnetic field structures are being formed in the remnant. Clearly, there is no significant difference (less than 7% inλ B ) between the simulations with different initial topologies considered here.
In Fig. 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. Notice 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 a (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.
Finally, notice 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.

IV. CONCLUSIONS
We have studied the influence of different initial magnetic configuration on the evolution of BNS mergers, using high-resolution large-eddy simulations employed also in the accompanying paper [39], for which we refer for further details on the methods and amplification mechanisms at work. In particular, we have considered ini- tial 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 12 G in only 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: presence of a toroidal field; nonaxisymmteric topology; magnetic field extended outside the stars; more complex, small-scale dominated configuration... 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 smallscale magnetic amplification, the physics involved in the post-merger phase (neutrino transport, temperaturedependent equation of state...).
This universality of the magnetic field outcome after the merger sets serious doubts on 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 kind of outflows that may appear during the first ms after the merger, which would have information mainly on the initial topology and intensity of the magnetic field [20,56,57].
The final message is that the commonly used simplification on the topology of the magnetic field is acceptable in BNS mergers, as far 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.