Gradient sub-grid-scale model for relativistic MHD Large Eddy Simulations

TL;DR

This paper introduces a novel gradient sub-grid-scale model for relativistic MHD large eddy simulations, enabling better representation of small-scale physics in astrophysical turbulence, especially in neutron star mergers.

Contribution

It extends the gradient model to relativistic MHD equations, providing the first comprehensive sub-grid-scale terms applicable to any equation of state.

Findings

Model accurately fits residuals from high-resolution simulations.

Effective in capturing small-scale dynamics in relativistic MHD turbulence.

Potential extension to general relativity for improved binary merger simulations.

Abstract

MHD turbulence is likely to play an important role in several astrophysical scenarios where the magnetic Reynolds is very large. Numerically, these cases can be studied efficiently by means of Large Eddy Simulations, in which the computational resources are used to evolve the system only up to a finite grid size. The resolution is not fine enough to capture all the relevant small-scale physics at play, which is instead effectively modeled by a set of additional terms in the evolution equations, dubbed as sub-grid-scale model. Here we extend such approach, commonly used in non-relativistic/non-magnetic/incompressible fluid dynamics, applying the so-called gradient model to a general set of balance-law equations, that includes the relevant case in which a non-trivial inversion of conserved to primitive fields is needed. In particular, we focus on the relativistic compressible ideal MHD scenario, providing for the first time (and for any equation of state) all the additional sub-grid-scale terms. As an application, we consider box simulations of the relativistic Kelvin-Helmholtz instability, which is also the first mechanism responsible for the magnetic field amplification in binary neutron star mergers and cannot yet be fully captured by the finest-grid and longest simulations available. The performance of our model is numerically assessed by comparing it to the residuals arising from the filtering of high-resolution simulations. We find that the model can fit very well those residuals from resolutions a few times higher. Although the application shown here explicitly considers the Minkowski metric, it can be directly extended to general relativity, thus settling the basis to implement the gradient sub-grid model in a GRMHD binary merger. Our results suggest that this approach will be potentially able to unveil much better the small-scale dynamics achievable in full GRMHD simulations.

Equations

$$f(\mathbf{x},t)=\overline{f}(\mathbf{x},t)+f^{\prime }(\mathbf{x},t).$$

$$\overline{f}(\mathbf{x},t)={\int }_{-\infty }^{\infty }G(\mathbf{x}-{\mathbf{x}}^{\prime })f({\mathbf{x}}^{\prime },t)d^3{x}^{\prime }.$$

$${\partial }_t{U}^a+{\partial }_k{F}^{ka}(U)=0,$$

$${\partial }_t{\overline{U}}^a+{\partial }_k{F}^{ka}(\overline{U})={\partial }_k{\overline{\tau }}^{ka}$$

$${\overline{\tau }}_F^{ka}:=F^{ka}(\overline{U})-\overline{F^{ka}(U)}.$$

$${\partial }_{t}u+\frac{1}{2}{\partial }_x{u}^2=0.$$

$${\partial }_t\overline{u}+\frac{1}{2}{\partial }_x{\overline{u}}^2={\partial }_x\overline{\tau },\overline{\tau }\equiv {\overline{u}}^2-\overline{u^2}.$$

$$\overline{\tau }={\overline{u}}^2-\overline{u^2}\approx \nu {\partial }_x\overline{u},$$

$$G(|\mathbf{x}-{\mathbf{x}}^{\prime }|)=\prod _{i=1}^D{G}_i(|x_i-x_i^{\prime }|),$$

