Extension of the sub-grid-scale gradient model for compressible magnetohydrodynamics turbulent instabilities

TL;DR

This paper extends the sub-grid-scale gradient model for compressible MHD turbulence simulations, improving the accuracy of energy evolution modeling and enabling better low-resolution dynamo effect predictions.

Contribution

The paper introduces an extension of the gradient model to include energy scales in compressible MHD, enhancing simulation accuracy in unresolved turbulent regimes.

Findings

Gradient model outperforms other models in correlation with high-resolution residuals.

Low-resolution runs with the gradient model accurately reproduce magnetic energy evolution.

Model extension is a step towards relativistic MHD applications.

Abstract

Performing accurate large eddy simulations in compressible, turbulent magnetohydrodynamics is more challenging than in non-magnetized fluids due to the complex interplay between kinetic, magnetic and internal energy at different scales. Here we extend the sub-grid-scale gradient model, so far used in the momentum and induction equations, to account also for the unresolved scales in the energy evolution equation of a compressible ideal MHD fluid with a generic equation of state. We assess the model by considering box simulations of the turbulence triggered across a shear layer by the Kelvin-Helmholtz instability, testing cases where the small-scale dynamics cannot be fully captured by the resolution considered, such that the efficiency of the simulated dynamo effect depends on the resolution employed. This lack of numerical convergence is actually a currently common issue in several astrophysical problems, where the integral and fastest-growing-instability scales are too far apart to be fully covered numerically. We perform a-priori and a-posteriori tests of the extended gradient model. In the former, we find that, for many different initial conditions and resolutions, the gradient model outperforms other commonly used models in terms of correlation with the residuals coming from the filtering of a high-resolution run. In the second test, we show how a low-resolution run with the gradient model is able to quantitatively reproduce the evolution of the magnetic energy (the integrated value and the spectral distribution) coming from higher-resolution runs. This extension is the first step towards the implementation in relativistic magnetohydrodynamics.

Tables

Table 1. Table 1: 3D KHI models: initial setup parameters (marked in italic the variations over the fiducial ones). See text for the other parameters, common to all simulations.

Model	$N^3$	$\gamma$	$v_{x0}$	$v_{z0}$	${\rho }_1$
KH3D250	${250}^3$	$4/3$	0.50	0	0.5
KH3D500	${500}^3$	$4/3$	0.50	0	0.5
KH3D500h	${500}^3$	$4/3$	0.50	1	0.5
KH3D500x	${500}^3$	$4/3$	0.25	1	0.5
KH3D500d	${500}^3$	$4/3$	0.50	0	1.0
KH3D500i	${500}^3$	1.0001	0.50	0	0.5
KH3D1000	${1000}^3$	$4/3$	0.50	0	0.5
KH3D1000h	${1000}^3$	$4/3$	0.50	1	0.5
KH3D2000	${2000}^3$	$4/3$	0.50	0	0.5

Table 2. Table 2: Table of 2D single layer models, with the theoretical and best-fit values ( 1 ​ σ 1 𝜎 1~{}\sigma error typically ≲ 2 % less-than-or-similar-to absent percent 2 \lesssim 2\% ) of the growth phase for the simulations run with WENO5Z.

model	N	$v_0$	$a_l$	$k_x$	$B_0$	${\alpha }_{\mathrm{an}}$	${\alpha }_{\mathrm{num}}$
grw1	50	0.645	0.05	1	0	1.73	1.53
grw2	100	0.645	0.05	1	0	1.73	1.70
grw3	200	0.645	0.05	1	0	1.73	1.75
grw4	400	0.645	0.05	1	0	1.73	1.75
grw5	200	0.645	0.025	1	0	2.4	2.37
grw6	200	0.645	0.1	1	0	0.66	0.68
grw7	200	0.3225	0.05	1	0	1.09	1.07
grw8	200	0.9215	0.05	1	0	1.77	1.79
grw9	200	0.645∘	0.05	2	0	1.36	1.32
grw10	200	0.645	0.05	1	0.129	1.69	1.70
grw11	200	0.645	0.05	1	0.258	1.56	1.54

Equations

$${\partial }_t\rho +{\partial }_k\left[\rho v^k\right]=0,$$

$${\partial }_t{S}^i+{\partial }_k\left[\rho v^k{v}^i-B^k{B}^i+{\delta }^{ki}\left(p+\frac{B^2}{2}\right)\right]=0,$$

$${\partial }_t{B}^i+{\partial }_k\left[v^k{B}^i-v^i{B}^k\right]=0,$$

$${\partial }_{t}U+{\partial }_k\left[\left(U+p+\frac{B^2}{2}\right)v^k-(v_j{B}^j)B^k\right]=0,$$

$$U=e+\frac{\rho v^2}{2}+\frac{B^2}{2},$$

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

$${\tilde{v}}^i:=\frac{{\overline{S}}^i}{\overline{\rho }}.$$

$${\partial }_t\overline{\rho }+{\partial }_k\left[\overline{\rho }{\tilde{v}}^k\right]=0.$$

$${\partial }_t(\overline{\rho }{\tilde{v}}^i)+{\partial }_k\left[\overline{\rho }{\tilde{v}}^k{\tilde{v}}^i-{\overline{B}}^k{\overline{B}}^i+{\delta }^{ki}\left(\tilde{p}+\frac{{\overline{B}}^2}{2}\right)\right]={\partial }_k{\overline{\tau }}_{\mathrm{mom}}^{ki},$$

$${\partial }_t{\overline{B}}^i+{\partial }_k\left[{\tilde{v}}^k{\overline{B}}^i-{\tilde{v}}^i{\overline{B}}^k\right]={\partial }_k{\overline{\tau }}_{\mathrm{ind}}^{ki},$$

$${\overline{\tau }}_{\mathrm{mom}}^{ki}:=$$

$${\overline{\tau }}_{\mathrm{kin}}^{ki}=$$

$${\overline{\tau }}_{\mathrm{mag}}^{ki}=$$

$${\overline{\tau }}_{\mathrm{ind}}^{ki}=$$

$${\overline{\tau }}_{\mathrm{pres}}=$$

$${\partial }_t\overline{U}+{\partial }_k\left[\tilde{\Theta }{\tilde{v}}^k-({\tilde{v}}_j{\overline{B}}^j){\overline{B}}^k\right]={\partial }_k[{\overline{\tau }}_{\mathrm{adv}}^k-{\overline{\tau }}_{\mathrm{hel}}^k],$$

$$\tilde{\Theta }=\overline{U}+\tilde{p}+\frac{{\overline{B}}^2}{2},$$

$${\overline{\tau }}_{\mathrm{adv}}^k$$

$${\overline{\tau }}_{\mathrm{hel}}^k=({\tilde{v}}_j{\overline{B}}^j){\overline{B}}^k-\overline{(v_j{B}^j)B^k}.$$

$${\tilde{p}}_{\mathrm{ideal}}=(\gamma -1)\left(\overline{U}-\frac{\overline{\rho }{\tilde{v}}^2}{2}-\frac{{\overline{B}}^2}{2}\right),$$

$${\tilde{\Theta }}_{\mathrm{ideal}}=\gamma \overline{U}+\frac{1-\gamma }{2}\overline{\rho }{\tilde{v}}^2+\frac{2-\gamma }{2}{\overline{B}}^2.$$

$${\overline{\tau }}_{\mathrm{pres},\mathrm{ideal}}=\left(\frac{1-\gamma }{2}{\overline{\tau }}_{\mathrm{kin}}^{lm}+\frac{2-\gamma }{2}{\overline{\tau }}_{\mathrm{mag}}^{lm}\right){\delta }_{lm}.$$

$$\overline{\tau }\approx C\tau .$$

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

$$\overline{fg}$$

$$\overline{vg}$$

$${\tau }_{\mathrm{kin}}^{ki}$$

$${\tau }_{\mathrm{mag}}^{ki}$$

$${\tau }_{\mathrm{ind}}^{ki}$$

$${\tau }_{\mathrm{adv}}^k$$

$${\tau }_{\mathrm{hel}}^k$$

$${\tau }_{\mathrm{pres},\mathrm{ideal}}^{ki}=\xi [(\gamma -1)\overline{\rho }{\partial }_j{\tilde{v}}^m{\partial }^j{\tilde{v}}_m+(\gamma -2){\partial }_j{\overline{B}}^m{\partial }^j{\overline{B}}_m].$$

$${\tau }_{\mathrm{kin}}^{ki}={\Delta }_f^2\overline{\rho }|\tilde{S}|{\tilde{S}}^{ki}$$

$${\tau }_{\mathrm{mag}}^{ki}={\Delta }_f^2|\overline{M}|{\overline{M}}^{ki}$$

$${\tau }_{\mathrm{ind}}^{ki}={\Delta }_f^2\frac{|\overline{J}|}{{\overline{\rho }}^{1/2}}$$

$${\tau }_{\mathrm{mom}}^{ki}={\Delta }_f^2|\overline{\rho }\tilde{S}:\overline{M}{|}^{1/2}{\tilde{S}}^{ki},$$

$${\tau }_{\mathrm{ind}}^{ki}={\Delta }_f^2\mathrm{sgn}(\overline{J}:\tilde{\Omega })\frac{|\overline{J}:\tilde{\Omega }{|}^{1/2}}{{\overline{\rho }}^{1/4}}{\overline{J}}^{ki},$$

$${\tilde{\Omega }}_{ki}:=\frac{1}{2}({\partial }_k{\tilde{v}}_i-{\partial }_i{\tilde{v}}_k),$$

Full text

I Introduction

In the context of computational fluid dynamics, direct numerical simulations (DNS) are employed to resolve the equations in all the dynamically relevant scales, including the dissipative ones. This means that the numerical grid size $\Delta$ has to be smaller than the dissipative scale. Since the computational cost of DNSs scales as the cube of the Reynolds number(Sagaut, Deck, and Terracol, 2006), in many cases it is unfeasible to have a small enough grid size and a large enough domain at the same time, thus representing all the scales of interest. This is the case for a variety of problems within different disciplines, ranging from industrial, laboratory-testable applications (such as gas turbines, steady turbulent flows in channels…), to astrophysical scenarios (such as convection zones of stars, planetary atmospheres, accretion disks, relativistic jets…). The increasing computational resources and the gradual sophistication of numerical techniques are constantly improving the achievable resolution and widening the number of scenarios where DNS will be feasible; however, currently there are many cases with very large Reynolds number, where simulations are still very far from capturing all the scales, and an alternative to the brute force of increasing resolution is then needed.

