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

# Universal small-scale structure in turbulence driven by magnetorotational instability

Vladimir Zhdankin,<fnref n="1" />★ Justin Walker,<fnref n="2" /> Stanislav Boldyrev<fnref n="2" />,<fnref n="3" /> and Geoffroy Lesur<fnref n="1" />,<fnref n="4" />,<fnref n="5" />†

Accepted 2017 February 9. Received 2017 February 7; in original form 2016 August 13

## ABSTRACT

The intermittent small-scale structure of turbulence governs energy dissipation in many astrophysical plasmas and is often believed to have universal properties for sufficiently large systems. In this work, we argue that small-scale turbulence in accretion discs is universal in the sense that it is insensitive to the magnetorotational instability (MRI) and background shear, and therefore indistinguishable from standard homogeneous magnetohydrodynamic (MHD) turbulence at small scales. We investigate the intermittency of current density, vorticity and energy dissipation in numerical simulations of incompressible MHD turbulence driven by the MRI in a shearing box. We find that the simulations exhibit a similar degree of intermittency as in standard MHD turbulence. We perform a statistical analysis of intermittent dissipative structures and find that energy dissipation is concentrated in thin sheet-like structures that span a wide range of scales up to the box size. We show that these structures exhibit strikingly similar statistical properties to those in standard MHD turbulence. Additionally, the structures are oriented in the toroidal direction with a characteristic tilt of approximately 17:5, implying an effective guide field in that direction.

**Key words:** accretion, accretion discs – instabilities – MHD – plasmas – turbulence.

## 1 INTRODUCTION

Turbulence plays a central role in accretion discs by governing the angular momentum transport responsible for mass accretion on to compact objects (Balbus & Hawley 1998), amplifying magnetic fields via the dynamo (Brandenburg et al. 1995), and intermittently dissipating magnetic energy, which may lead to non-thermal particle acceleration and X-ray flares (Baganoff et al. 2001; Markoff et al. 2001; McClintock & Remillard 2006). Many questions on the nature of turbulence in accretion discs remain, and it is important to understand to what extent insights from standard magnetohydrodynamic (MHD) turbulence can be adapted. In this work, we use the term standard MHD turbulence to refer to externally driven or decaying MHD turbulence in a periodic box in the absence of background gradients (and therefore statistically homogeneous). Standard MHD turbulence serves as a conventional testing ground for controlled numerical studies of turbulence in the inertial range.

Turbulence can be driven in ionized accretion discs by the magnetorotational instability (MRI) when a weak background magnetic field and Keplerian shear are present (Balbus & Hawley 1991). The dynamics are often studied in the shearing-box approximation, which provides a simple framework for describing turbulence in a small corotating patch of the accretion disc. Previous numerical studies of MRI-driven turbulence in a shearing box have primarily focused on large-scale quantities, such as the angular momentum transport coefficients, magnetic and kinetic energy spectra (Fromang 2010; Lesur & Longaretti 2011), anisotropy (Nauman & Blackman 2014; Murphy & Pessah 2015) and dynamo cycles (Bai & Stone 2013). Physical and numerical degrees of freedom such as the net flux (e.g. Salvesen et al. 2016), stratification (e.g. Simon, Hawley & Beckwith 2011), dissipation coefficients (Fleming, Stone & Hawley 2000; Fromang et al. 2007; Lesur & Longaretti 2007), box size (Bodo et al. 2008, 2011; Simon, Beckwith & Armitage 2012; Shi, Stone & Huang 2016) and boundary conditions were also explored. These numerical studies uncovered many important features of turbulence in a shearing box and reproduce aspects of global accretion disc simulations (Sorathia et al. 2012). Despite their relative physical simplicity, shearing-box simulations of MRI-driven turbulence have only recently identified a

<footnote n="1">¹JILA, NIST and University of Colorado, 440 UCB, Boulder, CO 80309, USA</footnote>

<footnote n="2">²Department of Physics, University of Wisconsin-Madison, 1150 University Avenue, Madison, WI 53706, USA</footnote>

<footnote n="3">³Space Science Institute, Boulder, CO 80301, USA</footnote>

<footnote n="4">⁴Université Grenoble Alpes, IPAG, F-38000 Grenoble, France</footnote>

<footnote n="5">⁵CNRS, IPAG, F-38000 Grenoble, France</footnote>

<footnote>* E-mail: zhdankin@jila.colorado.edu</footnote>

<footnote>† Visiting fellow.</footnote>
putative inertial range (Kunz, Stone & Quataert 2016; Walker, Lesur & Boldyrev 2016), with previous numerical simulations showing little quantitative resemblance to standard MHD turbulence (e.g. Fromang 2010; Lesur & Longaretti 2011).

