Energy flux decomposition in magnetohydrodynamic turbulence

Abstract

In hydrodynamic (HD) turbulence, an exact decomposition of the energy flux across scales has been derived that identifies the contributions associated with vortex stretching and strain self-amplification (Johnson, Phys. Rev. Lett., vol. 124, 2020 104501; J. Fluid Mech., vol. 922, 2021, A3) to the energy flux across scales. Here, we extend this methodology to general coupled advection–diffusion equations and, in particular, to homogeneous magnetohydrodynamic (MHD) turbulence. We show that several MHD subfluxes are related to each other by kinematic constraints akin to the Betchov relation in HD. Applied to data from direct numerical simulations, this decomposition allows for an identification of physical processes and for the quantification of their respective contributions to the energy cascade, as well as a quantitative assessment of their multi-scale nature through a further decomposition into single- and multi-scale terms. We find that vortex stretching is strongly depleted in MHD compared with HD, and the kinetic energy is transferred from large to small scales almost exclusively by the generation of regions of small-scale intense strain induced by the Lorentz force. In regions of large strain, current sheets are stretched by large-scale straining motion into regions of magnetic shear. This magnetic shear in turn drives extensional flows at smaller scales. Magnetic energy is transferred from large to small scales predominantly by the aforementioned current-sheet thinning in regions of high strain. The contributions from current-filament stretching – the analogue to vortex stretching – and from bending of magnetic field-lines into current filaments by vortical motion are both almost negligible, although the latter induces strong backscatter of magnetic energy. Consequences of these results for subgrid-scale turbulence modelling are discussed.

1 Introduction

Turbulence in electrically conducting fluids and plasmas is of relevance to a variety of processes in geophysical and astrophysical situations, as well as in industry (Weiss & Proctor 2014; Davidson 2016) and for nuclear fusion under magnetic confinement. For example, the solar wind is turbulent (Bruno & Carbone 2013), convection-driven turbulence occurs in planetary cores and in the outer layers of stars (Jones 2011), turbulence on ion and electron scales affects plasma confinement in magnetic confinement fusion reactors (Freidberg 2007), and the heat transfer in liquid metal cooling applications is dependent on the level of turbulence in the flow (Davidson 1999). Even though these systems are very different in terms of features like the presence of strong background magnetic fields, the level of magnetic field fluctuations, temperature gradients, density fluctuations, domain geometry or the level of collisionality, they nonetheless share fundamental nonlinear processes that define energy conversion and inter-scale energy transfers, at least on scales where the fluid approximation is applicable. Even in the simplest case of magnetohydrodynamic (MHD) turbulence, despite considerable theoretical progress (e.g. Goldreich & Sridhar 1995; Biskamp 2003; Zhou, Matthaeus & Dmitruk 2004; Petrosyan et al. 2010; Brandenburg, Sokoloff & Subramanian 2012; Tobias & Cattaneo 2013; Beresnyak 2019; Oughton & Matthaeus 2020; Schekochihin 2022), the physical nature of these processes remain opaque.

Moreover, the typical parameter ranges in which MHD turbulence develops in Nature are far from those attainable with direct numerical simulation (DNS) (Plunian, Stepanov & Frick 2013; Miesch et al. 2015; Schmidt 2015). As a consequence, the demand for approximations and subgrid-scale (SGS) models for large-eddy simulation (LES) of MHD turbulence that are able to capture the effects of unresolved small-scale fluctuations – that govern important processes such as magnetic reconnection and plasma heating – is increasing (Miesch et al. 2015). However, constructing such models is a challenge due to small-scale anisotropy (Shebalin, Matthaeus & Montgomery 1983; Oughton, Priest & Matthaeus 1994; Goldreich & Sridhar 1997; Tobias & Cattaneo 2013), strong intermittency in the magnetic fluctuations as observed in numerical simulations (e.g. Mininni & Pouquet 2009; Sahoo, Perlekar & Pandit 2011; Yoshimatsu et al. 2011; Rodriguez Imazio et al. 2013; Meyrand, Kiyani & Galtier 2015) and in the solar wind (e.g. Veltri 1999; Salem et al. 2009; Wan et al. 2012; Matthaeus et al. 2015), insufficient magnetic-field amplification and dynamo growth observed in LES of MHD turbulence (Haugen & Brandenburg 2006), and whether magnetic-field fluctuations are maintained by the flow or by an external electromagnetic force (Alexakis & Chibbaro 2022). For a summary of the SGS modelling effort and its challenges, we refer to the review articles by Miesch et al. (2015) and Schmidt (2015).

Here, we focus on energy transfer across scales in statistically stationary homogeneous MHD turbulence in a saturated nonlinear dynamo regime without a mean magnetic field, and with negligible levels or cross- and magnetic helicity. The total energy cascade in this case is direct (Aluie & Eyink 2010), transferring energy from the large to the small scales in a scale-local fashion. Herein, we use the term cascade to indicate a mean flux and refer to pointwise upscale and downscale energy transfers as inverse and direct transfer events, respectively, or in the former case also as backscatter. The aims of this paper are: (i) to understand the physical mechanisms that govern the MHD energy cascade and (ii) to quantify their importance and provide guidance for SGS modelling (Johnson 2022). In terms of turbulence theory, this corresponds to understanding physical properties of the SGS stresses as a function of scale. Eyink (2006) introduced a viable approach for this that involves expanding the SGS tensors in terms of vector field gradients. Here, we follow a filtering approach, generalising an exact gradient-based decomposition of the hydrodynamic (HD) energy fluxes (Johnson 2020, 2021) to coupled advection–diffusion equations and hence to the MHD equations. This methodology distinguishes between terms that are local in scale, corresponding to the first term in the gradient expansion and those which are truly multi-scale, providing a closed expression for the remainder of the series expansion. Expressing SGS stresses through vector-field gradients results in a decomposition of the energy fluxes in terms of different tensorial contractions between strain-rate, vorticity, current and magnetic strain, and as such facilitates the physical interpretation of such sub-fluxes. The provision of closed expressions allows for a quantification of the relative contribution of all terms to the energy cascade using data obtained by direct numerical simulation (DNS). In homogeneous and isotropic HD turbulence, where only the inertial term is present, the decomposition identifies three processes that transfer kinetic energy across scales, vortex stretching, strain self-amplification and strain-vorticity alignment, and quantifies their relative contribution to the energy cascade (Johnson 2020, 2021). Similarly, the direct cascade of kinetic helicity is carried by three different processes, vortex flattening, vortex twisting and vortex entanglement (Capocci et al. 2023).

In MHD, the total energy transfer can be split into four subfluxes, Inertial, Maxwell, Dynamo and Advection,¹ with the former two originating from the Reynolds and Maxwell stresses in the momentum equation, and the latter two from stresses in the induction equation that result in the advection and bending/stretching of magnetic field lines by the flow. As Dynamo and Advection terms have a common physical origin, the electric field, often only their sum is considered in a priori analyses of DNS data (Aluie 2017; Offermans et al. 2018; Alexakis & Chibbaro 2022) and a posteriori in LES (Zhou & Vahala 1991; Müller & Carati 2002; Grete et al. 2016; Kessar, Balarac & Plunian 2016; Vlaykov et al. 2016). However, in the present work, it will prove instructive to consider them separately. Here, we generalise and apply the aforementioned decomposition to each of these four subfluxes.

In doing so, we identify a single process, current-sheet thinning, to be the main contributor to the forward cascade of magnetic energy. Contributions from strain-induced current-filament stretching, the formal analogue to vortex stretching, are subdominant. Furthermore, we find that vortex-stretching and strain self-amplification are strongly suppressed at all length scales. Instead the back-reaction of current-sheet thinning on the flow through the Lorentz force constitutes the main transfer of kinetic energy from large to small scales, with contributions related to current-filament stretching turning out to be negligible.

The structure of the paper is as follows: we begin in § 2 with an outline of how the generalised method can be applied to obtain MHD energy subfluxes. In § 3, we discuss the numerical details and the associated datasets on which we performed the filtering analysis. In § 4, we consider each subflux decomposition, showing results for both mean terms and fluctuations. Ramifications of those results for SGS modelling are considered in § 5. In § 6, we discuss our main results and indicate future work directions. Several appendices flesh out some aspects of the derivations and analysis.

2 Theory

In this section, we begin by sketching the derivation of the coarse-grained energy equations for MHD and giving the definitions of the scale-space energy fluxes that appear in them. Subsequently, we show how each flux can be decomposed in terms of physically distinct contributions, and discuss their physical interpretations.

2.1 Coarse-grained energy equations

Our starting point is incompressible three-dimensional (3-D) homogeneous MHD turbulence. The primary dynamical variables are then the fluctuation velocity $\boldsymbol {u} (\boldsymbol {x}, t)$ and the fluctuation magnetic field $\boldsymbol {b} (\boldsymbol {x},t)$, where we measure the latter in Alfvén speed units: $\boldsymbol {b} / \sqrt {4{\rm \pi} \rho } \to \boldsymbol {b}$, with $\rho$ the uniform mass density. We consider situations with no mean magnetic field since these are more likely to exhibit global isotropy. The governing equations, with allowance for hyper-dissipation, are

(2.1)\begin{gather} \partial_t u_i+ \partial_j \left( u_i u_j \right) ={-} \partial_i \left( p + \frac{b^2}{2} \right) + \partial_j \left( b_i b_j \right) + \nu_\alpha ({-}1)^{\alpha +1} {\nabla}^{2\alpha} u_i + F_i , \end{gather}

(2.2)\begin{gather}\partial_t b_i + \partial_j \left (b_i u_j \right) = \partial_j \left( u_i b_j \right)+ \mu_\alpha ({-}1)^{\alpha +1} {\nabla}^{2\alpha} b_i, \end{gather}

(2.3)\begin{gather}\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u} = 0 , \quad \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{b} = 0 . \end{gather}

Here, $p$ is the pressure, $\boldsymbol {F}$ is a (large-scale) velocity forcing, $\nu _\alpha$ and $\mu _\alpha$ are the hyper-viscosity and hyper-resistivity, and $\alpha$ denotes the power of the Laplacian operator employed in the hyper-dissipation. Standard Laplacian dissipation corresponds to the case $\alpha =1$.

The MHD variables, and equations, may be spatially coarse-grained using a suitable filtering field, $G^\ell (\boldsymbol {r} )$ (Germano 1992; Aluie 2017). The role of $G^\ell (\boldsymbol {r} )$ is to strongly suppresses structure at scales less than the filtering scale $\ell$. For example, the filtered velocity field is