The two most widely used numerical approaches to bypass the resolution limitations in turbulent flows are the Reynolds Average Numerical Simulations (designed for steady-state flows and used especially in industrial applications) and the Large Eddy Simulations (LES). For both cases, the idea is to directly evolve the quantities on the resolved scales, i.e. from the integral scales (the largest ones, corresponding to dominion size or to the typical scale where energy is injected), down to the grid scales, and to describe the effect of the smaller ones through a suitable sub-grid-scale (SGS) model. LES have been developed and applied to many different computational fluid dynamics applications, with a plethora of variants, applicable to different numerical schemes. (seePiomelli (2014) for an historical and prospective point of view of LES development). While LES were originally designed for simple, steady flows, applications to decayingYang and Domaradzki (2004) transitional flows have become more popular (see recent reviewsZhiyin (2015) and references within). Due to the variety of turbulent scenarios, LES with SGS models represent an active field of inter-disciplinary research, including engineering, plasma physics and astrophysics. For the latter, which is also our ultimate aim of implementation, the use of LES is still not very commonMiesch et al. (2015); Schmidt (2015), compared to terrestrial scenarios.

Generally speaking, the purely hydrodynamic turbulent dynamics is relatively well understood due to the universality of the spectral properties, and an improvement of the proposed SGS models have allowed an impressive accuracy in reproducing experimental data or results coming from the DNSZhiyin (2015). On the other hand, the dynamics becomes more problem-dependent when magnetic fields are included, and prevent one from defining universal laws, valid for different scenarios. Magneto-hydrodynamic (MHD) turbulence consists of a non-linear interplay, with transfer from the kinetic to the magnetic (dynamo mechanism) and internal (dissipation) energies, and, at the same time, between different scales.Beresnyak and Lazarian (2015) In general, at the smallest scales the dynamo mechanism is more effective. A further non-trivial effect of the magnetic fields is the breaking of isotropy and self-similarity of the flow. Despite the huge theoretical and numerical progress in half a century, there are still many open questions. The main reason is that experimental data for high magnetic Reynolds number are limited, despite recent interesting resultsTzeferacos et al. (2018). A living reviewBeresnyak (2019) extensively summarizes the current state-of-the-art of MHD turbulence, focusing especially in astrophysical scenarios, with a focus on numerical methods, and linear and non-linear regimes of the dynamo effects.

In these scenarios, the smallest relevant dynamical scales cannot be captured by the available resolution, so that numerical simulations are still far from reaching the numerical convergence. In particular, LES have been very rarely applied to study the development of MHD instabilities. This is the case, for instance, of the Kelvin-Helmholtz instability (KHI), which growth rate is inversely proportional to the scale of the perturbation. The fastest-growing modes are characterized by a scale inversely proportional to the Reynolds number, so that if the latter is large enough, the KHI cannot be fully captured numerically. Among different astrophysical cases where KHI takes place, we mention the shear layer developed in a binary neutron star merger, source of the recently detected gravitational waves: the smallest relevant scales are likely a fraction of a meter, and the integral scales are thousands of kilometers, so that the state-of-the-art general relativistic MHD simulationsKiuchi et al. (2015, 2018) reach at its very best only 12.5 meters (at a price of about 60 million CPU-hours), thus not showing numerical convergence yet.

The arguably most used SGS model in LES is the Eddy-dissipative one, introduced by Smagorinsky half a century ago in the context of purely hydrodynamic turbulenceSmagorinsky (1963). It allows to mimic the fluid’s viscosity, and, in its magnetic extension Theobald, Fox, and Sofia (1994), the turbulent magnetic diffusivity. The multiple variants of the Eddy-dissipative model have allowed to set-up optimal coefficients in order to fit with an increasing accuracy the results coming from experiments or DNSs, especially for non-magnetic problemsSchumann (1975); Kraichnan (1976); Chollet and Lesieur (1981); Germano et al. (1991); Vreman, Geurts, and Kuerten (1994); Zhiyin (2015). Other SGS models have also been proposed relying mainly on dimensional analysis (cross-helicity model, for instance) or physical assumptions, like the dynamo mechanism. Comparisons between them have been studied in a variety of scenarios.Müller and Carati (2002a); Chernyshov, Karelsky, and Petrosyan (2007); Grete (2017)

In this paper, we focus on the gradient modelClark, Ferziger, and Reynolds (1979); Yeo (1987); Liu, Meneveau, and Katz (1994); Müller and Carati (2002a); Balarac et al. (2013); Grete (2017), based on a mathematical formulation which relates the finite resolution to an effective filtering of the equations. By expanding the filtered non-linear terms up to the first order of the Taylor series, one can approximate the residuals appearing in a given evolution equation with a functional form which includes the derivatives of the resolved fields, thus effectively extrapolating to the SGS the trends seen at the smallest resolved scales.

Recently, new advances in SGS modeling have been focused mostly in proposing advance methods for the dynamical evaluation (inspired by the pioneering work of GermanoGermano et al. (1991)) of the free coefficient in different SGS formulationsHendra and Bushe (2019); Zhou, Wang, and Jin (2018), including the gradient model for scenarios including, e.g., combustion flame Lar Kermani, Roohi, and Porté-Agel (2018), channel flow Ghaisas and Frankel (2016); Jiang et al. (2016, 2018); Dupuy, Toutant, and Bataille (2019); Lu et al. (2019), or a set of these standard computational fluid dynamics problems Vollant, Balarac, and Corre (2016).

In this paper, we propose an extension of the gradient model to include the residuals appearing in the energy evolution equation, for a generic equation of state, with a particular focus on the ideal gas.

In order to test the extended gradient model, we run box simulations of a non-forced magnetic KHI evolving according to the compressible MHD equations for an ideal gas. This problem, commonly used in several studies, allows us to test the different phases of a highly turbulent, magnetized flow: the development of the instability at small scales (the higher the resolution, the smaller such scales and the faster is the growth of the instability), and its decay after reaching isotropy and homogeneity. We test different initial conditions. We use high-order-accurate finite-difference methods, that are in general less suitable than spectral methods for turbulent flows (due to their higher dissipation), but are adequate in astrophysical simulations to deal with the development of strong shocks.

We will perform a-priori and a-posteriori tests. The former relies on the filtering of high-resolution runs, in order to fit the filtering residuals with an SGS model. We compare the performance of the gradient model with a few other models for the different non-linear residuals. The second kind of test is more challenging and consists in the effective implementation of the SGS model in the simulation, thus including its feedback in the evolution. One can compare the results of a low-resolution case with SGS model to a high-resolution one, without it.

Very detailed works about LES with SGS models have been performed in the context of either incompressibleMüller and Carati (2002a, b); Kessar, Balarac, and Plunian (2016); Vollant, Balarac, and Corre (2016) or compressible Grete et al. (2015, 2016); Vlaykov et al. (2016); Grete et al. (2017a); Grete (2017); Balarac et al. (2013) MHD, but not, to our knowledge, including an ideal gas equation of state. While the momentum and induction equations have been explored in many different works, the energy evolution has usually been neglected, or only partially considered, since they impose either incompressibility or a polytropic/isothermal equation of state Grete et al. (2017b).

In § II we summarize the MHD equations and their filtered versions implicitly evolved in a LES. In § III we introduce the extended gradient model, and a few more SGS models, for comparison. In § IV we describe the three-dimensional (3D) KHI simulations, together with some details about the platform and the numerical methods used. In § V we present the main highlights of the a-priori results obtained by filtering snapshots at different times of simulations with a variety of initial conditions, equations of state, and resolutions. In § VI we show the results coming from the implementation of the proposed model in low-resolution LES. Conclusions are drawn in § VII. The Appendices provide some details about the spectral and fitting calculations, and the numerical validation of our code with well-known, numerically-converging two-dimensional (2D) benchmark tests.

II LES and sub-filter-scale residuals in compressible ideal MHD

II.1 Evolution equations

The conservative formulation of the compressible ideal (non-relativistic) MHD equations consists of the continuity, momentum, induction and energy evolution equations, respectively:

[TABLE]

where $i$ and $k$ represent the spatial components, $\delta^{ki}$ is the standard Kronecker delta, $S^{i}=\rho v^{i}$ is the lineal momentum. The magnetized fluid is described by the mass density $\rho$, the fluid pressure $p$, the internal energy density $e$, the velocity field $v^{i}$, and the magnetic field $B^{i}$. The total energy density is defined as

[TABLE]

where for simplicity we have absorbed the $1/4\pi$ factor in $B^{2}$, and set the speed of light and the magnetic permeability $c=\mu_{0}=1$.

The closure of the set of equations is defined by an equation of state, which relates the pressure to the other physical variables, $p=p(\rho,e,B^{i})$. At each step of a simulation, the set of the discretized conserved quantities, ${\cal C}^{a}=\{\rho,S^{i},B^{i},U\}$, is numerically evolved, and the internal energy is recovered by inverting eq. (2). The physical variables are recovered accordingly: the velocity is simply given by dividing $S^{k}$ by $\rho$, while the fluid pressure is a equation-of-state-dependent function of the conserved variables, $p({\cal C}^{a})$.

II.2 Filtering and inversion of variables

In a LES, the effect of finite resolution can be thought as equivalent to a low-pass spatial filter applied to the evolution equations. The spatial scale cut-off, $\Delta_{f}$, is of the order of the grid cell, $\Delta$. Indicating the filter kernel with $G$, the filtering operator over a field $f$ can be written as

[TABLE]

An extended literature is dedicated to the mathematical representation of the filter.Yeo (1987); Pope (2000); Grete (2017) The functional form of the filter kernel is intrinsically given by the numerical method, and can be evaluated by specific studies of mode propagation, which can show both a damping and dispersive character for high wavenumbers $k\gtrsim 2\pi/\Delta$.Zhao, Lardjane, and Fedioun (2014) Generally speaking, finite-difference methods, here used, are much more dissipative (i.e., the effective filter is wider) than spectral methods, even when they are built for high-order accuracy, like the ones used here.

When we evolve any discretized set of equations with a conservative approach, we effectively solve for the filtered conserved fields, indicated hereafter by $\overline{{\cal C}}^{a}$. Any other primitive or auxiliary field $f$ (e.g., velocity, pressure or internal energy) can be recovered by expressing it as a function of the conserved fields, $f(\overline{{\cal C}}^{a})$. These implicit dependences will be indicated by $\widetilde{f}$. The most common example of this approach consists in expressing the velocities as a function of the filtered conserved density and momenta: the well-known mass-weighted or Favre-filtered velocities