One aspect of MRI-driven turbulence that was largely neglected to date is the small-scale structure of the turbulence, involving quantities such as the energy dissipation rate, current density and vorticity. Small-scale turbulence is often believed to be universal in the sense that it is insensitive to the driving mechanism and boundary conditions, although it may depend on the system size due to intermittency (scale-dependent inhomogeneity) and scale-dependent anisotropy (Goldreich & Sridhar 1995; Boldyrev 2006; Mason et al. 2012; Perez et al. 2012). Intermittency causes a large fraction of energy dissipation in MHD turbulence to be concentrated in coherent structures, such as current sheets and vorticity sheets, which occupy a small fraction of the volume (Wan et al. 2016; Zhdankin, Boldyrev & Uzdensky 2016a). These structures may determine the kinetic mechanisms of energy dissipation (Matthaeus et al. 2015), contribute to non-thermal particle acceleration, appear as intense observable flares (Zhdankin, Uzdensky & Boldyrev 2015a,b) and cause temperature inhomogeneity that leads to mineral formation (McNally et al. 2014).

A central open question is whether the character of MRI-driven turbulence transitions to standard MHD turbulence at sufficiently small scales (which requires large simulations to properly discern), or whether the presence of the MRI across many scales modifies the dynamics at all scales. It is likewise unknown whether the MRI modifies scale-dependent effects such as intermittency and anisotropy; even the large-scale variability of the MRI may non-trivially reduce or increase the degree of inhomogeneity and anisotropy at smaller scales. Furthermore, whereas numerical simulations established that standard MHD turbulence is characterized by thin intermittent structures that span large scales (e.g. Zhdankin et al. 2013, 2014), it is unclear whether such structures can persist in the background Keplerian shear of an accretion disc.

In this work, we take a step towards answering these questions by investigating the intermittency of small-scale quantities in numerical simulations of incompressible MHD turbulence driven by the MRI in an unstratified shearing box. We demonstrate that the system is classically intermittent by showing that the distribution of coarse-grained energy dissipation rates broadens with decreasing scale, at a rate comparable to that in hydrodynamic and MHD turbulence. We then perform a statistical analysis of intermittent dissipative structures and characterize their energetics, morphology and orientation. We show that the statistical properties of these structures are in strong quantitative agreement with similar measurements in standard MHD turbulence driven at large scales (Zhdankin et al. 2014, 2016a,b). These results support the picture that the small-scale dynamics of turbulence, involving the current density, vorticity, energy dissipation rate and possibly the asymptotic energy spectrum, are insensitive to the large-scale shear and MRI for sufficiently large systems.

## 2 METHODOLOGY

We consider a periodic shearing box with $x$ representing the radial direction, $y$ representing the toroidal direction and $z$ representing the vertical direction. The incompressible MHD equations in a shearing box read

$$ (\partial_t - q\Omega_0 x \partial_y + \mathbf{v} \cdot \nabla)\mathbf{v} - \mathbf{B} \cdot \nabla \mathbf{B} = -\nabla P + \nu \nabla^2 \mathbf{v} \\ - 2\Omega_0 \times \mathbf{v} + q\Omega_0 v_x \hat{\mathbf{y}}, $$

$$ (\partial_t - q\Omega_0 x \partial_y + \mathbf{v} \cdot \nabla)\mathbf{B} - \mathbf{B} \cdot \nabla \mathbf{v} = \eta \nabla^2 \mathbf{B} - q\Omega_0 B_x \hat{\mathbf{y}}, \\ \nabla \cdot \mathbf{v} = \nabla \cdot \mathbf{B} = 0, \quad (1) $$

where $\mathbf{v}$ is the velocity field, $\mathbf{B} = \mathbf{B}_0 + \mathbf{b}$ is the magnetic field (where $\mathbf{B}_0$ is the background field and $\mathbf{b}$ is the fluctuating part), $\Omega_0 = \Omega_0 \hat{\mathbf{z}}$ is the orbital frequency, $q = 3/2$ is the Keplerian shear parameter; we absorbed a factor of $1/\sqrt{4\pi\rho}$ in the magnetic field, where $\rho$ is the plasma density. We focus our primary analysis on the current density $j = \nabla \times \mathbf{B}$ and vorticity $\omega = \nabla \times \mathbf{v}$. The local resistive and viscous energy dissipation rates per unit volume are given by $\epsilon^\eta = \eta j^2$ and $\epsilon^\nu = 2\nu\sigma_{ij}\sigma_{ij}$, respectively, where $\eta$ is the resistivity, $\nu$ is the viscosity and $\sigma_{ij} = (\partial v_j / \partial x_i + \partial v_i / \partial x_j) / 2$ is the rate-of-strain tensor. We note that terms proportional to $\Omega_0$ are of a lower order in gradients than the other terms, and hence can be anticipated to become unimportant at small scales, leaving only the usual MHD terms.

