# Non-linear energy transfers in accretion discs MRI turbulence. I. Net vertical field case

G. Lesur, Pierre-Yves Longaretti

► To cite this version:

G. Lesur, Pierre-Yves Longaretti. Non-linear energy transfers in accretion discs MRI turbulence. I. Net vertical field case. Astronomy and Astrophysics - A&A, 2011, 528, 10.1051/0004-6361/201015740. insu-03625203

HAL Id: insu-03625203

https://insu.hal.science/insu-03625203

Submitted on 30 Mar 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Non-linear energy transfers in accretion discs MRI turbulence

## I. Net vertical field case

G. Lesur¹,² and P.-Y. Longaretti²

¹ Department of Applied Mathematics and Theoretical Physics, University of Cambridge, Centre for Mathematical Sciences, Wilberforce Road, Cambridge CB3 0WA, UK

² UJF-Grenoble 1 / CNRS-INSU, Institut de Planétologie et d'Astrophysique de Grenoble (IPAG) UMR 5274, 38041 Grenoble, France
e-mail: geoffroy.lesur@obs.ujf-grenoble.fr

Received 10 September 2010 / Accepted 30 November 2010

### ABSTRACT

The magnetorotational instability (MRI) is believed to be responsible for most of the angular momentum transport in accretion discs. However, molecular dissipation processes may drastically change the efficiency of MRI turbulence in realistic astrophysical situations. The physical origin of this dependency is still poorly understood as linear and quasi linear theories fail to explain it. In this paper, we look for the link between molecular dissipation processes and MRI transport of angular momentum in unstratified shearing box simulations, including a mean vertical field. We show that magnetic helicity is unimportant in the model we consider. We perform a spectral analysis on the simulations tracking energy exchanges in spectral space when turbulence is fully developed. We find that the energy exchanges are essentially direct (from large to small scale) whereas some non-linear interactions appear to be non-local in spectral space. We speculate that these non-local interactions are responsible for the correlation between turbulent transport and molecular dissipation. We argue that this correlation should then disappear when a significant scale separation is achieved, and we discuss several methods by which one can test this hypothesis.

**Key words.** accretion, accretion disks – instabilities – protoplanetary disks – turbulence

## 1. Introduction

The transport of angular momentum in astrophysical discs is a central problem in accretion theory. To explain disc lifetime and accretion rate, it is often assumed that these objects are turbulent. Turbulence is then included in global models using a turbulent viscosity prescription as in the $\alpha$ disc model (Shakura & Sunyaev 1973).

The origin of this turbulence has been the subject of many debates over the past decades. It is now generally assumed that the magnetorotational instability, or MRI for short (Velikhov 1959; Chandrasekhar 1960; Balbus & Hawley 1991), is responsible for disc turbulence, although hydrodynamic processes might also be at work (Lesur & Ogilvie 2010; Lesur & Papaloizou 2010). Although MRI generated turbulence is generally efficient at transporting angular momentum (Hawley et al. 1995), recent results have shown a strong sensitivity for MRI turbulence on small-scale dissipation processes (Lesur & Longaretti 2007; Fromang et al. 2007), and in particular on the magnetic Prandtl number $Pm$ (ratio of microscopic viscosity to resistivity). This effect, called the $\alpha - Pm$ correlation, casts doubts on the actual efficiency of the MRI in realistic situations since $Pm$ can vary by several orders of magnitude in discs (Balbus & Henri 2008). Several attempts have been made to explain this correlation, whether from the linear theory of dissipative MRI modes (Pessah & Chan 2008) or from the quasi-linear parasitic modes theory (Pessah & Goodman 2009). However, these approaches were shown to be unsuccessful when compared to high Reynolds number simulations (Longaretti & Lesur 2010). Instead, Longaretti & Lesur (2010) suggested that the

$\alpha - Pm$ correlation could be due to the nature of the MHD cascade in MRI generated turbulence, in which one might expect inverse cascades and/or non-local interaction in spectral space. This kind of process would allow for direct communication between the injection scales (transport scales) and the largest dissipation scale (either resistive or viscous).

The purpose of the present work is to investigate some of the conjectures presented by Longaretti & Lesur (2010) regarding the nature of the MHD turbulent cascade in accretion discs. To this end, we consider several of the high-resolution simulations presented by Longaretti & Lesur (2010), and we analyse the energy exchanges in spectral space. This paper is organised as follows. We describe our model, equations, and the spectral analysis we use in Sect. 2. Section 3 is the core of this paper and discusses our numerical results. The implications of these results are presented in the final section.

## 2. Model and spectral analysis

### 2.1. Equations

In the following, we adopt the shearing box model that accurately represents the local physics of an accretion disc (Hawley et al. 1995; Balbus 2003; Regev & Umurhan 2008). The adopted coordinate system is such that $x = (r - r_0)$ and $y = r_0\phi$, where $r_0$ is the fiducial radius of the shearing box in the disc and $\phi$ the azimuthal coordinate in the rotating frame. The velocity can be decomposed as a mean velocity plus a fluctuating part $U = -q\Omega xe_y + v$, where $\Omega$ is the local rotation rate and $q = 3/2$ for a Keplerian disc rotation profile. As a simplification,
we assume the flow is incompressible, which corresponds to the “small shearing box limit” of Umurhan & Regev (2004). The equations of motion then read as

$$ \partial_t \mathbf{v} = -\mathbf{v} \cdot \nabla \mathbf{v} - \nabla \Pi + \mathbf{B} \cdot \nabla \mathbf{B} + q\Omega x \partial_y \mathbf{v} \\ +2\Omega v_y \mathbf{e}_x - (2-q)\Omega v_x \mathbf{e}_y + v \nabla^2 \mathbf{v}, \quad (1) $$

$$ \partial_t \mathbf{B} = -\mathbf{v} \cdot \nabla \mathbf{B} + \mathbf{B} \cdot \nabla \mathbf{v} \\ +q\Omega x \partial_y \mathbf{B} - q\Omega B_x \mathbf{e}_y + \eta \nabla^2 \mathbf{B}, \quad (2) $$

$$ \nabla \cdot \mathbf{v} = 0, $$

$$ \nabla \cdot \mathbf{B} = 0, $$

where $\Pi$ is the total pressure, $v$ the viscosity, and $\eta$ the ohmic resistivity. In the following, we impose a mean vertical field $B_0$, which will be conserved during the evolution of the flow thanks to the shearing-sheet boundary conditions. It should be noted that the magnetic field strength is expressed in Alfvén speed for simplicity.

Several dimensionless numbers characterise the equations of motions. In this paper, we use:

- the amplitude of the imposed mean vertical field measured by

$$ \beta = \frac{(q\Omega)^2 L_z^2}{B_0^2} \quad (3) $$

where $L_z$ is the vertical box size. This definition mimics the usual plasma $\beta$ in vertically stratified discs obeying the vertical hydrostatic equilibrium constraint $c_s \sim \Omega L_z$;