(2.4)\begin{equation} \bar{\boldsymbol{u}}^\ell(\boldsymbol{x}) = \int \,\mathrm{d} \boldsymbol{r} \, G^\ell (\boldsymbol{r} ) \, \boldsymbol{u} ( \boldsymbol{x} + \boldsymbol{r}) . \end{equation}

This can be interpreted as a weighted average of $\boldsymbol {u}$ centred on the position $\boldsymbol {x}$. The weighting function $G^\ell$ decays very rapidly to zero at distances greater than a few $\ell$ from $\boldsymbol {x}$ and satisfies some other weak restrictions, such as smoothness and having a volume integral of unity. Filtering is a linear operation and commutes with differentiation, properties of which we will make considerable use below. From § 2.2 onwards, we will specialise to a Gaussian filter, but in this section, a specific choice of filter is not needed.

Coarse-graining of (2.1)–(2.2) introduces four SGS stress tensors, $\tau ^\ell (\cdot, \cdot )$, associated with the advective-type nonlinear terms (those containing a $\partial _j$). These each have the form

(2.5)\begin{equation} \tau^\ell (f_i, g_j) = \overline{f_i g_j}^\ell - \overline{f_i}^\ell \, \overline{g_j}^\ell , \end{equation}

where $\boldsymbol {f}$ and $\boldsymbol {g}$ are the solenoidal vectors appearing in a $g_j \partial _j \, f_i$ advective-type term. We remark that with this notation, the advecting field is the second argument in a $\tau ^\ell (\cdot,\cdot )$.

To obtain the equations governing the (pointwise) evolution of the coarse-grained kinetic energy $E^\ell _u (\boldsymbol {x},t) = \tfrac {1}{2} \bar {\boldsymbol {u}}^\ell \boldsymbol {\cdot } \bar {\boldsymbol {u}}^\ell$ and magnetic energy $E^\ell _b (\boldsymbol {x},t) = \tfrac {1}{2} \bar {\boldsymbol {b}}^\ell \boldsymbol {\cdot } \bar {\boldsymbol {b}}^\ell$, one filters (2.1)–(2.2) and then multiplies by $\bar {\boldsymbol {u}}^\ell$ and $\bar {\boldsymbol {b}}^\ell$, respectively (e.g. Zhou & Vahala 1991; Kessar et al. 2016; Aluie 2017; Offermans et al. 2018; Alexakis & Chibbaro 2022). The result can be written as

(2.6)\begin{gather} \partial_t {E}^\ell_u + \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{{\mathcal{J}}}_u ={-} \varPi^{I,\ell} - \varPi^{M,\ell} - {\mathcal{W}}^\ell - {\mathcal{D}}_{u} + \bar{\boldsymbol{u}}^\ell \boldsymbol{\cdot} \bar{\boldsymbol{F}}^\ell, \end{gather}

(2.7)\begin{gather}\partial_t {E}^\ell_b+ \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{{\mathcal{J}}}_b ={-} \varPi^{A,\ell} - \varPi^{D,\ell} + {\mathcal{W}}^\ell - {\mathcal{D}}_{b} , \end{gather}

where the $\boldsymbol {{\mathcal {J}}}$ terms account for the spatial transport of energy and the $\varPi^\ell$ terms embody energy fluxes (i.e. transfer across scale $\ell$). Our sign convention for the definitions of the $\varPi^\ell$ (see below) means that $\varPi^\ell > 0$ corresponds to forward transfer of energy, i.e. to scales smaller than $\ell$. The ${\mathcal {D}}$ terms represent (hyper-)dissipative effects. Also present is the resolved scale conversion (RSC) term, ${\mathcal {W}}^\ell = \bar {b}^\ell_i \bar {b}^\ell_j \partial _j \bar {u}^\ell_i$, here expressed as in Aluie (2017). This appears with opposite sign in each equation, and represents an exchange between large-scale kinetic and magnetic energies. An important point is that it is not an energy flux term since it does not involve energy transfer across scale $\ell$. It supports several interpretations including as: (i) the rate of work done on the large-scale flow by the large-scale Lorentz force and (ii) the energy gained by $\bar {\boldsymbol {b}}^\ell$ as it is distorted by $\bar {\boldsymbol {u}}^\ell$ or vice versa. Detailed forms for the spatial transport currents, $\boldsymbol {{\mathcal {J}}}_u, \boldsymbol {{\mathcal {J}}}_b$, depend on the form employed for ${\mathcal {W}}^\ell$ and are available elsewhere (e.g. Kessar et al. 2016; Aluie 2017; Offermans et al. 2018; Alexakis & Chibbaro 2022), while the problem of Galilean invariance has been addressed by Offermans et al. (2018).

Our primary interest herein centres on the pointwise energy flux (at scale $\ell$) terms, denoted by $\varPi ^\ell (\boldsymbol {x}, t)$, together with their volume averages, $\left \langle \varPi ^\ell \right \rangle$. For any choice of filter, these can be expressed in terms of contractions of filtered gradient tensors and SGS stress tensors:

(2.8)\begin{gather} \varPi^{I,\ell} ={-} \frac{\partial{\bar{u}_i^\ell}}{\partial{x_j}} \, \tau^\ell(u_i,u_j) , \end{gather}

(2.9)\begin{gather}\varPi^{M,\ell} = \frac{\partial{\bar{u}_i^\ell}}{\partial{x_j}} \, \tau^\ell(b_i,b_j) , \end{gather}

(2.10)\begin{gather}\varPi^{A,\ell} ={-} \frac{\partial{\bar{b}_i^\ell}}{\partial{x_j}} \, \tau^\ell(b_i,u_j) , \end{gather}

(2.11)\begin{gather}\varPi^{D,\ell} = \frac{\partial{\bar{b}_i^\ell}}{\partial{x_j}} \, \tau^\ell(u_i,b_j). \end{gather}

Like the $\tau ^\ell ( \cdot, \cdot )$, these arise in connection with the four advection type nonlinearities in (2.1)–(2.2) that we refer to as the Inertial, Maxwell (meaning from the Lorentz force), Advection and Dynamo terms. Note the capitalisation. Taken together with (2.6) and (2.7), these definitions of the $\varPi^\ell$ terms mean that the interpretation of the direction of an energy flux does not depend on which flux it is. This is why (2.9) and (2.11) lack a leading minus sign. Specifically, a positive value for any one of these fluxes corresponds to transfer of energy from scales greater than $\ell$ to scales smaller than $\ell$. Clearly, $\varPi ^{I,\ell } + \varPi ^{M,\ell }$ is the net flux of $E^\ell _u$, and $\varPi ^{A,\ell } + \varPi ^{D,\ell }$ that for $E^\ell _b$. As is well known, $\varPi ^{A,\ell }$ and $\varPi ^{D,\ell }$ have a common origin and may be readily combined to obtain the magnetic energy flux associated with the curl of the induced electric field.

We remark that energy transfer in MHD turbulence has traditionally been discussed in a Fourier-space approach (e.g. Dar, Verma & Eswaran 2001; Verma 2004; Alexakis, Mininni & Pouquet 2005; Mininni, Montgomery & Pouquet 2005; Teaca et al. 2009; Linkmann et al. 2017; Verma 2019). The filtering approach used here results in equivalent expressions for the mean inertial flux if the filter is chosen to be a Galerkin projector. However, in MHD, the two approaches differ in important details. In the Fourier-space approach, the transfer terms originating from the Lorentz force (momentum equation) and the field-line stretching term (induction equation) retain contributions that only involve a coupling among resolved scales. That is, the RSC term ${\mathcal {W}}^\ell$ is not separated out from the fluxes. This is the reason these terms do not vanish in the limit of $k\to \infty$ in the Fourier-based approach, in contrast to the terms defined in (2.8)–(2.11) using the SGS stresses, which therefore are genuine flux terms. As a consequence, the Fourier-based definitions can result in misinterpretations of the multi-scale nature of kinetic-to-magnetic energy conversion and of the role of the Lorentz force in inducing kinetic energy transfer across scales. For further details, see the works of Aluie (2017) and Offermans et al. (2018). An in-depth discussion of the SGS-based definition of energy flux and a comparison with the Fourier-based approach is supplied by Aluie & Eyink (2009).

2.2 Gaussian filter

Rather remarkably, the choice of a Gaussian filter enables an analytic determination of the SGS stresses and fluxes in terms of field gradients with contributions from the resolved scales and the subfilter scales. In the remainder of the paper, we therefore employ an isotropic Gaussian filter

(2.12)\begin{equation} G^\ell (\boldsymbol{r} ) = \frac{1}{(2{\rm \pi}\ell^2)^{3/2}} \exp \left( -\frac{|\boldsymbol{r}|^2}{2 \ell^2} \right) . \end{equation}

2.3 Exact expressions for the $\tau ^\ell$ and $\varPi ^\ell$

Employing the Fourier transform of the Gaussian filter (2.12), Johnson (2020, 2021) showed that the filtered version of an arbitrary field $\boldsymbol {u}(\boldsymbol {x},t)$ satisfies a diffusion equation,

(2.13)\begin{equation} \frac{\partial{ \bar{u}^\ell_j}}{\partial{(\ell^2)}} =\frac{1}{2} \nabla^2 \bar{u}^\ell_j, \quad \left.\bar{u}^\ell_j \right\vert_{\ell=0} = u_j( \boldsymbol{x}, t ), \end{equation}

where $\ell ^2$ is the time-like variable. It was further shown that the associated SGS stress tensor, $\tau ^\ell (u_i, u_j )$, obeys a forced version of this diffusion equation. The forcing term is $\bar {A}^\ell _{ik} \bar {A}^\ell _{jk}$, where $\bar {A}^\ell _{ik} = \partial _k \bar {u}^\ell _i$ is the gradient tensor for the filtered field. An exact solution for $\tau ^\ell ( u_i, u_j)$ was obtained that depends on $\bar {A}^\phi _{ik}$ for all scales $\phi \le \ell$. Substituting this into the equation for the SGS Inertial flux, (2.8), produces an exact solution for this flux corresponding to the exact summation of the perturbation series proposed by Eyink (2006).

Happily, this approach is readily extended to MHD and may be used to calculate the elements contained in (2.8)–(2.11). Below, we outline how to achieve this for the particular case of the magnetic energy subflux (2.10) that originates with the advection term in the induction equation (i.e. $\boldsymbol {u} \boldsymbol {\cdot } \boldsymbol {\nabla } \boldsymbol {b}$). More details, plus the general case of three distinct solenoidal fields, are available in Appendix A; see also Appendix C and Capocci et al. (2023).