We consider the suite of numerical simulations described in Walker et al. (2016). The dimensions of the box are chosen to be ($L_x, L_y, L_z$) = (2, 4, 1) and the simulation lattice is 1024 × 1024 × 512 cells. We fix the Reynolds number $Re \equiv \Omega_0 L_z^2 / \nu = 45000$ (which differs from the conventional turbulence Reynolds number, $Re_{turb} \equiv v_{rms} (L/2\pi) / \nu \sim 4000$) and magnetic Prandtl number $Pm \equiv \nu / \eta = 1$, while varying the superimposed vertical background field $\mathbf{B}_0 = B_0 \hat{\mathbf{z}}$, taking $B_0 \in \{0.03, 0.01, 0.005\}$ (in numerical units), which correspond to the initial fastest growing linear MRI modes at length-scales $L_{MRI} \equiv 2\pi/k_{MRI} = 2\pi(4/\sqrt{15})B_0 / \Omega_0 \in \{0.195, 0.065, 0.032\}L_z$, where $k_{MRI}$ is the wavevector for the fastest growing mode (Balbus & Hawley 1998). Hence, the simulations have varying dimensionless parameters $\beta = \Omega_0^2 L_z^2 / B_0^2$ and Elsässer number $\Lambda_\eta = B_0^2 / \Omega_0 \eta$. For each simulation, we consider a minimum of seven snapshots, each separated in time by $4\Omega_0^{-1}$, during a statistical steady state. We show results from $B_0 = 0.03$ unless otherwise stated; the results that we describe are very similar for all three simulations.

For reference, we now explicitly show that energy and cross-helicity injection occur mainly at large scales. The total energy injection rate is given by $\alpha_E = q\Omega_0 \int (v_x v_y - B_x B_y) d^3x$, while the total cross-helicity injection rate is given by $\alpha_C = (q - 2)\Omega_0 \int (v_x B_y - B_x v_y) d^3x$. The integrand can be evaluated in Fourier space to determine the injection rates at different wavenumbers, which we denote by $\tilde{\alpha}_E(k)$ and $\tilde{\alpha}_C(k)$. The (normalized) cumulative injection rate for energy is then given by $I_E(k) = \int_k^\infty \langle \tilde{\alpha}_E(k) \rangle dk / \int_k^\infty \langle \tilde{\alpha}_E(k) \rangle dk$, where $\langle \dots \rangle$ denotes ensemble average (i.e. time average). Similarly, the cumulative injection rate for cross-helicity can be characterized by $I_C(k) = \int_k^\infty \langle |\tilde{\alpha}_C(k)| \rangle dk / \int_k^\infty \langle |\tilde{\alpha}_C(k)| \rangle dk$, where the absolute values are taken since cross-helicity is not sign-definite and is zero on average. The injection rates $I_E$ and $I_C$ are shown in Fig. 1. We find that the injection rates strongly decrease with increasing wavenumber, becoming negligibly small compared to the cascade rate (which is unity in Fig. 1 across inertial-range scales).

## 3 RESULTS

We first consider the distribution of coarse-grained energy dissipation rates, which is a conventional method for characterizing the scale-dependent nature of intermittency (Biskamp 1995; Frisch 1995; Merrifield et al. 2005). The coarse-grained resistive (viscous) energy dissipation rates at the $n$th level, denoted by $\epsilon_n^\eta (\epsilon_n^\nu)$, are obtained by subdividing the simulation box $n$ times (in each direction) to obtain a set of subdomains of size
<figure><img src="image_2.png" /><figcaption>Figure 1. The cumulative injection rate of energy $I_E(k)$ (blue) and of cross helicity $I_C(k)$ (red) at wavenumbers greater than $k$.</figcaption></figure>

$(L_x, L_y, L_z) \cdot 2^{-n}$ and measuring the average $\epsilon^n (\epsilon^v)$ in each of these subvolumes. The (total) coarse-grained energy dissipation rate in each subdomain is given by $\epsilon_n = \epsilon_n^v + \epsilon_n^h$.