[TABLE]

Hereafter we will write the relevant expressions as functions of $\widetilde{v}^{k}$, but all equations could be easily rewritten in terms of the conserved filtered momentum, $\overline{S}^{k}$, simply using eq. (4) and employing the identity $\overline{\rho v}^{k}=\overline{\rho}~{}\widetilde{v}^{k}$.

Besides the velocity, we will make use of the same notation for other two non-conserved fields: the pressure, $\widetilde{p}:=p(\overline{\rho},\overline{S}^{i},\overline{B}^{i},\overline{E})$, for which the functional form depends on the equation of state, and a quantity closely related to it.

II.3 Filtered equations

The filtered continuity equation is simply

[TABLE]

The filtered momentum and induction equations are well known to be

[TABLE]

where

[TABLE]

The sub-filter scale (SFS) terms account for the non-commutativity of the filtering operator with the non-linear fluxes. There are then three SFS tensors related to the products $v^{k}v^{i}$, $B^{k}B^{i}$ and $v^{k}B^{i}$ in the fluxes,111In literature, different sign conventions are used. Hereafter we will always define the SFS terms as the product of the filter minus the filter of the product. and a scalar related to the pressure:

[TABLE]

Notice that, by construction, $\overline{\tau}_{\rm kin}$ and $\overline{\tau}_{\rm mag}$ are symmetric tensors, while $\overline{\tau}_{\rm ind}$ is anti-symmetric and its inclusion is equivalent to having a SFS electric field with components given by $\epsilon_{ijk}\overline{\tau}_{\rm ind}^{ki}$.

By following a similar approach, we can write the energy evolution equation as

[TABLE]

where we have introduced the advected quantity

[TABLE]

and the additional vectorial SFS terms on the right-hand side of eq. (13) arise from the non-linear terms in the flux of the filtered equation. The first SFS term arises from the product of two non-conserved fields, $\widetilde{v}^{k}$ and $\widetilde{\Theta}$:

[TABLE]

and is related to $\overline{\tau}_{\rm pres}$ by construction, since $\widetilde{\Theta}$ contains the same factor $\widetilde{p}+\overline{B}^{2}/2$, while the additional term $\overline{U}$ does not imply any SFS residual since it is conserved. The second SFS term in the right hand side of eq. (13) is the cross-helicity term, which is partially related to the induction equation tensor, $\overline{\tau}_{\rm ind}$, but cannot be written only as a function of it (or a simple contraction with another filtered quantity):

[TABLE]

We stress that the functional forms of the terms related to the pressure ($\overline{\tau}_{\rm pres}$, $\overline{\tau}_{\rm adv}^{k}$) depend on the equation of state and, like $\tau_{\rm adv}^{k}$, had been neglected or only partially considered before. The pressure term $\overline{\tau}_{\rm pres}$ has been studied for the non-magnetized fluid in a few scenariosBorghesi and Bellan (2014, 2015); Taskinoglu and Bellan (2010); Taşkinoğlu and Bellan (2011), for instance with a Peng-Robinson equation of stateMa, Korucu, and Miller (2015), considering the performances of the gradient modelSelle et al. (2007). Here we consistently include all the terms at play.

To summarize, we need to give the prescriptions for 22 independent components of: three tensors (two symmetric and one anti-symmetric), eqs. (9)-(11), the scalar term $(\widetilde{p}-\overline{p})$ appearing in eqs. (12) and (15), and two vectors, eqs. (15)-(16).

II.4 A specific equation of state: ideal gas

So far, our formalism is general and can extend to any equation of state, with the only difference arising from the dependence $p({\cal C}^{a})$. In the next sections, we will present numerical implementation and results for an ideal gas equation of state, namely $p_{\rm ideal}=(\gamma-1)e$, where $\gamma>1$ is the ideal gas adiabatic index. In this case, considering eqs. (2) and (14), we have:

[TABLE]

One of the 22 SFS independent functional forms (see previous sub-section), the pressure term, is automatically set by the traces of the tensors $\overline{\tau}_{\rm kin}$ and $\overline{\tau}_{\rm mag}$, since it reads:

[TABLE]

In the pressure residuals, besides the well-known magnetic pressure contribution, there is an additional (negative) term related to the internal energy (i.e., to the fluid pressure $p$), proportional to $(1-\gamma)/2$ times the SFS kinetic and magnetic energy. In the limit of isothermal gas ($\gamma\rightarrow 1^{+}$), such extra contributions are nullified: in that case, or if fluid incompressibility is assumed, the fluid pressure does not depend on the filtered conserved fields, so that $\widetilde{p}=\overline{p}$, and the only surviving term is the magnetic pressure, in agreement with the commonly used formalismGrete et al. (2016).

III Sub-grid-scales modeling

The explicit expressions of the entire set of SFS terms above (which by definition cannot be computed or known a-priori) represent the closure of the system of equations in a LES. The closures, or SGS models, considered hereafter are local, i.e. depend only on the derivatives of the fields at the smallest resolved scales: we do not consider the multi-scale modelsSondak and Oberai (2012); Sondak et al. (2015); Sagaut, Deck, and Terracol (2006), where two different resolved scales are compared in order to extrapolate, by self-similarity, the SGS model at the unresolved scales. There are virtually infinite possible combinations of mutual contractions and functional forms involving the first-derivatives of the resolved fields, which can be considered.(Speziale, 1991; Kosović, 1997; Lu and Rutland, 2016) Far from being exhaustive, below we list the ones we explicitly use in this work, as a comparison with the extended gradient model. Note that the comparison is possible only for some of the SGS terms appearing in the momentum and induction equations: $\tau_{\rm kin}$, $\tau_{\rm mag}$, (or their combination $\tau_{\rm mom}$), and $\tau_{\rm ind}$. Since the terms $\tau_{\rm hel}$ and $\tau_{\rm adv}$ were never formalized before to our knowledge, there is no available SGS models alternative to the one proposed in this paper to compare with. The physical meaning of these new terms is less intuitive, so it is not clear how to model them on a physically-based grounds.

Hereafter we will indicate with an overline the SFS terms $\overline{\tau}$ resulting from either the formal analysis (like in the previous section) or the post-process filtering of a numerical solution (see § V), and with a simple $\tau$ the analytical SGS models, which are implicitly defined proportional to a constant parameter $C$, that are meant to approximate the SFS terms:

[TABLE]

III.1 Extended gradient model

The gradient (or non-linear) model Leonard (1975); Clark, Ferziger, and Reynolds (1979); Yeo (1987); Liu, Meneveau, and Katz (1994); Müller and Carati (2002a); Balarac et al. (2013); Vollant, Balarac, and Corre (2016); Grete (2017) arises from the analytical Taylor expansion of the SFS terms appearing in the MHD equations, under the hypothesis of having a filter with a Gaussian kernel. In the space domain this kernel can be written as

[TABLE]

where $\xi=\Delta_{f}^{2}/24$ provides a filtering width equivalent (up to the second moment) to a box filter with size $\Delta_{f}$. When we apply such filter to a product of fields, $fg$, the Taylor expansion of its Fourier transform and its inverse in terms of the filter scale allow to write the following rulesYeo (1987); Grete (2017):

[TABLE]

which are accurate up to ${\cal O}(\xi^{2})$. The second equation is obtained by recurrence relations and is applicable to the product involving the mass-weighted velocity $\widetilde{v}$, eq. (4). Analogous (although possibly more elaborated) relations can be obtained for other non-conserved fields.

Applying these simple rules to the non-linear products of the fields appearing in the SFS tensors, we can obtain the expression for each tensor. The well-known expansions of the terms (9), (10), (11), at leading order in $\xi$, read:

[TABLE]

These expressions have been commonly derived and used, both in a-priori and a-posteriori studies cited above.

Here we propose the general extension to the SFS terms related to the energy evolution equation, eqs. (15)-(16), for a general equation of state in compressible ideal MHD:

[TABLE]

where we have used eq. (23) and noted that $\xi\partial_{j}\overline{\Theta}=\xi\partial_{j}\widetilde{\Theta}+{\cal O}(\xi^{2})$.

Last, we have to define the functional form of the pressure term $\widetilde{p}-\overline{p}$, which enters in eqs. (12) and (27). In general, one has to consider the functional form $p({\cal C}^{a})$, and apply the rules above to the non-linear products of conserved variables contained in it.

For an ideal gas equation of state, $\tau_{\rm pres}$ is given by construction by the traces of the kinetic and magnetic SFS tensors, according to eq. (19), so that:

[TABLE]

III.2 Eddy-dissipative models

As a comparison, we consider a simple version of the Eddy-dissipative SGS model (also known as Smagorinsky model).Smagorinsky (1963) It accounts for the fluid dissipation at small scales, introducing an effective viscosity tensor proportional to the strain tensor $\widetilde{S}^{ki}$. A subsequent workTheobald, Fox, and Sofia (1994) extended the Eddy-dissipative model to the MHD, introducing the analogous turbulent magnetic diffusivity in the induction equation, proportional to the current density tensor. In this paper, we will test the following tensors:

[TABLE]

where we have introduced the standard resolved fluid strain tensor $\widetilde{S}_{kj}$, the resolved magnetic strain tensor $\overline{M}_{kj}$, and the resolved current tensor $\overline{J}_{kj}$. The kinetic and induction SGS models are the ones commonly used in the literature, with the addition of the factor $1/\overline{\rho}^{1/2}$ in the last term, introduced on dimensional basis; we add a further term to account for the magnetic non-linear term.

By construction, the contractions of eq. (30) with the fluid strain tensor and of eq. (32) with the current tensors provide definite positive values of large-to-small-scale transfer energy rate: this shows their dissipative nature. In its original version, $\tau_{\rm kin}^{ki}=(c_{s}\Delta_{f})^{2}\sqrt{2}|\widetilde{S}|\widetilde{S}^{ki}$, the pre-factor is usually found to be $c_{s}\sim 0.1-0.2$(Lu and Rutland, 2016) for incompressible steady flows. In more elaborated versions, $c_{s}$ can also be determined with a dynamical procedure, instead of being fixed, as proposed originally for a channel flow problem.Germano et al. (1991); Oberai and Wanderer (2005) Besides the normalization, different recipes are available to treat differently the diagonal and off-diagonal components of $\tau_{\rm mom}$ (for instance, to see how the energy transfer is affected Grete et al. (2017b)). A plethora of detailed studiesbar (1980); Rogallo and Moin (1984); Moin and Kim (1984); Germano et al. (1991); Moin et al. (1991); Lilly (1992); Brun, Friedrich, and da Silva (2006); Lévêque et al. (2007) are dedicated to fine-tune the dissipation rate for specific problems (especially hydro-dynamical, but see a recent workGrete et al. (2017a) and the references within).