We seek an exact solution for $\tau ^\ell ( b_i, u_j) = \overline {b_i u_j}^\ell - \overline {b_i}^\ell \overline {u_j}^\ell$. Clearly, $\overline { b_i u_j }^\ell$ will also satisfy (2.13) since that equation holds for any (Gaussian) filtered field. Together with the product rule expansion of $\partial ( \bar {b}^\ell _i \bar {u}^\ell _j ) / \partial {\ell ^2}$, this yields

(2.14)\begin{equation} \left( \frac{\partial{}}{\partial{\ell^2}} - \frac{1}{2} \nabla^2 \right) \tau^\ell (b_i, u_j) = \bar{B}^\ell_{ik} \bar{A}^\ell_{jk} , \quad \tau^{\ell=0} (b_i, u_j) = 0 , \end{equation}

where $\bar {B}^\ell _{ik} = \partial _k \bar {b}^\ell _i$ is the gradient tensor for $\bar {b}^\ell _i$. The solution can be written in the form

(2.15)\begin{equation} \tau^\ell (b_i, u_j) =\ell^2 \, \bar{B}^\ell_{ik} \bar{A}^\ell_{jk} + \int_0^{\ell^2} \,\mathrm{d}\theta\left( \overline{\overline{B}_{ik}^{\sqrt{\theta}} \, \overline{A}_{jk}^{\sqrt{\theta}} }^\phi -\overline{\overline{B}_{ik}^{\sqrt{\theta}}}^\phi \, \overline{\overline{A}_{jk}^{\sqrt{\theta}}}^\phi\right), \end{equation}

where $\phi = \sqrt {\ell ^2 - \theta }$. The first term on the right-hand side is a ‘single-scale’ piece as it contains only resolved scale terms (cf. Clark, Ferziger & Reynolds 1979). In terms of other work, it corresponds to the nonlinear model employed by Leonard (1975), Borue & Orszag (1998) and Meneveau & Katz (2000), and is the first-order term in the expansion of Eyink (2006). It is also the leading-order term of the power law expansion in the filter limit going to zero relative of any filter kernel with finite moments (cf. § 13.4.4 of Pope 2000).

The second term, involving an integral over all scales smaller than $\ell$, is manifestly a multi-scale contribution that includes subfilter-scale field gradients. Note that the integrand in (2.15) can itself be written as an SGS stress, but one based on the field gradients rather than the fields themselves: $\tau ^\phi ( \bar {B}_{ik}^{\sqrt {\theta }}, \bar {A}_{jk}^{\sqrt {\theta }} )$.

Contracting with $\bar {B}_{ij}^\ell$ provides an exact expression for (2.10):

(2.16)\begin{align} \varPi^{A,\ell} & ={-} \ell^2 \, \bar{B}_{ij}^\ell \bar{B}_{ik}^\ell \bar{A}_{jk}^\ell - \bar{B}_{ij}^\ell \int_0^{\ell^2} \,\mathrm{d}\theta \left( \overline{\overline{B}_{ik}^{\sqrt{\theta}} \,\overline{A}_{jk}^{\sqrt{\theta}} }^\phi - \overline{\overline{B}_{ik}^{\sqrt{\theta}}}^\phi \, \overline{\overline{A}_{jk}^{\sqrt{\theta}}}^\phi \right) \end{align}

(2.17)\begin{align} & = \varPi^{A,\ell}_s + \varPi^{A,\ell}_m , \end{align}

where the subscripts $s$ and $m$ denote the single- and multi-scale contributions, respectively. It is evident that all the SGS energy fluxes, (2.8)–(2.11), and also other SGS fluxes (e.g. for helicities), can be written strictly in terms of (multi-scale) gradients of the velocity and magnetic vector fields. Appendix A contains further details. When discussing the individual SGS energy fluxes in § 4, we will make regular reference to (A3), the generalised form of (2.17).

The tensor contractions present in (2.17) can be expressed as the trace of the matrix products involved, after appropriate use of the transpose operation (superscript $t$). For example, $\bar {B}^\ell _{ij} \bar {B}^\ell _{ik} \bar {A}^\ell _{jk} = \mathrm {Tr}\left \{ { (\bar {\boldsymbol {B}}^\ell )^t \,\bar {\boldsymbol {B}}^\ell (\bar {\boldsymbol {A}}^\ell )^t} \right \}$.

Further insight into the physics of the scale-space flux $\varPi ^{A,\ell }$ may be extracted by expressing each gradient tensor as the sum of its index-symmetric and index-antisymmetric components. Let us write $S_{ij} = ( A_{ij} + A_{ji} ) /2$ and ${\varOmega }_{ij} = ( A_{ij} - A_{ji} ) / 2$, respectively the (velocity) rate-of-strain tensor and the rotation rate tensor, with ${\varOmega }_{ij} = -\epsilon _{ijk} {\omega }_k /2$ in terms of the vorticity ${\omega }_k = \epsilon _{ijk} \partial _i {u}_j$. Similarly, $B_{ij} = \varSigma _{ij} + {J}_{ij}$, where the non-zero elements of $J_{ij} = ( B_{ij} - B_{ji} ) / 2 = -\epsilon _{ijk} j_k/2$ are essentially the components of the electric current density $\boldsymbol {j} = \boldsymbol {\nabla } \times \boldsymbol {b}$ and the magnetic strain-rate tensor $\varSigma _{ij} = (B_{ij} + B_{ji})/2$. We will of course require the filtered versions of all of these quantities.

Ostensibly, this decomposition gives eight single-scale and eight multi-scale sub-fluxes; see (A7). However, properties of the trace of particular products of symmetric or antisymmetric matrices mean some of these may vanish, cancel or be equivalent. In the present case, one obtains, in connection with the single-scale contributions,

(2.18)\begin{align} \mathrm{Tr}\left\{ { \left(\bar{\boldsymbol{\mathsf{B}}}^\ell\right)^t \, \bar{\boldsymbol{\mathsf{B}}}^\ell \, \left(\bar{\boldsymbol{\mathsf{A}}}^\ell\right)^t} \right\} & = \mathrm{Tr}\left\{ { \left( \bar{\boldsymbol{\varSigma}}^\ell - \bar{\boldsymbol{\mathsf{J}}}^\ell \right) \left( \bar{\boldsymbol{\varSigma}}^\ell + \bar{\boldsymbol{\mathsf{J}}}^\ell \right) \left( \bar{\boldsymbol{\mathsf{S}}}^\ell - \bar{\boldsymbol{\varOmega}}^\ell \right) } \right\} \end{align}

(2.19)\begin{align} & = \mathrm{Tr}\left\{ { \bar{\boldsymbol{\varSigma}}^\ell \bar{\boldsymbol{\varSigma}}^\ell \bar{\boldsymbol{\mathsf{S}}}^\ell - \bar{\boldsymbol{\mathsf{J}}}^\ell \bar{\boldsymbol{\mathsf{J}}}^\ell \, \bar{\boldsymbol{\mathsf{S}}}^\ell + 2 \bar{\boldsymbol{\varSigma}}^\ell \bar{\boldsymbol{\mathsf{J}}}^\ell \bar{\boldsymbol{\mathsf{S}}}^\ell } \right\}, \end{align}

so that there are only three distinct subfluxes; see Appendix A for details. We write the single-scale flux for the Advection term as

(2.20)\begin{equation} \varPi^{A,\ell}_s =\varPi^{A, \ell}_{s, \varSigma \varSigma S}+ \varPi^{A, \ell}_{s, J J S} + 2\varPi^{A, \ell}_{s, \varSigma J S}, \end{equation}

where $\varPi ^{A,\ell }_{s, PQR} = -\ell ^2 \mathrm {Tr}\left \{ {(\bar {\boldsymbol {\mathsf{P}}}^\ell )^t \bar { \boldsymbol {\mathsf{Q}}}^\ell (\bar { \boldsymbol {\mathsf{R}}}^\ell )^t} \right \}$ and each of $\boldsymbol {\mathsf{P}}, \boldsymbol {\mathsf{Q}},\boldsymbol {\mathsf{R}}$ are either symmetric or antisymmetric tensors.

The multi-scale flux contributions which explicitly contain subfilter-scale fluctuations can also be so decomposed, with the details given in Appendix A. The particular forms for $\varPi ^{A,\ell }_m$ are discussed in § 4.3.

3 Methods and data

To quantify the four energy fluxes $\varPi ^{Y,\ell }$ present in (2.6)–(2.7), and their decompositions, we employ outputs from numerical simulations of the MHD equations (2.1)–(2.3). To obtain an appreciable inertial range, we use hyperdiffusion ($\alpha = 5$), a comparison with standard diffusive MHD ($\alpha = 1$) is provided in Appendix E, always with $\nu _\alpha = \mu _\alpha$. The fluctuation fields, $\boldsymbol {u}$ and $\boldsymbol {b}$, have zero means and there is no background magnetic field (i.e. $B_0 = 0$). The equations are solved using fully dealiased Fourier pseudospectral codes in a triply periodic $(2{\rm \pi} )^3$ domain (Patterson & Orszag 1971; Canuto et al. 1988). The time advancement is via a second-order Runge–Kutta scheme with dealiasing implemented using the two-thirds rule.

As table 1 indicates, we use up to $2048^3$ grid points. The spatial resolution of the simulations is quantified by both the grid spacing $\Delta {x} = 2{\rm \pi} / N$ and the hyper-diffusive Kolmogorov scales $\eta ^u_{\alpha } = (\nu _\alpha ^3 / \varepsilon _u)^{1/(6 \alpha -2)}$ and $\eta ^b_{\alpha } = (\mu _\alpha ^3 / \varepsilon _b)^{1/(6 \alpha -2)}$, where $\varepsilon _u$ and $\varepsilon _b$ are the mean kinetic and magnetic energy dissipation rates (Borue & Orszag 1995). For adequate resolution, we require $\eta ^u_\alpha / \Delta {x} \gtrsim 1.3$ and $\eta ^b_\alpha / \Delta {x} \gtrsim 1.3$ (e.g. Donzis, Yeung & Sreenivasan 2008; Wan et al. 2010).