As shown in the first panel of Fig. 2, the probability density function for coarse-grained energy dissipation rates normalized to the mean, $P(\epsilon_n/\epsilon_0)$, broadens with the number of subdivisions $n$, indicating that strong dissipative events are increasingly localized in space at smaller scales. We find that $P(\epsilon_n/\epsilon_0)$, as well as $P(\epsilon_n^v/\epsilon_0^v)$ and $P(\epsilon_n^h/\epsilon_0^h)$, are reasonably well fit for all levels $n$ by the lognormal distribution:

$$P(\epsilon_n/\epsilon_0) = \frac{1}{\sqrt{2\pi\sigma_n^2}} \frac{\epsilon_0}{\epsilon_n} \exp \left[ -\frac{1}{2\sigma_n^2} \left( \log \frac{\epsilon_n}{\epsilon_0} - \mu_n \right)^2 \right], \quad (2)$$

where $\mu_n$ and $\sigma_n$ are the location parameter and scale parameter, respectively, with $\mu_n = -\sigma_n^2/2$ to ensure a mean of unity. As shown in the second panel of Fig. 2, $\sigma_n^2/2 \propto n$ for a range of intermediate scales (specifically, $4 \lesssim n \lesssim 9$), consistent with the lognormal model of intermittency developed for hydrodynamic turbulence (Kolmogorov 1962). The coefficient of proportionality is measured to be near 0.1, which implies a non-zero intermittency parameter $(2/\ln 2)0.1 \approx 0.29$ (Frisch 1995; Biskamp 2003), which is comparable to that measured in hydrodynamic turbulence (Sreenivasan & Kailasnath 1993) and standard MHD turbulence (Zhdankin et al. 2016b).

<figure><img src="image_3.png" /><figcaption>Figure 3. Skewness $S_n$ (blue) and excess kurtosis $K_n - 3$ (red) for the logarithm of coarse-grained energy dissipation rate, $\xi_n = \log(\epsilon_n/\epsilon_0)$, versus level $n$.</figcaption></figure>

The robustness of the lognormal fits can be characterized by measuring the skewness $S_n$ and kurtosis $K_n$ of the variable $\xi_n = \log(\epsilon_n/\epsilon_0)$, defined by

$$S_n = \frac{\langle (\xi_n - \bar{\xi}_n)^3 \rangle}{\langle (\xi_n - \bar{\xi}_n)^2 \rangle^{3/2}} \\ K_n = \frac{\langle (\xi_n - \bar{\xi}_n)^4 \rangle}{\langle (\xi_n - \bar{\xi}_n)^2 \rangle^2}, \quad (3)$$

where $\bar{\xi}_n = (\xi_n)$ and angled brackets indicate an ensemble average. For a normal distribution in $\xi_n$, which corresponds to a lognormal distribution in $\epsilon_n/\epsilon_0$, we expect that $S_n = K_n - 3 = 0$. We show the skewness $S_n$ and the excess kurtosis $K_n - 3$ in Fig. 3; since both quantities are significantly less than unity, the lognormal distribution is indeed a good fit to the overall data for energy dissipation rates. We note that the skewness increases from negative to positive in the inertial range, while the excess kurtosis decreases from positive to negative. We find that the quantities averaged for $n \ge 3$ are $\langle S_n \rangle_{n \ge 3} \approx 0.02$ and $\langle K_n \rangle_{n \ge 3} - 3 \approx -0.06$. Deviations from the lognormal distribution could be attributed to a number of factors, including corrections from a relatively short inertial range and alternative models of intermittency in MHD turbulence (e.g. the log-Poisson models considered in Chandran, Schekochihin & Mallet 2015; Mallet & Schekochihin 2017).

<figure><img src="image_4.png" /><figcaption>Figure 2. Left-hand panel: probability density functions for coarse-grained energy dissipation rate $\epsilon_n/\epsilon_0$ for $n \in \{3, 6, 9\}$ (in green, red and blue, respectively). Best-fitting lognormals are shown in black. Right-hand panel: scaling of $\sigma_n^2/2$ for best-fitting lognormal distributions versus level $n$. Results for total dissipation are shown in blue, for resistive dissipation in red and for viscous dissipation in magenta.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Top panel: fraction of total resistive dissipation (solid lines) and volume (dashed lines) in current density structures versus threshold $j_{thr}/j_{rms}$ for varying background fields ($B_0 = 0.03$ in blue, $B_0 = 0.1$ in red and $B_0 = 0.005$ in green). Bottom panel: exponential fits (black) to fraction of resistive dissipation in current density structures (blue) and fraction of viscous dissipation in vorticity structures (red) versus threshold. For comparison, the case for current density structures in standard MHD turbulence is also shown (dashed magenta line).</figcaption></figure>