Since the fine-tuning of this model is not this purpose and we aim at the comparison of its performances with the gradient model in a turbulent MHD instability case, we will consider a constant parameter $C$, as in eq. (20).

III.3 Cross-helicity model

This model(Müller and Carati, 2002a) is a variant of the Eddy-dissipative one, with a different normalization, which include products of derivatives of the velocity and magnetic field:

[TABLE]

where we have corrected dimensionally the original version (by introducing $\overline{\rho}$ factors), and we have introduced the resolved vorticity $\widetilde{\Omega}_{ki}$. As for the previous one, this model also transfers kinetic energy from resolved scales to SGS. On the other hand, in the induction equation the signum function allows for a back-scatter of magnetic helicity and magnetic energy (inverse cascade) if ${\rm sgn}(\overline{J}:\widetilde{\Omega})<0$.

III.4 Vorticity and Alfvén velocity models (for induction only)

Finally, we propose two new further models for the induction equation only. The first one is the vorticity model, defined by:

[TABLE]

Qualitatively, similar models (in their co-variant version) have been proposed in general relativity for neutron star mergers, as extra-terms in the induction equation. They are based on the vorticity, assumed to be a physical tracer of the dynamo mechanism (see its use in a relativistic MHD scenarioPalenzuela et al. (2015); Giacomazzo et al. (2015)), but they consist in an arbitrary functional form.

Finally, we test a new model based on the current, proportional to the resolved Alfvén velocity:

[TABLE]

Based on dimensional grounds only, this tensor is proportional to $\Delta_{f}$, instead of $\Delta_{f}^{2}$.

IV Simulations

IV.1 Platform and numerical schemes

The code presented here has been generated by using Simflowny Arbona et al. (2013, 2018) to run under the SAMRAI infrastructure,Hornung and Kohn (2002); Gunney and Anderson (2016) which provides efficient parallelization and adaptive mesh refinement (not used in this work). Simflowny is an open-source and user-friendly platform developed by the IAC3 group since 2008 to facilitate the use of HPC infrastructures to non-specialist scientists. It allows to easily implement scientific dynamical models, by means of a domain specific language, based on MathML and SimML, and a browser-based integrated development environment, which automatically generates efficient parallel code for simulation frameworks. Simflowny splits the physical models and problems from the numerical techniques. The automatic generation of the simulating code allows to properly include the parallelization features, which in this case rely on the SAMRAI infrastructure. The combination of these two platforms provides a final code with a good balance of speed, accuracy, scalability, ability to switch physical models (flexibility), and the capacity to run in different infrastructures (portability).

The MHD system of equations, given by eq. (1), can be written formally in conservation law form, namely

[TABLE]

where ${\cal C}^{a}=\{\rho,\rho~{}v^{i},B^{i},U,\Phi\}$ is the list of evolved conserved fields and ${\cal F}^{k}({\cal C}^{a})$ their corresponding fluxes which might be non-linear but depend only on the fields and not on their derivatives. Here $\Phi$ is a scalar field introduced to dynamically control the divergence of the magnetic field by means of a divergence-cleaning schemeDedner et al. (2002). For instance, in the simulations presented below, the dimensionless ratio between the L2-norm of $\Delta^{2}(\vec{\nabla}\cdot\vec{B})$ and the volume-integrated magnetic energy rise at the beginning of the simulation, reach a maximum of less than $10^{-4}$ and then stably decrease to $\sim 10^{-5}$, thus ensuring a negligible violation of the constrain. For further details about the divergence-cleaning formalism and numerical results, see our previous works Palenzuela et al. (2018); Viganò et al. (2019).

The time integration of the discretized 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$. For the spatial reconstruction, we employ High-Resolution Shock-Capturing (HRSC) methods(Toro, 1997) based on high-order accurate finite-difference schemes. They are commonly used in astrophysical scenarios and allow to deal with the 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 with high-order non-oscillatory reconstruction schemes.Shu (1998)

In Appendix B, we show in detail the validation of the code through benchmark tests in 2D, comparing different reconstruction schemes, of third and fifth order: FDOC3, WENO3YC, FDOC5, MP5, and WENO5Z. The latter is our final choice for all the results here presented. We defer the reader to previous worksPalenzuela et al. (2018); Viganò et al. (2019) for further details on these numerical schemes and an extensive analysis of the performance with different discretization schemes for different problems, including MHD. Here we only stress that these schemes, although very accurate, are more dissipative than spectral methods (more commonly used in turbulence, see below for the implications).

IV.2 Initial conditions

The KHI has been studied by a plethora of works with many different numerical schemes, also in the compressible MHD caseKeppens et al. (1999). We employ a setup similar to what found in literature, Beckwith and Stone (2011) consisting in a periodic cubic 3D box of domain $[-L/2,L/2]^{3}$ (described by Cartesian $\{x,y,z\}$ coordinates), with shear layers at $y=\pm y_{l}$ planes, and initial conditions given by:

[TABLE]

where $a_{l}$ is the shear layer thickness and $\sigma_{z}$ is the extension scale the profile of $v_{z}$, respectively. The pressure and the magnetic fields are assumed constant, $p=p_{0}$, $\vec{B}=\vec{B}_{0}$.

The main flow is then initially given by $v_{x0}$. For our initial conditions, the vorticity, $\vec{\omega}=\vec{\nabla}\times\vec{v}$, is concentrated around the shear layers $y\pm y_{l}$, with maximum values of the order $|\omega_{z}|\sim v_{x0}/a_{l}$ and $|\omega_{x}|\sim v_{z0}/\sigma_{z}$. The integrated values in the volume are $\int_{V}~{}\vec{v}~{}dV={\cal O}(\delta v)$ and $\int_{V}\vec{\omega}~{}dV={\cal O}(k~{}\delta v)$ due to the symmetry of the main flow with respect to each layer, where $k$ is the inverse of the typical scale of the perturbation. Therefore, $v_{z0}$ is used to primarily control the initial net kinetic helicity ${\cal H}_{k}=\int_{V}~{}\vec{\omega}\cdot\vec{v}~{}dV\sim 2~{}v_{z0}v_{x0}a_{l}$.

We consider fiducial values of $L=1$, $\rho_{0}=1$, $\rho_{1}=0.5$, $a_{l}=0.01$, $y_{l}=0.25$, $v_{x0}=0.5$, $\sigma_{z}^{2}=0.01$, $\vec{B}_{0}=0.001~{}\hat{x}$, $p_{0}=1$, assuming in general an ideal gas equation of state with $\gamma=4/3$. We consider models with number of points between $250^{3}$ and $2000^{3}$, and a possible variation of the fiducial parameters as described in Table 1. The setup is initially dominated by the large-scale kinetic energy, similarly to what happens in astrophysical scenarios, where the instability is actually the trigger for the growth and amplification of magnetic fields. The initial value of the internal energy depends on the choice of the pressure, and, in our case, we have initially $E_{\rm int}\gg E_{\rm kin}\gg E_{\rm mag}$. Note that in the model KH3D500i, we consider an isothermal fluid, approximated by $\gamma=1.0001$, in order to study the effect of the evolved pressure-related contributions.

We evolve the system up to $t=20$, except for the KH3D2000 model, which required over two million CPU-hours to reach $t=10$. Similar 2D setups and related results are described and analyzed in the Appendix B.2.

IV.3 The initial perturbation and its implications

The components of the three initial velocity perturbations appearing in eqs. (38)-(40), are periodic, single modes with a wavenumber related to an integer $n_{i}\in[1,N/2]$, and a different propagation direction:

[TABLE]

Hereafter, we use $n_{x}=11$, $n_{y}=7$, $n_{z}=5$, $\delta v_{x0}=\delta v_{z0}=0.01$, $\delta v_{y0}=0.1$, $\sigma_{y}^{2}=0.1$. It is important to note that this multi-modal perturbation in different directions, with different prime numbers for the wavenumber values, allows us to excite virtually all the modes spectrum. The non-linear interaction of the three initial different modes is effectively similar to imposing a random perturbation. The specific forms (values of $\delta v_{i0}$, $\sigma_{y}^{2}$ and $n_{i}$) of the initial perturbations have no influence on the asymptotic turbulent behavior, as long as the spectrum is entirely excited, which is guaranteed if the wavenumbers are prime number between them.

This choice gives a qualitatively difference compared to the same case with a partial excitation of the modes (e.g., $\delta v_{i0}\neq 0$ only for one component, or $n_{x}=n_{y}=n_{z}$): in that case, the turbulence will be restrained to a single mode and its harmonics, and, if the mode wavelength is much larger than the grid resolution, numerical convergence can be easily reached. In this sense, see the 2D, single-mode perturbed benchmark tests, taken from Obergaulinger, Aloy, and Müller (2010), successfully replicated in Appendix B.1, where one can also get an insight on the growth rates of the KHI.

The excitement of the entire spectrum, together with the ideal flow assumption (absence of a viscous scale), has also important implications for the general behaviour and the numerical convergence, as explained in the next section.

IV.4 General behaviour and magnetic energy growth

In Fig. 1 we show the density distribution over the slice $z=0$, at $t=1,5,9$, for the highest resolution case, KH3D2000. The following qualitative behaviour is the same in all cases. 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 maintained fast up to $t\sim 2$ in all models (left panel). Then, a larger-scale mixing takes place (middle and right panel), during which the magnetic growth rate reduces. At $t\gtrsim 10-12$, the fluid approaches the isotropy and homogeneity.

In Fig. 2 we compare the evolution of the integrated kinetic, magnetic and internal energies for the different models. Notice that the initial kinetic energy varies by a factor of a few depending on the explored initial conditions (different colors in the figures), while the magnetic energy is the same in all cases (we always consider the same initial magnetic field). The initial value of the internal energy depends the initial pressure (i.e., on $\gamma$), so that there are large differences between the isothermal KH3D500i and the other cases. For comparison purposed, we show, for each case, the evolution of the internal energy relatively to its initial value.

For all models, the KHI makes the magnetic energy grow very fast, by orders of magnitude, at the expense of the kinetic budget. Part of the kinetic energy is constantly dissipated into internal energy (right panel).