Table 1. Simulation parameters and key observables, where $N$ is the number of collocation points in each coordinate, $\alpha$ is the power of $\nabla ^2$ used in the hyper-diffusion, $E_u$ the (mean) total kinetic energy, $\nu _\alpha$ the kinematic hyperdiffusivity, $\varepsilon _u$ and $\varepsilon _b$ are the kinetic and magnetic energy dissipation rates, $L_u = (3 {\rm \pi}/4 E_u) \int _0^{k_\text {max}} \,\mathrm {d}{k} E_u(k)/k$ the longitudinal integral scale, $\tau = L_u/\sqrt {2E_u/3}$ the large-scale eddy-turnover time, and ${Re}$ is the Reynolds number. Furthermore, $\eta _{\alpha } = (\nu _\alpha ^3 / \varepsilon )^{1 / (6 \alpha -2)}$, $\eta ^u_{\alpha } = (\nu _\alpha ^3 / \varepsilon _u)^{1 / (6 \alpha -2)}$ and $\eta ^b_{\alpha } = (\mu _\alpha ^3 / \varepsilon _b)^{1 / (6 \alpha -2)}$ are the hyperdiffusive Kolmogorov scales calculated with respect to the total, viscous and Joule dissipation rates, respectively, $k_\text {max}$ the largest retained wavenumber component after de-aliasing, $\Delta t$ the mean of the snapshots sampling intervals, and # indicates the number of snapshots used in the averaging. The magnetic Prandtl number, $Pm = \nu _\alpha / \mu _\alpha$, the ratio between the hyperviscosity and magnetic hyperdiffusivity, equals unity for A1.

The forcing $f$ applied to the system is a drag-free Ornstein–Uhlenbeck process, active in the wavenumber band $k \in [2.5,5.0]$ for the MHD simulations, while the hydrodynamics dataset H2 is forced in the band $k \in [0.5,1.5]$. The snapshots, consisting of instantaneous velocity and magnetic fields, have been sampled about once per large-scale turnover time after the simulations reach statistically stationary states.

Figure 1 shows the time-averaged omnidirectional kinetic and magnetic spectra. The peak in the kinetic spectrum is due to the activity of forcing for $k \in [2.5,5.0]$, while the magnetic spectrum is considerably lower over that interval since the induction equation is not forced. In the (approximate) inertial range, both spectra have approximately power law scaling, with $E_b(k)$ close to $k^{-5/3}$ and $E_u(k)$ significantly shallower. As is typically seen in MHD simulations with no mean field, the Alfvén ratio, $E_u(k)\,/\,E_b(k)$, is less than unity in the inertial range, i.e. magnetic energy predominates at these scales.

Figure 1. Time-averaged omnidirectional spectra for the velocity and magnetic field (run A1). The grey region indicates the wavenumber band where the velocity field is forced: $k \in [2.5,5.0]$.

At high $k$, there is a steep and roughly coincident decrease of both spectra. This is a consequence of two factors. First, the employment of hyperviscosity makes the dissipation range more concentrated in the small scales, leading to the sharp fall-off. Second, unit Prandtl number ensures that the dissipation wavenumber band is the same for the kinetic and magnetic spectra. Note the clear superequipartition in the entire scaling range of the spectra, with $E_b(k) > E_u(k)$ for $10 \leq k \leq 400$. Equipartition at large $k$ in the inertial range as previously predicted for decaying turbulence (Müller & Grappin 2004; Haugen & Brandenburg 2006) is not observed, however, as $E_b(k)$ is steeper than $E_u(k)$ it does approach $E_u(k)$ resulting in approximate equipartition in the dissipative range.

To calculate an effective Reynolds number for a hyperdissipative system, we follow the approach described by Buzzicotti et al. (2018). There, the standard (Laplacian dissipation) integral-scale Reynolds number $Re=U L_u \, / \,\nu \propto \left(L_u \, / \, \eta_2 \right)^{4/3}$ (e.g. Batchelor 1970; Pope 2000) is replaced with one based on the ratio between the integral scale $L_u$ and the effective dissipation range scale $I_d$. Specifically, we employ

(3.1)\begin{equation} Re = C \left( \dfrac{L_u}{I_d} \right)^{4/3} , \end{equation}

where $I_d = \pi\,/\,\text{argmax}\left( k^2 E_u(k) \right)$ is the scale where the dissipation spectrum $k^2 E_u(k)$ has a maximum. Here, $C$ is a fit parameter that has to be estimated by comparing (3.1) with the common definition of the Reynolds number in a standard-viscosity run. Making use of this procedure, we obtain $C=40$ and $Re=9931$ for run A1 (table 1). Generalisation of the Reynolds number for systems with hyper-dissipation has been discussed by Spyksma, Magcalas & Campbell (2012).

Figure 2 displays the four MHD energy subfluxes, introduced in (2.6)–(2.7), and their sum. Also shown is $\mathcal{W}^\ell$, the resolved scale conversion of kinetic energy to magnetic energy (recall this does not represent energy transfer across $\ell$). The two panels present fluxes obtained through different filters, results shown in figure 2(a) correspond to Galerkin truncation and those shown in figure 2(a) to the Gaussian filter of (2.12). The data shown in figure 2(a,b) are qualitatively similar but with quantitative differences. Focusing on the similarities, we see that the Inertial term $\langle \varPi^{I,\ell} \rangle$ is relatively weak and is the only one to exhibit inverse transfer regions, in the intervals $1.5 \times 10^2 \lesssim k \eta_\alpha \lesssim 2 \times 10^{-1}$ and $4\times 10^{-3} \lesssim k \eta_\alpha \lesssim 6 \times 10^{-3}$. All the other subfluxes are associated with a direct cascade from the large to the small scales. The energy transfer from the momentum equation of (2.6) is almost entirely dominated by the Maxwell subflux ($-\Pi^{M,\ell}$) whose peak occurs in proximity to the forcing region. In contrast, the advection term from the induction equation ($\Pi^{A,\ell}$) is peaked at the small scales, close to the dissipative range. The conversion term, $\mathcal{W}^\ell$, is positive and increases monotonically with $k \eta_\alpha$. Particularly for the Fourier filter, it is roughly constant (i.e. scale independent) in the region $0.2 \lesssim k \eta_\alpha \lesssim 0.9$ where kinetic and magnetic subfluxes are in equipartition (Bian & Aluie 2019); see figure 11 in Appendix D. Also evident is a sudden increase of $\mathcal{W}^\ell$ in the dissipative range ($k\eta_\alpha \gtrsim 1$), due to hyperdiffusion, where the conversion rate saturates to $\varepsilon_b$ as already pointed out by Bian & Aluie (2019).

Figure 2. Terms contributing to the MHD filtered energy flux across scale $\ell$, along with the resolved scale conversion (kinetic to magnetic) term, as a function of the adimensional parameter $k \eta_\alpha = \pi \eta_\alpha / \ell$: (a) Fourier filter; (b) Gaussian filter. All terms are normalised by the mean total energy dissipation rate $\varepsilon = \varepsilon_u + \varepsilon_b$. The dashed horizontal line indicates the normalised magnetic dissipation rate .. The error bars, although not fully visible, indicate one standard error. See (2.6)–(2.7).

Turning to the differences between the two kinds of filtering, we observe that the bandwidth of the inertial range plateau is narrower for the Gaussian filter case, roughly $k \eta_\alpha \in [0.04,0.4]$ versus $k\eta_\alpha \in [0.012,0.7]$. In general, the Gaussian filtering peaks are of lower amplitude and a little less localised. Linked to this is a more gradual roll-off of the fluxes at high $k\eta_\alpha$ and a slower convergence of $\mathcal{W}^\ell$ to $\varepsilon_b$ with increasing $k$. These effects arise because Gaussian filtering at scale $\ell$ retains some effects from scales $\leq \ell$, unlike the situation for the sharp Galerkin truncation of the Fourier filter. For $\Pi^{D,\ell}$, there is also a qualitative difference, with the high-$k$ local minimum and maximum seen with the Fourier filter essentially absent when the Gaussian filter is used.

4 Analysis and discussion

4.1 Inertial flux and comparison with hydrodynamics

The exact decomposition of the Inertial term $\Pi^{I,\ell}$, (2.8), is

(4.1)\begin{equation} \varPi^{I,\ell} = \varPi^{I,\ell}_{s, SSS} + \varPi^{I,\ell}_{m,SSS} + \varPi^{I,\ell}_{s, S \varOmega \varOmega} + \varPi^{I,\ell}_{m,S \varOmega \varOmega} + \varPi^{I,\ell}_{m, S \varOmega S} , \end{equation}

where $\Pi^{I,\ell}_{s,S \Omega S}$ is not included as it is identically zero. This is the special case of (A3) where all the fields are the velocity and naturally it coincides with the original decomposition provided by Johnson (2020) for Navier–Stokes turbulence. It is thus of interest to investigate whether there are differences between the HD and MHD instances of (4.1), and how these might arise.

Figure 3 compares these two cases. For the hydrodynamic case, panel (a), there is a relatively clear and extended plateau for the total flux (and some of the subfluxes) corresponding to an inertial range. In contrast, the MHD case shown in panel (b) lacks such a plateau and the individual subfluxes are mostly much smaller than their hydrodynamic counterparts. Indeed, only $\langle \Pi^{I,\ell}_{m,SS\Omega} \rangle$ has values comparable to its hydrodynamic counterpart, albeit with a different functional form, being negative for almost all $k$. Intriguingly, this term is the only one that does not vanish point-wise in two-dimensional (2-D) turbulence, as discussed by Johnson (2021).² A depletion of the inertial flux in MHD turbulence has been observed by Alexakis (2013) and Yang et al. (2021) for configurations with large-scale electromagnetic forcing and by Offermans et al. (2018) for a saturated dynamo at lower Reynolds number. We have verified that the relation of Betchov (1956), $\langle \Pi^{I,\ell}_{s,SSS} \rangle = 3 \langle \Pi^{I,\ell}_{s,S\Omega \Omega} \rangle$, holds for both datasets, as it must.

Figure 3. Contributions to the Inertial energy flux $\langle \Pi^{I,\ell} \rangle$ for (a) the HD dataset H2 and (b) the MHD dataset A1, as a function of the (non-dimensionalised) reciprocal scale $\ell$, i.e. $\pi \eta_\alpha / \ell$. All fluxes are normalised by the mean total energy dissipation rate . for the dataset. Filled symbols correspond to single-scale contributions while hollow markers indicate multi-scale contributions. Error bars are for one standard error. Both panels indicate that $\langle \Pi^{I,\ell}_{m,SSS} \rangle \approx \langle \Pi^{I,\ell}_{m,S \Omega \Omega} \rangle$. In panel (b), the purple arrow locates $k \eta_\alpha = 5.4 \times 10^{-2}$ (equivalent to $kL_u \approx 14$), the value used for the probability density functions (p.d.f.s) shown in figures 4, 7 and 10.

