Nonlinear closures for scale separation in supersonic magnetohydrodynamic turbulence

Turbulence in compressible plasma plays a key role in many areas of astrophysics and engineering. The extreme plasma parameters in these environments, e.g. high Reynolds numbers, supersonic and super-Alfvenic flows, however, make direct numerical simulations computationally intractable even for the simplest treatment -- magnetohydrodynamics (MHD). To overcome this problem one can use subgrid-scale (SGS) closures -- models for the influence of unresolved, subgrid-scales on the resolved ones. In this work we propose and validate a set of constant coefficient closures for the resolved, compressible, ideal MHD equations. The subgrid-scale energies are modeled by Smagorinsky-like equilibrium closures. The turbulent stresses and the electromotive force (EMF) are described by expressions that are nonlinear in terms of large scale velocity and magnetic field gradients. To verify the closures we conduct a priori tests over 137 simulation snapshots from two different codes with varying ratios of thermal to magnetic pressure ($\beta_\mathrm{p} = 0.25, 1, 2.5, 5, 25$) and sonic Mach numbers ($M_s = 2, 2.5, 4$). Furthermore, we make a comparison to traditional, phenomenological eddy-viscosity and $\alpha-\beta-\gamma$ closures. We find only mediocre performance of the kinetic eddy-viscosity and $\alpha-\beta-\gamma$ closures, and that the magnetic eddy-viscosity closure is poorly correlated with the simulation data. Moreover, three of five coefficients of the traditional closures exhibit a significant spread in values. In contrast, our new closures demonstrate consistently high correlation and constant coefficient values over time and and over the wide range of parameters tested. Important aspects in compressible MHD turbulence such as the bi-directional energy cascade, turbulent magnetic pressure and proper alignment of the EMF are well described by our new closures.


Introduction
Turbulence is ubiquitous in astrophysical plasmas, ranging from coronal mass ejections and stellar winds [1], through star formation in molecular clouds [2], to the gas in the interstellar [3] and intracluster medium. While experimental setups [4] become increasingly more realistic, they are still far away from the regime acting in such extreme conditions. For numerical simulations it is computationally too expensive (if even possible) to capture the entire range of physical processes from plasma kinetics to the integral scales of turbulence. In an astrophysical context, one has to further contend with the additional complications brought about by high compressibility and the accompanying supersonic and super-Alfvenic motion.
Possible ways to circumvent the infeasibility of direct numerical simulations are the use of calculations based on mean-field theories or large-eddy simulations [5]. These simulations only resolve the energy containing large scale dynamics and require a subgrid-scale (SGS) model to account for unresolved effects. While a lot of research has been successfully carried out in the realm of hydrodynamics [6], compressible magnetohydrodynamic SGS closures are essentially unexplored. Previous research is mainly based on the concept of turbulent dissipation in incompressible flows [7,8,9,10]. They expand the idea of a turbulent eddy-viscosity to an additional eddy-resistivity in the induction equation and propose different phenomenological models. Even though these models are then evaluated a posteriori , a general verification and justification a priori has so far only been considered for a single incompressible dataset [11]. Thus, our objective is to establish the validity of closures for the filtered, compressible MHD equations by coarse-graining multiple datasets from high-resolution simulations of statistically homogeneous, forced MHD turbulence.
In general, the effect of finite resolution in numerical simulations can be mimicked by applying a low-pass filter to the standard, ideal MHD equations. This is achieved by convolving the equations with a suitable filter kernel G. See e.g. Garnier et al [6] for details on the properties of low-pass filtering and the conditions that G needs to satisfy. For a homogeneous, isotropic, stationary kernel, under periodic boundary conditions [12] the equations take the following form ∂ρ ∂t The units of the magnetic field B incorporate 1/ √ 4π. An overbar denotes § filtered and a tilde mass-weighted filtered quantities [12]. For instance, the filtered density § Throughout the paper the symbol is used as a generic placeholder for variables. designates closure expressions. Furthermore, we employ Einstein summation convention and i,k is identified with the k-th partial derivative of the i-th component of . field is given by ρ = G * ρ, while the mass-weighed filtered velocity field is u = uρ/ρ. In this formalism, all filtered primary quantities, density ρ, velocity u, magnetic field B, and thermal pressure P are presumed to be known and directly accessible. Due to the introduction of mass-weighted filtering the only remaining terms that require closure are the SGS stress τ and the electromotive force (EMF), E. They are analytically expressed [10] as The SGS stress tensor can be decomposed into the well-known turbulent Reynolds stress τ u , a turbulent Maxwell stress τ b and a magnetic pressure term. Furthermore, the definitions of the SGS energies are obtained from applying the filter to the total filtered energy density E, which can be decomposed into the contribution due to resolved fields only and a remainder, designated as SGS energy .
It is important to point out that in general the filtering operator is not a Reynolds operator, in particular = . It follows that SGS terms, like E SGS , carry information not only about the interactions between unresolved fields but also about cross-scale interactions between unresolved and resolved fields. In addition to this, the turbulent magnetic pressure is identical to the magnetic SGS energy E b sgs and both kinetic and magnetic SGS energies are directly given by 2E u sgs = Tr (τ u ) and 2E b sgs = Tr τ b , i.e. they constitute the isotropic parts of the respective SGS tensors. Following the general tensor decomposition, the deviatoric, traceless parts are then given by τ * ij = τ ij − 1 3 δ ij τ kk .

