GRMHD large eddy simulations with gradient subgrid-scale model

The detection of binary neutron star mergers represents one of the most important astrophysical discoveries of the recent years. Due to the extreme matter and gravity conditions and the rich dynamics developed, it becomes a tremendous challenge to accurately simulate numerically all the scales present during the collision. Here we present how to study such systems by using large eddy simulations with a self-consistent subgrid-scale gradient model, that we generalized to the special relativistic case in a previous work and now extend to the general relativistic case. Adapted from nonrelativistic scenarios, the so-called gradient model allows to capture part of the effects of the hidden dynamics on the resolved scales, by means of a physically-agnostic, mathematically-based Taylor expansion of the nonlinear terms in the conservative evolution equations' fluxes. We assess the validity of this approach in bounding-box simulations of the magnetic Kelvin-Helmholtz instability. Several resolutions and a broad range of scenarios are considered in order to carefully test the performance of the model under three crucial aspects: (i) highly curved backgrounds, (ii) jumps on the fluid density profiles and (iii) strong shocks. The results suggest our extension of the gradient subgrid-scale model to general relativistic magnetohydrodynamics is a promising approach for studying binary neutron stars mergers, and potentially to other relevant astrophysical scenarios.


I. INTRODUCTION
One of the most important scientific discoveries of the last decade has been the first detection of gravitational waves from the coalescence of a binary neutron star (BNS) system by the LIGO/Virgo facilities [1,2]. The event GW170817 was followed by a broadband electromagnetic signal: a gamma-ray burst (GRB) after less than two seconds (GRB170817A) [3], a thermal infrared/optical radiation consistent with a kilonova [4,5] starting from hours after, and X-ray [6] and radio [7] afterglow, the latter been visible for many months. The observed signals are consistent with a BNS merging into a short-lived hypermassive neutron star (HMNS), which collapses to a black hole (BH) surrounded by an accretion disk. In particular, short GRBs are believed to be the high-energy manifestation of a relativistic jet, mainly powered by the rotational energy of the BH and the surrounding post-merger disk through the accretion process [8]. This single multi-messenger detection already led to important breakthroughs in our understanding of the universe, providing clues to fundamental physics such as the properties of matter at nuclear densities, emission mechanisms in extreme plasma conditions, and stringent tests of General Relativity (GR) (e.g. [9]) among others.
Magnetic fields are believed to play a fundamental role on the dynamics of these systems, especially in launching the magnetically dominated jets. However, the formation channels of BNS systems and their magnetic fields are quite uncertain. Observationally, we can infer the surface magnetic fields in old pulsars belonging to binary systems to be in the range of 10 8 −10 11 G [10]. Such magnetic fields are expected to be modified by local, small-scale amplification of the intensity by several orders of magnitude during the merger through several processes. The first one acting is the Kelvin-Helmholtz instability (KHI), which develops at the shear layer that forms between the external layers of the colliding NSs. The KHI acts only in the first few milliseconds of the merger and leads to a fast growth of small-scale magnetic fields, with a cut-off wavelength of the order of the size of the discontinuity layer. The simulations of BNS mergers with the current finest numerical resolutions, approaching ∆ ∼ 10 m (where ∆ is the grid size) [11,12], are still not able to fully capture the KHI: although they provide some notable amplification of small-scale magnetic structures, there is no yet sign of numerical convergence, meaning that dynamically relevant scales are still unresolved.
After the quick growth of the magnetic fields due to the KHI, there are two mechanisms associated to the differentially rotating HMNS which are expected to play an important role on longer timescales of 10-100 ms: magnetic winding, which creates and linearly amplifies a toroidal magnetic fields starting from a poloidal component, and the magneto-rotational instability (MRI). For the latter, the wavelength of the fastest growing modes is proportional to the magnetic fields. Some recent works manage to fully resolve the MRI by modelling only the post-merging phase and setting an initial largescale strong magnetic field of 10 15 −10 16 G (e.g., [12,13]). During and after these processes, the resulting magnetic fields intensity and configuration, the energy/helicity kinetic/magnetic spectra, and the formation of coherent structures are important features that affect the jet properties, the mass ejecta and the properties and fate of the newly formed HMNS [14]. Note however that the seed magnetic field configuration assumed in the simulations is usually a large-scale dipolar field, which artificially favours the appearance and amplification of an ordered magnetic field and the jet formation, but is inconsistent with the expected outcome turbulent spectra emerging from the KHI, which distributes magnetic energy over all scales. Therefore, there is great interest in understanding how these magnetic field amplification on the small scales can be transferred to the large ones, which might happen due to the combination of dynamo effect and inverse cascade of magnetic energy.
In brief, although the main magnetohydrodynamics (MHD) during and after the merger is qualitatively understood, only very accurate and demanding numerical simulations can model these processes and foresee their possible outcome. Such simulations are extremely difficult due to the intrinsic nonlinearity of the equations involved and on the wide range of spatial and temporal scales of MHD processes. To even further complicate the problem, strong gravity effects are important in NSs, implying that Einstein equations need also to be solved. Due to the unfeasible computational resources required, no study has so far consistently simulated all the phases of the magnetic dynamics described above with the required accuracy to address these issues. At contrast with the purely hydrodynamic case, the finite resolution hampers the ability of reproducing an efficient dynamo mechanism and therefore affects the large-scale dynamics, since MHD turbulence tends to distribute the magnetic energy (amplified at small scales) across the whole spectrum. In the absence of computationally viable direct numerical simulations with a spatial resolution able to capture all the relevant scales, the dynamo processes have been simulated in GRMHD by including additional terms in the induction equation [15][16][17]. While these approaches provide an effective growth of the magnetic fields, they have fundamental issues: they do not converge to the continuum version of the equations for very high resolution, and they need a tuning and switching on/off by hand of the extra terms, which is arbitrary and virtually impossible to be calibrated on physical grounds.
In order to find an inexpensive approximation able to capture the turbulent regime and allow back-scattering and dynamo mechanisms, in this work we borrow the explicit Large Eddy Simulation (LES) approach 1 , common in different areas like engineering and plasma physics (see, e.g., [18]), but scarcely utilized in relativistic MHD. In an explicit LES approach, the evolution equations for the resolved fields are supplemented with explicit subgrid-scale (SGS) terms which effectively model the small scales [19]. So far, to the best of our knowledge, only one work has so far introduced an SGS model into a GR simulation, limiting the study to the non-magnetic hydrodynamic case [20]. The author took into account the dissipative effects of magnetic fields on the post-merger differentially rotating fluid by implementing a purely viscous SGS model [21], originally designed for incompressible hydrodynamics. Note that, by construction, that particular SGS model only allows for an effective inclusion of dissipative mechanisms, and it cannot reproduce neither the relevant transfer of kinetic to magnetic energy (dynamo), nor the redistribution of the latter over a broad range of scales (back-scattering), relevant in BNS mergers and in most turbulent MHD scenarios. Therefore, in order to capture further effects arising from the SGS dynamics, one may need more sophisticated SGS models.
In this paper, we extend our previous SGS gradient models for nonrelativistic [22] and special relativistic [23] MHD, to full GR(MHD). The main advantage of the gradient SGS model (originally designed for incompressible fluids [24,25]) is that it does not rely on any physical assumption: it is based on the mathematical properties of the particular set of equations describing the dynamics. In this paper, we implement the proposed SGS model in a code using high-order accurate finite-difference numerical methods. We validate the approach by focusing on 3D bounding box simulations of the KHI. We explore a broad range of problem's parameters, including the most relevant aspects for the BNS scenario, namely: (i) spacetime curvature, (ii) sharp gradients in the initial fluid density profile, and (iii) strong shocks in the velocity field. We perform two different kind of tests: a-priori, comparing the residuals of the subfilter-scale (SFS) tensors from a high resolution run -without the SGS model-to the SGS tensors proposed by the model; and a-posteriori, in which the high-resolution run is now compared against a low-resolution run evolved using the SGS terms. The results of these tests support the present GR extension of our SGS model, and suggest it could be a promising tool to tackle the KHI developed during the merger in BNS. The model is general, and relevant also for the subsequent post-merger phase, as well as applicable to many other astrophysical scenarios.
This article is organized as follows. We summarize the formalism in §II and the numerical methods and analysis in §III. We then show the applications in bounding box KHI simulations in §IV, where we compare the results between a high-resolution, accurate simulation, with an explicit LES at a lower resolution. Finally, conclusions are drawn in §V.