At the end of our simulations the system does not reach yet a quasi-stationary state, and the dynamo mechanism is still acting, though at a slower pace. The magnetic energy has not saturated yet and will be, in the best case, limited by equipartition with kinetic energy, but our highest-resolution case shows is still far from it. Longer, possibly higher resolution runs would be needed to know when the magnetic energy saturates (and thus approaching a numerical convergence of the asymptotic solution), but this is not the aim of this paper.

A crucial thing to note is that we do not reach numerical convergence even with the high resolution employed. This is partially expected, since we employ no physical viscosity or resistivity (formally infinite fluid and magnetic Reynolds numbers), and we excite the entire spectrum of the instable modes (see § IV.3). The more resolution we have, the higher is the wavenumber of the excited modes. Since the dynamo mechanism is more efficient at small scales, the net effect is a much faster rise of the magnetic field for higher resolution. In order to have numerical convergence, one should either (i) excite a single (or a fixed amount of) mode, which wavelength is larger than the resolution; (ii) employ a finite Reynolds number, in order to approach convergence when the grid size becomes comparable to the dissipation scales; or (iii) have a much higher resolution, which would allow an even faster rise of the magnetic field to its saturation value.

Note also that initially larger values of kinetic helicity (blue curves) provide similar growth rate, but a larger asymptotic values of the integrated magnetic energies. On the other hand, changing the initial density jump (green) or the equation of state (cyan) has an important effect on the rise of the internal energy, but not on the kinetic or magnetic energy evolution. This is compatible with the relatively small contribution of the pressure-related terms to the energy interchanges found in the literature Salvesen et al. (2014).

IV.5 Spectra

Fig. 3 shows the kinetic (solid) and magnetic (dotted) radial-averaged spectra (see Appendix A.1 for details about their practical calculation), for different models, at $t=5,10$ and 20 (first three panels). They confirm the considerations done in the previous subsection, and add some insights related to the relevant scales and the dissipation of the numerical methods used.

As turbulence develops, the magnetic and kinetic spectral energy densities become comparable below a certain scale. The scale of equipartition evolves to larger values (i.e., smaller $k$) in time. At large scales, instead, the kinetic energy remains always dominant.

Spectra show an inertial range (between the integral scale and the scale at which the numerical dissipation becomes important). In this range, the kinetic spectrum shows the typical Kolmogorov slopeKolmogorov (1941) (thin solid black line $\propto k^{-5/3}$). The magnetic spectrum, instead, peaks at small $k$, being possibly compatible with the Kazantsev slopeKazantsev (1968) (dotted thin line $\propto k^{3/2}$). These results are largely independent on the specific perturbation or details of the initial data, as it is shown in Fig. 4 in the energy spectra at $t=10$ for different initial conditions.

The numerical dissipation of the finite-difference scheme is well visible as a change of slope (knee) in all the spectra, located at a scale a few times larger than the grid size. As the resolution increases, the knee moves to the right, but the overall shape is the same. We see the same knee in 2D, for all models and for other fifth-order methods, while the third-order WENO3 fails even to accurately reproduce the Kolmogorov slope at any range (see details in Appendix B). The dissipation visible in the spectra is quantitatively similar to what obtained in other finite-difference codes employing splitting schemesRadice and Rezzolla (2012, 2013), and contrasts with the low-dissipative spectral schemes often used in turbulence box simulations.

Spectra give a further insight on the lack of numerical convergence already noted in the integrated quantities. The more we refine the grid, the more fast-growing excited modes (high $k$) will be included, and the more effective is the kinetic-to-magnetic energy transfer. Since we have no viscosity in our equations, there is no limit to the value of $k$ at which the instability develops. For a given time, as the resolution increases, the magnetic spectra is roughly rescaled to larger values for all $k$, and the equipartition is reached at larger scales (smaller $k$). As mentioned in the previous subsection, this is due to the faster growth of the KHI at smaller scales, at which the dynamo mechanism is also more effective. This can be seen clearly in Fig. 5, where we show the representative images for the $B_{z}$ component in the $z=0$ plane at $t=9$, for the same model with different resolutions, KH3D500, KH3D1000, KH3D2000. It is clear how stronger small-scale structures are created in the best-resolved cases.

V A-priori results

After the previous qualitatively description, let us focus on the a-priori fitting of SFS residuals with SGS models. We consider a simulation with a certain grid step $\Delta=L/N$, and its snapshot at a given time. We define the filter as a simple average over groups of $S_{f}^{3}$ cells, where we define $S_{f}$ as the filter factor and the filtered grid size, $\Delta_{f}=S_{f}\Delta$. The a-priori test consists in comparing each of the SFS residuals, $\tau$ (see § II) to a given SGS model, $\bar{\tau}$ (see § III). We evaluate the correlation between them (each consisting of $N_{f}^{3}$ values) by means of the Pearson correlation coefficient ${\cal P}$, and obtain the best-fit parameter $C_{\rm best}$ which enters in eq. (20). Further details of this standard analysis are available in Appendix A.2.

Below, we will extensively compare the evolution in time of ${\cal P}$ and $C_{\rm best}$ (i.e., analyzing and fitting different snapshots), for different SGS models and/or for different simulations of Table 1. Initially we calculate the values of ${\cal P}$ and $C_{\rm best}$ for each of the 21 independent components of the tensors ($\tau_{\rm kin},\tau_{\rm mag},\tau_{\rm ind},\tau_{\rm adv},\tau_{\rm hel}$). However, after having checked that there are no significant differences among the different components, we will show only their average over all the components (for instance, the values for $\tau_{\rm ind}$ are the average among its three components).

Ideally, ${\cal P}$ should be as close to one as possible, and should not exhibit important variations in time: the model should be able to fit all snapshots. The values of $C_{\rm best}$ should be $\sim{\cal O}(1)$ for the gradient model. A limited variation over time supports the idea that it is possible to have an analytical model without a dynamical procedure to adjust the pre-coefficient. Besides the correlation properties, we will further compare the distributions of $\tau$ vs. $\bar{\tau}$ by looking at the 2D histograms and computing different high-order moments.

At the end of the section we have summarized the outcomes from these comparisons.

V.1 Representative case: comparison between SGS models

As a representative example, we consider the model KH3D1000, filtered with $S_{f}=2$ (thus, averaging the $N=1000^{3}$ values to $N_{f}=500^{3}$). In Fig. 6, we show the evolution in time of ${\cal P}$, and $C_{\rm best}$ for the different SGS tensors.

As a first check, we compare in Fig. 6 the averaged values of ${\cal P}$ (top panels), and $C_{\rm best}$ (bottom panels) for the diagonal (first column) and off-diagonal (second column) components of $\tau_{\rm kin}$ for the gradient, Eddy-dissipative and cross-helicity models. Differences between the values related to the diagonal and off-diagonal components are visible only at the very beginning, when the best-fit parameters for the kinetic SGS are positive for the diagonal terms (indicating dissipation), while the shear terms (off-diagonal) are negative. In this early stage of localized turbulence development (see left panel of Fig. 2), we register the lowest values of $\cal P$ for all models. This is due to the high degree of anisotropy and discontinuity in the velocity and density field at the beginning. After this phase, correlation further improves and the best-fit parameters tend instead to be very similar, for all models and snapshots studied. Therefore, hereafter we focus only on the values of $\cal P$ and $C_{\rm best}$ averaged over all components of a given tensor.

Coming to the comparison between SGS models, we find that the gradient model performs better than the others, since the associated Pearson correlation coefficients are very close to one for all tensor components (${\cal P}\gtrsim 0.9$) at all times. Among the other models, the Eddy-dissipative and cross-helicity model for $\tau_{\rm kin}$ and $\tau_{\rm ind}$ have a lower, but still non-negligible Pearson correlation coefficient (${\cal P}\sim 0.1-0.3$). The same holds for the $\tau_{\rm ind}$ Alfvén model. On the other hand, the Eddy-dissipative model for $\tau_{\rm mag}$, and both the cross-helicity and vorticity models for $\tau_{\rm ind}$ show basically no correlation.

Furthermore, the best-fit parameters for the gradient model are fairly constant in time and $C_{\rm best}\sim{\cal O}(1)$ (see bottom panels of Fig. 6). In comparison, the $C_{\rm best}$** values for the other SGS models** can change by orders of magnitude depending on the time considered (Eddy-dissipative model for $\tau_{\rm kin}$ and $\tau_{\rm ind}$, cross-helicity model for $\tau_{\rm mom}$ and Alfvén model for $\tau_{\rm ind}$). The best-fit values $C_{\rm best}$ of the Eddy-dissipative and cross-helicity model are often below unity, in line with what found in literature. For models unable to fit SFS (very low ${\cal P}$), the best-fit parameters can even change sign (vorticity models, Eddy-dissipative model $\tau_{\rm mag}$, cross-helicity models for $\tau_{\rm ind}$), confirming their unreliability.

A closer look to this correlation can be seen, for a given time, by looking at the 2D distribution of the $N_{f}^{3}$ pairs of local values $\tau$-$\bar{\tau}$, for a given tensor component. In Fig. 7 we show a representative case, the model KH3D500 at $t=10$. These plots are very similar for all times, components and models, with only differences in the range of the value. We plot with colors the 2D SFS and SGS probability distribution for different terms, $\tau_{\rm kin}^{zx}$ (top panels) and $\tau_{\rm ind}^{xy}$ (bottom) calculated with gradient (left), Eddy-dissipative (center) and cross-helicity (right, for which we show $\tau_{\rm mom}^{zx}$ instead of $\tau_{\rm kin}^{zx}$, since the latter is not available) models. For these components, the distributions are quite symmetrical around zero, as expected for shear terms under isotropic conditions. A perfect correlation would be indicated by a single line with a finite slope corresponding to $1/C_{\rm best}$. Only the gradient model (left panels) is approaching such distribution, indicating its ability to systematically reproduce quite well the lost information on average. All the others (center and right) visually show a very poor correlation, as the values of ${\cal P}$ show.

Another way to compare the statistics is to look at the moments of the SFS and SGS distributions. In Fig. 8 we compare the evolution of the mean, variance, skew and kurtosis of the SFS terms and the related different models $C_{\rm best}\tau_{\rm kin}$, for the simulation KH3D500, with filter factor $S_{f}=2$. For the mean and the variance, the gradient model adheres to the SFS case, with a good performance also for the higher moments (within a factor $\lesssim 2$). On the other hand, the Eddy-dissipative and the cross-helicity models show smaller values for all quantities. This is consistent with the distribution comparisons in literatureKessar, Balarac, and Plunian (2016). The same relative behaviour is valid for other tensors, resolutions, initial conditions and filter size.