Traditional closures
In hydrodynamics, the traceless part of the SGS stress tensor is commonly closed by means of the eddy-viscosity hypothesis τ u * ij = −2ν u ρ S * ij in analogy to the molecular viscosity term in the momentum equation, where S ij ≡ 1 2 ( u i,j + u j,i ) is the filtered kinetic rate-of-strain tensor. This introduces a turbulent kinetic eddy-viscosity ν u = C u ν ∆ E u sgs /ρ 1/2 which is proportional to a characteristic velocity, commonly given by the kinetic SGS energy, and a characteristic length scale ∆. This closure has already been applied directly to MHD [7,13] by neglecting the magnetic contribution τ b * ij in the momentum equation. A turbulent magnetic viscosity ν b = C b ν ∆ E b sgs 1/2 was used in [10] with the closure The SGS energies can either be determined by individual evolution equations, where several terms again require closure, or by an instantaneous closure. Smagorinsky [14] introduced such an instantaneous closure in pure incompressible hydrodynamics (B = 0) by assuming the SGS energy flux to be in equilibrium with the rate of dissipation Here, | S * | ≡ 2 S * ij S * ij denotes the rate-of-strain magnitude. Finally, the EMF is commonly modeled (e.g. [7,8,13,10]) by variations of [15] with resolved current J = ∇ × B and vorticity Ω = ∇ × u. The coefficients α, β, and γ are typically related to the α-effect, turbulent resistivity and turbulent cross helicity, respectively. The commonly used closures for these coefficients,

Nonlinear closures
In our new approach we adopt the compressible hydrodynamic nonlinear closure for the kinematic deviatoric stress tensor τ u * from [16], similar to the incompressible one from [17]. We propose the straightforward extension to MHD with The tensorial structure, e.g. u i,k u j,k , can be obtained by a Taylor expansion discarding terms with 2 nd and higher order gradients of the resolved fields. The overall normalisation with the subgrid-scale energies comes from the constraint that the SGS stresses vanish in laminar flows.
Applying the nonlinearity idea to the EMF generalizes the closure proposed by [11] to the compressible regime The closure explicitly preserves the anti-symmetry between velocity and magnetic field in E, which in turn helps in capturing their relative geometry. Finally, to complete the set of nonlinear closure equations, we use the Smagorinsky expression for the turbulent kinetic energy (3) and propose an analogous extension to the magnetic part Here, the turbulent magnetic energy is proportional to the magnetic rate-of-strain magnitude |M| ≡ 2M ij M ij . There are two advantages of closing E u sgs and E b sgs separately and not jointly via the total SGS energy. First, there is no additional need to close the often neglected turbulent magnetic pressure, as it is given by E b sgs . Second, the individual energies provide closures to the isotropic parts of the turbulent stress tensors τ u and τ b .