II. LARGE EDDY SIMULATIONS OF GRMHD
LES techniques are classically employed in scenarios where not all the dynamically relevant scales can be resolved simultaneously: therefore, the large scales are sim-ulated, while the small ones are only modeled by the inclusion of SGS terms, effectively capturing the effects of these small scales on the resolved ones. In any numerical evolution problem, the effect of discretization over a finite-resolution grid can be associated to an effective filtering of the continuum equations, with the filter size being of the order of the grid-cell. In other words, the information inside the cell is partially lost by the discrete representation of the numerical solution, regarded in this context as a weighted average of the values within the cell. If one formally applies a filtering operator to the set of evolution equations, new terms will appear due to the non-commutativity of the filtering operator acting on the nonlinear terms of the equations. These new, apriori unknown terms, represent the SFS residuals. The aim of an explicit LES approach is to model them, by using the available resolved fields, in order to close the system of evolution equations. For completeness, we will briefly summarize some of the results presented in our previous works [22,23]. Let us start from a generic set of conservation laws, where C a is a set of conserved evolved variables and F ka (P ) is the flux along the direction k of the field C a , which can be expressed in terms of the primitive fields P a . The connection between conserved and primitive fields can be formally written as 2 C a = f a (P ) , P a := (f −1 ) a (C) ≡ g a (C) .
For any given variable C a (x, t), the filtering operator separates the resolved part C a (x, t) from the unresolved one C a (x, t) as, with being G(x − x ) the filter's kernel function, with compact support. The filtered version of the system in eq. (1) then reads where the SFS residuals (related to the fluxes) are defined by 2 Although the relation f a between conserved and primitive fields is known explicitly, the inverse function g a is not analytical in relativistic ideal MHD, where it needs to be solved numerically.
Here we have defined P a := g a (C), which can be directly computed from the evolved C a , as opposed to P a ≡ g a (C), which are the ones directly appearing after the filtering, but not known in a finite-size-grid simulation due to the unavoidable SFS information loss. The second term in eq. (5) is, therefore, the piece that needs to be conveniently modeled by SGS terms, which may depend only on the known available fields, i.e. P a and C a .
After a detailed study in nonrelativistic KHI simulations [22], for which different available SGS models were compared, we have focused our attention on the gradient model. At contrast with other SGS models, this model is physically agnostic, in the sense that it does not depend on guessing which kind of physical effects dominate at small scales (dissipation at a certain scale, inverse cascade, etc.) in a specific scenario. Additionally, it can be applied to any problem and set of equations, since it relies solely on the mathematical structure of the system, extrapolating, cell by cell and term by term, the SGS dynamics by using the information contained in the gradients of the fields involved. In that sense, it resembles high-order reconstruction methods, effectively capturing small-scale effects with no significant extra computational cost.
Following Ref. [23], the gradient SGS model was applied to the general system (1) above, yielding to a formal general solution, first order-accurate in the gradient expansion, namely where ξ := ∆ 2 /24 for a flat Minkowski spacetime and ∆ represents the numerical grid spacing. Note that the gradient model (as most SGS models) is designed to converge to the continuous solution in the high resolution limit, since the SGS terms are proportional to ∆ 2 . Below, we apply this formulation to the GRMHD equations. After summarizing our formulation of the Einstein and the GRMHD equations, we extend our previous special relativistic work [23] to the full GR setting, under the assumption that the metric components are smooth and slowly varying, as compared to the matter fields.