V.2 Performance of the new gradient terms

The overall performance of the gradient model in the a-priori tests, as shown so far, are well known (Grete et al., 2015, 2016, 2017a, 2017b; Grete, 2017; Kessar, Balarac, and Plunian, 2016). Besides confirming those results for our compressible ideal MHD problem, we here show how the SFS tensors $\bar{\tau}_{\rm hel}$ and $\bar{\tau}_{\rm adv}$ can also be fit very well by the corresponding gradient models, with the values of ${\cal P}$ and $C_{\rm best}$ in line with the other tensors. This is shown in Fig. 9, where we show them for the simulations with $500^{3}$ points, and $S_{f}=2$. Fig. 10 compares the moments of their SFS and SGS distributions, for the KH3D500 case, $S_{f}=2$. It qualitatively behaves like the tensors shown above, with an important mismatch (still within a factor $<2$) only for the kurtosis.

Finally, we have studied how the compressibility terms $\propto\partial_{i}\rho/\rho$ affects the correlations, for this specific problem. The well-known kinetic, magnetic and induction tensors show a negligible change in the Pearson if we include the term or not. The same holds for $\tau_{\rm hel}$. On the other hand, if we neglect such correction, $\tau_{\rm adv}$ loses all the correlation with the SFS term $\bar{\tau}_{\rm adv}$, showing that the compressibility term is actually dominant. Notice however that the importance of each term in a SGS tensor generally depends on the problem. For instance, if we consider a case with much larger changes in density (which, again, is the standard case for astrophysical bodies), then the compressibility correction is fundamental and cannot be neglected.

Below, we explore the results for different initial conditions, resolutions and filter sizes.

V.3 Dependence with initial conditions and dimensionality

By comparing all the simulations with the same resolutions but different initial conditions, we notice that the obtained values of ${\cal P}$ and $C_{\rm best}$ are almost coincident with the previous case. An example is given in Fig. 11, where we show the results of the gradient models for $\tau_{\rm kin}$ and $\tau_{\rm ind}$ respectively, with filter factors $S_{f}=2$ and $S_{f}=4$. Differences with initial conditions are really minor, compared to the ones obtained by changing $S_{f}$, as we will see below. The same holds for the other SGS models, and other resolutions.

As a further check, we also changed the dimensionality of the problem. We consider a 2D case, with the periodic box setup described in Appendix B.2, with $N=500^{2},1000^{2},2000^{2}$. We obtain very similar results regarding the Pearson correlation coefficients and the best-fit parameters, for all models and tensors. As a representative example, we show the results for $\tau_{kin}$, for the gradient model, in Fig. 12. Actually, the best-fit parameters are even more constant than in the 3D case, and the Pearson correlation coefficient is slightly better. The degrading trend of the Pearson correlation coefficient for increasing filter factors $S_{f}$ (see below) is instead the same. The lack of dependence of the results on the dimensionality is remarkable, given the qualitative differences in the kinetic cascade in the two cases.

V.4 Dependence with filter size and resolution

Changing the resolution (i.e., $\Delta$) or the filter size $\Delta_{f}$ affects the scale range of the lost information. Therefore, the values of ${\cal P}$ and $C_{\rm best}$ can change. For all cases, the gradient model results to be the best one. However, there are some differences in the quantitative performance indicators.

First of all, in Fig. 13 we show that, for a given $S_{f}$ (in this case 2 or 4), the dependency on the initial resolution is mild. If we maintain the same filter factor, the results are very similar. This reflects the fact that the relative amount of information lost, for instance, between $N=2000$ and $N_{f}=500$, is similar to the one lost between $N=500$ and $N_{f}=125$. On the other hand, by increasing the filter size, instead, the relative amount of "hidden" information increases. Thus, the fit gets worse for higher filter size, as shown with more detail in Fig. 14 for the model KH3D1000 and $S_{f}=2,4,10$. Increasing the size of the filter, the values of ${\cal P}$ for the gradient model decrease from $\sim 0.95$ to $\sim 0.6-0.7$ for $\tau_{\rm kin}$ (top panels), and from $\gtrsim 0.9$ to $\sim 0.3-0.5$ for $\tau_{\rm ind}$ (middle). In Fig. 15 we show the 2D correlations of the gradient model $\tau_{\rm ind}^{xy}$ for the KH3D1000 simulation, comparing different size filters. It is evident how increasing the filter (i.e., the amount of information loss), the correlation degrades.

The values of $C_{\rm best}$ for the gradient model increase by up of a factor of a few as the filter factor increases, deviating from the theoretical unitary value. This can be understood taking into account that the higher-order ${\cal O}(\xi^{2})$ terms will become more important for larger $S_{f}$. The same trends (degrading of $\cal P$ and mild increasing of the best-fit parameter, for increasing $S_{f}$) are seen for the Eddy-dissipative model (see $\tau_{\rm kin}$ in bottom panels of Fig. 11, and $\tau_{\rm ind}$ in Figs. 13 and 14), and for the cross-helicity model for $\tau_{\rm kin}$. For the Alfvén model, instead, the best-fit parameter scales as $\Delta_{f}$, as expected, since we defined it to be proportional only to $\Delta_{f}$, instead than $\Delta_{f}^{2}$. On the other hand, the vorticity model, the Eddy-dissipative model for $\tau_{\rm mag}$, and the cross-helicity model for $\tau_{\rm ind}$ do not show any apparent trend with filter size and resolution, probably because of the absence of a true minimum in the least square residuals function.

The study of the trends for varying filter factors allows to understand that even the best SGS models are able to fit the unseen tensors only until an effective improvement of resolution of a factor of a few.

V.5 Summary of the a-priori results

We summarize the findings obtained in this section as follows:

• •

The gradient model is able to fit very well (${\cal P}\sim 0.9$) the SFS residuals, for a limited filter size $S_{f}=2$.

• •

The best-value pre-coefficients show little variations in time, and $C_{\rm best}\sim 1$, as expected by theory.

• •

The gradient model is able to fit the new proposed terms, $\tau_{\rm hel}$ and $\tau_{\rm adv}$, besides the other traditional terms.

• •

The fitting performance of the gradient model degrades with the increasing size of the filter, but is independent on the original resolution: what matters is $S_{f}$. Reasonable correlations are seen also for $S_{f}=4,10$, indicating the potential ability to partially reproduce information hidden well deep inside the SGS (i.e., scales up to one order of magnitude smaller than the grid size).

• •

When several models are available for the same SGS term ($\tau_{\rm kin},\tau_{\rm mag},\tau_{\rm ind}$), the comparison of their performance show that the gradient model is always the best one for this particular problem, for all SGS terms.

• •

The out-performance of the gradient model is confirmed for all times considered and for all different initial conditions, resolutions and filter sizes.

• •

The comparison of the higher-order statistics (variance, skew, kurtosis) and the 2D-correlation plots confirm that the gradient model is able to fit the SFS residuals with no important systematic biases.

VI A-posteriori test with LES

While the a-priori test is informative about how well the functional form of a SGS tensor fits SFS residuals, the non-linear dynamics and their feedback over time can be tested only by means of direct implementation of such tensors in a LES. We implement the SGS gradient tensors, eqs. (24)-(29), as an additional part of the fluxes in the filtered MHD equations, as in eqs. (6), (7), (13). Each SGS tensor includes the same free parameter $C$, eq. (20).

We consider the same 3D KHI setup and numerical methods as before, and use the fiducial models KH3D250 and KH3D500 from table§ 1. Then we compare the results with the previous runs with no SGS included (§ IV.4), at different resolutions. Fig. 16 compares the magnetic energy evolution obtained without or with SGS models (setting $C=\{1,4,16\}$), for different resolutions. The first evident effect is that, by increasing the value of $C$, the magnetic energy attained increases as well. Setting $C=1$ seems to have a negligible effect, which start to be visible for $C=4$. The magnetic energy value obtained by the gradient model for KH3D250 with $C=16$ is as high as the one obtained with KH3D500 without any SGS model. Similarly applying the SGS gradient model with $C=16$ to KH3D500 gives a growth of the magnetic field slightly smaller than KH3D1000.

More detailed information can be obtained by analyzing the spectra for the same cases, displayed in Fig. 17 at $t=\{10,20\}$. While the kinetic spectra is barely affected due to the dominance of the large scales, the magnetic spectra show some interesting features. At both times we can see that, including a SGS gradient model with $C=16$, we reach the same values of the spectra in the large scales (inertial range) as the highest-resolution cases without SGS models. Due to the intrinsic dissipation, the knee is still located at the same wavenumber, so that spectra unavoidably decays for high $k$. Since the magnetic spectra peak at small scales, a low resolution by definition cannot take into account an important part of the total magnetic energy. This is why the gain in total magnetic energy is limited, as it was shown in Fig. 16.

However, it is crucial to underline that the feedback on the large scales is properly captured by the SGS gradient model up to resolutions: the magnetic energy growth and spectra in the inertial range are comparable with what obtained using a resolution $\sim 4-8$ times higher.

Although these final results are encouraging, there is a mismatch between the a-priori values of $C_{\rm best}\sim 1-2$ and the ones that we need to implement in the LES (i.e., $C\gtrsim 10$) to have an important effect. A possible explanation might be that the Taylor expansion is not accurate enough and additional terms need to be included. In order to explore such possibility, we implemented the leading next-order terms in the $\tau_{\rm kin}$ and $\tau_{\rm ind}$ gradient models (following the already developed formalismVlaykov et al. (2016)), setting them proportional to $C^{2}\xi^{2}$. We do not see any relevant change in the overall evolution (integrated energies and spectra), compared to the standard first-order gradient models considered above (again, with different values of $C=1,4,16$). This excludes the second-order terms are the reason of the mismatch.

The high values of $C$ can be instead understood considering the scheme’s intrinsic dissipation. The gradient model is known generally to over predict back-scatter effects,Leonard (1975) leading to an over-estimation of the transfer from SGS to resolved scales. This is favored especially for spectral methods, where the very low numerical dissipation shows accumulation of energy at the smallest resolved scales. This is why many works try to improve the numerical stability considering a mixed model given by a superposition of both gradient and Eddy-dissipative models Müller and Carati (2002a); Kessar, Balarac, and Plunian (2016) . However, in our case the intrinsic dissipation of the scheme (visible from the steep spectral slope at high $k$) prevents this instability, unless we use high values of the parameter $C\gtrsim 20$.