Validation Method
In order to evaluate the proposed closures, we perform an a priori comparison using simulation data obtained from two, grid-based MHD codes (Enzo [18] and FLASHv4 [19]). This way the results are less likely to hinge on the particulars of the numerical implementation. In both cases we follow the evolution of a compressible, isothermal fluid in a cubic box with resolution of 512 3 grid cells and periodic boundary conditions, starting from uniform initial conditions. In Enzo we use an ideal equation of state with adiabatic exponent κ = 1.001 in order to approximate isothermal gas. Enzo is a finite-volume code, i.e. the evolution equations are evaluated in integral form by solving a Riemann problem for the mass, momentum and energy flux through cell walls. This allows for the conservation of MHD invariants (e.g. energy) to machine precision. We use a MUSCL-Hancock scheme [20] (a second-order accurate Godunov extension) with second-order Runge-Kutta time integration and Harten-Lax-van Leer (HLL) Riemann solver (a two-wave, three-state solver) to solve the ideal MHD equations. The FLASHv4 code is similar to Enzo (second-order accurate in space and time), but uses the positive-definite HLL3R Riemann solver [20]. Another difference is that FLASHv4 uses a polytropic equation of state to keep the gas exactly isothermal. Moreover, explicit kinematic viscosity and magnetic resistivity terms are included in the momentum, energy, and induction equations. We set the kinematic and magnetic Reynolds numbers to Re = Rm = 3780. Consequently, the magnetic Prandtl number Pm = Rm/Re is unity. For details on the numerical methods used in FLASHv4, including viscous and resistive dissipation, see [21] and [22]. Both codes employ divergence cleaning [23] to maintain ∇ · B = 0. A state of homogeneous and isotropic turbulence is reached by supersonic stochastic driving in the momentum equation (given by an Ornstein-Uhlenbeck process) at small wave-numbers, similar to [24,25]. Thus, the forcing field is evolving in time and space. The associated large auto-correlation time-scale T of the forcing translates to the eddy turnover time of the largest, energycontaining eddies. It is therefore the chosen unit of time in the following.
We explore a range of parameters. The initial strength of the magnetic field is set by the plasma β p -the ratio of thermal to magnetic pressure. The final sonic Mach number M s is determined by the forcing amplitude. For the Enzo simulations we have initial β p = 0.25, 2.5, 25, with M s ≈ 2.5 after t ≈ 2T turnover times. The FLASHv4 simulations reach M s ≈ 4, 2 for initial β p = 1, 5, keeping instead constant Alfvenic Mach number M a ≈ 3. We discard all initial data affected by transients (before a simulation time of t = 2T ) and analyze consequent snapshots taken approximately in intervals of 0.15T and 0.1T for the Enzo and FLASHv4 datasets, respectively.
The analysis begins with the application of a low-pass Gaussian (test) filter to the equations of motion. In the context of the closures we investigate, filtered quantities (i.e. density, velocity, and magnetic field) are interpreted as resolved, while the remainders represent the unresolved small scales. We can then compute τ and E both directly from (1) and (2), and from their respective traditional and nonlinear closures.
The determination of the length-scale of the filter bears some consideration. It needs to fall within an intermediate range of length scales, away from the particular effects of both the large-scale forcing and the small-scale dissipation. The largest scale of the system is the full box size L and corresponds to the wavenumber n = 1, while the smallest scale is given by the Nyquist wavenumber n Nyq = N/2 = 256 for the linear numerical resolution N = 512 grid cells. The turbulence injection wavenumber is n inj = 2 (corresponding to half the box size L/2) in both codes, which is why the energy spectra, as illustrated in figure 1, peak there. The figure shows the mean kinetic and total (kinetic plus magnetic) energy spectra as a function of wavenumber. Since the stochastic forcing is implemented only in the momentum equation both for Enzo and FLASHv4, the kinetic energy spectrum exhibits the most direct imprint of the forcing itself. Conversely, the total energy spectrum carries the overall effect of the smallscale dissipation through both kinetic and magnetic channels. Figure 1 demonstrates that our simulations produce approximate power-law scaling within a narrow range of wavenumbers, which is indicative of self-similar turbulent fluctuations [26]. This can be interpreted as inertial range dynamics, although the nature of the inertial range in compressible MHD turbulence is still not fully understood [27,28,29,30]. Furthermore, this range separates the forcing scale and the dissipation scales and is not affected by numerical diffusion in the absence of a bottleneck effect as demonstrated by [31]. The vertical dotted line in figure 1 indicates our chosen filter length scale, corresponding to ∆ = 16 grid cells or wavenumber n filter = N/(2∆) = 16. This filter scale falls within the range of the self-similar power-law range for both the kinetic and total energy spectra. This is why we use this ideal scale for our filter in the following analysis.
Additionally, this provides the motivation to treat data from both simulations on equal footing, even though FLASHv4 has explicit viscosity and diffusivity while Enzo solves the ideal MHD equation, subject to numerical dissipation only.
In order to incorporate coordinate independence, a scalar field is chosen for comparison, the SGS energy flux Σ, i.e. the term responsible for the transfer of SGS energy between resolved and unresolved scales. Its components associated with the Reynolds and Maxwell stresses and the EMF are Σ u = τ u ij S ij , Σ b = τ b ij S ij and Σ E = E ·J , respectively (see appendix of [10] for the detailed SGS energy equation). Here, we substitute (1) and (2) to obtain exact fluxes Σ and match these to model fluxes Σ that employ the corresponding closures. For example, in the case of the eddy-viscosity  closure for the deviatoric turbulent Reynolds stress tensor we compare the exact flux with the model flux On the one hand, the comparison involves the determination of the constant (in space and time), dimensionless closure coefficients C . They are computed individually for each snapshot by minimizing the error between Σ and Σ in the least-square sense. This allows to further test the constancy of the coefficients with respect to time and plasma parameters. On the other hand, the general performance of the closure is gauged by computing the Pearson correlation coefficient of Σ and Σ , where the obtained closure coefficients are substituted in. Several assumptions should be pointed out concerning this validation technique. Firstly, the simulation data we have available for comparison fall short of realistic astrophysical parameters, e.g. with regards to Reynolds numbers and resolution. In that sense, it would be interesting to use higher resolution direct numerical simulation data or three-dimensional observations or experimental results. The problem is that experimental data for supersonic compressible turbulent plasmas are not available and obtaining realistic Reynolds numbers is computationally challenging for astrophysical parameters. However, as seen from figure 1, the data we have are sufficiently well resolved for our analysis. Secondly, in choosing the SGS energy flux Σ, as a diagnostic variable, we implicitly assume that in the context of homogeneous and isotropic turbulence the turbulent transport (encoded by terms of the form ∇ · ( u · τ ) and ∇ · B × E ) averages out to zero on subgrid scales. This assumption can nevertheless be easily relaxed by incorporating further diagnostic variables. Finally, we have focused on the SGS energy since it increases monotonically with the strength of turbulence regardless of the type of turbulence (e.g. compressive or solenoidal, weak or strong, etc.).
As an extension, the other two quadratic MHD invariants -the magnetic helicity and cross-helicity, may further highlight distinct turbulence properties present in particular flow configurations. These should be kept in mind as further avenues of investigation, once a preferred closure has been identified by the described validation technique.