A. Einstein equations
The spacetime geometry is described by the Einstein equations, which can be extended in a convenient way by using the covariant conformal Z4 formulation (CCZ4) [26], namely where R ab is the Ricci tensor associated to the spacetime metric g ab and T ab is the total stress-energy tensor with trace trT ≡ g ab T ab . The Z a four-vector measures deviations from Einstein's solutions [27,28] (i.e., those satisfying the constraints). Notice that suitable damping terms, proportional to the parameter κ z > 0, have been included in order to enforce the dynamical decay of the constraint violations associated with Z a [29].
The covariant equations can be written as an evolution formalism by performing the 3 + 1 decomposition [30]. The line element can be decomposed as where α is the lapse function, β i is the shift vector, and γ ij is the induced metric on each spatial foliation. In this foliation, the normal to the hypersurfaces can be defined as n a = (−α, 0) and the extrinsic curvature as where L n is the Lie derivative along n a . Subsequently, a conformal decomposition is applied to the evolved fields, which consists of performing a conformal transformation to the metric and the extrinsic curvature, i.e., γ ij =γ ij /χ with χ being the conformal factor andγ ij having a unit determinant, and K ij into a trace trK and a trace-less partÃ ij . This transformation leads to two new constraints, which can also be enforced dynamically by including additional damping terms to the evolution equations [29]. The final set of evolution fields, together with the gauge conditions setting the choice of coordinates, can be found in [31].

B. GRMHD equations
Here we consider the usual stress-energy tensor associated to a magnetized, non-viscous and perfectly conducting fluid (e.g., [17]), and we use the usual relativistic units. The set of primitives variables describing this perfect magnetized fluid is given by P a = ρ, v i , , B i : ρ is the rest-mass density of the fluid, its specific internal energy, v i is the velocity vector and B i is the magnetic field. However, evolution equations are expressed in terms of a different set of fields, known commonly as the (densitized) conserved ones, given by The relation between these two sets of fields is given by where W = (1 − v 2 ) −1/2 is the Lorentz factor, the pressure p is defined through an Equation of State (EoS) and the enthalpy is defined as h = ρ(1 + ) + p. The evolution equations for conserved fields can be written as follows: where we have introduced a new field φ associated to the hyperbolic divergence cleaning of the solenoidal constraint (as in previous works, e.g. [31]). The fluxes have been written explicitly in connection with their special relativistic counterparts: Note that the ideal MHD condition E i = − ijk v j B k allows to write easily the previous fluxes and sources in terms of evolved and primitive fields. For numerical accuracy reasons, in our simulations we actually recombine two equations in order to evolve the conserved field The source terms, written already as a function of conformal variables, are