We next describe the fraction of volume occupied by intermittent structures and their contribution to the overall energy dissipation, which is determined from cumulative distributions conditioned on thresholds in the current density and vorticity (denoted $j_{thr}$ and $\omega_{thr}$, respectively). In the first panel of Fig. 4, we show the resistive energy dissipation and volume filling fraction versus the normalized current density threshold $j_{thr}/j_{rms}$ (where $j_{rms}$ is the time-averaged rms current density). We find that the cumulative distributions are nearly identical for all three values of $B_0$. For large thresholds ($j_{thr}/j_{rms} \gg 1$), which are associated with intermittent current density structures, the fraction of energy dissipation declines exponentially, as $\exp(-\lambda j_{thr}/j_{rms})$ with $\lambda \approx 1/1.42$. We find that roughly 17 per cent of the resistive dissipation occurs in 1 per cent of the volume, and 50 per cent of the resistive dissipation in 7.5 per cent of the volume, indicating a high degree of inhomogeneity. In the second panel of Fig. 4, we show that the fraction of viscous dissipation (estimated by integrating $\nu\omega^2$) in vorticity structures at a threshold $\omega_{thr}/\omega_{rms}$ is also well fit by an exponential, with a slightly stronger decay constant, $\lambda \approx 1/1.17$. For comparison, we also show the resistive energy dissipation versus $j_{thr}/j_{rms}$ in simulations of standard homogeneous MHD turbulence with a background magnetic field

that is comparable to turbulent fluctuations (taken from the data set in Zhdankin et al. (2016b)). We find that the two cases agree very well. In contrast, we note that standard MHD turbulence with a strong guide field is significantly more intermittent, with a smaller decay constant $\lambda \approx 1/3$.3 (see, e.g. Zhdankin et al. 2016a).

We now arrive at our statistical analysis of intermittent dissipative structures. We apply the methodology described in Zhdankin et al. (2014) and mainly show the results for structures in the current density (i.e. current sheets). Each structure is defined as a cluster of connected points with current density magnitudes exceeding a fixed threshold, $|j| > j_{thr}$. For each current sheet, we measure the resistive energy dissipation rate $\mathcal{E} = \int dV \epsilon^\eta$ (integrated across the constituent points) and three characteristic scales: the length $L$, width $W$ and thickness $T$. The characteristic scales are measured in three orthogonal directions such that the length $L$ is the maximum distance between any two points in the structure, $W$ is the maximum distance between any two points in the plane orthogonal to $\hat{L}$ (the direction of $L$) coinciding with the point of maximum current density inside the structure, and $T$ is the distance across the structure in the direction $\hat{L} \times \hat{W}$, passing through the point of maximum current density inside the structure. The statistical results are insensitive to the threshold as long as $j_{thr} \gg j_{rms}$; we choose $j_{thr} = 4j_{rms}$. We show the current density profile in representative $yx$ and $xz$ slices of the simulation in Fig. 5, with current sheets exceeding a threshold of $2.5j_{rms}$ marked in red.

We show probability distributions for the current sheet measurements in Fig. 6; for comparison, we also show distributions for the same measurements in a reduced MHD (RMHD) simulation with $Re = 3200$ and strong guide field $B_0 \approx 5b_{rms}$ (taken from Zhdankin et al. 2014, 2016a). The distribution for energy dissipation rates $\mathcal{E}$ is consistent with a power law with a critical index of $-2.0$ across more than three decades, which implies that the populations of weak structures and strong structures both contribute equally to the overall energy dissipation. This critical index is also observed in the RMHD case. The distributions for lengths $L$ and widths $W$ also show power laws agreeing with RMHD, which has indices near $-3.2$ (the distribution for $W$ is possibly shallower, with index $-3.0$). Lengths and widths extend to scales comparable to the box size in the toroidal direction and the vertical direction, respectively. Thicknesses $T$ are peaked at small scales and do not exhibit a clear power law (although possibly a very steep one with index $\sim -5$). The distributions are very similar in all three simulations; in particular, the cases with $B_0 = 0.03$ and $0.01$ are nearly indistinguishable, while differences in the $B_0 = 0.005$ case (such as the shorter power law in the distribution for $W$) can be attributed to a shorter inertial range.