Figure 4. Standardised p.d.f.s of the MHD Inertial subfluxes $\Pi^{I,\ell}_X$ at $k\eta_\alpha = 5.4 \times 10^{-2}$ (equivalent to $kL_u =14$) from dataset A1, where $X$ identifies the specific subflux. (a) Single-scale fluxes; (b) multi-scale fluxes.

Also of interest is that there appears to be an approximate multi-scale analogue of the Betchov relation, with $\langle \Pi^{I,\ell}_{m,SSS} \approx \langle \Pi^{I,\ell}_{m,S\Omega \Omega} \rangle$. For hydrodynamics, this was already noted by Johnson (2021) and justified by Yang et al. (2023). Evidently, this approximate degeneracy also holds in this MHD situation, although the smallness of the terms makes this difficult to appreciate from figure 3(b).

Having discussed mean fluxes, we now examine some statistical properties of their pointwise contributions. Figure 4 presents standardised probability density functions (p.d.f.s) of the MHD Inertial subfluxes at $k\eta_\alpha = 5.4 \times 10^{-2}$. It is evident that the distributions are strongly non-Gaussian and exhibit very wide tails, with fluctuations at tens of standard deviations. (As we shall see, this is a common characteristic for all the MHD energy fluxes and subfluxes.) Most importantly, we observe strong backscatter in all subfluxes. That is, the depletion of the mean Inertial flux results from cancellations between forward and inverse transfer events made possible by an increase in backscatter events. Similar observations for the total Inertial term have been reported by Offermans et al. (2018), where a comparison between the kinematic, nonlinear and saturated stages of the dynamo had been carried out. The variance, skewness and kurtosis for each p.d.f. are reported in table 2. Although the p.d.f. of $\Pi^{I,\ell}_{s,S\Omega \Omega}$, which corresponds to the scale-local vortex stretching, has more asymmetrical tails than those of $\Pi^{I,\ell}_{s,SSS}$, the skewness for the latter is nonetheless larger. The term $\Pi^{I,\ell}_{s,S \Omega S}$ does not show up in figure 4 because it is identically zero due to the symmetries of the tensors involved in the corresponding trace of (A3). In contrast, panel (b) shows that . is the most negatively skewed p.d.f. among the Inertial ones.

Table 2. Values of variance $\left(\sigma_X^Y\right)^2$, skewness $S_X^Y= \left\langle \left( \Pi^{Y,\ell}_X - \langle \Pi^{Y,\ell}_X\rangle \right)^3 \right\rangle/\left(\sigma_X^Y\right)^3$ and kurtosis $K_X^Y= \left\langle \left( \Pi^{Y,\ell}_X - \langle \Pi^{Y,\ell}_X\rangle \right)^4 \right\rangle/\left(\sigma_X^Y\right)^4$ for the subflux p.d.f.s shown in figures 4–10, where $X$ indicates the subflux identifier and $Y$ denotes the term identifier.