C. Filtered GRMHD equations
We now apply the filtering operator to the equations. First, notice that it is transparent to the partial derivatives, so it is straightforward to calculate the filtered MHD equations. Second, as it was mentioned at the beginning of the section, during a numerical simulation one can not calculate the filtered nonlinear , S k ( P )}, that is, the fluxes calculated with the primitive fields P a , obtained from the inversion of the filtered conserved fields C a . The difference results in SFS terms, which are corrections to the fluxes. Third, we assume that SFS terms arising from the metric functions are negligible compared to the ones arising from the turbulent fluids (see a more detailed discussion at the end of this section). Under these considerations, the filtered version of the GRMHD system can be written as: where we have introduced the SFS corrections inside each flux. The SGS gradient terms approximating such corrections can be written considering the non-trivial nonlinear dependencies of fluxes on the conserved variables, and following the general rule (6): where ξ = γ 1/3 ∆ 2 /24 scales appropriately to be consistent with the volume element of the spacetime. The parameter C is theoretically meant to be one, but, as we will see below and as shown in [22], the optimum value (i.e., the one best mimicking the SFS residuals in a LES) can differ depending on the numerical methods employed, being in general C 1. Therefore, it is advisable to leave C as the only free parameter of the SGS model. The set of the H tensors, after some algebraic manipulations, can be obtained explicitly by computing the set of equations written below, in the order in which they appear, where the quantities Ψ denote auxiliary fields which are used to simplify the implementation: where we have used the shortcuts 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 met-ric γ ij . Notice that, in order to include generic EoS, we also need the derivatives of pressure with respect to other thermodynamic fields. Therefore, the set of auxiliary fields in our implementation is given by P a = ρ, , p, v i , B i , φ, dp/dρ, dp/d . The expressions above extend our previous calculations [23], which assumed a flat Minkowski spacetime with √ γ = α = 1 and β i = 0, to the general relativistic case. Here, the metric is instead arbitrary, but is considered transparent to the gradient operators, which implicitly assumes that the metric is much smoother than the fluid fields. This approximation has two consequences: (i) we neglect any SFS correction to the Einstein equations, (ii) we deal with the SFS terms of the MHD equations (within the 3 + 1 formalism) essentially as in the nonrelativistic case, accommodating the curved spacetime elements (lapse function, shift vector and spatial metric) which appear in the flux terms.
In other words, we neglect the metric gradients in the fluxes C a = √ γ D, S i , U , B i , φ , and we neglect the SFS terms in the source terms. This is justified by the fact that the fields prone to turbulence, to which the evolution of resolved scales is sensitive, are the MHD ones, and not the metric. In other words, derivatives of the metric components are negligible compared to the ones of velocity, magnetic field, density and pressure (and their combinations). These are expected to be reasonable approximations when MHD turbulence is developed, even in the presence of strong curvature like in a BNS merger, because the large variations of the turbulent fields at small scales would be the dominant effects we want to capture with the SGS terms.
Finally, we also neglect the SFS terms in the solenoidal magnetic field constraint equation, since it is designed to keep φ under control but does not represent a meaningful physical turbulent field.

A. Numerical methods
The generation of our computational code is performed by using the in-house-built and publicly available platform Simflowny [32], developed to solve generic partial differential equations. It has a modular design which allows an independent and user-friendly implementation of the physical model, problem setup, and numerical schemes (time advance, space discretization and reconstruction methods). It generates efficient codes for the SAMRAI infrastructure [33,34], which allows an excellent scaling of parallel computation performances over thousands of processors.
The code has been deeply tested for different scenarios [31,35,36], including basic tests of MHD and GR. As in our previous works [22,23], we use the Method of Lines, which allows to address separately the time and the space discretization. The space-time evolution equations are discretized in space by using centered finite-difference, fourth-order-accurate operators for the derivatives, and sixth-order Kreiss-Oliger dissipation to filter the high-frequency modes unresolved in our grids. For the fluid part, we employ High-Resolution Shock-Capturing (HRSC) methods [37] to deal with the possible appearance of shocks and to take advantage of the existence of weak solutions in the equations. The fluxes at the cell interfaces are calculated by combining the Lax-Friedrich flux splitting formula [38], with the fifth-order, monotonicity-preserving reconstruction method MP5 (which showed better performances for turbulent fluids, compared to other methods [31]). The SGS terms are computed with fourth-order centered finitedifference operators. The time integration of the resulting semi-discrete equations is performed by using a fourth-order Runge-Kutta scheme, which ensures the stability and convergence of the solution for a small enough time step ∆t ≤ 0.4 ∆. We employ a divergence-cleaning scheme, which ensures that the deviations from the constraint, quantified by the dimensionless ratio between the volume-integrated values of the L2-norm of ∆( ∇ · B) and the magnetic energy, is always (typically) less 10 −4 . For further details about the divergence-cleaning formalism and numerical results, see our previous works [31,35].