Results
The fitting results for the SGS stress tensors' energy flux are given in figure 2 for the isotropic components and figure 3 for the deviatoric components. The isotropic parts of τ u and τ b are given by the SGS energies τ ii = 2 3 E sgs from (3) and (7). Both, the coefficient values of kinetic part (figure 2(a) top panel) and the magnetic part ( figure 2(b) top panel), have a small spread within a factor of two across time and all simulations. Furthermore, closure and data are highly correlated (bottom panels) with a median correlation coefficient of 0.90 and 0.91, respectively. More detailed numerical values of these and all following coefficients and correlations are listed in table 1.
The differences in the deviatoric parts τ u * in figure 3(a) and τ b * in figure 3(b) between the nonlinear and the eddy-viscosity closures are apparent. While our nonlinear closure exhibits approximately constant coefficient values and correlations over time in all simulations, the kinetic eddy-viscosity closure shows a correlation weaker by ≈ 0.2 and bigger spread in coefficient values. Moreover, the magnetic eddy-viscosity closure is effectively uncorrelated with the simulation data and the coefficients can even switch sign at different times. The performance of the different closures can be understood from figure 4, where we plot the energy flux distributions Σ u * and Σ b * for a single snapshot. A negative flux Σ u * < 0 corresponds to a forward energy cascade -the transfer of energy     from resolved to subgrid scales, because Σ u * appears as a sink term in the SGS kinetic and magnetic energy evolution equations and as a source term in the respective resolved energy equations. Conversely, a positive flux corresponds to an inverse energy cascade, i.e. transport of energy from subgrid to resolved spatial scales. The general distribution of the actual fluxes in figure 4 is representative for all snapshots. The kinetic SGS energy fluxes are globally almost 1:1 in both directions of the turbulent cascade with a slight tendency towards the forward cascade. However, the forward cascade is about 3-10 times stronger depending on the parameters as indicated by the position of the peaks in the distribution. For this reason, the kinetic eddy-viscosity closure shows a moderate correlation even though it captures only the forward energy cascade -from large to small scales. In fact, since under the eddy-viscosity hypothesis the kinetic SGS energy flux has the form Σ u * = −ν u ρ| S * | 2 , see (9), any model in which the eddy-viscosity ν u has a definite signature with respect to space cannot reproduce a bi-directional energy cascade that is well represented by the nonlinear closure. In contrast to the kinetic SGS energy flux, the global magnetic flux clearly has a preferred direction. Depending on the parameters, between 60% and 80% of cells have a positive SGS magnetic flux indicating energy transfer from large to small scales. Nevertheless, the difference in strength is less pronounced as the overall forward flux is about 2 times stronger than the backward one. Again, these properties are well captured by the proposed nonlinear closure whereas the magnetic eddy-viscosity closure is poorly correlated in both strength and magnitude. Finally, it should be noted, that the nonlinear closures also work very well with the Smagorinsky energy closure. Exchanging the exact expressions E sgs in (4) and (5) with E sgs only slightly reduces the correlations (max. 5%) and the coefficients remain constant up to the second significant figure (not plotted here). Moving on to the EMF, the nonlinear closure outperforms the traditional α − β − γ closure in almost all datasets, maintaining a constant coefficient with a median correlation of 0.79 (figure 5). The traditional closure exhibits consistently weaker correlations, despite the increased flexibility of three free coefficients. Only C E β , related to the turbulent resistivity term in the EMF, is approximately constant, whereas C E α and C E γ fail to maintain steady values or consistent signature. The reason for the consistently better correlations of the nonlinear closure is hinted at in figure 6. This probability density plot of the local alignment between E and E demonstrates that the traditional closure is almost randomly aligned (flat distribution) whereas the nonlinear closure approaches the desired δ-distribution at 0 • .

