Turbulent magnetic-field amplification in the first 10 milliseconds after a binary neutron star merger: comparing high-resolution and large eddy simulations

The detection of binary neutron star mergers represents one of the most important and complex astrophysical discoveries of the recent years. One of the unclear aspects of the problem is the turbulent magnetic field amplification, initially triggered by the Kelvin-Helmholtz instability at much smaller scales than any reachable numerical resolution nowadays. Here we present numerical simulations of the first ten milliseconds of a binary neutron star merger. First, we confirm in detail how the simulated amplification depends on the numerical resolution and is distributed on a broad range of scales, as expected from turbulent MHD theory. We find that an initial large-scale magnetic field of $10^{11}\,$G inside each star is amplified in the remnant to root-mean-square values above $10^{16}\,$G within the first $5$ milliseconds for our highest-resolution run. Then, we run large eddy simulations, exploring the performance of the subgrid-scale gradient model, already tested successfully in previous turbulent box simulations. We show that the addition of this model is especially important in the induction equation, since it leads to an amplification of the magnetic field comparable to a higher-resolution run, but with a greatly reduced computational cost. In the first 10 milliseconds, there is no clear hint for an ordered, large-scale magnetic field, which should indeed occur in longer timescales through magnetic winding and the magneto-rotational instability.

Although the central aspects of BNS systems are qualitatively understood, the details of the merger and postmerger dynamics remain only poorly constrained, with many important questions still open. In this paper we are mainly concerned with one of such issues: the amplification and large-scale (re-)organization of the magnetic field, arguably required to launch the successful jet outflows associated to the short gamma-ray burst (SGRB). Despite the recent progress of general-relativistic magnetohydrodynamics (GRMHD) simulations [15][16][17][18][19][20][21][22][23][24][25][26], the impact of magnetic turbulence on the evolution of the hypermassive neutron star (HMNS) remnant is highly uncertain, mostly due to the lack of a spatial resolution able to capture all the relevant scales. It has been recognized that the effects of turbulent viscosity and dynamo (so far numerically under-resolved), along with neutrino transport, can be crucial for the redistribution of angular momentum, mass ejecta, lifetime of the remnant and production of the jet (e.g., [27]).
Observationally, the typical range of magnetic field strengths characterizing Gyr-old NSs (typical age at which binaries can merge) is 10 8−11 G [28]. 1 Magnetic field amplification occurs during and after merger through a number of distinct MHD mechanisms, channeling a fraction of the abundant orbital kinetic energy (∼ 10 53 erg) of the system. The Kelvin-Helmholtz instability (KHI), originated in the shearing layer at the collision interface, drastically enhances the magnetic field by stretching and folding embedded field lines in a process known as small-scale turbulent dynamo. Local specialrelativistic MHD simulations have shown that the development of the KHI at merger can generate magnetar-level magnetic field strengths within the first few milliseconds [31,32]. Later, GRMHD simulations of BNS mergers of unprecedented high-resolution (grid-spacing of 17.5 m) [18] showed that an initial magnetic field of moderate strength 10 13 G can be amplified up to ∼ 10 16 G within ∼ 5 ms after merger, reaching magnetic saturation levels at energies E B 10 50 erg. However, no sign of numerical convergence was found, meaning that the KHI is not yet fully resolved even at those resolutions.
After the quick growth of the magnetic field due to the KHI, there are two other mechanisms associated to the differentially rotating HMNS that dominates on longer timescales 10 ms: magnetic winding, which linearly amplifies the toroidal components of the field from the poloidal ones, and the magneto-rotational instability (MRI). For the latter, the wavelengths of the fastest growing modes are proportional to the magnetic fields. Therefore, even the highest-resolution GRMHD simulations to date cannot resolve the MRI, unless artificially large initial magnetic fields above 10 13 G are adopted as to increase the associated cutoff length scales. Even in this way, simulations are far from capturing the turbulent cascade all the way down to the viscous scale (determined by neutrino viscosity [33]), as it would be required for a direct numerical simulation (DNS). Finally, efficient MRI amplification is expected to continue acting inside the accretion disk after the remnant collapses to a black hole.
In the absence of computationally viable DNS to consistently evolve all the phases of the magnetic dynamics described above, different approaches were considered. Many studies have imposed rather large initial pre-merger (e.g., [21][22][23][24][25]) or post-merger ( [26]) magnetic field strengths ∼ 10 14−16 G, to compensate the inability to capture the KHI amplification. However, the quantitative results may not be fully reliable, since the amplification via KHI happens over a broad range of scales and does not preserve a large-scale ordered field.
One of the most promising alternatives is performing large eddy simulations (LESs), in which the evolution equations are modified in order to account for the unresolved subgrid-scale (SGS) dynamics [34]. This method was applied, in the present context, by including new terms (chosen proportional to the fluid vorticity) into the induction equation [19,20,35]. While the results of these studies show an effective growth of the magnetic field, they do not match the physical MHD dynamics and rely on arbitrarily tuning and switching "by hand" of the extra terms. Other approaches have, instead, centered their attention on the turbulent viscous effect during the post-merger phase, evolving viscous hydrodynamics (HD) in substitution of the MHD equations [36][37][38][39]. These models are however unable, by construction, to capture the dynamo mechanism and depend on parameters to be calibrated via very high resolution GRMHD simulations (e.g., [40]).
A more sophisticated alternative, based on the socalled gradient SGS model [41,42], has been proposed recently for Newtonian, special and general relativistic MHD, respectively in [43][44][45]. It was proven to have very good performance (in terms of capturing the magnetic amplification especially) in box simulations of the KHI, for a variety of initial conditions and resolutions, but it was not yet implemented in BNS mergers. The advantage of this approach is that it relies on the mathematical expansion of the fields involved in the dynamics, with no a-priori physical assumptions. In that sence, this SGS model is conceptually similar to high-order reconstruction methods used in finite-volume numerical schemes.
In this paper we perform BNS merger simulations, focusing on the magnetic field amplification due to the KHI during the first ∼ 10 ms after merger. In contrast to previous studies, our simulations begin with each star having realistic magnetic field strength values of about 10 11 G.
We use high-order numerical methods and the elaborated gradient SGS model already presented for GRMHD box simulations of the KHI [45], which is applied for the first time to the BNS merger scenario.
This article is organized as follows: our LES approach for GRMHD is briefly revisited on §II. The general setup, as well as the numerical methods, is described on §III. The results of the simulations are presented and analyzed in §IV. Conclusions are drawn on §V.