One implication of the above distributions is that energy dissipation is spread across current sheets spanning a broad range of scales. This is evident when plotting compensated dissipation-weighted distributions $E(X)X$ for the current sheet sizes $X \in \{L, W, T\}$, i.e. the distribution weighted by the energy dissipation rates and compensated by the given scale, which is shown in Fig. 7. The broad peaks of $E(L)L$ and $E(W)W$ reveal that energy dissipation is dominated by current sheets with lengths in the range of $0.1L_z \lesssim L \lesssim 2L_z$ and widths in the range of $0.02L_z \lesssim W \lesssim 0.3L_z$. Meanwhile, thicknesses are peaked at small scales, $T \sim 3 \times 10^{-3}L_z$. We interpret these lengths and widths to be the inertial range in the corresponding directions, and the thickness to be the dissipation scale. We note that although the typical thicknesses are only marginally resolved in this case, a similar distribution was found in the much better resolved RMHD simulations (Zhdankin et al. 2016a).
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Images of current density magnitude *j* in *yx* (azimuthal/radial) and *xz* (radial/vertical) slices of the simulations. Points exceeding a threshold of 2.5*j*<sub>rms</sub> are marked in red (note that a higher threshold is used in the analysis). Note the characteristic tilt of ~17.5 for structures in the *x*-*y* plane.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 6.** Probability distributions for resistive energy dissipation rates *E*, lengths *L*, widths *W* and thicknesses *T* of current sheets for $B_0 = 0.03$ (blue) and 0.005 (red). For comparison, distributions from RMHD simulations are shown in magenta.</figcaption></figure>

The various measurements performed on the current sheets often exhibit correlations. These can be characterized from 2D distributions, as shown in Fig. 8 for *L* versus *W* and for *E* versus *L*. For clarity, these 2D distributions are weighted by the dissipation rates. We find that *L* ~ *W* and *E* ~ *L*², similar to the scalings measured for RMHD (Zhdankin et al. 2016a). The thickness varies over a much smaller range and does not exhibit a clear scaling.

Unlike standard MHD turbulence, where structures are elongated along the mean magnetic field, the typical orientation of structures in MRI-driven turbulence is not *a priori* obvious. Therefore, we next