Conclusions and outlook
In summary, we have proposed a set of constant coefficient closures for the SGS stress and EMF in the filtered MHD equations and conducted a priori tests. The tests we performed do show that the new nonlinear closures perform significantly better than traditional, phenomenological closures with respect to both structural and functional diagnostics. The tests consist of filtering Enzo and FLASHv4 simulations of homogeneous, isotropic turbulence and comparing the resulting SGS terms to their respective closures (dependent only on the filtered fields). All quantities are compared via their contributions to the SGS energy flux Σ , where the closure coefficients are computed by individual least-square fitting. In addition, the alignment for the EMF vector is investigated. All new coefficients correlate well with the data. They are constant over time and as a direct consequence the proposed closures may be implemented in large-eddy simulations without the need for a computationally expensive dynamical procedure which computes the coefficient values at run time. In addition, the coefficients remain constant across simulation runs from two different codes and a wide range of plasma parameters, suggesting that the proposed closures capture an underlying physical mechanism at work in highly compressible turbulent plasma flows. Moreover, the new closures successfully represent the turbulent magnetic pressure, reproduce the bi-directional energy cascade and are well aligned with the EMF. We recognize the slightly lower correlation of the nonlinear closure in the EMF than in the SGS stress counterpart, suggesting small room for improvement.
Nevertheless, the performance improvement over the traditional closures already supports the implementation and validation of the new closures in an SGS model for large-eddy simulations of compressible turbulent plasma flows. These simulations would then allow us to infer the effect of the proposed model on the large scale flow in practice. Potential applications include accretion disks [32], star-forming magnetized clouds [33,34] and plasmas on cosmological scales [35,36,37,38,39,40].