B. A-priori tests
A first assessment of the performance of the SGS model, including its correct implementation, is the apriori test. In order to do it, one has to run simulations with a certain grid-spacing ∆ and consider a snapshot at a given time. Then, one has to apply a spatial filter to all the conserved fields; the simplest recipe is to use a simple average groups of S 3 f cells, where we define S f as the filter factor. The information lost in the filtering process contained in the scales ∈ [∆, S f ∆] (therefore, it is not the entire loss of information, which is by definition unknown), for a given nonlinear term, is the difference between the filter of the entire term and the nonlinear combination of the filtered values of the conserved fields. Such partial SFS information can be numerically evaluated as τ ( x f ) at each location x f of the filtered grid, and compared to the corresponding SGS model τ ( x f , C). The performance is then evaluated by the classical Pearson correlation coefficient P between the SFS and SGS quantities, together with the best-fit value of the free parameter C in the SGS terms, as defined by eq. (22). This best-fit value, which theoretically should be close to 1, can be easily calculated by using the classical analytical formula for a linear regression fit. The calculation of P and C is done for each SFS/SGS term. In our previous works [22,23], the gradient model showed very good performance (P 0.8) for all terms (better than other SGS models available in the nonrelativistic case), and degraded to P 0.5 only for quite large filter factors S f 8. Below we will show similar results for relativistic KHI in a curved background.

C. A-posteriori tests
Although informative, the a-priori tests only provide how well a SGS model fits the instantaneous SFS residuals between ∆ and ∆ f in a snapshot, over a certain range of parameters. The second class of tests is instead more challenging, since it considers the feedback of the included SGS terms, accumulated over time. It consists in comparing a low-resolution+SGS simulation with reference results. The latter are ideally represented by an analytical solution, or by numerical solutions to which simulations show numerical convergence. However, most of the times they are not available, so one has to take a high-resolution simulation as a reference.
The comparison between simulations can be done at different levels. A first one is represented by the volumeintegrated quantities, especially the magnetic energy in our case. At a more detailed level, very informative whenever there is turbulence, it is illustrative to compute also the radially-averaged spectrum [22,39]. 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 radial wave-number. This defines the kinetic and magnetic spectra, that we define for simplicity as in the nonrelativistic case. Below, we will do a-posteriori tests with different resolutions, activating the SGS terms with different values of the free parameter C.

IV. NUMERICAL SIMULATIONS
We test our approach in bounding box simulations, analogous and generalized compared to what can be found in the literature for the nonrelativistic cases [22,40,41], and in our previous special relativistic work [23]. Here, we consider a fixed metric, setting to zero the time evolution of the spacetime fields, evolving only the MHD part. This allows us to test the effect of a curved spacetime on the turbulence, but without the additional complication of a dynamical metric. In all the following simulations, we evolve the system at least to the saturation time, after which the solution approaches a stationary regime, with the turbulence completely developed and slowly decaying due to numerical dissipation.