$$G_i(|x_i-x_i^{\prime }|)=\{\begin{array}{l}1/{\Delta }_f\mathrm{if}|x_i-x_i^{\prime }|\le {\Delta }_f/2\\ 0\mathrm{otherwise}\end{array}$$

$$G_i(|x_i-x_i^{\prime }|)={\left(\frac{1}{4\pi \xi }\right)}^{1/2}\exp \left(\frac{-|x_i-x_i^{\prime }{|}^2}{4\xi }\right),$$

$$\hat{G}(\mathbf{k})=\exp (-\xi {\mathbf{k}}^2),$$

$$\hat{G}(\mathbf{k})$$

$$\frac{1}{\hat{G}(\mathbf{k})}$$

$$\overline{f}\equiv G\ast f$$

$$f\equiv G^{-1}\ast \overline{f}$$

$$\overline{fg}$$

$$\overline{fgh}$$

$$\overline{f(g)}$$

$$\overline{\tau }\equiv {\overline{u}}^2-\overline{u^2}\backsimeq -2\xi {\partial }_x\overline{u}{\partial }_x\overline{u}.$$

$${\partial }_t{C}^a+{\partial }_k{F}^{ka}(P)=0,$$

$$C^a=f^a(P),P^a:=(f^{-1}{)}^a(C)\equiv g^a(C).$$

$${\partial }_t{\overline{C}}^a+{\partial }_k{F}^{ka}(\tilde{P})={\partial }_k{\overline{\tau }}^{ka},$$

$${\overline{\tau }}_F^{ka}:=F^{ka}(\tilde{P})-\overline{F^{ka}(P)}.$$

$${\overline{\tau }}_P^a:=g^a(\overline{C})-\overline{g^a(C)}.$$

$$\overline{F^{ka}(P)}\backsimeq F^{ka}(\overline{P})+\xi \left({\nabla }^2{F}^{ka}(\overline{P})-\frac{d{F}^{ka}}{d{\overline{P}}^b}{\nabla }^2{\overline{P}}^b\right),$$

$$F^{ka}(\overline{P})\backsimeq F^{ka}(\tilde{P})+\xi \frac{d{F}^{ka}}{d{\tilde{P}}^b}\left({\nabla }^2{\tilde{P}}^b-\frac{d{\tilde{P}}^b}{d{\overline{C}}^e}{\nabla }^2{\overline{C}}^e\right),$$

$$\frac{d\widetilde{P}}{d\overline{C}}\equiv\left(\frac{dC}{dP}\right)^{-1}\bigg{\rvert}_{\widetilde{P},\overline{C}}$$

$${\tau }_F^{ka}=\xi \left(\frac{d{F}^{ka}}{d{\tilde{P}}^b}\frac{d{\tilde{P}}^b}{d{\overline{C}}^e}{\nabla }^2{\overline{C}}^e-{\nabla }^2{F}^{ka}(\tilde{P})\right)$$

$${\tau }_F^{ka}=-\xi \nabla \frac{d{F}^{ka}}{d{\overline{C}}^b}\cdot \nabla {\overline{C}}^b$$

$$C^a=f^a(P)=\left\{\rho ,\rho v^i,\rho ϵ+\rho \frac{v^2}{2}+\frac{B^2}{2},B^i\right\}$$

$${\partial }_t\rho +{\partial }_k{N}^k=0,$$

Full text

I Introduction

The first long-awaited simultaneous detection of gravitational waves (GW) and electromagnetic radiation from a binary neutron star (BNS) merger [1] contains a wealth of data, including multi-wavelength electromagnetic radiation ranging from radio to gamma [2]. The existing models of BNS mergers allow to compare the observed GW signal with the theoretical waveform, thus inferring the chirp mass, and, secondarily, constraining the mass ratio and the deformability of nuclear matter. This analysis, together with the associated short gamma-ray burst and the kilonova emission [3], served to constrain some unanswered questions, regarding for instance the equation of state of nuclear matter [4], the maximum mass supported by a neutron star [5, 6] and its radius [7], the amount of ejecta [8, 9] and the contribution of BNS mergers to heavy elements production via r-processes [10].

Given the intrinsic difficulties of numerical general relativistic magneto-hydrodynamics (GRMHD) simulations, a few works have consistently included the full set of combined Einstein equations and MHD equations to solve such scenarios. Although the main features of the dynamics are quite clear from these simulations, some interesting and more subtle questions regarding the role of magnetic fields and the process they undergo after the merger are still unclear. In particular, during the merger, the Kevin-Helmholtz instability (KHI) is triggered at the shear layer between the colliding stars, and induces a fast growth of any seed magnetic field [11, 12, 13, 14]. The instability develops faster for small-scale perturbation, with a cut-off wavelength of the order of the shear layer thickness (likely being $\sim$ meter or less). Additionally, other two mechanisms related to the remnant’s rotation are supposed to play an important role at longer timescales: the magnetic winding, which creates and amplifies toroidal magnetic field starting from a radial component, and the magneto-rotational instability (MRI) [15].

A fundamental component of most theoretical models of binary NS mergers involved the formation of a jet perpendicular to the orbital plane, which is expected to be produced due to the accretion of the disk onto a central black hole. This jet is supposed to power the short gamma ray burst [16, 17], and it is often thought to require or be favoured by the presence of a strong, large-scale (dipolar) magnetic field. The creation of the jet is not trivial, with some simulations succeeding in generating it when starting with a strong magnetic field [18], but most of them not yet able to see it (see for instance [19]). In any case, the growth and creation of a large scale magnetic field has not yet been observed, not even in simulations with the highest numerical resolutions of O(10m) [14]. Even in that case, they are not yet able to fully capture the KHI, whose smallest and fastest-growing scales are of the order of the discontinuity layer (probably a few meters or smaller).

Limitations are being slowly overcome thanks to the growing computational resources. However, the full capture of all scales by means of direct numerical simulation (DNS) is still relatively far out of reach even for the most advanced codes and infrastructures. Thus, it is well worth to study alternative approaches that try to simulate the physical processes at play with a much lower computational cost. There have been several attempts to effectively model the dynamics occurring in the small scales. For instance, Ref. [20, 21] employs the viscous hydrodynamic equations to model the differentially rotating remnants of binary neutron stars mergers. Other works concern dynamo mechanisms, in which the magnetic field growth can be modeled by a simple algebraic modification of the induction equation [22, 13].

An alternative approach to the problem comes from computational fluid dynamics and it can be summarized in the following statements. First, any numerical simulation can be seen as a filtered version of the evolution equations. Secondly, the filtering of any non-linear combination of the evolved fields (e.g., in the fluxes of the MHD equations) implies the appearance of residual terms in the evolution equations, corresponding to the unresolved dynamics in the small scales. In this approach, commonly known as explicit Large Eddy Simulations (LES), relies first on this separation of the scales of the solution (i.e., resolved and unresolved), combined with an explicit sub-grid scale (SGS) model to account for the dynamics occurring at the unresolved small scales [29, 30, 31, 32]. Although LES are commonly used in many engineering applications, including combustion, acoustics, and simulations of the atmospheric boundary layer and other fluids, its extension to astrophysical magnetized plasmas is less spread [23]. Its application to strongly self-gravitating fluids is even more limited. To our knowledge, only Ref. [24] has so far explicitly developed explicit LES for the relativistic hydrodynamic equations with a covariant form of the viscous SGS model, proposed by Smagorinsky in the sixties [25] and commonly used in other contexts. Since the Smagorinsky model is by definition purely dissipative, it easily allows numerical stability, but it cannot capture any possible backscatter towards larger scales. This is one of the reason why more elaborated models have been considered in the non-relativistic case.

In this paper we propose a more general formulation of the so-called gradient model [26, 27], that we already assessed and implemented in the non-relativistic case [28]. It allows to capture part of the unresolved small-scale dynamics, regardless of the specificity of the problem. The main advantage of this model is that it does not rely on any physical assumption in the functional form of these residual terms. First we have to extend the LES formalism to generic conservative evolution equations like the relativistic MHD ones. This system has some particularities; first, the evolution equations of the evolved fields depend on some other quantities, usually called primitive fields. Second, the relation between the conserved and the primitive fields is explicit and analytical, but the inverse is usually not. After writing the LES for generalized conservation equations, we will extend the gradient model to this system.

Let us emphasize again our aim: assuming that the relativistic MHD equations can be written as an evolution system by using the 3+1 decomposition, how can we efficiently mimic the effect of the small-scale dynamics, which are missed due to the finite numerical resolution? The goal of this work is, hence, to extend the use of a well established numerical technique (i.e. LES and a particular SGS model) to relativistic MHD. It is worth to discuss here some possible limitations of our approach. First, our proposed SGS model for relativistic MHD is not gauge invariant. The construction of the SGS model (and, more generally, any LES) begins from the filtering operation, which is to be associated in our context with the discretization of the equations at the 3+1 level. Thus, the filtering operator is not invariant but depends on the foliation, and therefore, the SGS tensors arising from such filtering are not gauge invariant neither. That means that different gauges will have different averages and, as a consequence, different SFS tensors. Since the SGS terms are intended to effectively capture part of the missing small scale dynamics (contained in the SFS tensors), it is quite natural to expect a gauge dependent SGS model, adapted instead to the particular discretization of the system at the 3+1 level. It is worth to reiterate that the gradient model 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 commonly used in the HRSC schemes, with the important difference that a careful analysis of the equations needs to be done in order to derive the functional forms of the SGS terms. Once the latter is done, it can be applied to any system of hyperbolic PDEs after the 3+1 decomposition. Another important point that we want to stress is that the LES non-relativistic MHD model has the mathematical properties (conservation of total energy, magnetic helicity, approximate conservation of the cross helicity) expected of a model derived from the MHD equations by an averaging operation [33]. Similar extensions would be expected to hold in the relativistic case. Finally, probably the strongest limitation which could affect the effectiveness of SGS model is the presence of strong shock in the simulations. Although the High-Resolution-Shock-Capturing methods that we use can deal with them by reducing the order of accuracy, this introduces an additional dissipation in the numerical scheme that overcomes the effects of the SGS model, as it is discussed for instance in [34].

In § II, we summarize the formalism of the filtering process in LES, very common in a plethora of hydro-dynamical simulations related to different fields, but less known in relativistic hydrodynamics. § III summarizes the basics of the gradient SGS model. In § IV, we extend the LES formalism and the SGS gradient model to a general system of conservation laws, applied to non-relativistic and relativistic compressible ideal MHD in § V. We validate the proposed model in a set of box simulations of the relativistic KHI in § VI, by comparing the SGS model to the residuals coming from filtering high-resolution simulations in § VII. Finally, we summarize our results in § VIII.

II Large Eddy Simulations and sub-filter-scale residuals

From a formal point of view, the effect of finite resolution in a numerical simulation is equivalent to a low-pass spatial filter, with a filter size of the order of the grid cell. Applied to a field $f({\bf x},t)$, the filtering operation separates the resolved part from the sub-filter scale (SFS) residuals:

[TABLE]

Indicating the filter kernel with $G$, the filtering operator over a field $f$ can be written generically as

[TABLE]

Large-Eddie-Simulations consist on applying the filter to the evolution equations under consideration and write them down as a function only of the resolved fields.

Let us consider a generic evolution system written as

[TABLE]

where $U^{a}$ is a set of evolved fields and $F^{ka}(U)$ is the flux along the direction $k$ of the field $U^{a}$. The equation for the filtered fields can be obtained by applying the filtering operator to the equations. Since the filtering operator commutes with both spatial and time derivatives, we have

[TABLE]

where we have defined the SFS tensors related to the fluxes as

[TABLE]

In order to see the effect of this scale-separation in a simple application, let us consider the well-known Burgers equation, a longstanding field of research in LES [35]. In the one-dimensional case, given the velocity field $u$, the equation reads:

[TABLE]

The equation for the filtered fields can be obtained by applying the filtering operator to the equations, namely

[TABLE]

Notice that the filtered equation is equivalent to the original one except by a term proportional to $\overline{\tau}$ which represents the SFS residuals, that is, the loss of small-scale information due to the filtering process related to the non-linear terms. If in the region over which we average $u$ has a definite sign (which is usually the case except close to the shock), then $\overline{u}^{2}<\overline{u^{2}}$ and $\overline{\tau}$ is negative-definite, so that $-\overline{\tau}$ represents the kinetic energy hidden in the SFS.

Regardless of the specific system considered, these new terms always appear due to the non-commutativity of the filtering operator with the non-linear terms of the equations. In a LES, which effectively only evolves $\overline{u}$ and cannot simulate the unresolved scales, the explicit expression of $\overline{\tau}$ is a-priori unknown, and needs to be written as a function of the filtered fields in order to close the system. This implies in practice to approximate the SFS tensors with a sub-grid scale (SGS) model.

There are several ways of doing that, based mostly either on physical arguments or on an expected self-similarity of the solution. For instance, a popular, historical approach to fluid dynamics [25] is the setup of an artificial viscosity $\nu$, namely

[TABLE]

The viscosity parameter is often taken to be proportional to $\Delta^{2}|\partial_{x}{\overline{u}}|$, due to dimensional reasons [25] and in order to ensure that this term vanishes in the continuum limit $\Delta\rightarrow 0$, thus guaranteeing numerical convergence to the continuum (or DNS) solution. The numerical value of the proportionality coefficient can be fixed by hand or estimated by means of dynamical procedures which assume self-similarity (thus similar to a multiscale/multivariational approach) [36]. Regardless of the value of $\nu$, this functional form of $\tau$ is equivalent to the viscous Burger’s equation.

However, in many cases the physics involved does not consist only in dissipation, but can involve, for instance, inverse cascade and scale-dependent transfers of energy between the fluid kinetic energy and the magnetic one. Therefore, in those cases any result coming from the LES might be biased by the a-priori choice of a given SGS model. Hereafter, we will focus on the so-called gradient model, which has already been shown to capture the main features of the turbulent dynamics in our previous study for non-relativistic MHD turbulence box simulations [28].

III The SGS gradient model

Although many of the SFS terms are modeled based on some physical properties, it is possible to compute them relying only on mathematical arguments by considering the analytical Taylor expansion of the SFS terms[26, 27], which rely on the properties of the filtering operator $G$.

An homogeneous isotropic low-pass filter is independent on the direction (i.e., $G({\bf x}-{\bf x^{\prime}})=G(|{\bf x}-{\bf x^{\prime}}|)$) and only smooths out fluctuations on length scales smaller than the filter size, leaving unchanged the variations of the solution at larger length scales. In addition, the filter operator is linear and commutes with spatial derivatives. Generically, it can be written for any dimension $D$ as

[TABLE]

where $G_{i}(|x_{i}-x^{\prime}_{i}|)$ is just the one-dimensional kernel function.

The simplest low-pass filter is the mean value in a cubic domain with size $\Delta_{f}$ in each Cartesian direction $\{x_{i}\}$, described by the normalized kernel

[TABLE]

Despite the appealing simplicity of the box filter, which makes it very useful to perform numerical calculations, we will see below that it is not suitable for analytical calculations involving its derivatives, since they are not continuous. Therefore, at a formal level, it is more practical to introduce the normalized Gaussian kernel, which in the space domain can be written as

[TABLE]

where $\xi$ defines the effective filtering width. Besides having the same zeroth and first moments, Gaussian and box filters have the same second moment if we set $\xi=\Delta_{f}^{2}/24$. The filter is a useful mathematical tool to analyze not only the different scales on the solution, but also the structure of the evolution equations, allowing us to distinguish between the resolved or filtered fields and the unresolved ones on the SFS.

The main hypothesis is that the discretization of the equations over a finite-size grid is approximated by a filtering operator acting on the equations, with a Gaussian kernel equivalent to a box filter with width $\Delta_{f}$, eq. (11). Its $D$-dimensional Fourier transformed function is

[TABLE]

where $\bf k$ is the wavenumber. The main idea is to compute an approximation of the inverse filtering operator based on gradient expansion of the filter kernel $G$, that is, an approximation of the inverse Fourier transform of $1/{\hat{G}}$. The first step is to perform a Taylor expansion of the transformed function and its inverse in terms of the filter scale, that is,

[TABLE]

Considering the expansions to the fields ${\hat{f}}$ and ${\hat{\overline{f}}}$, the application of the inverse Fourier transformation yields to an infinite series representation of the filter operator and its inverse in terms of gradient operators acting on the fields, namely[37]

[TABLE]

These expressions are absolutely convergent and formally accurate at all orders, since the Gaussian kernel is infinitely differentiable and with unbound support. In fact, it was found that these series converges for all canonical filters.

The main idea is to filter the evolution system by means of a spatial filtering with Gaussian kernel (e.g. [38, 39]). Then, there appear unknown SFS tensors, which can be computed explicitly by using the Taylor expansion Eqn. (16). Therefore, we have obtained the following useful relations to first order in $\xi$:

[TABLE]

where “$\backsimeq$” means (here and thereafter) “accurate up to $\mathcal{O}(\xi^{2})$ terms” and $\nabla()\cdot\nabla()$ denotes contractions of spatial derivatives. Also, note that by $\frac{df}{d{\overline{g}}}$ we really mean $\frac{df}{dg}({\overline{g}})$.

Going back to our previous Burger’s equation example, it is easy to model the SFS terms appearing in Eq. (7) using the gradient expansion above, namely

[TABLE]

In principle, corrections to these expression are expected due to the form factor (due to the kernel shape) and the contribution from higher orders. In other words, this approximated expression can help as long as we can capture most of the dynamics with our LES, leaving to the SGS the task of mimicking the small-scale contributions. Notice that the SGS terms scale with $\xi\propto\Delta_{f}^{2}$, thus ensuring the numerical convergence (i.e., it vanishes in the continuous limit $\Delta_{f}\rightarrow 0$).

Regardless on these quite obvious caveats, notice that this prescription is very different from the viscous one given by eqn. (8): one arises from physical considerations, while the other just from mathematical ones. The first one is very useful in scenarios where the physics involved is well known and present universal behaviour, like in pure hydro-dynamical, non-relativistic cases. However, in the case of MHD, the non-linear interplay between kinetic and magnetic energy is rich of different physical mechanisms (dynamo effect, dissipation, helicity transfer…), at different scales. In these cases, it is more useful to have a SGS model which relies only on the mathematical features (basically, the gradients) of the involved field.

IV Gradient model for general system of conservation laws

The gradient model described before has been applied mainly to the non-relativistic HD or MHD equations, often in their incompressible version. In order to generalize it to the relativistic case, we tackle the problem from a broader perspective by considering a general system of conservation laws,

[TABLE]

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 transformations between conserved and primitive fields can be expressed formally as

[TABLE]

Notice that, although the relation $f^{a}$ between conserved and primitive fields is always known explicitly, the inverse function $g^{a}$ is not analytical in the relativistic case, where it usually needs to be solved numerically.

Conservative-scheme simulations are then equivalent to a filtered version of the evolution equations, such that the values of the filtered conserved variables $\overline{C}^{a}$ are numerically known. The filtering operation, implicitly contained in a LES, allows us to evaluate $\widetilde{P}^{a}:=g^{a}(\overline{C})$ but not $\overline{P}^{a}\equiv\overline{g^{a}(C)}$, which is what really appears in the filtered version of the equations. This distinction suggests a different definition of the SFS tensor in the filtered equations, which can be obtained again by applying the filter to eqs. (18), namely

[TABLE]

where hereafter the indexes “$a,b,e$” represent the different fields of the set of primitive and conserved quantities, while “$i,j,k$” represent the spatial directions. Here we have defined the SFS tensors related to the fluxes as

[TABLE]

that is, in terms of $\widetilde{P}^{a}$ instead of $\overline{P}^{a}$, which are the fields that we can compute from the evolved $\overline{C}^{a}$. Notice that these SFS tensors for the fluxes are not special, and similar ones can be calculated for any function. In particular, such residuals can be computed for the primitive fields $P^{a}$ (or for any other non-conserved field), namely,

[TABLE]

where the labels on $\overline{\tau}$ indicates to which field (and specific component of the field) the SFS residuals are referred.

In general, neither the SFS tensors associated to the primitive fields nor to the fluxes are known, and they can be important whenever a non-negligible part of the dynamics occurs in the small scales. This is the typical scenario in turbulence, where the fluctuations are relevant. The task is now to exploit the above gradient expansion of the filter kernel to provide a closed expression for $\overline{\tau}^{ka}_{F}$ (20) in terms only of the filtered variables and its derivatives. As we will see, this process essentially involves the filtering of composed functions; the Taylor expansion; and the inverse function theorem.

Let us start by expanding the expression

[TABLE]

where repeated index “$a,b,..$” denotes summation over the space of fields. Then, we can Taylor expand $F^{ka}(\overline{P})$ around $\widetilde{P}$, to first order in $\xi$, as

[TABLE]

and by the inverse function theorem we can express (locally) the Jacobian of the inverse variable transformation as the inverse of the Jacobian. That is,

[TABLE]

meaning the matrix inversion of the Jacobian $\frac{dC^{a}}{dP^{b}}$, evaluated at the filtered variables (namely, either $\overline{C}^{a}$ or $\widetilde{P}^{a}$).

Finally, by noticing that we can interchange $\overline{P}\rightleftarrows\widetilde{P}$ (at first order in $\xi$) in the last term of (22), and then combining it with (23) we finally get the functional form of the SGS tensor $\tau^{ka}_{F}$, namely

[TABLE]

which approximates the SFS residuals of the fluxes (i.e., $\tau^{ka}_{F}\backsimeq\overline{\tau}^{ka}_{F}$). In most of the cases it is useful to re-express the SGS tensor $\tau^{ka}$ above in the following equivalent form

[TABLE]

Due to its simplicity we will use this relation in the following, although in some specific evolution systems other equivalent forms can be preferred111 For instance, another interesting equivalent expression is given by $\tau^{ka}=-\xi\,\frac{d^{2}F^{ka}}{d\overline{C}^{b}d\overline{C}^{c}}\,\nabla\overline{C}^{b}\cdot\nabla\overline{C}^{c}$, which may be useful in situations where the fluxes can be expressed in terms of conserved quantities in a close form. .

V Applications to compressible ideal MHD

We will check the validity of our approach by applying first to the non-relativistic MHD equations, for which the LES equations with the gradient model has been studied for decades and were derived in the compressible case for a generic equation of state [28]. After that, we will consider the unexplored MHD relativistic case.

V.1 Non-relativistic case

The set of primitives fields for the non-relativistic MHD is given by $P^{a}=\left\{\rho,v^{i},\epsilon,B^{i}\right\}$ (density, velocity, specific internal energy and magnetic field, respectively). The conserved ones are $C^{a}=\left\{\rho,N^{i},U,B^{i}\right\}$ (density, momentum density, energy density and magnetic field), which can be written explicitly as a function of the primitive ones as:

[TABLE]

The evolution system is thus written as

[TABLE]

where the fluxes read

[TABLE]

with $2v^{[i}B^{k]}:=v^{i}B^{k}-B^{i}v^{k}$ denoting the usual anti-symmetrization operation. Notice that these evolution equations involve also the pressure $p(\rho,\epsilon)$, so an equation of state is required in order to close the system. The filtered version of the system, including the SFS terms in the right hand sides, is

[TABLE]

In this problem the relation between conserved and primitives is easily invertible, allowing to express all the fluxes as explicit functions of the conserved quantities. Thus, it is possible to compute $\frac{dF^{ka}}{d\overline{C}^{b}}\equiv\left\{\frac{dN^{k}}{d\overline{C}^{b}},\frac{dT^{ki}}{d\overline{C}^{b}},\frac{dS^{k}}{d\overline{C}^{b}},\frac{dM^{ki}}{d\overline{C}^{b}}\right\}$, and then use equation (26) to solve for the SGS tensors $\tau_{F}^{ka}$. The former can be written explicitly as:

[TABLE]

where $\widetilde{p}\equiv p(\overline{\rho},\widetilde{\epsilon})$, $\widetilde{S}^{k}\equiv S^{k}(\widetilde{P})=\widetilde{\Theta}\,\widetilde{v}^{k}-(\widetilde{v}\cdot\overline{B})\overline{B}^{k}$ and $\widetilde{\Theta}:=\overline{U}+\widetilde{p}+\frac{\overline{B}^{2}}{2}$. Also we will need the following derivative,

[TABLE]

Finally, we obtain

[TABLE]

where the SGS tensors have been splatted conveniently in order to connect with previous works, being:

[TABLE]

Notice that these results, obtained starting from a general formulation and applying it to the non-relativistic MHD case, agree with previous derivations calculated for this specific case [28]. The results here were derived for a general equation of state, which is reflected on the expression for $\tau_{\rm pres}$. In particular, the form of $\tau_{\rm pres}$ is greatly simplified for an ideal gas equation of state.

V.2 Relativistic case

Let us now consider the special relativistic MHD, for which a gradient SGS model has not been calculated so far. In this case, the set of primitives variables is given by $P^{a}=\left\{\rho,v^{i},\epsilon,B^{i}\right\}$, the conserved ones are $C^{a}=\left\{D,S^{i},U,B^{i}\right\}$, and the relations among them are given by

[TABLE]

where $W=(1-v^{2})^{-1/2}$ is the Lorentz factor, the pressure $p$ is defined through an equation of state and the enthalpy is defined by, $h:=\rho(1+\epsilon)+p$ . The evolution system is written,

[TABLE]

where

[TABLE]

Since the electric field can be obtained from the ideal MHD condition $E^{i}=-\epsilon^{ijk}\,v_{j}\,B_{k}$, all the previous fluxes can be easily written in terms of the primitive fields. Finally, we shall define $\mathcal{E}:=hW^{2}$ and $\Theta:=\mathcal{E}+B^{2}$, in order to simplify our following calculations.

The filtered version of the system can be written as:

[TABLE]

where on the right-hand side we have introduced the SFS tensors associated to each flux, as defined above.

As in the non-relativistic case, we propose to model the filtered SFS terms appearing in the equations (41)-(44) by means of a compact application of the gradient model. First of all, we define the double gradient operator $H$, acting on any given field $X$, as

[TABLE]

which satisfies a sort of generalized Leibniz’s rule,

[TABLE]

Notice that, when acting on any conserved field, it vanishes by its definition (i.e., $H_{D}=H_{S}^{i}=H_{U}=H_{B}^{i}=0$). On the other hand, when the operator applies to a non-conserved variable, the quantity is non-zero. This holds in particular for $\{\widetilde{p},\widetilde{\Theta},\widetilde{v}^{k}\}$ as we shall see below.

By applying these rules and using eq. (26), the SGS gradient tensors approximating the SFS terms of eqs. (41)-(44) read:

[TABLE]

where the set of the $H$ tensors, after some algebraic manipulations, can be written as:222We use a mixed notation for scalar products when it comes about the gradients, i.e. $\nabla X\cdot\nabla Y$ instead of $\nabla_{i}X\nabla^{i}Y$, in order to make well visible the gradient terms, which are the core of the SGS model and always appear contracted to each other.

[TABLE]

where $H_{E}^{i}$ is just the Hodge dual of $H^{ij}_{M}$, i.e. $H_{E}^{i}=\frac{1}{2}\epsilon^{i}_{\phantom{i}jk}H_{M}^{jk}$. Notice the values of the double gradient appearing above, $\{H_{p},H_{\Theta},H_{v}^{k},H_{E}^{k},H_{\mathcal{E}}\}$, are meant to approximate the SFS residuals related to the non-conserved fields, $\{\bar{\tau}_{p},\bar{\tau}_{\Theta},\bar{\tau}_{v}^{k},\bar{\tau}_{E}^{k},\bar{\tau}_{\mathcal{E}}\}$, defined according to equation (21) as $\bar{\tau}_{X}\backsimeq-\xi H_{X}$ (exactly like for the conserved field SFS residuals). Their explicit expressions are obtained by computing the following set of equations in the order in which they appear, where the quantities $\widetilde{\Psi}$ denote the auxiliary fields which are used to simplify the presentation (and also to facilitate their implementation):

[TABLE]

[TABLE]

Note that these equations are remarkably much more involved than in the well-known non-relativistic case [28]. The main difference being, of course, the more complicate relationship between the conserved and primitive variables in the relativistic setting. Notice that the continuity equation now acquires a SFS term, while the energy equation does not contain any SFS residual. Such “exchanged roles” of the continuity and energy equations comes from the non-relativistic limit and can be seen already at the level of the equations.

Notice also that, depending on the equation of state, some terms can be considerably simplified. For instance, $\widetilde{\Psi}_{A}=W^{2}\Gamma p\rho$ for an ideal gas EoS with coefficient $\Gamma$.

VI Application to magnetic relativistic simulations

One way to check the validity of the relativistic SGS gradient model is by performing a detailed analysis of numerical simulations displaying turbulence. Therefore, we will consider a decaying (i.e., non-forced) turbulent dynamics in relativistic MHD, triggered by the KHI, a scenario already studied especially in the non-relativistic case [40, 41]. Notice that the KHI provides different stages of the turbulent flow: the initial development at small scales, the transfer of kinetic to magnetic energy, the saturation and mixing, and the final slow decay [40, 28]. Moreover, the KHI is thought to take place in binary neutron star mergers, and will be the natural mechanism to be studied in the forthcoming general relativistic simulations.

VI.1 Setup

The initial setup is just an extension from the non-relativistic [41, 28] to the special relativistic case. We evolve eqs. (37)-(40) by using the SAMRAI infrastructure [42, 43] with a code generated by the platform Simflowny [44, 45]. The numerical schemes are the same used in Ref. [28], which were described in detailed in [46, 47]. The discretization of the continuum equations is performed by using the Method of Lines, which allows to address separately the time and the space discretization. We employ High-Resolution Shock-Capturing (HRSC) methods [48] 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-Friedrichs splitting [49] with the WENO5Z [50] high-order non-oscillatory reconstruction scheme. 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 $\Delta t\leq 0.4~{}\Delta$.

We set our problem in Cartesian coordinates, considering a periodic box $[-L/2,L/2]^{3}$. We shall consider different resolutions between $128^{3}$ and $1024^{3}$, and evolve the system up to $t=20$. The primitive fields read initially:

[TABLE]

where $a_{l}$ is the mixing layer scale, $y_{l}$ is the distance of the shear layers to the plane $y=0$, $\sigma_{y}$ and $\sigma_{z}$ are the extension scale of the initial perturbation in the $y$-direction and the profile of $v_{z}$, respectively. The main flow is initially given by $v_{x0}$. The standard values that we consider are $L=1$, $y_{l}=1/4$, $\rho_{0}=1.5$, $\rho_{1}=0.5$, $a_{l}=0.01$, $v_{x0}=0.5$, $B_{x0}=0.001$, $B_{y0}=B_{z0}=0$, $p_{0}=1$ and $\sigma_{z}^{2}=0.1$. We consider an ideal gas equation of state, $p=(\gamma-1)\rho\epsilon$, with $\gamma=4/3$. The purpose of this paper is not to explore the dynamics for different parameters (see for instance [40] for a discussion about the role of the initial values of Mach number and magnetic field).

The initial perturbation, $\delta v_{i}$, is a superposition of single-mode perturbation with a number of nodes $n_{i}\in[1,N/2]$, periodic in the boundary box:

[TABLE]

We underline that the specific form of the initial perturbation has no influence on the asymptotic turbulent behavior, as long as we excite the entire spectrum of modes, which can be achieved easily if the modes are not the same (or multiple) to each other, $n_{i}\gg 1$ and $\delta v_{i}\ll v_{0x}$. Hereafter, we use $n_{x}=11$, $n_{y}=7$, $n_{z}=5$, $\delta v_{x}=\delta v_{z}=0.01$, $\delta v_{y}=0.1$, $\sigma_{y}^{2}=0.01$. Since the KHI is known to grow faster for smaller scales, and due to the absence of physical viscosity in this test, we do not expect a numerical convergence, at least in the growth phase, since the more we refine the grid, the more fast-growing excited mode will be included. This actually reproduces the scenario of the binary neutron star mergers, where the finest available resolutions are likely still very far from being able to capture all the relevant modes, and the simulations do not show numerical convergence, in terms of total magnetic energy and its spectra. See also our previous work [28] for more details and for the general, similar behavior in the non-relativistic case.

VI.2 General behaviour

The following qualitative behavior is similar for all the resolutions considered. The development of the turbulent dynamics is illustrated in Fig. 1, where the density distribution over the slice $z=0$ is displayed at times $t=\{2,10,20\}$ for the highest resolution case with $N=1024^{3}$. At the beginning, the instability develops as small-scale structures with modes given by the initial perturbation. The development starts at the shear layer, where the transfer of kinetic to magnetic energy (dynamo effect) is fast up to $t\sim 2$ for all the resolutions (left panel), extending up to $t\sim 4$ for the highest one. Then, a larger scale mixing takes place (middle and right panel), during which the effectiveness of the dynamo mechanism reduces. After that period the mixing is completed and the fluid looks isotropic and homogeneous.

In Fig. 2 we compare the evolution of the integrated kinetic and magnetic energy for different resolutions. In our initially kinetic-dominated setup, both the internal and magnetic energies grows during the evolution, at the expense of the kinetic energy. For this particular setup, the internal energy is quantitatively the dominant one (due to the chosen values of the initial pressure), while the magnetic energy is always quite smaller than the kinetic energy.

As it has been observed already by several works, the transfer from kinetic to magnetic energy is more prominent at small scales. Therefore, since the initial perturbation spans all the scales and we have no viscosity included, when we consider higher resolution we capture more unstable modes, which in turn enhance further the magnetic growth. For the explored resolutions, we do not observe yet a saturation of the magnetic field in the homogeneous phase. Actually, in all cases, the dynamo mechanism slowly goes on, allowing in time to reduce the difference between the total kinetic and magnetic energy (possibly showing a slow approach to equipartition at the end of our $N=1024^{3}$ simulation). Similar behaviours are observed in the non-relativistic case for this problem. The saturation level will arguably be set by the strong feedback of large magnetic scales, but longer, and possible even more accurate simulations are needed to have a definitive answer.

VI.3 Spectra

In addition to volume-integrated quantities, whenever there is turbulence it is illustrative to compute also the radially-averaged spectrum [51, 28]. For a given field $f$ defined in a periodic box of side $L$, we use common python functions to calculate its discrete fast Fourier transform $\hat{f}(\vec{k})=\Sigma_{\vec{x}}f(\vec{x})e^{-i\vec{k}\cdot\vec{x}}$, where the sum is performed over the $N^{3}$ spatial points equally spaced in each direction, with $k_{j}=n~{}\Delta k$, where $\Delta k=\frac{2\pi}{L}$ and $n\in[0,N/2]$ is an integer.

We consider the radial coordinates of the Fourier space, describing it with $\Delta k$-wide radial bins also centered on $k_{r}=\{n~{}\Delta k\}$. Then, we calculate the spectra $\mathcal{E}(k)$ as averages $<\cdot>_{k_{r}}$ over the annular bins of the power density per unit of radial wavenumber in 3D:

[TABLE]

where $k^{2}=k_{x}^{2}+k_{y}^{2}+k_{z}^{2}$ and the normalization arises from the Parseval identity. For further technical considerations about normalizations, caveats (e.g. the systematic noise introduced by the conversion to radial coordinates in the Fourier space) and possible corrective factors, we refer to a recent dedicated paper [51]. The analysis of our simulations with a large number of points have been parallelized by using the python package [52] due to the large memory required.

The obtained spectra for $t=\{2,10,20\}$ are displayed in Fig. 3. First of all, note that all spectra show a change of slope at $k\sim$ few times $2\pi/\Delta$, due to numerical dissipation. This depends on the scheme, and it is a feature that already appeared in other relativistic hydrodynamical turbulence works employing finite-differences [53]. A possible cure could be to use spectral methods, which are suitable for bounding box simulations, but not for a complex astrophysical scenario. The change of slope means that the dynamics of the smallest resolvable scale is partially numerically damped. In the literature, this is actually a known issue that go under the name of implicit LES, meaning that the numerical dissipation effectively acts as an (uncontrollable) SGS Smagorinsky-like model within the simulation.

With this caveat in mind, we can analyze the spectra, focusing mostly on the large and intermediate scales. As time goes on, there is a direct cascade transferring kinetic energy from large to small scales (i.e., low to high wavenumbers). Therefore, the kinetic energy spectra extends to larger wavenumbers (i.e., smaller scales) with a Kolmogorov’s slope $\propto k^{-5/3}$, as the resolution is increased.

When one looks at the magnetic spectra, one can see that most energy tends to be stored at small scales, because these are where they are injected via dynamo process. The magnetic energy is spread also at larger scales (what is known as inverse cascade), so that we obtain a similar shape for the magnetic energy (roughly compatible with the analytical Kazantsev dependence $\propto k^{3/2}$ [54]). The more we rise the resolution, the more effective is this dynamo mechanism. Since the physical limit to the small scale is set by the viscosity (which in turns set the thickness of the shear layer), which is neglected here, we can understand the lack of numerical convergence, at least in the growth phase. This explains the evolution of the total magnetic energy seen above.

Ideally, the implementation of an effective SGS model should properly include the feedback of the small-scale dynamics on the magnetic energy distribution and its growth. That should provide a magnetic spectrum which, at the intermediate and large scales, should be similar to a case with a much higher resolution (see our non-relativistic results in [28]).

VII A-priori fitting

VII.1 Methodology

The most important analysis of our simulations is the a-priori test of the gradient SGS model. We run simulations with a certain grid step $\Delta=L/N$. Then, we consider a snapshot at a given time, and spatially filter all the evolved fields. The simplest recipe is to use a simple average groups of $S_{f}^{3}$ cells, where we define $S_{f}$ as the filter factor. This corresponds to apply a filter in the real space, with a box kernel of size $\Delta_{f}=S_{f}\Delta$, obtaining filtered fields evaluated over $N_{f}^{3}=(N/S_{f})^{3}$ points. The information lost in the filtering process is the field variation contained between the scales represented by $N_{f}$ and $N$. The solution in these SFS can be quantitatively evaluated by the explicit formal definitions of $\overline{\tau}$ as defined above in Eq. (20). Let us consider a simple case as an illustrative example, with a SFS tensor defined as $\overline{\tau}(\vec{x}_{f})=\overline{f}(\vec{x}_{f})\overline{g}(\vec{x}_{f})-\overline{fg}(\vec{x}_{f})$. We can evaluate it at each of the $N_{f}^{3}$ positions of the filtered mesh $\{\vec{x}_{f}\}$, as

[TABLE]

where $i$ indicate each of the $S_{f}^{3}$ discrete positions considered inside the cell centered in $\vec{x}_{f}$. Note that this estimation is not an exact evaluation of the loss information, since, by construction, it can only include the range of scales $[\Delta,S_{f}\Delta]$. The information for scales $<\Delta$ cannot be evaluated.

Once built each component of each SFS tensor, one can consider a given SGS model $\tau$. A measurement of the linear correlation between the numerical data and the different models can be estimated with the Pearson correlation coefficient between the SFS and SGS quantities,

[TABLE]

While the Pearson coefficient tests the functional form, one can also consider each SGS component with a pre-coefficient $C$ to be adjusted. Its best-fit value can be calculated by the minimizing the L2-norm, $\Sigma\,[\overline{\tau}(\vec{x}_{f})-C\tau(\vec{x}_{f})]^{2}$, where the sum is performed over all the positions $\{\vec{x}_{f}\}$. The minimization gives simply:

[TABLE]

This procedure can be repeated independently for each SGS component, for each tensor $\tau$. As we showed in our previous work [28], when one compares the performance of the gradient model with other SGS models available for the non-relativistic case, the Pearson correlation of the gradient model stands out, being always much closer to one than the others, and degrading to $\lesssim 0.5$ only for quite large filter sizes.

Due to these previous findings and the absence of consistent models to compare with in the relativistic case, we assess the gradient model by a-priori fits for different times, resolutions, and filter size. Below, we report our main findings, usually represented by averaging the different independent components of a given tensor. Differences between different components are statistically negligible, and can arise only temporarily during the initial phase of the development when fields are not homogeneous (see [28] for more details).

VII.2 Results

As a representative example, we consider the different resolutions filtered with $S_{f}=2$ (thus, averaging to $N_{f}=512^{3}$). In Fig. 4, we show the correlation Pearson values and the best-fit coefficients between the main SFS tensors, eqs. (41)-(44) and the corresponding SGS components tensors, eqs. (47), as a function of time. The Pearson coefficients are very close to one for all tensor components (${\cal P}\gtrsim 0.8$) at all times, indicating a very good correlation between the SFS and the SGS tensors even during the transitional development of the instability (where we have the minimum value of the correlation). Moreover, the best-fit coefficients are fairly constant in time and $\sim{\cal O}(1)$. This is the most important result, which confirms that the proposed model is actually fitting well for a variety of time (i.e., MHD configurations) and resolutions. As a comparison, other SGS models in the non-relativstic cases studied in [28] showed little correlation, or even not at all. In particular, the Smagorinsky model, the only one tested so far in GRMHD [24], showed a non-zero correlation (${\cal P}\lesssim 0.3$ at its best, for $S_{f}=2$, with a largely-varying best pre-coefficient) due to the ability to catch the transfer from resolved to SFS scales. Notice however that this is only a part of the non-linear dynamics involved in MHD turbulence.

The overall great performance of the gradient model are well known [55, 38, 39, 56, 57, 58, 28], due to the strong mathematical basis on which it relies. Let us stress that we have explored different initial conditions, finding mainly the same results. Another interesting result can be obtained by studying the simulation with $N=1024^{3}$ points with different filter size $\Delta_{f}$, as displayed in Fig. 5. As the filter size increases, there is more information lost by averaging in the cells. This leads to a degradation of the Pearson coefficient, although it is still above ${\cal P}\gtrsim 0.5$ even for $\Delta_{f}\geq 16$, indicating that the LES gradient terms are able to fit the functional form of the non-linear terms corresponding to a quite higher resolution. The variation over time of the best-fit coefficients are limited, which legitimates one to consider a constant pre-coefficient in the LES implementation, without any dynamical procedure (like in other non-relativistic models [36]) for its estimation.

One can also assess the correlations between the auxiliary double gradient terms given by eqs. (V.2)-(54) and the corresponding SFS terms. In Fig. 6, we show as an illustrative example the good agreement obtained by comparing the SGS term $(-\xi H_{v})$ with $(\widetilde{v}-\bar{v})$ (left panel), and $(-\xi H_{\Theta})$ with $(\widetilde{\Theta}-\overline{\Theta})$ (right). This is an important check, since the contribution to the overall SGS tensors of Fig. 4 could be given by a dominant component, thus “hiding” the others. Our results show that, taken one by one, each testable term (i.e., the $\propto H$ terms) correlates well with the corresponding SFS value.

VIII Conclusions

In this paper we have presented a formalism to study LES for a generic system of conservation equations, in which the evolved fields appear through non-linear (and possibly not analytical) relations in the fluxes. We have also extended the gradient SGS model for this generic system of conservation laws.

We have verified that, by applying this formalism to the non-relativistic MHD equations, one recovers known results available in the literature. Furthermore, we have applied the formalism to the relativistic MHD system, extending for the very first time the gradient SGS model to the relativistic case.

We have performed 3D numerical simulations of the KHI in a bounding box at different resolutions with high-order numerical schemes. Within these simulations we have been able to compare the residual SFS tensors with the gradient SGS tensors (i.e., a priori tests), showing a high correlation between these two quantities, as indicated by a Pearson number close to 1. Notably, the resulting best fit-parameters are also close to the expected value $C_{best}\approx 1$. Finally, we have seen that the Pearson value of these tensors is still significant (i.e., larger than 0.5) with a filter size up to $\Delta_{f}=16\Delta$, indicating that this approach can effectively account for a resolution between a factor of 4 and one order of magnitude times larger.

Reproducing small-scale dynamics by means of a validated LES can be seen as a computationally efficient way to solve the equations, equivalent to consider an effective higher resolution. When implemented in a LES, this approach can guarantee a considerable saving of computational time, by using a relatively low resolution, or if combined with high resolution, a feasible way to capture previously unaccessible small-scale dynamics. Quantifying the gain is not trivial, since it depends on the numerical scheme used (see [28] for considerations about the intrinsic dissipation in commonly-used finite-difference schemes). Our previous results in the non-relativistic case and the a-priori assessment here presented, allow us to estimate a gain in the effective resolution of at least a factor 4, and possibly up to a factor 8.

Finally, let us mention that although the LES and gradient SGS model presented here are valid only for the special relativistic case since we have used the Minkowski metric, the extension to General Relativity should be straightforward. As a matter of fact, although Einstein equations are highly non-linear, they do not show a “turbulent” regime: the metric tends to be smooth and slowly varying, compared to the matter fields. Therefore, the metric functions entering in the filtered MHD equations vary much more smoothly than the fields involved in the turbulent dynamics. Such extension, together with its application to binary neutron star simulations, will be presented in a forthcoming paper.