- the viscous Elsasser number:

$$ \Lambda_v = \frac{B_0^2}{\Omega v}, \quad (4) $$

which is related to the Reynolds number used in Longaretti & Lesur (2010) by $\Lambda_v = qRe/\beta$;

- the resistive Elsasser number:

$$ \Lambda_{\eta} = \frac{B_0^2}{\Omega \eta}, \quad (5) $$

connected to the magnetic Reynolds number by a similar relation;

- the magnetic Prandtl number:

$$ Pm = \frac{v}{\eta} = \frac{\Lambda_{\eta}}{\Lambda_{v}}, \quad (6) $$

which compares the amount of viscous and resistive dissipation.

## 2.2. Fourier transform in sheared flows

It is convenient to introduce the shearing frame $(x', y', z')$:

$$ x = x', $$

$$ y = y' - q\Omega t x', $$

$$ z = z'. $$

Writing the equations of motions in the sheared frame allows us to eliminate the explicit spatial dependency:

$$ \partial_t \mathbf{v} = -\mathbf{v} \cdot \nabla \mathbf{v} - \nabla \Pi + \mathbf{B} \cdot \nabla \mathbf{B} \\ +2\Omega v_y \mathbf{e}_x - (2-q)\Omega v_x \mathbf{e}_y + v \nabla^2 \mathbf{v}, \quad (7) $$

$$ \partial_t \mathbf{B} = -\mathbf{v} \cdot \nabla \mathbf{B} + \mathbf{B} \cdot \nabla \mathbf{v} \\ -q\Omega B_x \mathbf{e}_y + \eta \nabla^2 \mathbf{B}, \quad (8) $$

where $\mathbf{v}$ and $\mathbf{B}$ are now assumed to be functions of $x'$ so that the nabla operator expression becomes

$$ \nabla = \mathbf{e}_x (\partial_{x'} + q\Omega t \partial_{y'}) + \mathbf{e}_y \partial_{y'} + \mathbf{e}_z \partial_{z'}. \quad (9) $$

In the sheared frame, the shearing sheet boundary conditions make every physical quantity $X(x')$ periodic so that $X$ can be expanded in Fourier series:

$$ X(x') = \sum_{k'} X_{k'} \exp(i\mathbf{k}' \cdot \mathbf{x}') = \sum_{k(t)} X_{k(t)} \exp(i\mathbf{k}(t) \cdot \mathbf{x}). \quad (10) $$

This relation defines the time dependent unsheared wave vectors:

$$ k_x = k_{x'} + q\Omega k_{y'} t, \quad (11) $$

$$ k_y = k_{y'}, \quad (12) $$

$$ k_z = k_{z'}. \quad (13) $$

As expected, the application of the $\nabla$ operator to $X(x')$ corresponds to a multiplication of its Fourier components by $i\mathbf{k}(t)$.

This definition of the unsheared wave vectors with the Fourier decomposition (10) is usually referred to as a “shearing wave” decomposition of a sheared flow. It was first used by Lord Kelvin to study the stability of sheared flows (Thomson 1887) and later in the astrophysical context by Goldreich & Lynden-Bell (1965) for spiral arms of galaxies.

## 2.3. Shell filter decomposition

Following Frisch (1995) and Alexakis et al. (2007), one defines shell-filtered quantities in the unsheared Fourier space. At any given time, a series of linearly spaced shell sizes $K$ is defined from $K_1$ to $K_{max}$; by construction $\delta K = K_j - K_{j-1}$ is the shell width (for any $j$). We define the shell-filtered field $X_{K_j}$ in shell $K_j$ by

$$ X_{K_j} = \sum_{K_j - \delta K/2 < k \le K_j + \delta K/2} X_k \exp(i\mathbf{k} \cdot \mathbf{x}). \quad (14) $$

It should be noted that, since $\mathbf{k}$ depends on time, the exact number and the distribution of the modes entering the above formula for any given $K$ might vary in time, adding an extra complication compared to the homogeneous case of Alexakis et al. (2007). This point is discussed in the appendices.

## 2.4. Energy transfer equations

The transfers we are interested in relate to the equations of the kinetic and magnetic energies. These shell-restricted, box-averaged equations involve a number of transfer functions that are introduced along with the related equations. We follow here the logic of Alexakis et al. (2007) and extend it to shear flows. The analysis of transfers gives indications about the locality of interactions in Fourier space.

Because of the incompressibility condition, energy transfers (but also stress and magnetic helicity) involve only at most the product of three components of the velocity and magnetic fields and their derivatives. Therefore, couplings in Fourier space only depend on triads of wave vectors, noted $\mathbf{k}_1, \mathbf{k}_2, \mathbf{k}_3$. The closing condition $(\mathbf{k}_1 + \mathbf{k}_2 = \mathbf{k}_3)$ furthermore imposes that at least two of the wave vectors $\mathbf{k}_i$ are of the same magnitude; the third one can be either much smaller (implying non-local couplings through large scales) or of the same magnitude as the other two (implying locality of couplings in Fourier space).
Using the equations of motion in the sheared frame (7–8), one can derive the equation for the shell-filtered energy density:

$$ \partial_t E_K = \sum_Q [T_{vv}(Q, K) + T_{bv}(Q, K)] + S_{v,K} + q\Omega I_{v,K} - vD_{v,K} \quad (15) $$

$$ \partial_t M_K = \sum_Q [T_{bb}(Q, K) + T_{vb}(Q, K)] + S_{b,K} + q\Omega I_{b,K} - \eta D_{b,K} \quad (16) $$

where $E_K = \langle v_K^2/2 \rangle$ and $M_K = \langle B_K^2/2 \rangle$. In the above expression, we have defined the transfer functions by

$$ \begin{aligned} T_{vv}(Q, K) &= -\langle \boldsymbol{v}_K \cdot (\boldsymbol{v} \cdot \boldsymbol{\nabla}) \boldsymbol{v}_Q \rangle, \\ T_{bb}(Q, K) &= -\langle \boldsymbol{B}_K \cdot (\boldsymbol{v} \cdot \boldsymbol{\nabla}) \boldsymbol{B}_Q \rangle, \\ T_{bv}(Q, K) &= \langle \boldsymbol{v}_K \cdot (\boldsymbol{B} \cdot \boldsymbol{\nabla}) \boldsymbol{B}_Q \rangle, \\ T_{vb}(Q, K) &= \langle \boldsymbol{B}_K \cdot (\boldsymbol{B} \cdot \boldsymbol{\nabla}) \boldsymbol{v}_Q \rangle, \end{aligned} $$

where $\langle \cdot \rangle$ denotes an spatial average on the shearing box volume, and $T_{ij}(Q, K)$ represents the transfer from energy “i” (kinetic or magnetic) from shell $Q$ to energy “j” (kinetic or magnetic) in shell $K$. We note that $T_{ij}(Q, K) = -T_{ji}(K, Q)$, so that whatever is taken from one shell of one type is totally transferred to the other shell. Similarly, $T_{ii}(K, K) = 0$: there is no effective transfer from a shell to itself, as should be. This justifies the identification of these quantities as shell-to-shell energy transfer functions; in fact, the third member of a triad is only a relay in effective energy exchanges between shells $Q$ and $K$ (for a more detailed discussion, see Verma 2004). In these expressions, $\boldsymbol{v}_K \cdot \boldsymbol{v}_K = \boldsymbol{v}_K \cdot \boldsymbol{v}$ has been used, as well as $\boldsymbol{v} = \sum_Q \boldsymbol{v}_Q$.

The next terms in the shell energy budget involve energy transfers due to the mean shear $S_{v,K}$ and $S_{b,K}$. These terms are singular in time as they correspond to energy fluctuations due to wave entering or leaving the shell $K$ as $k(t)$ evolves. They can be formally defined by

$$ \begin{aligned} S_{x,K} &= \sum_{k'} \frac{X_{k'}^* X_{k'}}{2} \delta(t - t_{k'}) \epsilon_{k'} \\ &= \sum_{k'} \frac{q \Omega k_y k_x(t) X_{k'}^* X_{k'}}{|k(t)|} \\ &\quad \times [\delta(|k(t)| - K + \delta K/2) - \delta(-|k(t)| + K + \delta K/2)], \end{aligned} $$

where $t_{k'}$ is the instant when the wave $k'$ enters or exits the shell $K$ and $\epsilon_{k'} = \pm 1$ for an entering/exiting wave (see Appendices A, B).

The remaining terms in the shell kinetic and magnetic energy budgets,

$$ \begin{aligned} I_{v,K} &= \langle v_{x,K} v_{y,K} \rangle, \\ I_{b,K} &= -\langle B_{x,K} B_{y,K} \rangle, \\ D_{v,K} &= \langle (\boldsymbol{\nabla} \times \boldsymbol{v}_K)^2 \rangle, \\ D_{b,K} &= \langle (\boldsymbol{\nabla} \times \boldsymbol{B}_K)^2 \rangle, \end{aligned} $$

represent the energy injection through the shear and dissipation through viscosity and resistivity in shell $K$.

## 2.5. Energy fluxes in spectral space

Using the transfer function defined above, it is possible to introduce energy flux in Fourier space. In this work, we will use the

following fluxes:

$$ \mathcal{F}_v(K_0, t) = \sum_{K=K_0}^{K_{\max}} \sum_Q T_{vv}(Q, K) \quad (17) $$

$$ \mathcal{F}_b(K_0, t) = \sum_{K=K_0}^{K_{\max}} \sum_Q T_{bb}(Q, K) \quad (18) $$

$$ \mathcal{F}_x(K_0, t) = \sum_{K=K_0}^{K_{\max}} \sum_Q T_{vb}(Q, K) + T_{bv}(Q, K) \quad (19) $$

$$ \mathcal{F}_{s,v}(K, t) = \sum_{k'} \frac{q \Omega k_y k_x(t) \mathbf{V}_{k'}^* \cdot \mathbf{V}_{k'}}{|k(t)|} \delta(|k(t)| - K) \quad (20) $$

$$ \mathcal{F}_{s,b}(K, t) = \sum_{k'} \frac{q \Omega k_y k_x(t) \mathbf{B}_{k'}^* \cdot \mathbf{B}_{k'}}{|k(t)|} \delta(|k(t)| - K) \quad (21) $$

where $\mathcal{F}_v$, $\mathcal{F}_b$, $\mathcal{F}_x$, and $\mathcal{F}_s$ are respectively the kinetic, magnetic, exchange, and shear fluxes. The shear fluxes are evaluated in $K = K_0 + \delta K/2$. Each of these fluxes computes the amount of the energy transferred from shells $K \le K_0$ to shells $K > K_0$ (i.e. the flux of energy “through” the outer boundary of shell $K_0$). The kinetic (respectively magnetic) fluxes compute transfers of kinetic to kinetic (respectively magnetic to magnetic) energy, whereas the exchange flux is a flux of total energy (magnetic plus kinetic) in which magnetic and kinetic energy are constantly transformed into one another. Shear fluxes are singular in time as they are non zero only when a wave enters or leave shells $K > K_0$. It should be noted that shear fluxes are statistically non zero only for anisotropic turbulence, as the amplitude of any mode ($k_x, k_y, k_z$) should be statistically equal to the amplitude of the mode ($-k_x, k_y, k_z$) in isotropic turbulence.

These fluxes allow one to check the direction of the energy transfers due to the non-linear terms. Indeed, a direct energy transfer (large to small scale) implies a positive flux with the above definition, whereas an inverse cascade can be characterised by a negative flux.

## 2.6. Numerical method

Equations (1) and (2) are solved using the Snoopy code. Snoopy is a 3D spectral (Fourier) method based on the shearing-wave decomposition (11)–(13). The Fourier transforms are evaluated using the FFTW 3 library, with both MPI and OpenMP parallelisation techniques. Non-linear terms are computed using a pseudo-spectral algorithm (Canuto et al. 1988) and antialiasing is enforced using the “3/2” rule. Time integration is performed by a third-order, low-storage Runge-Kutta scheme for non-linear terms, whereas an implicit scheme is used for viscous and resistive terms. This spectral scheme uses a periodic remap algorithm in order to continuously follow the smallest wave number of the system in the sheared frame $|k(t)| < k_{\max}$ (see Umurhan & Regev 2004, Appendix C for a complete description of the periodic remap algorithm). The periodic remap method used in this code is different from the continuous remap method used by Lithwick (2007). Our main motivation in implementing a periodic remap is the possibility of using the 3/2 antialiasing rule and power of 2 grid sizes<fnref n="1" /> for which Fourier transform and parallelization methods are more efficient. This code or its variant

<footnote n="1">¹ If the number of point in one direction is a multiple of 2, waves at the Nyquist frequency do not have any properly defined phase. This is not a problem for classical spectral methods or for the periodic remap method, since the Nyquist frequency is either in the dissipation range or in the antialiasing dump zone. However, when using a continuous</footnote>
<figure><img src="image_2.png" /><figcaption>Fig. 1. Energy spectrum at $Pm = 0.0625$ (left) and $Pm = 0.25$ (right). In the $Pm = 0.25$ case, a power-law spectrum is observed for the kinetic energy corresponding to a $k^{-3/2}$ spectrum.</figcaption></figure>

has been used in several context, including the MRI (Longaretti & Lesur 2010) and the subcritical baroclinic instability (Lesur & Papalozou 2010). It is available for download on the author's website.

## 3. Results

### 3.1. Simulations parameters and averaging procedure

The spectra and transfers presented in this section were all derived from two simulations of an MRI saturated state. These runs correspond to the $Pm = 0.25$ and $Pm = 0.0625$ high-resolution runs discussed in Longaretti & Lesur (2010). Both runs have a resolution<fnref n="2" /> of $N_x \times N_y \times N_z = 768 \times 384 \times 192$ with a box aspect-ratio $L_x \times L_y \times L_z = 4 \times 4 \times 1$. We imposed a mean vertical field in the box with $\beta = 10^3$ and $\Lambda_v = 30$ (Re = $2 \times 10^4$). Each simulation was integrated for 50 orbits starting from random noise, and the spectra were averaged from the last 40 orbits to remove any influence of the initial conditions. The two simulations considered in this section only differ by their ohmic resistivity, where the $Pm = 0.25$ run has $\Lambda_\eta = 7.5$ ($Rm = 5000$) and the $Pm = 0.0625$ run has $\Lambda_\eta = 1.87$ ($Rm = 1250$).

The statistical average of any quantity $X$ of interest $\langle X \rangle_{\text{stat}}$ should in principle be computed on different realisations but are evaluated in practice as usual via an ergodic hypothesis:

$$ \langle X \rangle_{\text{stat}} = \lim_{T \to \infty} \frac{1}{T} \int_0^T X(t) dt \approx \frac{1}{T} \int_0^T X(t) dt \approx \sum_i X(T_i), \quad (22) $$

where $0 \le T_i \le T$ are a sufficiently large number of instants of flow snapshots.

The spectra were averaged in the spherical shells introduced in (14). The shells were defined so that $K_n = 2\pi n / L_z$ and $\delta K = 2\pi / L_z$. This means that some power is present in the shell $K = 0$, as it contains large-scale horizontal waves with no vertical structure. The shell-integrated spectra and transfers obtained by this procedure are then averaged in time over 40 instantaneous

snapshots (1 snapshot per orbit). For simplicity, we renormalized the wavevectors $K$ so that $K' = K/2\pi$ on all the plots in this section. Shells $K > 32$ are incomplete in the $y$ direction since the resolution per scale height is lower in that direction. This is not a problem since these shells are in the dissipative range, and high $k_y$ modes are weaker than the equivalent high $k_x/k_z$ modes due to the anisotropy of MRI turbulence (see Sect. 3.2). With our procedure, one can reconstruct the box-averaged quantities by summing the spectra over the integers $K$.

The shear transfer terms $F_{s,v}$ and $F_{s,b}$ are computed in a special way. Indeed, one cannot compute $F_s$ for a given shell and snapshot time numerically because of the $\delta$ functions. Instead, we introduce a shell-averaged flux:

$$ F_s^c(K_0, t) = \frac{1}{\delta K} \int_{K_0 - \delta K/2}^{K_0 + \delta K/2} dK \mathcal{F}_s(K, t). \quad (23) $$

As $\langle \mathcal{F}_s \rangle_{\text{stat}}$ depends only on $K$ (the turbulence is statistically stationary) and varies little with $K$ on scales of the order of $\delta K$, one has $\langle \mathcal{F}_s \rangle_{\text{stat}} \approx \langle F_s^c \rangle_{\text{stat}}$. One can therefore use $F_s^c$ in the averaging procedure described above to estimate $\langle \mathcal{F}_s \rangle_{\text{stat}}$. The numerical flux we obtain is then averaged over time following the procedure described above.

### 3.2. Spectra and energy injection

We first present the energy spectra in Fig. 1 for $Pm = 0.0625$ and $Pm = 0.25$. The standard deviation, as measured from 40 instantaneous snapshots, is shown as a shaded region on these spectra. This dispersion stems from temporal fluctuations of the turbulence intensity. The most obvious feature observed in these spectra is the presence of a $k^{-3/2}$ spectrum for the kinetic energy; the traditional Kolmogorov scaling $k^{-5/3}$ appears to be excluded in the $Pm = 0.25$ run, but it cannot be strictly excluded in the $Pm = 0.0625$ run. A $k^{-3/2}$ power law was also found in zero net-flux MRI turbulence (Fromang 2010), and although the spectrum shape differs, our spectra being exempt of any “bump” at intermediate scale. As our runs do not resolve the inertial range of the turbulent cascade yet, these apparent spectral shapes require some comment. The presence of a $k^{-3/2}$ spectrum is usually related to the theoretical argument of Iroshnikov (1963) and Kraichnan (1965) (or IK). However, MRI turbulence is not strongly magnetized, so it falls outside the domain of validity of

<footnote>remap method, the Nyquist frequency waves are remapped to large-scale waves in physical space, which might lead to unphysical behaviours.</footnote>

<footnote n="2">² The quoted resolution now includes the aliasing domain. This convention differs from the one adopted in our previous papers, where only the “useful” domain was accounted for when quoting resolutions.</footnote>
<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Energy Injection spectrum at $Pm = 0.0625$ (**left**) and $Pm = 0.25$ (**right**). Although the injection is significantly reduced at small $Pm$, shape of the spectrum is similar and dominated by the largest scale.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Bidimensional magnetic energy spectra at $Pm = 0.25$. **Left**: $(k_x, k_z)$ spectrum, averaged in $k_y$. **Centre**: $(k_x, k_y)$ spectrum, averaged in $k_z$. **Right**: $(k_y, k_z)$ spectrum averaged in $k_x$. Each contour correspond to a factor 10 in magnetic energy.</figcaption></figure>

the IK phenomenology. Moreover, the magnetic field spectrum does not follow any well-defined power law, as expected from the wide and overlapping injection (see below) and dissipation spectra, indicating that the spectrum we get is *not* an inertial spectrum. We are therefore forced to conclude that although the kinetic spectrum looks like an IK or Kolmogorov spectrum, it is described neither by the IK or Kolmogorov phenomenologies nor by recent extensions (Boldyrev 2005).

Changing the magnetic Prandtl number does not change the power-law index for the kinetic energy. We note, however, two major effects: the overall spectra amplitudes are reduced and the dissipation scales move to larger scale as one reduces $Pm$. These two effects are expected since it is known that smaller $Pm$ turbulence is associated with lower transport efficiency hence a weaker injection of energy in the cascade. This effect is confirmed by the injection spectra (Fig. 2), which are significantly reduced at smaller $Pm$.

We note that the energy injection peaks at the largest scale of the box, although injection still exists at $k \sim 10$. Therefore, although a power-law spectrum is found for $2 < k < 10$, this spectrum cannot be described as an “inertial range” since energy is still injected at these intermediate scales.

We present in Fig. 3 bidimensional spectra of magnetic energy for $Pm = 0.25$. Kinetic spectra are not shown as they

share essentially the same properties. These spectra were obtained by averaging 3D energy spectra over 40 orbits and taking the average in the $k_x$, $k_y$, or $k_z$ directions. We first note a strong anisotropy in the $(k_x, k_y)$ plane which indicates that trailing shearing waves ($k_x k_y > 0$) have more energy than leading shearing waves ($k_x k_y < 0$). As we see below, this results in non-zero shear transfer terms.

Looking at the aspect ratio of the energy contours, we see that turbulence is slightly less anisotropic at large $k$ than at small $k$ (the contours are less “elongated” at large $k$), although complete isotropy is not yet reached in this simulation. Let us, however, point out that the spectral truncation (due to the finite resolution) tends to deform the contours at large $k$, which might accelerate the return to isotropy. One should therefore perform higher resolution runs (or at least double $N_y$) in order to confirm this return to isotropy. In principle, one would expect a return to isotropy at small scales if the non-linear transfer terms dominate all the other terms (injection, body forces) at large enough $k$. However, this is not always the case (e.g. in the presence of a strong mean magnetic field).

The $(k_x, k_z)$ spectrum shows that turbulence is essentially isotropic at large $k$ in that plane. For $k \sim 1$, we find a slight anisotropy where modes with $k_z \neq 0$ are favoured. This is probably a result of large-scale MRI unstable modes, which all have
<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Average relative magnetic helicity (**left**) and cross helicity (**right**) spectra in the $Pm = 0.25$ case (black lines). Instantaneous spectra are represented in light blue. The absolute value of relative helicities is plotted here, since the helicity sign is constantly changing.</figcaption></figure>

$k_z \neq 0$ in the presence of a mean vertical field. Note that this anisotropy disappears very quickly as one moves to larger $k$.

### 3.3. Magnetic helicity and cross helicity

The presence of kinetic and/or magnetic helicity in MHD turbulence is often invoked to explain large-scale dynamo action. Indeed, it is known that an inverse cascade of magnetic helicity can appear in fully developed helical MHD turbulence (Frisch et al. 1975), potentially leading to a buildup of magnetic energy on a large scale. Some kinematic effects often used in mean field dynamos, like the $\alpha$ effect (Moffatt 1978), also lead to the generation to large-scale helical fields (Brandenburg & Subramanian 2005). Magnetic helicity has therefore been suggested as a possible driving mechanism (or at least a tracer) of disc dynamos (Blackman 2010). Moreover, several authors have tried to link magnetic helicity conservation and magnetic helicity flux to the saturation properties of the MRI (Vishniac 2009; Käpylä & Korpi 2010).

We define the magnetic helicity $\mathcal{H}^M = \langle A \cdot (B - \langle B \rangle) \rangle$, where $A$ is the vector potential of the fluctuation (this expression is gauge-invariant in the shearing box), and $\langle \cdot \rangle$ denotes a volume average. We show in Fig. 4 (left) the spectrum of relative helicity $K|H_K|/2M_K$ for the $Pm = 0.25$ run. As can be seen, the relative helicity is less than 1% for all scales of these simulations. Moreover, this quantity is strongly fluctuating and its sign is not well defined<fnref n="3" />.

These results tend to indicate that magnetic helicity is dynamically unimportant in the unstratified simulations presented here and that MRI saturation is not related to a magnetic quenching effect due to magnetic helicity accumulated at large scale. This was to be expected in the first place as unstratified shearing boxes (both with and without mean field) are mirror symmetric; however, this picture might change when stratification is included.

Another quantity that might play a role in the MHD turbulence cascade is the cross helicity $\mathcal{H}^C = \langle v \cdot b \rangle$ (see e.g. Perez & Boldyrev 2010). When cross helicity is non zero, the energy of Alfvén waves travelling along and against the mean field are not equal. For this reason, turbulence with cross helicity is often called *imbalanced* turbulence. Locally imbalanced turbulence is

often observed in strong MHD turbulence, where the guide field is weaker than the turbulent fluctuations. To check whether unstratified MRI turbulence was imbalanced, we computed cross-helicity spectra of our simulations (Fig. 4, right). As for the magnetic helicity, we find that the relative cross helicity is small ($<10^{-1}$) and highly fluctuating at all scales. This tends to indicate that MRI turbulence is not imbalanced in our setup. As for magnetic helicity, this result was to be expected because of the mirror symmetry properties of the unstratified shearing box. The absence of any significant cross-helicity also shows that energy spectra in Elsässer variables $z^{\pm} = v \pm b$ are equal and proportional to the kinetic plus magnetic energy spectrum.

### 3.4. Energy fluxes

The energy fluxes (17)–(21) allow one to check the average direction of the energy flux in spectral space. To explain the dependence of the turbulent transport on $Pm$, several authors (e.g. Lesur & Longaretti 2007; Fromang et al. 2007) have suggested that an inverse cascade driven by resistive and viscous scales might be at work. Since magnetic helicity is irrelevant to this problem, only the kinetic, magnetic, exchange, and shear fluxes are important for the unstratified shearing box and should be checked for an inverse cascade.

We present the energy fluxes at $Pm = 0.25$ and $Pm = 0.0625$ in Fig. 5. Standard deviations are shown for kinetic and magnetic fluxes as shaded regions. These deviations are computed following the procedure described in Sect. 3.2. We always find positive fluxes, meaning that the non-linear transfers are *forward* or direct (from large to small scales), on average. However, at larger scale, the standard deviation may allow for an inverse cascade of kinetic energy in the $Pm = 0.25$ run. This indicates that, although the kinetic cascade is direct on (time) average over most of the spectrum, inverse cascades can sometimes be observed on the largest scales. This inverse cascade of kinetic energy could be an explanation for the large-scale hydrodynamic structures that are observed in several MRI turbulence simulations, such as vortices (Fromang & Nelson 2005) and zonal flows (Johansen et al. 2009).

We also observe that the energy cascade is dominated by the magnetic and exchange fluxes down to the resistive scale, since the kinetic flux and shear fluxes are almost negligible. The shear fluxes are also always positive. This is due to the anisotropy

<footnote n="3">³ The absolute value of the relative magnetic helicity is plotted in Fig. 4.</footnote>
<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Energy fluxes as computed from Eqs. (17)–(21) at $Pm = 0.0625$ (**left**) and $Pm = 0.25$ (**right**). Energy fluxes are always direct (from large to small scales) and dominated by magnetic and exchange fluxes.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Zoom on the energy fluxes in the dissipative range for $Pm = 0.0625$ (log-log representation with $K > 10$).</figcaption></figure>

of MRI turbulence in which shearing waves statistically have a larger amplitude than leading waves (see Sect. 3.2).

In the $Pm = 0.0625$ case, the kinetic flux becomes dominant at subresistive scale ($K \ge 10$), indicating that the cascade becomes essentially hydrodynamic below the resistive scale, as expected for low $Pm$ MHD turbulence (see also Fig. 6). Moreover, the kinetic flux dominates the kinetic shear transfer term at least for the larger $k$, which indicates that, as far as the non-linear transfers are concerned, the cascade is close to isotropic at small scales, as noted in Sect. 3.2. Finally, none of the fluxes reaches a plateau at intermediate scales, which would be expected in the presence of an inertial range. This indicates that the $k^{-3/2}$ kinetic spectrum found in Fig. 1 is not, properly speaking, an inertial spectrum.

### 3.5. Energy transfer locality

To test the locality of the energy transfer in spectral space, we plotted the transfer functions in the $Pm = 0.25$ case for several values of $K$: at the injection scale ($K = 1$), in the intermediate range ($K = 5$) and in the resistive range ($K = 20$). We first plot the kinetic to kinetic and magnetic to magnetic transfers in Fig. 7. The transfers $T_{vv}$ and $T_{bb}$ obtained at all scales perfectly illustrate local energy exchanges. Energy is

taken from wavenumbers slightly smaller than $K$ and is transferred to wavenumbers slightly larger than $K$, except (not surprisingly) for $K = 1$. As expected from the energy flux, we also find that the cascade is direct, with energy going from small to large wavenumbers. Finally, the $T_{vv}$ transfers are always much smaller than the $T_{bb}$ transfers above the dissipation range, illustrating the magnetically dominated energy transfer described above.

We next plot the exchange transfers $T_{vb}$ and $T_{bv}$ in Fig. 8. We note in this case that the scales involved in each transfer are much broader. In particular, $T_{vb}$ measured at the resistive scale ($K = 20$) has contributions coming from all scales, including the largest injection scales. This effect can also be seen in $T_{bv}(Q, K = 1)$, which exhibits a very long tail toward large $k$, down to the resistive scale. Comparing these transfers to $T_{vv}$ and $T_{bb}$ directly show that these terms are highly non-local. In turns, this indicates that the exchange flux computed in the previous section is non-local. The results in the $Pm = 0.0625$ case are not shown here as they are very similar to the $Pm = 0.25$ case.

We note that despite the non-locality of the energy exchanges, the overall cascade direction is still forward, confirming our previous interpretation regarding the exchange flux. We also note that the shear transfer terms are local by definition, since they transfer energy to neighbouring shells.

## 4. Summary and discussion

In this paper, we have described some properties of the turbulent cascade found in incompressible MRI turbulence. We showed that compared to isotropic MHD turbulence, the presence of a mean shear led to several new transfer terms and introduced a source of anisotropy. We computed the effect of each non-linear term and found that all the terms contribute to a direct cascade of energy (from large to small scales), but some terms involved non-local transfers in Fourier space. This lack of locality is due to the Lorentz force and to the magnetic stretching term of the induction equation (combined here in the exchange transfer term). We also showed that magnetic helicity, although non zero, was totally negligible and should not play any role in the behaviour of MRI turbulence.

The presence of non-local transfer terms in the MRI turbulent cascade is the most important finding of this work. It indicates that, in principle, the large scales – responsible for the transport – can directly interact with the small dissipative ones
<figure><img src="image_8.png" /><figcaption>Fig. 7. Transfers function $T_{vv}(Q, K)$ and $T_{bb}(Q, K)$ in the $Pm = 0.25$ run for $K = 1; 5; 20$. These transfers are local in Fourier space (see text).</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 8. Transfers function $T_{vb}(Q, K)$ and $T_{bu}(Q, K)$ in the $Pm = 0.25$ run for $K = 1; 5; 20$. These transfers are non-local in Fourier space (see text).</figcaption></figure>

through non-linear terms. This direct interaction could of course explain the correlation observed between $Pm$ and the turbulent transport of angular momentum $\alpha$ (Longaretti & Lesur 2010). However, it should be pointed out that non-local transfers have already been found in isotropic MHD turbulence (e.g. Alexakis et al. 2007). Therefore, MRI turbulence is not special regarding the nature of these non-linear transfers.

Although some non-linear terms are found to directly connect injection and dissipation scales in current simulations, one might wonder if this could be true in a more realistic setup where the injection and dissipation scales are separated by a wide range of scales (typically $10^{10}$). In other words, what is the maximum scale separation these terms can connect? A partial answer to this question is given by Aluie & Eyink (2010). To describe their result, let us define the structure functions:

$$ \delta v_{l,p} = \langle |v(\mathbf{x} + \mathbf{I}) - v(\mathbf{x})|^p \rangle. \quad (24) $$

In the inertial range the structure function depends only on $|\mathbf{I}|$ and $\delta v_{l,p} \propto |\xi_p|$, where $\xi_p$ is the structure function index of order $p$. It is then possible to derive an upper bound to the non-linear transfer terms thanks to the Hölder inequality. Applying this procedure to the non-local transfer $T_{ub}$, Aluie & Eyink (2010) found

$$ |T_{ub}(Q, K)| < (\text{const.}) Q^{1-\xi_3^u/3} K^{-2\xi_3^b/3}, \quad (25) $$

where $Q$ and $K$ are dyadic (octave) wavenumbers and $K > Q/2$. Similar terms can be obtained for $T_{bu}$ and $K < Q/2$. If one assumes Iroshnikov-Kraichnan theory, one has $\xi_3^u = \xi_3^b = 3/4$. In constrast, considering Goldreich-Sridhar (GS) phenomenology, which should be valid for MRI turbulence, one gets $\xi_{3\parallel}^u = \xi_{3\parallel}^b = 3/2$ and $\xi_{3\perp}^u = \xi_{3\perp}^b = 1$. In all these cases, (25) indicates that the non-locality of these transfer terms cannot extend over several decades, with a typical scaling $T_{ub}(Q, K) \sim \epsilon(K/Q)^{-2/3}$ for GS turbulence (with $\epsilon$ the usual turbulence energy injection rate).

We therefore conclude that the non-locality in Fourier space is somewhat relative. Although $T_{ub}$ and $T_{bu}$ are non-local compared to $T_{uu}$ and $T_{bb}$, these terms should be local when one

considers transfers over several decades. Unfortunately, separating the injection scale from the dissipative scales by several decades is numerically difficult. It is even harder for MRI turbulence since the injection term is rather broad in spectral space compared to forced turbulence. Assuming the injection and dissipation scales, both spread over one decade in Fourier space, one typically needs $20000^3$ simulations to get a 2-decade inertial range in which non-local transfers are significantly reduced. This kind of resolution is for the moment out of reach of the best computational facilities.

Nevertheless, we can conjecture that if the $Pm-\alpha$ correlation is actually due to the non-local transfers, then it should vanish when the injection and dissipation scales are well separated, as is the case in some accretion discs. Although this conclusion looks rather reassuring for the relevance of today simulations regarding small-scale dissipation, it tells us neither what the asymptotic value of $\alpha$ is in this limit nor how MRI turbulence behaves when the scale separation is not achieved, a situation that probably occurs in the inner regions of protoplanetary discs where $\Lambda_\eta$ is not very large.

## Appendix A: Shearing-wave approach to the shear transfer term

The shell-filter decomposition can be properly defined using a projector operator $\Pi$ on the field $F$ in the sheared frame:

$$ [\Pi_{K_j}(F)](x', t) = \sum_{k' \in \Sigma_j(t)} F_{k'}(t) \exp(ik' \cdot x') \quad (A.1) $$

where $\Sigma_j$ is the shell containing all the shearing waves with a norm between $K_j - \delta K/2$ and $K_j + \delta K/2$:

$$ \Sigma_j(t) = \{ K_j - \delta K/2 < |k' + q\Omega k'_j t e_x| \le K_j + \delta K/2 \}. \quad (A.2) $$

Our notation indicates that the projected $F$ is function of space and time. Also, $[\Pi_{K_j}(F)]$ is real for real fields $F$.
<figure><img src="image_10.png" /><figcaption>**Fig. A.1.** Evolution of shearing waves in the presence of a fixed shell in Fourier space. Some waves can either enter or exit the shell as time evolves (see text).</figcaption></figure>

As can be seen, the waves included in $\Sigma_j(t)$ change with time. This is to be expected since shearing waves should enter and exit shells as does the one defined above (see Fig. A.1). As a result, the projector operator $\Pi_{K_j}$ has an explicit time dependence that leads to non-trivial transfer effects. The above projector operator can be written using Heaviside function $\Theta$:

$$[\Pi_{K_j}(F)](x', t) = \sum_{k'} \Theta(|k(t)| - K_j + \delta K/2) \times \Theta(-|k(t)| + K_j + \delta K/2) \times F_{k'}(t) \exp(ik' \cdot x'), \quad (A.3)$$

where we have defined $k(t)$ as a function of $k'$ as in Eqs. (11)–(13).

One next defines the energy within a shell $E_{K_j}^F = \langle \Pi_{K_j}^2(F) \rangle/2$ where $\langle \cdot \rangle$ denotes a volume average. The energy equation then reads as

$$\frac{dE_{K_j}^F}{dt} = \langle \Pi_{K_j}(F) \frac{\partial \Pi_{K_j}(F)}{\partial t} \rangle = \langle \Pi_{K_j}(F) \Pi_{K_j}(\partial_t F) \rangle + \langle \Pi_{K_j}(F) [\partial_t \Pi_{K_j}](F) \rangle. \quad (A.4)$$

The first term on the righthand side leads to the terms obtained in isotropic turbulence (Alexakis et al. 2007) and introduced in Eqs. (15) and (16) along with the injection terms $I_{v,K}$ and $I_{b,K}$. The second one, however, comes from shearing waves entering and exiting the shells. Using (A.3), it is possible to obtain an exact (though singular) expression of the operator time derivative:

$$\partial_t[\Pi_{K_j}](F) = \sum_{k'} \frac{q\Omega k_y k_x(t)}{|k(t)|} \times \left[ \delta(|k(t)| - K_j + \delta K/2) - \delta(-|k(t)| + K_j + \delta K/2) \right] \times F_{k'}(t) \exp(ik' \cdot x'). \quad (A.5)$$

This expression can be interpreted easily. As an example, let us consider waves with $k_y k_x(t) > 0$. Then, the first $\delta$ function represents waves entering the shell, the second delta represents waves leaving the same shell, and the factor in front of the delta functions quantifies the “flux” of waves going through a shell boundary. This interpretation is similar to the phenomenological picture one can have of waves travelling through a fixed shell in the unsheared Fourier space (Fig. A.1).

We then deduce from (A.4)

$$\frac{dE_{K_j}^X}{dt} = \langle \Pi_{K_j}(F) \Pi_{K_j}(\partial_t F) \rangle + \sum_{k'} \frac{q\Omega k_y k_x(t) F_{k'}^* F_{k'}}{|k(t)|^2} \times \left[ \delta(|k(t)| - K_j + \delta K/2) - \delta(-|k(t)| + K_j + \delta K/2) \right], \quad (A.6)$$

where we have used the property $\delta(x)\theta(x) = \delta(x)/2$ (Fisher 1971). As expected, this expression shows two contributions to the energy evolution inside a shell: a volume contribution and a surface contribution equal to the energy of the waves entering and leaving the shell. Introducing the equation of motions in sheared space (7)–(8) in the above relation leads to the energy equations (15)–(16).

Finally, using the relation

$$\partial_t |k(t)| = \frac{q \Omega k_y k_x(t)}{|k(t)|}, \quad (A.7)$$

one can write (A.6) as

$$\frac{dE_{K_j}^X}{dt} = \langle \Pi_{K_j}(F) \Pi_{K_j}(\partial_t F) \rangle + \sum_{k'} \frac{F_{k'}^* F_{k'}}{2} \delta(t - t_{k'}) \epsilon_{k'} \quad (A.8)$$

where $t_{k'}$ is the instant when the wave $k'$ enters or exits the shell $K_j$ and $\epsilon_{k'} = \pm 1$ for an entering/exiting wave. This somewhat simpler expression has the same interpretation as (A.6).

## Appendix B: Unsheared Fourier transform approach to the shear transfer term

It is possible to understand the origin of the shear transfer term (A.5) starting from the equations of motion in unsheared coordinates (1), (2) with an appropriate use of continuous Fourier transforms. Shear-periodic functions are not absolutely integrable, but this difficulty can be circumvented because their continuous Fourier transform is well-defined as a distribution. To demonstrate this point, we consider a 2D infinite medium in which a field $F$ obeys the model equation

$$\partial_t F(\mathbf{x}, t) - q\Omega x \partial_y F(\mathbf{x}, t) = 0. \quad (B.1)$$

Let us introduce the *unsheared* Fourier transform:

$$\tilde{F}(\mathbf{k}, t) = \frac{1}{(2\pi)^3} \iiint \mathrm{d}\mathbf{x} \, F(\mathbf{x}, t) \exp(-i\mathbf{k} \cdot \mathbf{x}) \\ F(\mathbf{x}, t) = \iiint \mathrm{d}\mathbf{k} \, \tilde{F}(\mathbf{k}, t) \exp(i\mathbf{k} \cdot \mathbf{x}). \quad (B.2)$$
If $F$ obeys the shearing sheet boundary conditions, the solution of the equation is a Fourier series $F_p$ of the form (10) with Fourier coefficients $F_{k'}$. Its Fourier transform in the unsheared spectral space is then

$$ \tilde{F}_p(\mathbf{k}, t) = \sum_{k'} \delta(\mathbf{k} - \mathbf{k}_0(t)) F_{k'}, \quad (B.3) $$

where $\mathbf{k}_0(t) = \mathbf{k}' + \mathbf{V}_{k'}t$ and where we have defined the “Fourier velocity” $\mathbf{V}_k = q\Omega k_y \mathbf{e}_x$, evaluated in $\mathbf{k}'$. By construction $\tilde{F}_p(\mathbf{k}, t)$ is solution of

$$ \partial_t \tilde{F}(\mathbf{k}) + q\Omega k_y \partial_{k_x} \tilde{F}(\mathbf{k}) = 0, \quad (B.4) $$

which is the Fourier transform of our model equation$^4$. The left-hand side of this equation can be interpreted as a comoving derivative of $\tilde{F}(\mathbf{k}, t)$ in Fourier space with the Fourier velocity $\mathbf{V}_k$. This equation tells us that the amplitude of the waves is constant when one moves at velocity $\mathbf{V}_k$ in Fourier space consistently with the form of the solution $\tilde{F}_p$. It can also be interpreted as constant amplitude shearing waves, as expected.

It is then possible to introduce the projector operator, now time-independent because it has been defined in unsheared coordinates:

$$ [\Pi_{K_j}(F)](x, t) = \iiint d\mathbf{k} \Theta(|\mathbf{k}| - K_j + \delta K/2) \\ \times \Theta(-|\mathbf{k}| + K_j + \delta K/2) \\ \times \tilde{F}(\mathbf{k}, t) \exp(i\mathbf{k} \cdot \mathbf{x}). \quad (B.5) $$

As previously, the shell energy ($E_{K_j}^F = \langle \Pi_{K_j}^2(F) \rangle/2$) time variation follows from

$$ \frac{dE_{K_j}^F}{dt} = \langle \Pi_{K_j}(F) \partial_t [\Pi_{K_j}(F)] \rangle \\ = -\langle \Pi_{K_j}(F) \Pi_{K_j}(\nabla_k \cdot (\mathbf{V}_k \tilde{F})) \rangle \quad (B.6) $$

where (B.4) has been used in the last equality along with $\nabla_k \cdot \mathbf{V}_k = 0$. Because the volume average selects the zero-frequency contributions, this leads us back to (A.6) with the help of the relation

$$ \int dk_x \Theta(|k| - C) \partial_{k_x} \tilde{F}(\mathbf{k}) = - \int dk_x \frac{k_x}{|k|} \delta(|k| - C) \tilde{F}(\mathbf{k}), $$

and with the use of $\delta(x)\theta(x) = \delta(x)/2$ and (B.3).

In the spirit of the integral expressions used in this appendix, Eq. (A.6) can be recast in integral form by introducing the energy density in Fourier space $\mathcal{E}^F(\mathbf{k}, t)$ defined by

$$ \mathcal{E}^F(\mathbf{k}, t) = \sum_{k'} \delta(\mathbf{k} - \mathbf{k}_0(t)) \frac{F_{k'} F_{k'}^*}{2}. \quad (B.7) $$

With this definition,

$$ E_{K_j}^F = \iiint d\mathbf{k} \Theta(|\mathbf{k}| - K^-) \Theta(-|\mathbf{k}| + K^+) \mathcal{E}^F(\mathbf{k}, t) \quad (B.8) $$

while

$$ \frac{dE_{K_j}^F}{dt} = \iiint d\mathbf{k} \mathbf{V}_k \cdot \mathbf{n} \mathcal{E}^F(\mathbf{k}, t) [\delta(|\mathbf{k}| - K^-) - \delta(-|\mathbf{k}| + K^+)] \quad (B.9) $$

$$ = - \iint_{\partial K_j} d\mathbf{k} \mathbf{V}_k \cdot \mathbf{n} \mathcal{E}^F(\mathbf{k}, t). \quad (B.10) $$

In these relations, $K^\pm = K \pm \delta K/2$ has been defined, the second integration is performed on the surface of the shell $K_j \equiv \partial K_j$, and $\mathbf{n}$ is the normal to this surface. The last expression has an explicit flux form.

This approach can be applied to the original shearing sheet MHD equations. The time dependence due to the shear term will produce the just computed shear flux contribution.

**Acknowledgements.** G.L. acknowledges support by STFC. P.Y.L. acknowledges the hospitality of the Isaac Newton Institute and of the DAMTP in Cambridge, where parts of this research has been conducted. G.L. thanks S. Fromang for his comments on the initial version of the manuscript. This work was granted access to the HPC resources of IDRIS under allocation x2009042231 made by GENCI (Grand Equipement National de Calcul Intensif).

## References

Alexakis, A., Mininni, P. D., & Pouquet, A. 2007, New J. Phys., 9, 298

Aluie, H., & Eyink, G. L. 2010, Phys. Rev. Lett., 104, 081101

Balbus, S. A. 2003, ARA&A, 41, 555

Balbus, S. A., & Hawley, J. F. 1991, ApJ, 376, 214

Balbus, S. A., & Henri, P. 2008, ApJ, 674, 408

Blackman, E. G. 2010, Astron. Nachr., 331, 101

Boldyrev, S. 2005, ApJ, 626, L37

Brandenburg, A., & Subramanian, K. 2005, Phys. Rep., 417, 1

Canuto, C., Hussaini, M. Y., Quarteroni, A., & Zang, T. A. 1988, Spectral Methods in Fluid Dynamics (Springer)

Chandrasekhar, S. 1960, Proc. Nat. Acad. Sci., 46, 253

Fisher, B. 1971, The Quarterly J. Math., 22, 291

Frisch, U. 1995, Turbulence (Cambridge University Press)

Frisch, U., Pouquet, A., Leorat, J., & Mazure, A. 1975, J. Fluid Mech., 68, 769

Fromang, S. 2010, A&A, 514, L5

Fromang, S., & Nelson, R. P. 2005, MNRAS, 364, L81

Fromang, S., Papaloizou, J., Lesur, G., & Heinemann, T. 2007, A&A, 476, 1123

Goldreich, P., & Lynden-Bell, D. 1965, MNRAS, 130, 125

Hawley, J. F., Gammie, C. F., & Balbus, S. A. 1995, ApJ, 440, 742

Iroshnikov, P. S. 1963, AZh, 40, 742

Johansen, A., Youdin, A., & Klahr, H. 2009, ApJ, 697, 1269

Käpylä, P. J., & Korpi, M. J. 2010, MNRAS, submitted [arXiv:1004.2417]

Kraichnan, R. H. 1965, Phys. Fluids, 8, 1385

Lesur, G., & Longaretti, P.-Y. 2007, MNRAS, 378, 1471

Lesur, G., & Ogilvie, G. I. 2010, MNRAS, 404, L64

Lesur, G., & Papaloizou, J. C. B. 2010, A&A, 513, A60

Lithwick, Y. 2007, ApJ, 670, 789

Longaretti, P., & Lesur, G. 2010, A&A, 516, A51

Moffatt, H. K. 1978, Magnetic Field Generation in Electrically Conducting Fluids (Cambridge University Press)

Perez, J. C., & Boldyrev, S. 2010, Phys. Plasmas, 17, 055903

Pessah, M. E., & Chan, C.-K. 2008, ApJ, 684, 498

Pessah, M. E., & Goodman, J. 2009, ApJ, 698, L72

Regev, O., & Umurhan, O. M. 2008, A&A, 481, 21

Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337

Thomson, W. 1887, Phil. Mag., 24, 188

Umurhan, O. M., & Regev, O. 2004, A&A, 427, 855

Velikhov, E. P. 1959, Sov. Phys.-JETP, 36, 1398

Verma, M. K. 2004, Phys. Rep., 401, 229

Vishniac, E. T. 2009, ApJ, 696, 1021

<footnote n="4">⁴ The Fourier transform of (B.1) gives (B.4) except for an extra term that vanishes once the Fourier series expression of $F$ is used.</footnote>