A. 2D single-vortex KHI
Since turbulent MHD is intrinsically a nonlinear and complex phenomenon, our first benchmark test considers only a simplified case with the creation of a single vortex via KHI in a flat metric, which is just the relativistic analogous of previous studies [22,40].
We set a rectangular domain in the x − y plane, centered in the origin, with side length L x = 1, L y = 2, with a number of points (N, 2N ) equally spaced in the two directions. Boundary conditions are periodic in x-direction and radiative in the y-direction. The initial conditions consist of ρ = 1 and p = 0.24 everywhere, with a shear layer of thickness a l representing the initial discontinuity at y = 0, and the following velocity and magnetic fields: where v 0 is the shear velocity, and δv y v 0 represents the amplitude of the initial perturbation, which consists of a single mode k x . We also consider an ideal EoS p = (Γ − 1)ρ with Γ = 5/3. This basic case has the advantage that there is only one initial excited mode, so that the complexity is greatly reduced and we can isolate better the effects of the SGS model implementation. The outer boundary conditions allow perturbations perpendicular to the shear layer to go out of the system, thus helping the creation of a single vortex. This corresponds to a spectrum power distributed uniquely among k x and, at much lesser extent, its higher harmonics. In other words, turbulence does not develop completely and the relevant scales to be resolved are basically set by k x . Therefore, the numerical convergence is ensured for high enough resolutions ∆ 2/k x . Indeed, in previous works [22,40] it is evident how the numerical resolution affects the growth rate.
In order to evaluate convergence, we consider the evolution of the volume-integrated y-component of the kinetic and magnetic energies, tracers of the turbulence, and defined as in the nonrelativistic case for simplicity (relativistic corrections would not change the comparison between different resolutions): After the initial triggering stage, both quantities tend to grow exponentially until a saturation level. The growth rate and the saturation value of E ky are physically controlled by the values of v 0 , a l , k x and B 0 , but not by the amplitude of the perturbation, as long as δv y v 0 . Analytical values of the growth rate of E ky are known for the nonrelativistic case only [42]. At saturation, the solution approaches a stationary regime, with a fullydeveloped vortex showing periodic oscillations in both directions, but maintaining the shape. The full development of the vortex occurs only under certain conditions of Mach number (roughly of order 1) and Alfven number (i.e., the initial large-scale magnetic field has to be weak enough to leave the fluid free to develop small structures), as explained in detail for the nonrelativistic case in [40].
We set up a test with parameters δv y = 0.01, a l = 0.05, k x = 1, v 0 = 0.5, and B 0 = 0.0005. With such a choice, numerical convergence of the solution is reached with only few hundreds of points in each direction, allowing fast tests (see [22,40] for an extended exploration of parameters). In both case, we employ 3 resolutions: N = {25, 50, 100}. In the case of the lowest resolutions, we apply the SGS model with different values of the free parameter C = {0, 1, 2, 4, 8}.
We monitor both E ky and E my as a function of time to calculate their corresponding growing rates, obtained by fitting these curves to a functional ∝ e λ k t and ∝ e λmt , respectively, where λ parameters are the growth rate values. After the initial transient, the growth rate show convergence to a well-defined value as the resolution is increased (i.e., λ k = 0.97 and λ m = 0.81 for N = 25, converging to λ k = 1.22 and λ m = 1.09 for N = 100). Interestingly, a similar growth rate and evolution can be achieved with very low resolution N = 25, using a subgrid model with C ≈ 8 (i.e.,λ k = 1.10 and λ m = 0.97). In order to see some effects compared to the case without SGS, we need at least C = 4. The evolution of E ky is explicitly displayed in Fig. 1 only for the extreme cases C = {0, 8}. We also repeated the same tests in the nonrelativistic case (the models grw1-2-3-4 and 10 of [22,40]), finding an optimal value of C ∼ 3. This stresses the fact that the exact value of C would need a finetuned calibration, which unfortunately depends on both the physical model and initial conditions and the intrinsic dissipation of the numerical methods used. However, generally speaking, the inclusion of the model (typically with C ∼ 1 − 10) will anyway improve the accuracy in reproducing the effects of the SGS dynamics.

Setup
Next, we consider a fully turbulent 3D problem. We set our problem in Cartesian coordinates considering a periodic box [−L/2, L/2] 3 . The primitive fields read initially: where a l is the mixing layer scale, y l is the distance of the shear layers to the plane y = 0, σ y and σ z are the extension scale of the initial perturbation in the ydirection and the profile of v z , respectively. The main flow is initially given by v x0 and v y0 , having opposite signs across the layer. The standard values that we consider are L = 1, y l = 1/4, a l = 0.01, v x0 = 0.5, v y0 = 0, p 0 = 1, B x0 = 0.001 and σ 2 z = 0.1. The initial perturbation, δv i , is a superposition of single-mode perturbation with a number of nodes n i ∈ [1, N/2], periodic in the boundary box, δv i = δv 0 sin(2πx i n i /L), with n x = 11, n y = 7, n z = 5 (different prime numbers-based initial modes in each direction ensure the excitation of the entire spectrum of modes, see [22] for a discussion), δv x = δv z = 0.01, δv y = 0.1, σ 2 y = 0.01. We employ an ideal gas EoS with Γ = 4/3.
As noted in previous works, the fully-developed (i.e., coming from a random/multimode initial perturbation) the KHI is known to grow faster for smaller scales, and due to the absence of physical viscosity in this benchmark test, there is no physical lower limit to the scale of the excitable modes. Since our initial perturbation is designed to excite the entire spectrum of modes, we do not expect a numerical convergence in the growth phase, since the more we refine the grid, the more fast-growing excited modes will be included.
We validate the gradient SGS model for GRMHD under a broad range of physical conditions, extending some of the tests already performed in previous works [22,23]. First, we vary different values of the maximum density ρ 0 + ρ 1 ∈ [2 : 50] (keeping the minimum density ρ 0 − ρ 1 = 1), representative for small to extreme initial density jumps. Secondly, we consider a fixed, nonflat metric by setting the conformal factor to a Gaussian shape, with given depth χ 0 and width σ χ , being r is the radial distance to the center of the box. We consider depth up to χ 0 = 0.8, which is roughly the value reached in realistic BNS merger simulations. We set the width either to σ χ = 0.5 or σ χ = 3.0, leading to a total variation of χ between χ = [0.20 − 0.90] and χ = [0.20 − 0.24], respectively. Thirdly, we test different values of the initial vertical velocity v y0 , such that strong shocks, traveling along the y-direction, are initially induced in the fluid.