This is probably the same reason why we need high values of $C$ to see some effects in the LES, much higher than the a-priori best-fit ones. As a matter of fact, the intrinsic dissipation partially damps the small-scales structures, i.e., the local gradients. Since the gradient model relies on them, we need a higher value of $C$ to compensate the damping. Note that this effect in the a-priori test is not visible, since we filter the simulation, reaching a scale $S_{f}\Delta$, where the intrinsic numerical dissipation is not that significant. Said in other words, if one looks at the spectra, the steep decay is especially severe at the smallest scales $\sim\Delta$, so that applying the SGS model there (like in the a-posteriori LES tests) require a higher factor than if we apply it at $S_{f}\Delta$ (like in the a-priori test). We then expect that the $C$-mismatch should decrease by using a less dissipative numerical scheme, probably at the cost of a less stable solution.

VII Conclusions

This work generalizes the SGS gradient model for the compressible ideal MHD with a generic equation of state $p(\rho,e)$, although we focus here on the specific case of the ideal gas. We have extended the formalism of filtering, including also the energy evolution equation.

We have considered a turbulent scenario with a very large Reynolds number (formally, infinite), which is common in astrophysics. In those cases, the finest available numerical resolutions are often very far from being able to capture all the relevant small scales and the integral ones. With this in mind, we have considered box simulations of magnetic KHI, with the initial excitation of perturbations which extend to the smallest possible wavelengths. If the modes are well below the best achievable resolution, the simulations do not show numerical convergence, in terms of total magnetic energy and its spectra. The absence of physical viscosity and/or resistivity in our runs, and the specific problem, imply not having a reference solution to compare with. This actually represents the main motivation for the use of SGS in astrophysical LES: to explore unaccessible scales at an affordable computational cost.

We have generated our computational code by using the platform SimflownyArbona et al. (2013, 2018), endowed with HRSC high-order finite-difference schemes commonly used in numerical relativity applications, already described in depth in our previous workPalenzuela et al. (2018). We stress that we use methods that may not be optimal for periodic box turbulence, but are adequate in astrophysical simulations with the development of strong shocks. As a matter of fact, many works in turbulence use spectral methods, which are particularly suitable for studies of smooth solutions (or with weak shocks) in periodic boxes due to their intrinsic low numerical dissipation. However, in a real problem, where the dominion is not periodical and the dynamics includes strong shocks, codes employing finite-difference or finite-volume schemes are preferred due to their robustness. With this in mind, we have considered box simulations of magnetic KHI, with the initial excitation of perturbations which extend to the smallest possible wavelengths.

Note that the gradient model does not assume or mimic any physical behavior like dissipation. Instead, it extrapolates the non-linearity of the flux terms to the SGSs. Therefore, it can be applied in principle to any set of equations. This is a key issue for our specific problem, for which the efficiency of the dynamo effect (conversion of kinetic into magnetic energy) strongly depends on the smallest resolved scales, where the non-linear coupling between velocity and magnetic fields needs to be captured.

We have assessed in detail the gradient model, and compared it with a few more SGS models. Our results agree with what found in literature: the gradient model outperforms all the others in the a-priori tests (fitting SFS residuals), with best-fit parameters that do not deviate from the expected unity more than a factor of 2. We obtained our results for a remarkable variety of: dimensionality, initial conditions, resolution, time within a simulation and filter factor. The Pearson correlation coefficient is excellent (${\cal P}\gtrsim 0.9$) for a filter factor $S_{f}=2$, and gradually degrades, being still decent (${\cal P}\gtrsim 0.5$) for a filter factor 10. The degradation of the performance is understandable: if the LES resolution is too coarse compared to the dynamically interesting scales, then the inclusion of the SGS model will not be able to represent all the missing scales.

Then we include the gradient model in the LES to perform a posteriori tests, less frequent in the literature, and compare it to the integrated energies and spectra coming from the runs with a higher resolution. The results can be summarized as follows: (i) we obtain a more effective dynamo mechanism, with a growth by up to a factor of a few of the magnetic energy, (ii) the gain of the magnetic energy comes from the large scales, where by adding the gradient with a large value of the parameter (i.e. $C=16$) we can reproduce a very similar magnetic spectrum to the one given by a much higher resolution. By definition, and due to the intrinsic dissipation, the magnetic energy stored at the small scales is not completely provided, but its feedback on the large scales (i.e., the overall large-scale dynamo) is effectively taken into account.

The use of such a large value of the parameter needed to reproduce the high-resolution results suffers a mismatch with the a-priori test. We think that this is due to the intrinsic dissipation of the finite-difference code, and less dissipative methods should need a smaller value of $C$. The intrinsic dissipation of our scheme also explains why our simulations are numerically stable for any value $\lesssim 20$, while spectral codes typically need to introduce an artificial viscosity (or use a mixed SGS model) for $C\sim 1$.

Although the precise value of $C$ to be used can be improved, here we stress how the implementation of the extended SGS gradient proposed manages to obtain our main aim for the considered problem: the growth to intermediate and large-scale magnetic energy, which cannot otherwise be captured. The ability of the gradient model to capture the SGS dynamics is limited to a factor of a few times the resolution used. Therefore, it can be seen as a computational-cost-effective (the extra cost is to add terms to the flux) way to capture the small-scale dynamics, equivalent to employing an effective higher resolution. This is especially important in order to improve the simulations, and possibly get a little bit closer to the numerical convergence, in the astrophysical turbulent scenarios mentioned above.

The main intrinsic limitation of the SGS gradient model arises from the mathematical assumption of having smooth fields: in the presence of strongly varying fields, the truncation to first-order in $\xi$ of the Taylor expansion becomes a poor approximation. Although this effect is not observed in our simulations, in other scenarios with strong shocks/discontinuities the gradient SGS model might loses part of its accuracy in describing the unresolved scales. A possible solution might be to include higher-order terms in the Taylor series. Nevertheless, we stress again that the LES combined with the gradient SGS model allows to include only part of the effects coming from the unseen dynamics, effectively digging information inside the numerical cell, but only to a certain effective depth (which is difficult to evaluate, since the performance of the gradient degrades gradually with the size filter).

This paper is the first step towards the generalization of SGS-gradient models in LES to the framework of General Relativity. An astrophysical case where it will be applied is the binary neutron star merger, which requires a dominion of thousands of km, and, at the same time, should ideally simulate the KHI happening in the shear layer forming during the collision, probably below the meter.

Acknowledgements.

We acknowledge support from the Spanish Ministry of Economy, Industry and Competitiveness grants AYA2016-80289-P and AYA2017-82089-ERC (AEI/FEDER, UE). CP also acknowledges support from the Spanish Ministry of Education and Science through a Ramon y Cajal grant. The computational time in the Barcelona Supercomputing Center’s MareNostrum has been granted by the $17^{th}$ PRACE regular call (project Tier-0 GEEFBNSM, P.I. CP), and by the RES calls AECT-2018-1-0005 (P.I. DV), AECT-2019-1-0007 (P.I. DV). We thank B. Miñano for the technical support, S. Liebling and F. Carrasco for useful comments.

Appendix A Evaluation of spectra and fitting

A.1 Spectra

For a given field $f$ defined in a periodic box $[0,L]^{3}$, 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 $x_{i}\in[0:L)$ 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 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}$. The normalization comes from the Parseval identity, so that $\int_{0}^{N/2}\mathcal{E}(k_{r})dk_{r}=\int_{V}f^{2}(\vec{x})dV$ where $f^{2}=\rho v^{2}/2,B^{2}/2$ are the local energy densities in the real space. Note that this identity, expressed in terms of the integral in $k_{r}$, strictly holds if the fluid is isotropic, so that the average over the annular bin does not lose statistically relevant information (for instance, a strong dependence on power on the direction of $\vec{k}$). 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 paperDurran, Weyn, and Menchaca (2017). Note that the calculation of spectra for simulations with a large number of points ($N\gtrsim 1000^{3}$) require a large computational memory. We use a parallelization by parallel slabs, with a dedicated python package, recently providedMortensen (2016). This calculation is different and more precise, compared with simplified methods used in other works, e.g. the spatial average over one dimension of the spectra obtained from the Fourier transform over the other two dimensions.Beckwith and Stone (2011)

A.2 A-priori fitting

For a simulation with a certain grid step $\Delta=L/N$, 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

[TABLE]

obtaining filtered fields evaluated over $N_{f}^{3}=(N/S_{f})^{3}$ points. When one considers a non-linear combination of conserved fields, the filtering process allows one to evaluate numerically the residuals contained between the scales represented by $N_{f}$ and $N$, which can be seen as an approximation of the formal definitions of $\overline{\tau}$, eqs. (9)-(12) and (15)-(16). For instance, if a SFS tensor is formally 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 by the Pearson correlation coefficient,

[TABLE]

Note that, due to the usually very large of degrees of freedom ($N_{f}^{3}$), the associated $p$-value is almost always very small, and will not be taken into account as a useful indicator.