II. LARGE EDDY SIMULATIONS IN GRMHD
The concept and the mathematical foundations behind the explicit LES with a gradient SGS approach have been extensively explained in our previous works (and references within) in the context of Newtonian [43] and relativistic MHD [44,45], to which we refer for details and further previous references. In brief, the space discretization in any numerical simulation can be seen as a filtering of the continuous solution, with an implicit kernel (numerical-method-dependent) having the size of the numerical grid, ∆. The evolved numerical values of the fields can be then formally interpreted as weighted averages (or filtered) over the numerical cell. Seen in this way, the subgrid deviations of the field values from their averages causes a loss of information at small scales, for those terms which are nonlinear functions of the evolved variables. SGS terms obtained from the gradient model are added to the equations in order to partially compensate such loss.
Under the 3+1 decomposition framework [46], the line element can be written as where α is the lapse function, β i is the shift vector, and γ ij is the induced metric on each spatial foliation, with determinant √ γ. We use the covariant conformal Z4 formulation [47,48] to evolve the Einstein equations. A summary of the final set of evolution equation for the spacetime fields, together with the gauge conditions setting the choice of coordinates, can be found e.g. in [49]. The GRMHD equations for a magnetized, non-viscous and perfectly conducting fluid [20] (in units G = c = M = 1) consider the set of conserved variables √ γD, √ γS i , √ γU, √ γB i . They are functions of the rest-mass density ρ, the specific internal energy , the velocity vector v i and the magnetic field B i (primitive fields), as follows: where W = (1 − v 2 ) −1/2 is the Lorentz factor. The pressure p is defined through the EoS detailed in §III. The discretized evolution equations, including the hyperbolic divergence cleaning via damping of the field φ [49], can be written as follows: 2 The fluxes consist of the following standard terms: , and of the additional SGS terms: The cumbersome expressions of the tensors H have been obtained in detail for the special [44] and general relativistic [45] cases. Here we apply the latter, following the expressions reported in the Appendix A. The coefficient ξ = γ 1/3 ∆ 2 /24 has the proportionality to the spatial grid squared, which is typical of SGS models and ensures by construction the convergence to the continuous limit (vanishing SGS terms for an infinite resolution). Importantly, for each equation there is a precoefficient C i , which is meant to be of order one for a numerical scheme having a mathematically ideal Gaussian filter kernel and neglecting higher-order corrections. However, finite-difference numerical methods are usually more dissipative (and dispersive). Therefore, as shown in [43][44][45], the value that best mimics the feedback of small scales onto the large scales in a LES can differ depending partially on the numerical methods employed and on the specific problem. In practice, one needs a calibration of the different SGS parameters to maximize the effectiveness of the gradient model. Finally, the set of source terms in (5), written already as a function of conformal variables, can be found explicitly in [45]. SGS terms are applied to the fluid equation only, considering the full general relativistic setting, with the assumption that the metric components are smooth and slowly varying, as compared to the turbulent and shocked matter fields (see [45] for a discussion).