A-priori tests
The results for the a-priori tests are summarized in Fig. 2, displaying P (top) and the best-fit C (bottom). In all these cases, for any given time the a-priori tests give C ∼ 1 − 2 and P 0.8 for the minimum filter factor S f = 2. Performances gradually degrade for increasing loss of information, e.g. higher S f (see [22,23] for a detailed exploration of other values of S f and N ).
This first test further confirms that the gradient SGS model is still valid when a high density jump is included, when shocks are present on the solution, or when nontrivial metric factors are incorporated. More importantly, it validates the smooth-metric assumption and reinforces that the presence of a curved spacetime background does not compromise the good performance of the model.

A-posteriori tests
We analyze now the simulations performed by including the gradient SGS model in the evolution equations, that is, with C ≥ 1 (if not specified, the same indicated value of C is applied to all SGS terms). For concreteness, we focus on the case with a flat metric background and ρ 0 = 1.5, ρ 1 = 0.5. We perform simulations with no SGS modeling (C = 0) for different resolutions N = {128, 256, 512, 1024}, and four LES+SGS with a low resolution, N = 128, for C = {1, 2, 4, 8}. We also perform a high-resolution simulation N = 512 with C = 8 to check the convergence of solutions when LES and SGS modeling are included. Note that higher values of C result in numerical instabilities, besides arguably becoming physically and mathematically not consistent. From now on, we will label the different simulations through the values of N and C, e.g. KH128C4 for N = 128 and C = 4.
In Fig. 3 we show the density distribution and the zcomponent of magnetic field B z , which develops due to the turbulent instability, at t = 6 for KH128C0, KH128C8 and KH1024C0. As expected, the intensity of the magnetic field rises as the value of C is increased (i.e, compare the left and middle panels). The same effect is achieved by increasing the resolution (right panel), getting up to an order of magnitude higher intensities than in the lower resolution. Although the global effects are similar, notice that the KH128C8 case cannot show fine details due to its intrinsic low resolution: this does not prevent it from reproducing the feedback unresolved dynamics, thanks to the SGS terms.
In Fig. 4 we can see this trend quantitatively: the volume-integrated magnetic energy rises more if we increase either the resolution or the value of C. A quasistationary state, indicated by saturation of the integrated magnetic energy, is only achieved at the end of the simulation in the high-resolution cases and in the KH128C8. This means the turbulence has been completely developed at this time. Nota that in the KH128C8 case we get a magnetic energy roughly an order of magnitude above the KH128C0 one and it is similar to KH256C0. This means we are capturing the same magnetic energy with half of the resolution. The magnetic energy difference between the KH512C0 and KH512C8 simulations is less than the difference between the KH128C0 and KH128C8 ones, showing two important results: (i) the SGS effects actually converge to zero as the resolution increases, and (ii) the solutions of the simulations with the SGS modeling actually converge to those obtained with higher resolutions, but using much less computational resources.
In Fig. 5 we show the spectra of both magnetic (dashed line) and kinetic (solid line) energy at different time snapshots corresponding to t = {6, 10, 16, 20}. To improve the visualization of the results, we do not include the C = {1, 2, 4} cases with N = 128: they lie in the middle between C = 0 and 8. The maximum of the energy spectra is reached before the spectral knee (located at k a few times smaller than k max ≡ π/∆), after which the slope steepens due to the intrinsic dissipation of the finite-difference method. At t = 6 we are still at the beginning of the simulation and eddies have not been developed completely yet. For this reason, al-most all magnetic energy spectra have similar values. At t = 10 we can clearly see the effects of the SGS model in the KH128C8 simulation, whose magnetic energy spectra arises above the KH256C0 case and reaches similar values to the KH512C0 one. At t = 16, the spectra of the high-resolution case KH512C0 and the KH128C8 are getting closer at large scales (low wave-numbers). At t = 20 the system has reached a quasi-stationary state and the energy spectra is comparable for the highest resolution simulations KH512C0, KH512C8 and the lowest resolution one with SGS modeling KH128C8. Notice that at all times the spectra of KH512C8 has still the highest magnetic energy, as it combines the small scale dynamics (i.e., accessible because of the high resolution) and the SGS effects.
In Fig. 6 we compare the energy spectra obtained from simulations with two resolutions, N = {128, 256}, for different metric backgrounds. First, it is represented the spectra energy with χ = 1 (left), then a curved metric with a range of χ = [0.20 − 0.90] (middle), and finally an extreme case with χ = [0.20 − 0.24] (right panel). With χ = 1, the magnetic energy spectra for KH128C8 is above the case KH256C0, as it was already shown previously. The same behavior, significantly moderated, is also observed when there is a large variation of the curved metric in the simulation domain (middle panel). In this case, the magnetic spectra for KH128C8 also rises with respect to KH128C0, but it is only slightly above the KH256C0 case. Finally, when the metric curvature is really strong and homogeneous (right panel), turbulence is strongly suppressed, since the fluid can not move freely on such strong gravitational field. Therefore, all the curves in the magnetic energy spectra are comparable. The fact that the KH128C8 and KH256C0 cases also coincide here reinforces the idea that LES with the SGS gradient model is actually working as expected, not producing artificial enhancements of magnetic energy. case we additionally set H v = 0. Clearly, almost all the rising of the magnetic energy spectra of KH128C8 is due to the Newtonian analogous of the magnetic SGS terms, achieved with H v = 0 and C M ≥ 1. This means that this contribution, which corresponds to the nonrelativistic limit dominates over the other terms (at least for these physical conditions), indicating an easy way to model the feedback of SGS dynamo at small scales.