While the Pearson correlation coefficient tests the functional form, one can also consider each SGS component with a parameter $C$ to be adjusted. Its best-fit value can be calculated by the minimizing the L2-norm, $\Sigma\,[\overline{\tau}^{ki}(\vec{x}_{f})-C\tau^{ki}(\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$.

Appendix B 2D tests

B.1 Open boundaries Kelvin Helmholtz test

As a validation test for the methods employed, we run a well-known 2D KHI test set-upObergaulinger, Aloy, and Müller (2010), with open boundaries in the vertical direction. We set a rectangular domain in the $x-y$ plane, with side length $L_{x}=1$, $L_{y}=2$, centered in $(0,0)$, with a number of points $(N,2N)$. The initial conditions consist of $P=\rho=1$ everywhere, with an ideal equation of state with $\gamma=5/3$, with a $a_{l}$-thick shear layer along the $x$ direction, with the following initial conditions:

[TABLE]

where $v_{0}$ is the shear velocity, and $\delta v_{y}\ll v_{0}$ represents the initial perturbation. We impose periodic boundaries in the $x$-direction, and open conditions (i.e., copying the values of the last cell to the ghost cell) in the $y$-direction (see a previous detailed studyObergaulinger, Aloy, and Müller (2010) for the effects of changing the boundaries in the $y$-direction from open to reflecting, but this is not the focus of our work).

We set $\delta v_{y}=10^{-6}$ and we considered different cases, as listed in Table 2, similar to Table A.1 in the reference paper.Obergaulinger, Aloy, and Müller (2010) We use different reconstruction methods: WENO3YC, FDOC3, FDOC5, MP5, WENO5Z (for details about the parameters of these models and further tests, see our dedicated paperPalenzuela et al. (2018)). We follow the simulations well beyond the saturation time, with a time-step $\Delta t=0.25~{}\Delta$, where $\Delta=L_{x}/N$ is the homogeneous resolution of the numerical grid.

As a first example, in Fig. 18 we show the plots of the density and velocity fields (white arrows), at $t=10$ and $t=15$ (left/right panels, respectively), for the reference model grw3 (top panels) and its magnetic version, grw10 (middle panels). For the latter, we also show in the bottom panel the plots of vorticity (with colors) and magnetic field (black arrows). After an initial phase, the instability is triggered and develops vortices, until a quasi-stationary configuration is reached.

A quantitative tracer of the degree of the turbulence developed at a given time during the development stage is the $y$-component of the kinetic energy integrated in the volume:

[TABLE]

After an initial time, $E_{ky}$ tends to grow exponentially until it saturates. We fit the exponential growth rate phase with a function $E_{fit}=E_{0}e^{2\alpha_{\rm kh}t}$, and compare the best-fit values of $\alpha_{\rm kh}$ with the numerical values obtained by the numerical simulations available in literatureObergaulinger, Aloy, and Müller (2010), and the ones theoretically obtained by means of a theoretical stability analysisMiura and Pritchett (1982).

The growth rate and the saturation value of $E_{ky}$ are physically controlled by different values of $v_{0}$, $a_{l}$, $k_{x}$ and $B_{0}$, but not by the amplitude of the perturbation $\delta v_{y}$, which only influence the time at which the instability starts to develop. Furthermore, the numerical methods and resolution are crucial to determine the growth rate and the instability trigger time.

In Fig. 19 we show the $E_{ky}$ growth for two of the eleven cases analyzed and listed in Table 2: grw1 (low resolution) and grw10 (magnetic case). Inspecting the first, low-resolution case, grw1, we underline that the third-order methods, FDOC3 and WENO3YC, are not able to capture the instability in the low resolution case, while the fifth-order methods do, even though the growth rate is different. In particular, WENO5Z gives the higher value, slightly larger than MP5 and and FDOC5, and slightly smaller than the theoretical reported valuesObergaulinger, Aloy, and Müller (2010) (see last two columns of Table 2). Such differences disappear already with intermediate resolutions, $N=200$ (model grw3, right panel). In this case, the third-order method still show a value of $\alpha_{\rm kh}$ slightly inferior than the theoretical value, which is instead perfectly matched by the simulations with the three fifth-order methods. The three of them are basically indistinguishable.

We then test all the cases of the benchmark paperObergaulinger, Aloy, and Müller (2010) with the method WENO5Z, which was the one with the best behavior. For all of them, we find an excellent agreement with the theoretical expectationsMiura and Pritchett (1982), as shown in Table 2. One on side, the trends between the growth rate and the physical parameters $v_{0}$ (compare grw3, grw7, grw8), $a_{l}$ (grw3, grw5, grw6), $B_{0}$ (grw3, grw10, grw11) and $k_{x}$ (grw3, grw9) are correctly recovered. On the other hand, the differences among the models grw1, grw2, grw3 and grw4 show how the resolution allows one to better capture the instability. For this particular case, having $N\gtrsim 200$ seems to be enough to capture the instability. However, we underline that the importance of the numerical resolution is amplified if the initial perturbation mode ($k_{x}$) is larger: the smaller the scale, the more important the resolution and methods become.

As a final preliminary test with this setup, we examine the possible issue considered in literatureHopkins (2015): in presence of a uniform, Galilean boost in one direction (i.e., an additional advection velocity, constant everywhere), the development of the KHI can be numerically inhibited, due to the numerical dissipation which is enhanced by the value of the advecting velocity. Thus, we compare the same model (grw3 in this case), with and without an additional advection in the direction $x$ with a value comparable to the difference of velocities, $v_{x,\rm adv}=1$, such that, initially, eq. (49) becomes $v_{x}=1+\frac{v_{0}}{2}\tanh{\frac{y}{a_{l}}}$. Fig. 21 shows the growth of the kinetic energy with and without the advection, with WENO5Z, for the problems grw1, grw3, grw10. No important differences are noticed. We repeated the test with for the cases with lower-order methods (WENO3YC, FDOC3) and/or the low-resolution problems, and in any case we found an appreciable difference.

This first test shows that the numerical implementation is correct and able to reproduce the expected growth rate of a single-mode KHI, as long as at least one among the resolution and the accuracy order is high enough.

B.2 Periodic 2D Kelvin Helmholtz tests

We also tested the KHI in a periodic 2D box with a double mixing layer, similar to the 3D case presented in the main text. This has allowed us to compare the available results for 2D turbulent periodic boxes(Beckwith and Stone, 2011; Stone et al., 2008; Lecoanet et al., 2016; McNally, Lyra, and Passy, 2012).222See also Athena webpage:

http://www.astro.princeton.edu/j̃stone/Athena/tests/ We consider a square domain with ${x,y}$ both in the range $[-L/2,L/2]$ with $N$ points in each direction. The initial positions of the mixing layers are $y=\pm L/4$, as the initial set-up describes:

[TABLE]

where the sign $\pm$ indicate the region $|y|>L/4$ and $|y|\leq L/4$, respectively. The perturbation in each $i$ direction is given either by an uniform random distribution within $[-\delta v_{0},\delta v_{0}]$, or by $\delta v_{i}=\delta v_{0}\cos(2\pi xn_{i}/L)$, where $n_{i}\in[0,N/2]$ are integers and they can be different in the two directions.

Hereafter we set $L=1$, $v_{0}=1$, $\rho=1$, $\delta\rho=0.5$, $p_{0}=2.5$, and we use an ideal equation of state with $\gamma=1.4$, with a sinusoidal perturbation with $\delta v_{0}=0.01$, $n_{x}=4$ and $n_{y}=7$, unless specified otherwise. We impose periodic boundary conditions in all directions and we run with a time-step $\Delta t=0.25~{}\Delta$, in order to be sure that the numerical errors are not dominated by the time discretization. We run up to $t=10$, a time long enough to reach the quasi-stationary state, and we run the same simulation with different resolutions, $N\in[100^{2},2000^{2}]$ with WENO5Z.

We have run different values of the initial magnetic field. In general, the initially constant magnetic field partially inhibits the KHI. This feature is well known and is quantified by the Alfvén factor, defined as the ratio between hydro and Alfvén velocities, $A=v_{0}/\sqrt{B^{2}/\rho}$. When $A\lesssim 1$, the KHI is totally inhibited, as explained in detailed in previous works.Obergaulinger, Aloy, and Müller (2010) In our case, we start from a value $A\sim 10^{3}$, thus the KHI takes place.

Fig. 22 shows the evolution for a $N=1000^{2}$ simulation with WENO5Z, including magnetic field with $B_{0}=0.001$. We show the density map (color scales) and the velocity and magnetic fields with arrows, at different times.

The stages of the dynamics are similar to the 3D case: at the beginning, vortical structures are formed with a pattern given by the $k_{x}$ wave-numbers of the initial perturbation ($t\lesssim 1$ for the chosen values of $\delta v_{y}$). The symmetry of the patterns is broken when the structures at the two mixing layers start to interact ($t\gtrsim 1$). Then, the turbulence, initially concentrated around the mixing layer, tends to homogeneously spread over the entire dominion, thanks to the periodic boundary conditions. Meanwhile, the vortexes tend to merge (the well-known inverse kinetic cascade in 2D), and at late time ($t\gtrsim 10$ in this case) only two large vortexes of opposite vorticity sign survive. The turbulent dynamics reaches a saturated, fully developed state, dominated by these large vortexes, and, given the absence of any forcing, it decays, transforming kinetic energy into internal energy.

For our setup, when the system reaches a quasi-stationary state, the magnetic energy remains still significantly smaller then the kinetic energy, and is significantly stored in small scales. Kinetic and magnetic spectra show the same characteristics as the 3D case (main text), with the expected 2D Kolmorogov and Kazantsev slopes in the inertial range.

The timescales and the resolved structures mainly depend on the problem parameters ($B_{0}$, $\delta v_{y}$), and on the resolution, respectively. At this respect, Fig. 23 shows the plot of density with colors, and the velocity and magnetic fields at $t=10$, for decreasing resolution $N=1000^{2},400^{2},100^{2}$. It is evident how the the resolution affects the capability to capture the turbulent details.

We have also performed the a-priori fitting for different resolutions and filter sizes. We report a representative case in Fig. 12, briefly mentioned in the main text. Overall, for all tensors and SGS modles results are really similar to the 3D ones.

We checked that the growth of $E_{ky}$, indicating the development of the turbulent state, is enhanced if high resolution or high-order methods are used, and MP5 and WENO5Z provide statistically compatible results, if the same resolution is used. The difference between two different numerical schemes are only stochastic: the particular position of the vortexes can be different, but the spectra and the integrated quantities are statistically indistinguishable. In the main text, all simulations are run with WENO5Z.

The same holds if we compare the fully developed turbulence obtained with different initial perturbations: the specific time of growth can change, but at saturation the runs are statistically equivalent. In this sense, we checked that the values of $\delta v_{0}$ or the specific form of the perturbation can affect the delay of the beginning of the exponential growth, but not the growth rate and the asymptotic value of $E_{ky}$. This holds always as long as the initial values of the perturbation is random, or it is sinusoidal but at least one among $n_{x}$ and $n_{y}$ of the initial perturbation has no common integer factors with $N$. When the three numbers have a common integer factor $C$, then the initial data are effectively periodic in the $x$-direction on a scale $L/C$. Therefore, in these cases, a mode with wavenumber $C$ and its higher harmonics are maintained during the entire simulation, since the $x$-scales larger than $L/C$ are not explored, thus preventing a proper spectral cascade (the spectrum is constituted by spike corresponding to the chosen modes and its harmonics) and the related full development of the turbulence. For instance, we numerically proved that, starting with $n_{x}=n_{y}=4$ or $n_{x}=n_{y}=8$, and $N=500$ (for which $C=4$), the simulations show a periodicity in the $x$-direction given by $L/4$, being 4 the common factor between $N$, $n_{x}$ and $n_{y}$. This is a further test for the correct numerical implementation, that shows no sign of numerical break of symmetry.

Summarizing, the 2D tests has allowed us to validate the capability of the code to simulate turbulent dynamics, and to identify WENO5Z as the best numerical method among the tested ones. That is what we use in the main text.