characterize the orientation of current sheets by measuring direction vectors for the sizes: $\hat{L}$, $\hat{W}$ and $\hat{T}$. We show scatterplots of the components of $\hat{L}$, $\hat{W}$ and $\hat{T}$ in the *x*-*y* plane in Fig. 9; for clarity, only intense current sheets exceeding an energy dissipation rate *E*/*E*<sub>tot</sub> > 4 × 10<sup>-5</sup> are shown. We find that the current sheets are typically extended in the azimuthal direction, but with a tilt of approximately *θ*<sub>tilt</sub> = 17.5. To better characterize the orientations, we utilize the tilted coordinates (*x'*, *y'*, *z*') obtained by rotating the original (*x*, *y*, *z*) coordinate system by the angle *θ*<sub>tilt</sub> counterclockwise around the *z*-axis. We show the components of $\hat{L}$, $\hat{W}$ and $\hat{T}$ in the tilted *x'*-*z*'
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Compensated, dissipation-weighted distributions $E(X)X$ for current sheet scales $X \in \{L, W, T\}$ (in blue, red and green, respectively).</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 8.** Top panel: 2D dissipation-weighted distributions showing current sheet lengths $L$ versus widths $W$. Bottom panel: similar for energy dissipation rates $\mathcal{E}$ versus $L$.</figcaption></figure>

plane in Fig. 10. Isotropy around the tilted y'-axis is evident from these scatterplots. Most of the current sheets have $\hat{L}_{y'} \approx \pm 1$, i.e. they are aligned with the tilt direction $\hat{y}' \approx \hat{y} \cos \theta_{\text{tilt}} - \hat{x} \sin \theta_{\text{tilt}}$, while $\hat{W}$ and $\hat{T}$ are statistically isotropic in the plane perpendicular to $\hat{y}'$. These orientations indicate that there may be an effective guide field in the $\hat{y}'$-direction. This tilt angle is somewhat larger than the angle associated with the angular momentum transport coefficient noted in previous studies (e.g. Simon et al. 2012). It, however, agrees well with distribution of the magnetic inclination

angle, $\theta_{\text{inc}} = \tan^{-1}(B_x/B_y)$, which we find to be strongly peaked near $\theta_{\text{inc}} \sim -18^\circ$ in our simulations.

In addition to the above analysis of current sheets, we performed a similar analysis for vorticity sheets by identifying structures that exceed a vorticity threshold, $|\omega| > \omega_{\text{thr}}$. The statistical results are essentially the same as for the current sheets, in agreement with the standard case (Zhdankin et al. 2016a). For example, we compare the distribution of viscous energy dissipation rates (estimated by integrating $\nu\omega^2$ across the structure) for vorticity structures to resistive dissipation rates of current sheets in Fig. 11. Vorticity structures tend to be somewhat weaker and smaller (with distributions for length and width being shifted to smaller scales by less than a factor of 2), but otherwise have quantitatively similar distributions. The orientations of the vorticity structures are also very similar to current sheets, being aligned with the same tilt axis.

## 4 CONCLUSIONS

In this work, we investigated the small-scale structure and intermittency of energy dissipation in a set of MHD simulations of MRI-driven turbulence. We found that MRI-driven turbulence leaves an imprint on small-scale quantities (e.g. current density, vorticity, dissipation rates) that is nearly indistinguishable from simulations of standard homogeneous MHD turbulence with a uniform guide field. This is in contrast to previous studies that focused on inertial-range quantities, such as the energy spectrum, and found little resemblance to standard turbulence (Fromang 2010; Lesur & Longaretti 2011), likely due to the effects of large-scale shear and broad-band energy injection from the MRI in a relatively small domain. Our results indicate that the small-scale structure is insensitive to these effects, and can be understood as standard MHD turbulence with a self-consistently generated guide field that is slightly tilted from the toroidal direction. Surprisingly, even the inertial-range lengths and widths of intermittent structures agree with the standard case. Other quantities such as the energy spectra may duplicate the results from standard MHD turbulence once simulations become sufficiently large.

If the above picture is correct, then the small-scale turbulence in accretion discs can be modelled with standard MHD turbulence phenomenology (Goldreich & Sridhar 1995; Boldyrev 2005). One implication is that the dynamics at kinetic scales, responsible for plasma heating and particle acceleration, will be insensitive to the large-scale accretion disc physics. For practical applications, future work is required to verify that these conclusions are insensitive to stratification, radiative effects and global geometry. It also remains to be shown that the results are unaffected by compressive fluctuations (Federrath et al. 2011, 2014), which may form shocks as an alternative route to intermittent dissipation (Schmidt, Federrath & Klessen 2008). For instance, it is known that the excitation of mildly subsonic MHD turbulence in Keplerian shear flows excites non-axisymmetric spiral density waves that are ultimately dissipated as shocks (Heinemann & Papaloizou 2009a,b).

In this study, we used numerical simulations corresponding to magnetic Prandtl number $Pm = 1$. This was done to maximize the inertial range of turbulence for the available numerical resolution. In our future work, we will address the cases of $Pm < 1$ and $Pm > 1$. In the latter case, MRI-driven turbulence may be excited solely by the non-linear magnetic dynamo action, i.e. without an imposed magnetic flux (e.g. Lesur & Longaretti 2007). Whether this affects the properties of MRI-driven turbulence is an interesting question.
<figure><img src="image_10.png" /><figcaption>**Figure 9.** Scatterplots of the components of $\hat{\mathcal{L}}$, $\hat{\mathcal{W}}$ and $\hat{\mathcal{T}}$ in the $x-y$ plane (for intense current sheets with $\mathcal{E}/\mathcal{E}_{\text{tot}} > 4 \times 10^{-5}$). A characteristic tilt of about 17.5, counterclockwise, is evident.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 10.** Scatterplots of the components of $\hat{\mathcal{L}}$, $\hat{\mathcal{W}}$ and $\hat{\mathcal{T}}$ in the tilted $x'-z'$ plane (for intense current sheets with $\mathcal{E}/\mathcal{E}_{\text{tot}} > 4 \times 10^{-5}$). Isotropy around the (out-of-plane) $y'$-axis is evident.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 11.** Distribution of viscous energy dissipation rates for vorticity sheets compared to resistive energy dissipation rates for current sheets.</figcaption></figure>

## ACKNOWLEDGEMENTS

We thank Dmitri Uzdensky and Mitch Begelman for fruitful discussions and for comments on the manuscript. VZ acknowledges support from NSF grant AST-1411879. JW is supported by the Wisconsin Alumni Research Foundation at the University of Wisconsin-Madison. SB is partly supported by NSF grant AGS-1261659, and by the Vilas Associates Award from the University of Wisconsin-Madison. Simulations were performed at the Texas Advanced Computing Center (TACC) at the University of Texas at Austin under the NSF-Teragrid Project TG-PHY110016.

## REFERENCES

* Baganoff F. et al., 2001, Nature, 413, 45

* Bai X.-N., Stone J. M., 2013, ApJ, 767, 30

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

* Balbus S. A., Hawley J. F., 1998, Rev. Mod. Phys., 70, 1

* Biskamp D., 1995, Chaos Solitons Fractals, 5, 1779

* Biskamp D., 2003, Magnetohydrodynamic Turbulence. Cambridge Univ. Press, Cambridge

* Bodo G., Mignone A., Cattaneo F., Rossi P., Ferrari A., 2008, A&A, 487, 1

* Bodo G., Cattaneo F., Ferrari A., Mignone A., Rossi P., 2011, ApJ, 739, 82

* Boldyrev S., 2005, ApJ, 626, L37

* Boldyrev S., 2006, Phys. Rev. Lett., 96, 115002

* Brandenburg A., Nordlund A., Stein R. F., Torkelsson U., 1995, ApJ, 446, 741

* Chandran B. D., Schekochihin A. A., Mallet A., 2015, ApJ, 807, 39

* Federrath C., Chabrier G., Schober J., Banerjee R., Klessen R. S., Schleicher D. R., 2011, Phys. Rev. Lett., 107, 114504

* Federrath C., Schober J., Bovino S., Schleicher D. R., 2014, ApJ, 797, L19

* Fleming T. P., Stone J. M., Hawley J. F., 2000, ApJ, 530, 464

* Frisch U., 1995, Turbulence: The Legacy of AN Kolmogorov. Cambridge Univ. Press, Cambridge

* Fromang S., 2010, A&A, 514, L5
Fromang S., Papaloizou J., Lesur G., Heinemann T., 2007, A&A, 476, 1123

Goldreich P., Sridhar S., 1995, ApJ, 438, 763

Heinemann T., Papaloizou J., 2009a, MNRAS, 397, 52

Heinemann T., Papaloizou J., 2009b, MNRAS, 397, 64

Kolmogorov A. N., 1962, J. Fluid Mech., 13, 82

Kunz M. W., Stone J. M., Quataert E., 2016, Phys. Rev. Lett., 117, 235101

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

Lesur G., Longaretti P.-Y., 2011, A&A, 528, A17

McClintock J. E., Remillard R. A., 2006, in Lewin W., van der Klis M., eds, Compact Stellar X-ray Sources. Cambridge Univ. Press, Cambridge, p. 157

McNally C. P., Hubbard A., Yang C.-C., Mac Low M.-M., 2014, ApJ, 791, 62

Mallet A., Schekochihin A., 2017, MNRAS, 466, 3918

Markoff S., Falcke H., Yuan F., Biermann P. L., 2001, A&A, 379, L13

Mason J., Perez J. C., Boldyrev S., Cattaneo F., 2012, Phys. Plasmas, 19, 055902

Matthaeus W., Wan M., Servidio S., Greco A., Osman K., Oughton S., Dmitruk P., 2015, Phil. Trans. R. Soc. A, 373, 20140154

Merrifield J. A., Müller W.-C., Chapman S. C., Dendy R. O., 2005, Phys. Plasmas, 12, 022301

Murphy G. C., Pessah M. E., 2015, ApJ, 802, 139

Nauman F., Blackman E. G., 2014, MNRAS, 441, 1855

Perez J. C., Mason J., Boldyrev S., Cattaneo F., 2012, Phys. Rev. X, 2, 041005

Salvesen G., Simon J. B., Armitage P. J., Begelman M. C., 2016, MNRAS, 457, 857

Schmidt W., Federrath C., Klessen R., 2008, Phys. Rev. Lett., 101, 194505

Shi J.-M., Stone J. M., Huang C. X., 2016, MNRAS, 456, 2273

Simon J. B., Hawley J. F., Beckwith K., 2011, ApJ, 730, 94

Simon J. B., Beckwith K., Armitage P. J., 2012, MNRAS, 422, 2685

Sorathia K. A., Reynolds C. S., Stone J. M., Beckwith K., 2012, ApJ, 749, 189

Sreenivasan K., Kailasnath P., 1993, Phys. Fluids A, 5, 512

Walker J., Lesur G., Boldyrev S., 2016, MNRAS, 457, L39

Wan M., Matthaeus W., Roytershteyn V., Parashar T., Wu P., Karimabadi H., 2016, Phys. Plasmas, 23, 042307

Zhdankin V., Uzdensky D. A., Perez J. C., Boldyrev S., 2013, ApJ, 771, 124

Zhdankin V., Boldyrev S., Perez J. C., Tobias S. M., 2014, ApJ, 795, 127

Zhdankin V., Uzdensky D. A., Boldyrev S., 2015a, Phys. Rev. Lett., 114, 065002

Zhdankin V., Uzdensky D. A., Boldyrev S., 2015b, ApJ, 811, 6

Zhdankin V., Boldyrev S., Uzdensky D. A., 2016a, Phys. Plasmas, 23, 055705

Zhdankin V., Boldyrev S., Chen C. H. K., 2016b, MNRAS, 457, L69

This paper has been typeset from a TeX/Latex file prepared by the author.