V. DISCUSSION
In this paper, we have introduced the first GRMHD explicit LES with the physically-agnostic gradient SGS model. The aim is to better capturing small-scale effects of MHD turbulence at relevant astrophysical scenarios. Here we have focused on validating our numerical code, which combines our proposed SGS model for GRMHD with the use of high-order accurate numerical schemes, in box simulations developing the KHI with different parameters.
The extension from our previous studies of nonrelativistic [22] and special relativistic [23] MHD, to full GR, was carried under the following assumptions: (i) the SGS model is constructed starting from the filtering operation (associated to the numerical discretization) at the 3 + 1 decomposition level; (ii) when filtering the MHD equations, the metric components are considered "transparent" to the gradient operators involved in the SGS tensors; and (iii) all SGS corrections arising from the Einstein equations are neglected. The assumption (i) means the construction of our model is not fully co-variant. Indeed, we argue that the need for a SGS model (to effectively capture part of the missing small scale dynamics) arises from the discretization itself, which is not an invariant operation. Thus, we find quite natural the SGS model adapts instead to the particular discretization at the 3 + 1 level. It is worth to emphasize that the gradient model (by construction, as most SGS models) does not modify the continuous limit, nor the principal part of the evolution system. In this sense, the inclusion of the SGS terms is analogous to the numerical reconstruction methods, which also introduce a violation of the GR covariance. The approximation (ii) has to do with the fact that the metric components can be considered smooth in comparison with the turbulent MHD fields, and thus, the contribution from their gradients in the SGS tensors should be sub-dominant. The last assumption (iii), linked to the previous one, states the dominant SGS corrections on the turbulent MHD variables are not expected to produce significant deviations on the spacetime geometry. We expect these assumptions to be fairy reasonable in the context of BNS mergers, where the large variations of the turbulent fields at small scales should contribute the dominant effects we aim to capture with the model. We test our approach using 2D and 3D bounding box simulations of the relativistic KHI in a curved background, finding essentially the same results as in their nonrelativistic and special relativistic counterparts [22,23] and supporting the present GR extension of our SGS model. We have considered a variety of different scenarios with different resolutions, and compared via apriori tests the SFS residuals within a certain scale range [S f ∆, ∆] with the proposed SGS model. We obtain bestfit values C ∼ 1−2 (as expected) and high values P 0.8 for S f = 2. This is consistent with our previous studies and indicate that the SGS is suitable to fit well the SFS residuals down to S f 4 at least (gradually degrading for smaller scales). The most important novelty of the a-priori tests presented here are the inclusion in the problem of nontrivial metric components (i.e., curvature) and strong shocks in the fluid initial conditions. Moreover, we have also performed LES with the proposed SGS models, thus allowing a-posteriori tests, in which a high-resolution run is compared against lowresolution explicit LES with different values of C ∼ O(1 − 10). The integrated magnetic energy in the KHI rises as one increases the resolution, while similar behaviour is shown to occur when increasing C at a lower resolution run including SGS terms. When spectra are looked, the N = 128 run with C = 8 can achieve comparable magnetic energy spectrum with respect to highresolution cases. Note that C significantly larger than 1 are needed. In order to explain this, note that in our previous work [22], we saw that including the next-leading order terms in the Taylor expansion of the SGS tensors give negligible differences in both magnetic energy and spectra, as compared with the first order approximation. We can conclude the required relatively high values of C 1 can be associated to the intrinsic dissipation of the numerical finite-difference scheme employed (being likely smaller in case of using the intrinsically less dissipative spectral methods).
With the tests here presented, we have consistently extended to GRMHD the good performance of the gradient SGS model already found in nonrelativistic MHD [22,25] and special relativistic MHD [23]. Overall, we have shown how the implementation of the gradient SGS terms in GRMHD LES can indeed reproduce part of the feedback given by unresolved dynamics over the large scales.
The main applicability relies on the fact that the mere inclusion of the gradient SGS terms allows the saving of at least one order of magnitude in computational time, providing similar numerical results. This approach is then promising and supports our goal of applying the proposed gradient SGS model for GRMHD to realistic astrophysical scenarios where small scale dynamics can be crucial, like the supposedly spectacular growth of the magnetic field occurring during the merger of BNS, and the consequences it has for the production of observable jets.