A. Numerical methods
As in our previous works, we use the code MHDuet, generated by the platform Simflowny [50,51] and based on the SAMRAI infrastructure [52,53], which provides the parallelization and the mesh refinement.
The code has been deeply tested for different scenarios [45,49,54,55], including basic tests of MHD and GR. Briefly, it uses: fourth-order-accurate operators for the spatial derivatives in the SGS terms and in the Einstein equations (the latter are supplemented with sixth-order Kreiss-Oliger dissipation); a high-resolution shock-capturing method for the fluid, based on the Lax-Friedrich flux splitting formula [56] and the fifth-order reconstruction method MP5 [57]; a fourth-order Runge-Kutta scheme with a small enough time step ∆t ≤ 0.4. ∆; and an efficient and accurate treatment of the refinement boundaries when sub-cycling in time [58,59]. A complete assessment of the implemented numerical methods can be found in [49,54].

B. EoS and conversion to primitive variables
We consider a hybrid EoS during the evolution, with two contribution to the pressure. On one side, we use the piecewise polytrope fit to the SLy zero-temperature EoS [60], defined by p = K i ρ Γi , where i = 0, 1, 2, 3 indicates each of the four segments delineated by the transition density values log ρ = {14.165, 14.7, 15.0}, Γ i = {1.35692, 3.005, 2.988, 2.851} and K 0 = 3.59389 × 10 13 (all in cgs units). On the other hand, thermal effects are modeled by an additional pressure contribution given by the ideal gas EoS, with adiabatic index Γ th = 1.75 [61].
The conversion from the evolved or conserved fields to the primitive or physical ones is performed by using the procedure described in our previous works [45,55]. An exception is the highest resolution simulation, for which the strong magnetic fields developed in low-density regions forced us to use a more robust procedure [62]. To minimize further failures on the recovery procedure outside the dense regions, we impose a minimum density of 6.1 × 10 7 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.1 × 10 11 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 accounted for only in the most dense regions of the star.