Figure 5. Two-dimensional sketch of current-sheet thinning and strain rate amplification by the latter. A current sheet, $\boldsymbol{\mathsf{J}}$, is stretched by (a) large-scale strain $\boldsymbol{\mathsf{S}}$, into a magnetic shear layer (b, red arrows). This induces a stretching of the magnetic flux tubes in the sheet. By conservation of magnetic flux, the magnetic field strength at the thereby generated smaller scales increases. That is, magnetic energy is transferred from large to small scales. The resulting magnetic shear layer has an associated magnetic strain rate field, $\textsf{$\mathit{\Sigma}$}$, whose principal axes (solid blue arrows) are at $.$ to those of the large-scale strain rate tensor (straight black arrows). As the magnetic shear will align with the extensional direction of the (velocity) strain rate tensor, this causes the fluid to be accelerated along these extensional directions and slowed down in the compressional directions, thereby generating a stronger rate-of-strain field across smaller scales, $\boldsymbol{\mathsf{S}'}$, indicated by the dashed arrows. The principal axes of the large-scale strain rate tensor are denoted by $\boldsymbol{\mathsf{e}}_1$ in one extensional direction and $\boldsymbol{\mathsf{e}}_3$ in one compressional direction and analogously for the small-scale strain rate.

In connection with the approximate degeneracy between $\Pi^{I,\ell}_{m,SSS}$ and $\Pi^{I,\ell}_{m,S\Omega \Omega}$ discussed above, figure 4(b) reveals that their p.d.f.s coincide up to events with a standardised probability density of $10^{-7}$, potentially indicating that the approximate identity is true not just on average but also for higher-order moments, although further analysis is required.

4.2 Maxwell flux

The decomposition of the energy flux associated with the Lorentz force, which we refer to as the Maxwell term, $\Pi^{M,\ell}$ of (2.9), contains an extra (single-scale) term with respect to the Inertial flux. Specifically, from (A3), we obtain

(4.2)\begin{equation} \varPi^{M,\ell} = \varPi^{M,\ell}_{s,S \varSigma \varSigma} + \varPi^{M,\ell}_{m,S \varSigma \varSigma} + \varPi^{M,\ell}_{s, S J J} + \varPi^{M,\ell}_{m, S J J} + \varPi^{M,\ell}_{s, S J \varSigma} + \varPi^{M,\ell}_{m, S J \varSigma} . \end{equation}

Terms of type $S \Sigma \Sigma$ can be associated with strain rate amplification by magnetic shear, while terms of type $SJJ$ correspond to current-filament stretching that is analogous to vortex stretching in HD. The last two terms are of type $SJ\Sigma$ and describe the back-reaction of the magnetic field on the flow or, more specifically, how the velocity strain rate is modified (typically amplified) in connection with a current-sheet thinning process. As we shall see, this is by far the dominant process. It proceeds as follows (figure 5). First, a current sheet is stretched by large-scale straining motions into a magnetic shear layer, in a process similar to vortex thinning in HD (Kraichnan 1976; Chen et al. 2006; Johnson 2021). This results in a stretching of the magnetic flux tubes in the sheet. By conservation of magnetic flux, the magnetic field strength at the thereby generated smaller scales must increase. That is, magnetic energy is transferred from large to small scales. (We will revisit this process in § 4.3 in the context of the inter-scale transfer of magnetic energy.) The magnetic rate-of-strain field associated with the resulting magnetic shear layer now accelerates fluid along its extensional directions and slows it down in the compressional directions, thereby generating a stronger rate-of-strain field across smaller scales.

It is instructive to consider the process in two dimensions, in analogy to the vortex thinning of 2-D HD. In the reference frame of the rate-of-strain tensor at scale $\ell$, the associated terms are

(4.3)\begin{gather} \varPi^{M,\ell}_{s,SJ\varSigma}(\boldsymbol{x}) = 2 \lambda_S^\ell(\boldsymbol{x}) \, \lambda^{\ell}_\varSigma (\boldsymbol{x}) \, \bar{j}^\ell (\boldsymbol{x}) \sin{2\psi(\boldsymbol{x})}, \end{gather}

(4.4)\begin{gather}\varPi^{M,\ell}_{m,SJ\varSigma}(\boldsymbol{x}) = 2 \int_0^{\ell^2} \,\mathrm{d}\theta \int_{\mathbb{R}^3} G^\phi(\boldsymbol{r}) \, \lambda_S^\ell(\boldsymbol{x}) \, \lambda^{\sqrt{\theta}}_\varSigma (\boldsymbol{x} + \boldsymbol{r}) \, \bar{j}^{\sqrt{\theta}} (\boldsymbol{x} + \boldsymbol{r}) \sin{2\psi(\boldsymbol{x} + \boldsymbol{r})} \,\mathrm{d} \boldsymbol{r} , \end{gather}

where $\pm \lambda_S$ and $\pm \lambda_\Sigma$ are the eigenvalues of the velocity and magnetic rate-of-strain tensors, respectively, $\psi$ the angle between the respective eigenvectors, and $j$ the out-of-plane component of the current density. As can be seen from these formulae, a maximum energy transfer occurs when the principal axes of the magnetic rate-of strain tensor have a $\pm 45^{\circ}$ angle to those of the velocity rate-of-strain tensor. Depending on the sign of the out-of-plane current density and $\sin(2\psi)$, (4.3) and (4.4) result in a direct or an inverse energy transfer. Figure 5 presents a schematic depiction of the process. A direct energy transfer occurs if the angle between the principal axes of velocity and magnetic strain-rate tensors is in the same rotational direction as the out-of-plane current. A similar, albeit less straightforward, assessment is possible in three dimensions, where

(4.5)\begin{equation} \varPi^{M,\ell}_{s,SJ\varSigma}(\boldsymbol{x}) = 2 \ell^2 \bar{S}_{ij}^\ell \bar{J}^\ell_{ik}\bar{\varSigma}^\ell_{kj} = 2 \ell^2 \sum_{i=1}^3 \sum_{j=1}^3 \lambda^\ell_i \mu_j^\ell\cos^2{\psi_{ij}} , \end{equation}

and similarly for the multi-scale term. Here, $\mu_j^\ell$ is the $j-$th eigenvalue of the symmetric part of the product matrix $\overline{{J}}^\ell_{ik}\overline{\Sigma}^\ell_{kj}$ (a contribution to the Maxwell SGS stress tensor) and $\psi_{ij}$ the angle between the $i-$th eigenvector of $\bar{\boldsymbol{\mathsf{S}}}^\ell$ and the $j-$th eigenvector of the symmetric part of $\overline{J}^\ell_{ik}\overline{\Sigma}^\ell_{kj}$. For a forward cascade of kinetic energy, $\langle \Pi^{M,\ell} \rangle > 0$, which implies that $\overline{S}^\ell_{ij} \, \overline{J}^\ell_{ik} \overline{\Sigma}^\ell_{kj}$ should be preferentially positive. Due to the presence of the cosine squared factor, this implies that the principal axes of the rate-of-strain tensor and those of the subscale stress must preferentially align, resulting in a stretching of the magnetic flux tubes along the extensional directions of the strain-rate tensor, as discussed above.

The subfluxes on the right-hand side of (4.2) and the total Maxwell flux are shown in figure 6, as a function of (non-dimensionalised) reciprocal $\ell$. We see immediately that the net energy transfer proceeds from large scales to small scales with the total Maxwell flux $\langle \Pi^{M,\ell} \rangle$ being the dominant energy subflux for MHD, carrying approximately $80\%$ of the total energy dissipation rate at its peak. At large scales, the major contribution is from the multi-scale term $\langle \Pi^{M,\ell}_{m,SJ\Sigma}\rangle$, switching to its single-scale partner, $\langle \Pi^{M,\ell}_{s,SJ\Sigma}\rangle$, as the dissipation scale is approached. All remaining terms in (4.2) are negligible. Summarising the mean Maxwell flux behaviour, we may say that the net kinetic energy transfer in MHD proceeds by the back-reaction of the magnetic field on the flow during the aforementioned current-sheet thinning process, while the contribution from current filament stretching and strain-amplification by magnetic shear are negligible.

Figure 6. Contributions to the Maxwell energy flux .. Data are from dataset A1 and normalised by the mean total energy dissipation rate .. Error bars indicate one standard error.

The p.d.f.s for the Maxwell energy fluxes are shown in figure 7. The predominance of $\langle \Pi^{M,\ell}_{m,SJ\Sigma} \rangle$ in the direct cascade can be also appreciated by examining its p.d.f., which is the most positively skewed among the multi-scale terms of figure 7; see also table 2 for p.d.f. moments. As can be seen from the data shown in figure 7, while all terms of type $S\Sigma \Sigma$ and $SJJ$ nearly vanish on average, fluctuations of approximately 60 standard deviations are not uncommon, and their multi-scale contributions show considerable backscatter. That is, the mean effects of current-filament stretching and strain-amplification by magnetic shear on the flow are negligible because of cancellations between pointwise forward and inverse transfers, and both these processes result in extreme backscatter events.

Figure 7. Standardised p.d.f.s for Maxwell subfluxes $\Pi_X^{M,\ell}$, at $k\eta_\alpha = 5.4 \times 10^{-2}$, from dataset A1, where $X$ represents the subflux identifier. (a) Single-scale fluxes; (b) multi-scale fluxes. Note that the p.d.f.s for $\Pi^{M,\ell}_{m,S\Sigma \Sigma}$ and $\Pi^{M,\ell}_{m,SJJ}$ are approximately coincident.

In Appendix B, we show that the averages of the subfluxes $\Pi_{m,S\Sigma \Sigma}^{M,\ell}$ and $\Pi_{m,S J J}^{M,\ell}$ can be connected via an exact Betchov-like relation that holds for all homogeneous flows,

(4.6)\begin{equation} \langle \varPi^{M,\ell}_{m,S\varSigma \varSigma} \rangle = \langle \varPi^{M,\ell}_{m,S J J }\rangle + 2\langle\varPi_{m,\varOmega \varSigma J}^{M,\ell} \rangle . \end{equation}

Note that the final term, $\Pi_{m,\Omega \Sigma J}^{M,\ell}$, does not appear in the Maxwell flux, (4.2), and our numerical results indicate $\Pi_{m,S\Sigma \Sigma}^{M,\ell} \approx \Pi_{m,SJJ}^{M,\ell}$, see figure 6. Physically, we may interpret this approximate identity as indicating that the net ‘strain-production’ by magnetic shear is almost equal to the net strain production by current-filament stretching. However, we stress again that these contributions to the interscale kinetic energy transfer are negligible. Further discussion on terms associated with strain production and current-filament stretching is provided in Appendix B.

4.3 Advection and dynamo fluxes

In this section, we focus on the decomposition of both the Advection term, $\Pi^{A,\ell}$, and the Dynamo term, $\Pi^{D,\ell}$, as defined in (2.10)–(2.11). As is well known, these two SGS fluxes share the same physical origin, namely the induced (fluctuation) electric field, and this is associated with certain symmetries and equivalences between the Advection and Dynamo subfluxes.

From the application of (A3), we find

(4.7)\begin{align} \varPi^{A,\ell} & = \varPi^{A,\ell}_{s, \varSigma \varSigma S} + \varPi^{A,\ell}_{m, \varSigma \varSigma S} + \varPi^{A,\ell}_{m, \varSigma \varSigma \varOmega} + \varPi^{A,\ell}_{s, \varSigma J S} + \varPi^{A,\ell}_{m, \varSigma J S} + \varPi^{A,\ell}_{s, \varSigma J \varOmega}\nonumber\\ & \quad + \varPi^{A,\ell}_{m, \varSigma J \varOmega} + \varPi^{A,\ell}_{s, J \varSigma S} + \varPi^{A,\ell}_{m, J \varSigma S} + \varPi^{A,\ell}_{s, J \varSigma \varOmega} + \varPi^{A,\ell}_{m, J \varSigma \varOmega} + \varPi^{A,\ell}_{s, J J S} + \varPi^{A,\ell}_{m, J J S} + \varPi^{A,\ell}_{m, J J\varOmega} , \end{align}

(4.8)\begin{align} \varPi^{D,\ell} & = \varPi^{D,\ell}_{s,\varSigma S \varSigma} + \varPi^{D,\ell}_{m,\varSigma S \varSigma} + \varPi^{D,\ell}_{s,\varSigma S J} + \varPi^{D,\ell}_{m,\varSigma S J} + \varPi^{D,\ell}_{m,\varSigma \varOmega \varSigma} + \varPi^{D,\ell}_{s,\varSigma \varOmega J} + \varPi^{D,\ell}_{m,\varSigma \varOmega J}\nonumber\\& \quad + \varPi^{D,\ell}_{s,J S \varSigma} + \varPi^{D,\ell}_{m,J S \varSigma} + \varPi^{D,\ell}_{s, J S J} + \varPi^{D,\ell}_{m, J S J} + \varPi^{D,\ell}_{s, J \varOmega \varSigma} + \varPi^{D,\ell}_{m, J \varOmega \varSigma} + \varPi^{D,\ell}_{m,J \varOmega J} , \end{align}

where we do not list terms that vanish identically, see Appendices A and F.

Figure 8 displays most of the Advection and Dynamo (sub)fluxes in separate panels. For clarity, only the subfluxes relevant to the discussion below and to the net energy flux are shown. The cyclic property of the trace can be used to show that some terms vanish identically and that each Advection subflux (both single-scale and multi-scale) is equal to (plus or minus) a partner Dynamo subflux; see the subfluxes expressions in Appendix F. For this reason, in figure 8, the following subflux pairs are not displayed since they are opposite in sign, $\Pi_{s,\Sigma \Sigma S}^{A,\ell}$ and $\Pi_{s,\Sigma S \Sigma}^{A,\ell}$ as well as $\Pi_{s,\Sigma J \Omega}^{A,\ell}$ and $\Pi_{s,\Sigma \Omega J}^{A,\ell}$ together with their multi-scale counterparts. Hence, these cancel pairwise and make no contribution to the net magnetic energy flux $\Pi^{A,\ell} + \Pi^{D,\ell}$; see Appendix F. In contrast, $\langle \Pi^{A,\ell}_{s,JJS} \rangle$ and $\Pi^{D,\ell}_{s,JSJ}$ and the related multi-scale terms, are equal and thus do contribute to the net flux.

Figure 8. Decomposed fluxes for the (a) Advection term $\langle \Pi^{A,\ell}\rangle$ and (b) Dynamo term $\langle \Pi^{D,\ell}\rangle$. Fluxes are normalised by the mean total energy dissipation rate $\varepsilon$ for dataset A1. In panel (a), $\langle \Pi^{A,\ell}_{s,\Sigma J S} \rangle$ and $\langle \Pi^{A,\ell}_{s,J \Sigma S} \rangle$ perfectly coincide. The error bars indicate one standard error.

Note, too, that the two Betchov relations (B6)–(B10) can provide another source of symmetry, or approximate symmetry.

The physical interpretation of the respective terms is very similar to what has been discussed for the Maxwell flux, except that now the effect of the flow on the magnetic field must be considered. Recall from § 4.2 that the Maxwell flux terms of type $SJ\Sigma$ correspond to the stretching and, by incompressibility, thinning of current sheets into magnetic shear layers, and that the back-reaction on the flow induced by this process is responsible for the bulk of the kinetic energy transfer to smaller scales. As we shall see, and as expected from the discussion of current-sheet thinning in § 4.2, this process also transfers most magnetic energy from large to small scales. However, in contrast to the Maxwell flux situation (dominated by a multi-scale term), here it is two of the Advection single-scale terms, $\Sigma J S$ and $J \Sigma S$, that carry most of the magnetic energy flux.

Focusing on the Advection term, from figure 8(a), we observe that the net flux is everywhere positive and peaked at smaller scales, roughly at the end of the inertial range. Recall that the Maxwell flux is peaked at larger scales (figure 6). Due to the cyclic property of the trace, the terms $\Pi^{A,\ell}_{s,\Sigma J S}$ and $\Pi^{A,\ell}_{s,J\Sigma S}$ are equal while, because of the symmetry of the tensors involved, the subfluxes $\Pi^{A,\ell}_{s,JJ\Omega}$ and $\Pi^{A,\ell}_{s,\Sigma \Sigma \Omega}$ together with $\Pi^{D,\ell}_{s,J\Omega J}$ and $\Pi^{D,\ell}_{s, \Sigma \Omega \Sigma}$ are identically zero. It is also apparent from figure 8(a) that the remaining Advection terms make negligible contributions to the net flux. We note that $\Sigma \Sigma S$ type terms correspond to magnetic shear amplification due to straining motions and those of type $JJS$ to current filament stretching – the analogue to vortex stretching in HD – as depicted in figure 9. Terms of type $\Sigma J \Omega$ encode bending of magnetic field lines into current filaments by rotational flow, i.e. a change in the magnetic-field geometry induced by vortical motion. Thus, the net Advection term is primarily due to single-scale contributions, being approximately equal to $2\, \Pi^{A,\ell}_{s,\Sigma J S}$, and it carries approximately $40\%$ of the total energy flux at its peak.

Figure 9. Sketch of current-filament stretching. A current filament $\boldsymbol{\mathsf{J}}$, with associated magnetic field $\mathbf{b}$, is (a) stretched by large-scale strain, $\boldsymbol{\mathsf{S}}$ with two contractional and one extensional direction into (b) a longer and thinner filament along the extensional direction. This induces a stretching of the magnetic flux tubes in the filament. By conservation of magnetic flux, the magnetic field strength at the thereby generated smaller scales increases (red arrows). That is, magnetic energy is transferred from large to small scales. This process is analogous to vortex stretching in hydrodynamic turbulence. The principal axes of the large-scale strain rate tensor are denoted by ${\boldsymbol{\mathsf{e}}}_1$ in the extensional direction and ${\boldsymbol{\mathsf{e}}}_2$ and ${\boldsymbol{\mathsf{e}}}_3$ in the contractional directions.

For the Dynamo term, $\Pi^{D,\ell}$, we observe that the net flux is almost flat in the inertial range, substantially positive definite for these scales (figure 8b) and responsible for approximately $15\%$ of the total energy flux. All subfluxes, except $\Pi^{D,\ell}_{s,\Sigma SJ}$ shown in yellow and $\Pi^{D,\ell}_{s,JS\Sigma}$ indicated by the filled purple symbols, are negligible. However, using the cyclic property once again, one can show that $\Pi^{D,\ell}_{s,\Sigma SJ} = -\Pi^{D,\ell}_{s,J S \Sigma}$; hence, these two terms cancel out and do not contribute to the net flux. In contrast to the Advection term, the major contributions to the net Dynamo term are from subleading single and multi-scale subfluxes of various types, with each contributing only 1–2 % to the total energy flux adding up to a total of approximately $15\%$ of the total energy flux. In summary, single-scale current-sheet thinning is the dominant process transferring magnetic energy across scales. It solely originates from the advective term in the induction equation.

Consider now the p.d.f.s. As a consequence of the symmetries between Dynamo and Advection subfluxes, only the Advection p.d.f.s are shown in figures 10(a) and 10(c), respectively for the single-scale terms and the multi-scale terms. It is striking that the p.d.f. of $\Pi^{A,\ell}_{s,\Sigma J \Omega}$, which describes the bending of magnetic field lines into current filaments by rotational flow, is by far the most strongly fluctuating with huge fluctuations of more than 100 standard deviations, also showing strong magnetic backscatter. In contrast, the aforementioned current-sheet thinning process, that is mostly responsible for the forward cascade of magnetic energy, shows weak magnetic backscatter. The other Advection subflux p.d.f.s, both single-scale and multi-scale, span a range comparable with those associated with the p.d.f.s for the Inertial and Maxwell fluxes. Moreover, all the multi-scale fluxes (Advection and Dynamo) have quite similar p.d.f.s, and thus so do the net multi-scale flux p.d.f.s (figure 10d). In the case of the net single-scale p.d.f.s, the Advection–Dynamo agreement is still good in the cores of the distributions, but there is a significant difference at larger negative fluctuations (figure 10b).

Figure 10. Standardised p.d.f.s for dataset A1 at $k \eta_\alpha = 5.4 \times 10^{-2}$. (a) Single-scale Advection subfluxes $\Pi^{A,\ell}_X$, where $X$ represents the subflux identifier, and (b) net single-scale fluxes for the Dynamo ($\Pi_s^{D,\ell}$) and Advection terms ($\Pi^{A,\ell}_s$). In the $x$-axis titles for panels (b) and (d), $Y=A$ or $D$, as appropriate. Note that $\Pi^{A,\ell}_{s,\Sigma J S} = \Pi^{A,\ell}_{s,J S \Sigma}$ (see figure 8). (c,d) As for panels (a,b) except for the multi-scale subfluxes of $\Pi^{A,\ell}$ and $\Pi^{D,\ell}$. Note that panels (c,d) have the same $x$-axis range as figures 4 and 7.

4.4 Total energy flux, .

In the preceding subsections, we analysed the four contributions to the incompressible MHD energy flux, finding that each of them may be reasonably well approximated using just some of the subfluxes. These approximations may now be assembled to give an approximate form for the mean total MHD energy flux. Specifically, we suggest that

(4.9)\begin{align} \left\langle \varPi^{\ell} \right\rangle & = \left\langle \varPi^{I, \ell} \right\rangle + \left\langle \varPi^{M, \ell} \right\rangle + \left\langle \varPi^{A, \ell} \right\rangle + \left\langle \varPi^{D, \ell} \right\rangle \approx{-} 2 \langle \varPi^{M,\ell}_{s,SJ\varSigma} \rangle - \langle \varPi^{M,\ell}_{m,SJ\varSigma} \rangle \end{align}

(4.10)\begin{align} & = 4 \ell^2 \left\langle \mathrm{Tr}\left\{ { (\bar{\boldsymbol{S}}^\ell)^t \bar{\boldsymbol{J}}^\ell (\bar{\boldsymbol{\varSigma}}^\ell)^t} \right\} \right\rangle \end{align}

(4.11)\begin{align} & \quad + 2\left\langle \int_0^{\ell^2} \,\mathrm{d}\theta \, \mathrm{Tr}\left\{ {\left( \bar{\boldsymbol{S}}^\ell \right)^t \left( \overline{ \bar{\boldsymbol{J}}^{\sqrt{\theta}} \left(\bar{\boldsymbol{\varSigma}}^{\sqrt{\theta}} \right)^{t} }^\phi - \overline{\bar{\boldsymbol{J}}^{\sqrt{\theta}}}^\phi\, \overline{ \left( \bar{\boldsymbol{\varSigma}}^{\sqrt{\theta}} \right)^{t} }^\phi \right)} \right\} \right\rangle \end{align}

is a suitable expression. Note that, after using cyclic properties of the trace, this is expressed only in terms of two of the Maxwell subfluxes, one single-scale and one multi-scale. These are discussed in § 4.2 and we remind the reader that all subflux definitions can be found in Appendix F. Interestingly, the terms on the right-hand side of (4.11) are part of a group of terms that remain non-zero in 2-D MHD; see Appendix C. We intend to explore this intriguing feature in future work.

Equation (4.11) demonstrates that the mean total energy flux in MHD turbulence is largely given by the stretching and thinning of current-sheets into magnetic shear layers by large-scale strain, resulting in a transfer of magnetic energy from large to small scales. In addition, there is a back-reaction of this process on the flow, whereby the ensuing magnetic strain-rate field accelerates fluid along its extensional directions and slows it down in the compressional directions, thereby generating a stronger strain-rate field across smaller scales, as shown schematically in figure 5.

4.5 Current-sheet thinning and magnetic reconnection

As a means of inter-scale energy transfer, the current-sheet thinning process can only proceed as described in an inertial range, i.e. at scales where Joule and viscous dissipation are negligible. As the dissipative scales are approached, Alfvén's theorem ceases to hold and magnetic reconnection can occur, with associated changes to the topology of the magnetic field. Interestingly, models of MHD reconnection are geometrically very similar to the inertial-range energy transfer process described here. For example in the Sweet–Parker model (Parker 1957; Sweet 1958), a current sheet is thinned by a flow that pushes magnetic field lines closer and closer together. Eventually, when the distance between the field lines approaches scales where Joule dissipation becomes important, the topological conservation of the magnetic field is broken and magnetic field-lines reconnect. That is, the continuation of the current-sheet thinning process described herein, and shown conceptually in figure 5, to smaller and smaller scales can lead naturally to magnetic reconnection.

5 Consequences for MHD subgrid-scale modelling

MHD LES modelling usually proceeds through variations on the Clark and Smagorinsky models, here shown for simplicity for the HD case,

(5.1)\begin{equation} \tau^{\rm CLARK}_{ij} = \ell^2 \partial_k \bar{u}^\ell_i \partial_k \bar{u}^\ell_j \end{equation}

for the Clark model (Clark et al. 1979; Pope 2000),³ where we note that this model corresponds to what we have herein called single-scale flux contributions (see the first term on the right-hand side of (2.15)), and the Smagorinsky model

(5.2)\begin{equation} \tau^{\rm SMAG}_{ij} ={-}2\nu_e\bar{S}^\ell_{ij} , \end{equation}

where . is the turbulent or eddy viscosity given by

(5.3)\begin{equation} \nu_e = (C_S \ell)^2 \sqrt{\bar{S}^\ell_{ij} \bar{S}^\ell_{ij}} , \end{equation}

with . being a free parameter, the Smagorinsky constant (Smagorinsky 1963). This is the case for incompressible MHD (Zhou & Vahala 1991; Müller & Carati 2002; Kessar et al. 2016). Moreover, adaptations have been made for the compressible case (Chernyshov, Karelsky & Petrosyan 2010; Grete et al. 2016; Vlaykov et al. 2016), for channel flow at high magnetic Reynolds number (Hamba & Tsuchiya 2010; Jadhav & Chandy 2021, 2023) and for extensions of MHD taking various levels of microphysics terms into account, such as Hall MHD (Miura, Araki & Hamba 2016) or Braginskii-extended (two-fluid) MHD (Miura, Hamba & Ito 2017), with the latter specifically focussed on the ballooning instability in stellarator devices. Most such approaches use the same SGS model for the inertial and Maxwell stresses based on velocity-field gradients with eddy viscosities involving either only the strain rate tensor or a weighted sum of the squared strain rate and the squared current, and similarly structured SGS models for the induction equation, with either heuristics or trial-and-error approaches to find suitable values for model constants. In what follows, we briefly discuss how the present results can be used to construct suitably structured SGS models for each SGS-stress in the MHD equations.

In terms of SGS modelling, our simulation-based results suggest that the Inertial terms can be neglected and a dissipative model, such as the Smagorinsky model, for the Maxwell stresses should suffice to capture the (leading-order) mean effects. In terms of fluctuations, we find that the observed mean Inertial flux depletion is caused by considerable backscatter in all Inertial subfluxes (see figure 4). This could suggest that a more sophisticated model would be required for the Inertial term. However, in 3-D HD turbulence, backscatter-free SGS models such as the standard static Smagorinsky closure perform well in capturing high-order statistics, that is, anomalous exponents and multifractal predictions for the correlation between velocity-field increments and SGS stresses (Linkmann, Buzzicotti & Biferale 2018). Furthermore, as the SGS stresses enter the filtered Navier–Stokes equation only through their divergence, this results in a degree of (gauge) freedom to determine model stress tensors that produce much less backscatter than those constructed using the standard definition (Vela-Martín 2022). In fact, backscatter can be traced back to spatial fluxes disconnected from the (scale-space) energy transfer and as such does not require modelling (Vela-Martín 2022).

For the magnetic energy transfers, we observe the net transfers are from large to small scales, suggesting again that dissipative models should suffice. Indeed, as multi-scale terms in our flux decomposition are negligible, a Clark-type model involving only the coupling between current and strain-rate tensors may work well for a nonlinear saturated non-helical (small-scale) dynamo. However, additional stabilising terms may be required as the Clark (or gradient) model is know to result in numerically unstable LES of a mixing layer in HD (Vreman, Geurts & Kuerten 1996, 1997) and for MHD (Müller & Carati 2002; Kessar et al. 2016), as discussed in further detail below. Similar to the Inertial term results, we observe considerable backscatter in the magnetic energy fluxes (see figure 10), and it remains to be seen if the aforementioned results by Linkmann et al. (2018) on the effect of SGS closures on high-order statistics carry over from HD to MHD.

For a non-helical saturated dynamo, as is the case here, Müller & Carati (2002) and Kessar et al. (2016) carried out MHD LES with Clark-type models constructed from the full velocity and magnetic-field gradients for the sum of the Reynolds and Maxwell SGS stress in the momentum equation and for the magnetic stresses, resulting in unstable simulations as in the HD case. Using a priori analyses of DNS data, Kessar et al. (2016) trace the instabilities back to the Clark terms transferring an insufficient amount of kinetic and magnetic energy to small scales, and to a production of backscatter. According to our analysis, for the momentum equation, the former effect is due to the Maxwell stress having a significant multi-scale component which is not captured in the Clark model. Concerning the latter, we also observe strong backscatter from the kinetic and magnetic single-scale subfluxes except for terms connected with current-sheet thinning, see figures 4, 7 and 10. Furthermore, a Clark model based on the full gradients by generalisation of (5.1) to MHD will introduce effects that are not present in the full MHD dynamics especially concerning the Maxwell and magnetic stresses, as all combinations of vorticity, strain, current and magnetic strain are included in the model and equally weighted,

(5.4)\begin{gather} \tau^{{\rm CLARK},M}_{ij} = \ell^2 \partial_k \bar{b}^\ell_i \partial_k \bar{b}^\ell_j = \ell^2 \left (\bar{\varSigma}^{\ell}_{ik} \bar{\varSigma}^{\ell}_{jk} + \bar{\varSigma}^{\ell}_{ik} \bar{J}^{\ell}_{jk} + \bar{J}^{\ell}_{ik}\bar{\varSigma}^{\ell}_{jk} + \bar{J}^{\ell}_{ik} \bar{J}^{\ell}_{jk}\right) , \end{gather}

(5.5)\begin{gather}\tau^{{\rm CLARK},A}_{ij} = \ell^2 \partial_k \bar{b}^\ell_i \partial_k \bar{u}^\ell_j = \ell^2 \left (\bar{\varSigma}^{\ell}_{ik} \bar{S}^\ell_{jk} + \bar{J}^{\ell}_{ik} \bar{S}^\ell_{jk} + \bar{\varSigma}^{\ell}_{ik}\bar{\varOmega}^{\ell}_{jk} + \bar{J}^{\ell}_{ik} \bar{\varOmega}^{\ell}_{jk}\right), \end{gather}

(5.6)\begin{gather}\tau^{{\rm CLARK},D}_{ij} = \ell^2 \partial_k \bar{u}^\ell_i \partial_k \bar{b}^\ell_j = \ell^2 \left (\bar{S}^\ell_{ik} \bar{\varSigma}^{\ell}_{jk} + \bar{\varOmega}^{\ell}_{ik} \bar{\varSigma}^{\ell}_{jk} + \bar{S}^\ell_{ik}\bar{J}^{\ell}_{jk} + \bar{\varOmega}^{\ell}_{ik} \bar{J}^{\ell}_{jk}\right) . \end{gather}

Our analysis, however, shows that only terms stemming from the coupling between current and strain-rate tensors are significant, with all remaining contributions to the net Maxwell and magnetic fluxes being negligible. Thus, retaining these contributions in an SGS model may substantially (and inappropriately) affect the small-scale structure of the flow and the magnetic field, and lead to an overestimation of backscatter. It remains to be seen if backscatter-related instabilities can be suppressed if only the current-sheet thinning (CST) terms,

(5.7)\begin{gather} \tau^{{\rm CST},M}_{ij} = \ell^2 \left ( \bar{\varSigma}^{\ell}_{ik} \bar{J}^{\ell}_{jk} + \bar{J}^{\ell}_{ik}\bar{\varSigma}^{\ell}_{jk} \right) , \end{gather}

(5.8)\begin{gather}\tau^{{\rm CST},A}_{ij} = \ell^2 \left (\bar{\varSigma}^{\ell}_{ik} \bar{S}^\ell_{jk} + \bar{J}^{\ell}_{ik} \bar{S}^\ell_{jk} \right) , \end{gather}

(5.9)\begin{gather}\tau^{{\rm CST},D}_{ij} = \ell^2 \left (\bar{S}^\ell_{ik} \bar{\varSigma}^{\ell}_{jk} + \bar{S}^\ell_{ik}\bar{J}^{\ell}_{jk}\right) , \end{gather}

are retained. That is, both magnetic strain and current density tensors must be included in SGS models as recently discussed by Alexakis & Chibbaro (2022).

A further challenge for LES modelling of MHD turbulence is to accurately capture the transfer of magnetic to kinetic energy (Haugen & Brandenburg 2006) and vice versa. Using standard Smagorinsky models for both the momentum and induction equations, Haugen & Brandenburg (2006) showed that the Smagorinsky constant can be fine-tuned to obtain a good agreement between filtered DNS and LES for kinetic energy spectra, but only at the expense of a strong suppression of the (nonlinear) dynamo. Moreover, as discussed by Offermans et al. (2018), the resolved-scale conversion term must either be fully accounted for in LES, resulting in the need to resolve all scales where this term is active, or in the present case of a saturated dynamo, a model including an extra term accounting for the under-resolved dynamo effect must be provided.

In this paper, we have only considered the no mean magnetic field situation. The presence of a strong background magnetic field is likely to require an SGS modelling approach that differs from those just discussed, as the ensuing anisotropy and two-dimensionalisation of magnetic- and velocity-field fluctuations may result in partial inverse fluxes. We will report results on configurations with strong background magnetic fields in due course. Similarly, the large-scale (helical) dynamo requires an investigation in its own right, and the magnetic-field growth at large scales is likely to require a different type of SGS modelling approach, as Kessar et al. (2016) report that the Clark and even a standard static Smagorinsky model result in unstable LES.

The method discussed here can be readily extended to Hall and two-fluid MHD and other fluid models for plasma turbulence. For instance, the applicability of the Smagorinsky closure to Hall MHD has been assessed by an a priori analysis using sharp filtering (Miura & Araki 2012) prior to the deployment of said closure (Miura et al. 2016). An a priori analysis and decomposition of the Hall flux in analogy to the results presented here could lead to a better understanding of the physics of the interscale magnetic energy transfer induced by the Hall effect and thus to a refinement of Hall-MHD SGS models. Finally, we point out that a new LES method, so-called physics-inspired coarsening (PIC), has been devised recently (Johnson 2022). In the homogeneous case, this approach reduces to Gaussian filtering and the representation of SGS stresses in terms of field gradients as generalised herein. In PIC, the velocity field advanced in LES is formally obtained by artificial viscous smoothing, with the required pseudo-diffusion being introduced through an auxiliary Stokes equation. This approach may be generalisable to MHD and more complex fluid models applicable to plasma turbulence.

6 Conclusions

Generalising a method introduced by Johnson (2020) for Navier–Stokes turbulence, we have presented a general analytical method for obtaining exact forms for inter-scale fluxes in advection–diffusion equations through products of vector-field gradients, and applied it to kinetic and magnetic energy fluxes in homogeneous MHD turbulence. The aim was to provide expressions for subfluxes that are physically interpretable in terms of the action of the magnetic field on the flow and vice versa. A quantification thereof is of interest for the fundamental understanding of cascade processes in MHD turbulence, and also provides guidance as to what physics needs to be captured in subgrid-scale models and how such models should be constructed so that they preserve, at least approximately, empirical features of the mean energy fluxes and their fluctuations.

In MHD, scale-space energy fluxes are defined as the contraction of velocity- or magnetic-field gradients with the appropriate subgrid-scale stresses. Rewriting these in terms of symmetric and antisymmetric components of field-gradients tensors yields terms with clear physical meanings. For example, strain and vorticity in the case of the velocity field, and current and magnetic strain/shear for the magnetic field. Expressing the MHD SGS stresses in terms of vorticity, rate-of-strain, current and magnetic rate-of-strain results in an exact decomposition of magnetic and kinetic energy fluxes in terms of interactions between the symmetric and antisymmetric components of velocity- and magnetic-field gradients.

The kinetic energy flux comprises two terms, the Inertial flux (as in hydrodynamics) and a flux term associated with the action of the Lorentz force on the flow. The former is decomposed into terms associated with vortex stretching, strain self-amplification and strain-vorticity alignment. A term-by-term comparison between the Inertial fluxes in HD and MHD turbulence shows that all Inertial subfluxes are depleted due to cancellations between forward scatter and backscatter events, and are indeed almost negligible in MHD turbulence. That is, the physics of the kinetic energy cascade is very different in statistically steady MHD turbulence as compared with HD turbulence, as vortex stretching and strain self-amplification have, on average, very little effect. In MHD turbulence, almost all kinetic energy is transferred downscale by a current-sheet thinning process: in regions of large strain, current sheets are stretched by large-scale straining motion into regions of magnetic shear. This magnetic shear in turn drives extensional flows at smaller scales. The magnetic energy is mainly transferred from large to small scales by the aforementioned current-sheet thinning in regions of high strain, while the mean contributions from current-filament stretching – the analogue to vortex stretching – and bending of magnetic field lines in high magnetic strain regions into current filament by vortical motion are almost negligible. The latter effect, which results in a change in the magnetic field geometry at small scales is associated with strong magnetic backscatter. A decomposition into single- and multi-scale components of the subfluxes shows that the mean kinetic energy flux induced by the back-reaction of current-sheet thinning on the flow has a strong multi-scale component especially at large scales where the magnetic field is weak, while the multi-scale component of the magnetic energy flux is almost negligible. We consistently observe that the multi-scale components of the respective flux terms fluctuate less than the single-scale components.

Finally, we note that the method can be further expanded in various directions. Within MHD turbulence, for instance, decompositions of the magnetic helicity and cross-helicity fluxes could identify the physics driving these respective cascades with potential implication for selective decay for decaying MHD turbulence. Kinetic energy transfer across scales and its conversion to magnetic energy, and vice versa, depend on the value of the magnetic Prandtl number (Brandenburg & Rempel 2019). The methodology introduced here could be used to quantify the .-dependence of the different physical processes involved in the energy cascade. Another possibility would be to include temperature fluctuations. An extension or application to compressible flows would be of interest especially for astrophysical plasmas. As flux terms associated with any advective nonlinearity can be analysed by this method, a decomposition of Hall MHD and of fluid models of ion- or electron-temperature-gradient turbulence (Ivanov et al. 2020; Adkins et al. 2022; Ivanov, Schekochihin & Dorland 2022) may be of interest for the magnetic confinement fusion community.