C. Initial conditions
The initial data is created with the Lorene package [63], using the same piecewise polytropic EoS described above. We consider an equal-mass BNS in quasicircular orbit, with an irrotational configuration having a separation of 37.7 km and an angular frequency of 2254 rad s −1 . The total mass of the system is M = 2.67 M .
Each star initially has a purely poloidal magnetic field confined in its interior, calculated from a vector potential A φ ∝ r 2 (P − P cut ), where P cut is a hundred times the pressure of the atmosphere and r the distance to the axis perpendicular to the orbital plane passing through the centre of each star. The maximum magnetic intensity (at the centres) is 5×10 11 G, orders of magnitude lower than the large initial fields of other simulations (e.g., [18,[21][22][23][24][25]) and compatible with the upper range of the expected realistic intensities for old NSs. Such values are also at the lower border of the computational feasibility, since the accurate evolution for too small ratios of magneticto-kinetic pressure is hampered by round-off errors.

IV. RESULTS
We consider a numerical cubic domain, ranging from [−384, 384] km along each direction, large enough to reduce contamination from the boundaries. Our initial binary system evolves for 2-3 orbits before merging and forming a differentially rotating remnant that relaxes to an hypermassive-neutron star (HMNS) in a few milliseconds. We follow such inspiral with five nested levels of Fixed Mesh Refinement (FMR), each being a cube doubling the resolution of the previous one. The smallest and finest of them is 70 km wide, thus it encloses the stars during the inspiral and the forming remnant. At the merger time (hereafter, t = 0), we have then considered different simulations, summarized in Table I.
First, we present standard simulation without SGS terms (C i = 0), also called implicit LES 3 (iLES, denoted by C0 hereafter), with grid spacing corresponding to low (LR, finest level: 147 m), medium (MR, 74 m) and high resolution (HR, 37 m). The LR case has five FMR levels, like in the inspiral. In the MR and HR cases, we activate one and two additional Adaptive Mesh Refinement (AMR) levels (again doubling the resolution of the previous level), respectively, describing the regions exceeding certain density thresholds properly set, in order to better resolve the remnant.
Secondly, we perform LES with LR including the SGS models. Here we report the cases with a fixed C M = 8, spanning C T = C N = {0, 1, 2, 4, 8}. Other combinations of parameters with C M = C T = C N > 2 have been tested, but they produced an excessive dissipation in the momentum equation, leading to unrealistic results. Also, we have considered the non-relativistic limit of the SGS term in the induction equation proposed in [45] (i.e. neglecting the H k v contribution in Eq.(A3)). In contrast to the box simulations results in [45], we see that the relativistic corrections on these SGS terms produce here a significant increase of the magnetic field amplification, so we have kept the full expressions for our simulations. From now on, we will refer to the LES simulations by labeling them in a schematic way according to the C i values, as indicated in Table I.

A. Results with different resolutions
First, we consider the three iLES cases. In Fig. 1 we show the density and the magnitude of the magnetic field, in the equatorial plane z = 0, for the three resolutions (different rows) at t = {2.5, 5, 10} ms (different columns) after the merger. In agreement with previous results [22], the magnetic field grows on small structures especially in the outermost, less dense layers of the remnant, where plasma is closer to an equipartition between magnetic and kinetic energy. As expected, this amplification is en-hanced by a finer grid, since smaller wavelengths grow faster in the KHI. At t = 2.5 ms, the two cores are still clearly distinguishable, indicating that the remnant has not relaxed to a HMNS yet. At this early stage, magnetic fields locally exceed ∼ 10 17 G only in the HR simulation, with fine structures clearly visible. At t = 5 ms, the remnant and surrounding disk are forming and turbulence drives the magnetic field amplification to maximum values of ∼ 5 × 10 17 G in the HR, dropping one order of magnitude in the MR simulation and another one for the LR case (see a quantitative comparison of magnetic energy evolution below). It can be seen how the magnetic field is mostly confined to the outermost layers of the remnant, since the dense core is less prone to turbulent motions. At t = 10 ms, the strong magnetic fields has started to penetrate into the dense core of the remnant in the HR run, while a significant overall increase in the field strength is also noticeable for the lower resolutions. At this time, although small-scale structures still dominate, the rotation has acquired a visible imprint on the magnetic field distribution, developing spiral-like filaments at the outermost layers.
This visual inspection can be quantified by the study of the energy spectral distribution (see Appendix B for definitions and calculations), as shown in Fig. 2 for t = {5, 10} ms. Note that, in general, we can identify the inertial range between scales much larger than ∆ (around which numerical dissipation acts) and smaller than the energy-injection scales (set in this case by the rotation). In such range, the kinetic and magnetic spectra approximately follow the Kolmogorov (k −5/3 ) and Kazantsev (k 3/2 ) slopes (dotted black lines in the figures hereafter), as expected in turbulent MHD scenarios.
The kinetic energy distribution is dominated by large scales, so that the resolution has a lower impact on it. Instead, the absence of a peak in the magnetic energy at low k (at least until 10 ms) means that there is no hint, at these times, for the creation of a strong, largescale, ordered magnetic field. Small scales are the main form of magnetic energy storage, hence the importance of the numerical resolution. This can be clearly observed, especially at early times (5 ms, left panel): the higher the resolution, the larger the growth of the magnetic en-ergy is, even though the spectra have the same profiles. At later times (10 ms, right panel), the difference between different resolutions greatly decreases, especially at large scales. Thus, pointing to a saturation of the KHI, achieved by all the three resolutions.
The magnetic amplification here illustrated presents the typical dynamical stages of the KHI, described e.g. in [32] as follows: an initial startup transient associated to the full development of the turbulent cascade (triggered at the merger); the kinematic phase, in which the magnetic fields are still sub-dominant but grow exponentially, driven by an essentially hydrodynamical turbulent mechanism as in Kazantzev's theory [64]; the approach to saturation when the magnetic field becomes strong enough as to back-react on the fluid motion and establish a dynamical balance signaled by kinetic/magnetic spectral equipartition at small scales. Generally speaking, the magnetic saturation levels are expected to converge (at least above certain threshold resolution [32]), while the growth rates and timescales of each dynamical phase are highly sensitive to the numerical resolution.
Note that the drop in the spectra for high k (approaching the upper limit set by π/∆) is due to the intrinsic numerical dissipation of the finite-volume scheme (spectral methods should not show it). Overall, the same behaviour was observed in our box simulations of the KHI [43,45]: a rising of the magnetic energy as smaller and smaller eddies develop, finally reaching equipartition of the spectral distribution at small scales (high k), while at large scales (low k) the kinetic energy always dominates.
It is also interesting to look separately at the magnetic components. In Fig. 3 of a developing turbulence that stretches and twists the initially weak large-scale magnetic seed. At later times, when differential rotation starts to dominate the kinematics, the magnetic structures tend to follow the rotation, partially losing the isotropy. In particular, at t = 10 ms, the toroidal field is slightly predominant over the other components, although still highly turbulent.
More quantitatively, Fig. 4 shows the toroidal magnetic spectra, defined as the azimuthal component of the field (in simple words, the direction identified by the remnant's bulk rotation), and the poloidal one, defined by the remaining directions (the naming of such decomposition is strictly correct only in axial symmetry, but we adopt it here for simplicity). At t = 5 ms these two components are very similar for all resolutions. However, at later times, the toroidal field starts to grow more in the HR case, in agreement with the above-mentioned intuition from Fig. 3. The isotropy of the KHI-induced turbulence at early stages and the hints for a gradual ordering in the azimuthal direction at t = 10 ms is consistent with previous results by [22], who showed how magnetic winding start to play an important role at this stage.  Table I. Let us first focus on the iLESs (solid lines). The HR case grows much faster, since smaller scales are excited by the KHI. This qualitatively agrees with the exponential growth rate ∝ 1/∆ predicted by the KHI analytical theory [65] and seen in previous GRMHD results [18]. At t = 5 ms after merger, only the C0 HR run has reached magnetic saturation, approximately at 2 × 10 50 erg. At this time, the magnetic energy of the three iLESs are separated by more than an order of magnitude among them. Later on, at t = 10 ms, the difference on the magnetic energy between the three resolutions is reduced almost by half, suggesting similar saturation levels of the magnetic field at late times.
The bottom panel displays the root mean square (r.m.s.) magnetic field for the three resolution iLES and for the LES simulation with the optimal parameters, CM8, that will be discussed in the next subsection. The r.m.s. is computed on regions with ρ > ρ X g cm −3 , being ρ A = 6×10 9 , ρ B = 6×10 10 and ρ C = 6×10 11 , and taking ρ B for the magnetic energy in the top panel. The run C0 HR shows mean values of 10 15 G when considering only the most dense part of the star (i.e., ρ > ρ C ), but increases to 10 16 G when also the outer envelope is taken into account (i.e., ρ > ρ A ).

B. LESs with gradient SGS model
Let us now turn to the effects of including the SGS model (LESs), and continue our analysis of Fig. 5. We shall stress that the aim of the SGS model, applied on the LR setup (at this particular stage of the merger evolution), is to reproduce the magnetic field amplification observed on the higher resolution simulations C0 MR/HR.
All the LR LESs with C M = 8 show an enhanced growth in magnetic energy compared to the C0 LR. However, we find that increasing the value of C T tends to reduce the observed magnetic field amplification, presumably due to an additional effective viscosity included in the momentum equation. The magnetic growth is thus more prominent in the CM8 (i.e., when SGS terms are included only in the induction equation), with its integrated energy being as high as the C0 MR run at t = 5 ms. The r.m.s. magnetic field for the CM8 case is comparable to C0 MR for most of the times, and significantly larger than the C0 LR at t = 5 ms, although all the simulations seems to reach comparable values at t = 10 ms, as it occurred with the total integrated magnetic energy. Notice that the effect of the gradient SGS model is most pronounced on the initial startup stage of the KHI, whereas the magnetic growth-rates on the kinematic phase does not seem to deviate much from the C0 LR, at it can be observed at t = 10 ms.
In the simulations CM8, CM8C1 and CM8C2, the remnant approaches a quasi-stationary stage at late times. Instead, CM8C4 and C8 show a different qualitative behavior and collapse to a black hole only after few milliseconds after the merger. This dependence on to the parameters of the SGS model is analogous to the sensitivity of the collapse time for short-lived HMNSs with numerical resolution, which has been observed previously both in HD [66,67] and MHD [68] simulations. Notice also that increasing these parameters C T and C N above 2 reduce the growth of the magnetic field energy.
Analyzing more in depth what the SGS model actually does, Fig. 6 displays the density and the magnetic field magnitude in the orbital plane z = 0, for the CM8 (top) and CM8C1 (bottom) cases, at t = {2.5, 5, 10} ms (from left to right). In both cases, the LR by construction does not allow the formation of very fine structures like the ones of HR (see bottom panels of Fig. 1). However, despite the lack of resolution, the SGS model is able to provide a growth of magnetic field up to local maximum values of ∼ 10 16 G at t = 5 ms, earlier than in the C0 LR. Also for these cases, filamentary structures start to appear at about t = 10 ms.
A comparison among the spectra is shown in Fig. 7, for CM8, CM8C1 and C0 LR, at t = {5, 10} ms after merger. Overall, these profiles are similar to those of iLESs, with the main difference given by their integrated values (i.e., the total magnetic energy). This again shows that at t = 5 ms CM8 is the most amplified one among the LR cases, between two and three orders of magnitude higher than the others for all wavenumbers (except the very high ones, which are dominated by numerical dissipation). The CM8C1 case exhibit a moderate growth of the magnetic energy spectra with respect to the C0 LR run, but considerably smaller than CM8. At t = 10 ms, the spectral distribution for these three cases is quite similar, and very close to equipartition at large wavenumbers. This is again consistent with Fig. 5, where these low-resolution simulations reach nearly the same magnetic energy values at late times. This behaviour on the magnetic energy spectra of LES was also found in our bounding-box simulations [43][44][45].
In summary, the LR LES that have a closer resemblance to the higher-resolution iLES (i.e., C0 MR/HR), at least at these early times, is CM8.

V. CONCLUSIONS
In this article, we showed the first results from LESs of BNS with the extended gradient model, already pre-sented for non-relativistic [43], special [44] and general relativistic [45] MHD box simulations of the KHI; here it has been implemented in a full GRMHD code in order to study the BNS merger scenario. Moreover, our code implements overall fourth-order accurate numerical schemes, while most existing GRMHD simulations rely on second-order accurate approaches (see advantages in the use of fourth-order schemes in [69]).
We have focused on the the magnetic field amplification within the first ∼ 10 ms after the BNS merger. And analyzed the role of numerical resolution in capturing these MHD turbulent-dynamo effects. With our bestresolved run reaching a grid-spacing of ∆ ∼ 37 m (in the finest level) and relying on the use of high-order numerical methods, for our highest resolution simulation we were able to demonstrate an amplification of r.m.s. values between 10 15 G, in the densest regions of the remnant, and 10 16 G, when the less dense outer envelope is also considered.
We have tested the gradient SGS model, by studying different values of its pre-coefficients C i , one in each equation. We show that the SGS terms on the induction equation acting in a moderate resolution (with ∆ ∼ 147 m) are able to mimic at least the magnetic growth of a betterresolved simulation (with ∆ ∼ 74 m).
In our previous works [43,45] we observed that the best results, with the gradient SGS model and our numerical schemes, were achieved by setting the constants C i approximately up to one order of magnitude larger than their theoretical values C i = 1. We concluded this was due to the intrinsic dissipation of our numerical scheme. However, in the present context, we found that if we set all C i = 8 there is an excessive dissipation in the momentum equation which prevents a rapid growth of the magnetic field during the turbulent regime, whereas it accelerates the collapse to a black hole of the remnant. Taking this into account, we have found that the best calibration (in terms of reproducing the higher-resolution magnetic field amplification) consists in rather high values C M = 8 but C T = C N ∼ 0-2. The reason for this remain to be clarified, and probably lies in the presence of two scales in the system. First, a fairly well resolved hydrodynamical one which includes the differential rotation and convection within the remnant. Second, a much smaller MHD scale, involving the turbulent dynamics originated by the KHI, which is still far from being resolved. As a consequence, only the coefficients for the magnetic field evolution needs to be artificially enhanced from their theoretical values, in order to maximize the effect of the SGS model to approach the results obtained with higher resolution simulations. Introducing too large contributions of the SGS model in the momentum equation leads to a higher effective viscosity which finally hampers the turbulence, partially suppressing the magnetic field amplification. In order to understand better the details and to disentangle the numerical and physical reasons for this, it would be helpful to implement and test our SGS model in other codes, for different numerical methods and/or

scenarios.
Regardless on the details of the SGS modeling, we have shown how the energy and magnetic spectra follow, respectively, the expected Kolmogorov (k −5/3 ) and Kazantsev (k 3/2 ) slopes, as in our iLES simulations. The magnetic spectra have a peak at small k, very different from a large-scale ordered field. Therefore, we warn against the widespread argument that an initially strong large-scale magnetic field can compensate the lack of ability to follow the KHI growth: the latter is intrinsically turbulent and can easily provide local maximum values exceeding 10 17 G, but contained in very small structures. This is also consistent with the fact that at early times where the kinetic dynamics deriving from the collision of the two cores is still dominating. The KHI triggers a quite isotropic turbulence, which destroys any largescale weak field. Only at later times, the dominating differential rotation should provide (via winding and MRI) the necessary energy injection at large scales which could partially order such strong but finely structured magnetic field. This can happen via inverse cascade and isotropy breaking favouring in particular the stretching of magnetic field lines in the azimuthal direction. We hope that our approach, which is applicable to any GRMHD problem, can be used in the near future to explore the postmerger phase dynamics. Given its potential in capturing the turbulent-dynamo effect at a much lower computational costs, it could be useful to further assess the production of large-scale fields that is required for the jet formation associated to the SGRB.  The explicit expression of the H tensors appearing in the SGS gradient terms were first obtained in Ref. [44], and then extended to GR in Ref. [45]. They can be written as , Ψ A = W 2 p dp d + ρ 2 dp dρ , where h = ρ(1 + ) + p is the enthalpy, E = hW 2 , Θ = E +B 2 , and the two gradients ∇ (on each term) symbolize spatial partial derivatives ∂ i (and ∂ j ), with "·" indicating contraction among them with the spatial metric γ ij .

Appendix B: Spectra calculation
It is illustrative to compute the radially-averaged spectrum of the kinetic and magnetic energy [43,70]. For a given field f defined in a periodic box of side L, we use common python functions to calculate its discrete fast Fourier transformf ( k) = Σ x f ( x)e −i k· x , where the sum is performed over the N 3 points equally spaced in each direction, with k j = n ∆k, where ∆k = 2π L and n ∈ [0, N/2] is an integer. Then, we calculate the solidangle-averaged values 4π < k 2 |f | 2 > k over the radial bins in the Fourier space, centered at k = {n ∆k}, which represent the power density per unit of angular wavenumber. This defines the kinetic and magnetic spectra, that we define for simplicity as in the non-relativistic case.
The calculation of the spectra is done by choosing the same 70-km-wide cube of the fifth FMR level, which encloses the remnant (i.e., almost the totality of the kinetic and magnetic energy of the system). Within this domain, the information analyzed is the one of the finest grid available. Since the domains of the AMR levels in the MR and HR cases are smaller than such domain over which spectra are calculated, we interpolate the values of the fields from the coarser levels, filling a regular mesh with the same grid spacing as the highest level present in that simulation. By construction, such interpolation have effects on the spectra limited to the smallest scales (highest k), which are not resolved outside the AMR levels. In addition, and in order to reduce the contamination from the rarefied atmosphere, the spectra is computed by considering only the velocity and magnetic fields in the regions with density larger than 6 × 10 9 g cm −3 , setting them to zero otherwise.