ON THE ANISOTROPIC NATURE OF MRI-DRIVEN TURBULENCE IN ASTROPHYSICAL DISKS

The equations governing the local dynamics of a magnetized disk in the shearing box framework are given by (Hawley et al. 1995; Brandenburg et al. 1996)

$$\begin{eqnarray}\begin{array}{rcl} \frac{\partial {\boldsymbol{v}} }{\partial t}+{\boldsymbol{v}} \cdot \nabla {\boldsymbol{v}} & = & -2{\Omega }\times {\boldsymbol{v}} +2q{{{\Omega }}^{2}}x \\ {} & {} & -\frac{\nabla p}{\rho }+\frac{{\boldsymbol{J}} \times {\boldsymbol{B}} }{\rho }+\nu {{\nabla }^{2}}v, \\ \end{array}\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&\frac{\partial {\boldsymbol{B}} }{\partial t}=\nabla \times \left( {\boldsymbol{v}} \times {\boldsymbol{B}} \right)+\eta {{\nabla }^{2}}{\boldsymbol{B}} .\end{eqnarray} \tag{ 2 }$$

Here ρ is the mass density, ${\boldsymbol{v}}$ is the fluid velocity in the rotating frame, with angular frequency ${\Omega }={\Omega }{\boldsymbol{\hat{z}}}$, p is the pressure determined through an equation of state, and ${\boldsymbol{B}}$ is the magnetic field, with $\nabla \cdot {\boldsymbol{B}} =0$. The current density is ${\boldsymbol{J}} \equiv \nabla \times {\boldsymbol{B}} /{{\mu }_{0}}$, with ${{\mu }_{0}}$ a constant dependent on the unit system adopted. The first and second terms on the right-hand side of Equation (1) account for the Coriolis acceleration and the first-order expansion of the tidal potential, respectively. The shear parameter is defined as

$$\begin{eqnarray}&&q=-\frac{d{\rm ln} {\Omega }}{d{\rm ln} r},\end{eqnarray} \tag{ 3 }$$

so that q = 3/2 for a Keplerian disk. The coefficients ν and η stand for the viscosity and resistivity, assumed to be constant.

Let us consider an incompressible flow, which is a good approximation provided that the magnetic field involved is very weak, i.e., the magnetic energy density is small compared to its thermal counterpart. The set of Equations (1)–(2) admits exact, MRI-unstable solutions of the form

$$\begin{eqnarray}&&{\boldsymbol{v}} =-q{\Omega }x\check{{\boldsymbol{y}} }+{{{\boldsymbol{V}} }_{0}}{\rm sin} (Kz)\;{{e}^{{\Gamma }t}},\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&{\boldsymbol{B}} ={{\bar{B}}_{z}}\check{{\boldsymbol{z}} }+{{{\boldsymbol{B}} }_{0}}{\rm cos} (Kz)\;{{e}^{{\Gamma }t}}.\end{eqnarray} \tag{ 5 }$$

The ratio of the amplitudes, i.e., ${{B}_{0}}/{{V}_{0}}$, sets the initial amplitude of the MRI mode. In ideal MHD, the wavelength of the fastest mode is ${{\lambda }_{{\rm max} }}=2\pi /{{K}_{{\rm max} }}=\sqrt{15/16}\;{{\bar{v}}_{A,z}}/{\Omega }$, where ${{\bar{v}}_{A,z}}={{\bar{B}}_{z}}/\sqrt{4\pi \rho }$ is the Alfvén speed associated with the background magnetic field, ${{\bar{B}}_{z}}$. This mode grows at the rate ${{{\Gamma }}_{{\rm max} }}=3{\Omega }/4$, with ${{V}_{A,0}}/{{V}_{0}}=\sqrt{5/3}$, where ${{V}_{A,0}}={{B}_{0}}/\sqrt{4\pi \rho }$ is the Alfvén speed associated with the MRI-driven perturbation (Pessah et al. 2006).

These exact primary MRI modes are themselves, in turn, subject to secondary, parasitic instabilities. The fastest-growing parasitic modes, which are associated with KH modes, have the same vertical periodicity as the primary MRI mode itself (Goodman & Xu 1994; Pessah & Goodman 2009). As a prelude to analyzing in detail the nonlinear evolution of the MRI modes as they are affected by secondary instabilities, we perform a series of simpler numerical experiments to study the dynamics of the KH instabilities resulting from a periodic velocity profile with constant amplitude.

We begin by examining the development of KH modes that results from a background flow with $\delta {{v}_{x}}={{V}_{0}}{\rm sin} (2\pi z/{{L}_{z}})$. This velocity profile is akin to the velocity profile induced by an MRI mode, except for the fact that its amplitude is time independent. Incidentally, note that a time-independent MRI amplitude is usually invoked when calculating secondary instabilities affecting an MRI mode (Goodman & Xu 1994; Latter et al. 2009; Pessah 2010). The argument behind this approximation is that, provided that the amplitude of the MRI is large enough, the growth rate of the secondary instability should be much larger than the MRI growth rate. This preliminary exercise is useful for two reasons. It provides an independent check for the semi-analytical results predicting the parasitic growth rates (Goodman & Xu 1994) and mode structure (Pessah 2010). It allows us to get acquainted with the type of mode structure that emerges as a natural consequence of the action of parasitic modes in 3D shearing-box simulations of the MRI discussed below.

We solve the set of equations describing a hydrodynamic flow starting from the initial conditions

$$\begin{eqnarray}&&{\boldsymbol{v}} (x,z)=\left( \begin{array}{ccccccccccccccc} {{V}_{0}}{\rm sin} (Kz) \\ 0 \\ \epsilon {{V}_{0}}{\rm sin} ({{k}_{h}}{{x}_{h}}) \\ \end{array} \right),\end{eqnarray} \tag{ 6 }$$

where $\epsilon ={{10}^{-2}}$. The stability of this flow can be analyzed semi-analytically, as described in Goodman & Xu (1994) and Pessah (2010). The flow is unstable with a dimensionless growth rate $s/K{{V}_{0}}$ that depends only on ${{k}_{h}}/K$, i.e., the ratio of the wavenumbers associated with the periodic background shear flow and the excited perturbation. This growth rate attains a maximum at ${{k}_{h}}/K\simeq 0.6$ and vanishes beyond ${{k}_{h}}/K=1$. This is illustrated in Figure 1, where we compare the results obtained with linear theory and numerical codes. In order to gauge the importance of compressibility in the growth of the KH instabilities, we performed two separate sets of simulations in compressible and incompressible gas dynamics using the publicly available codes PLUTO (Mignone et al. 2012; where the equation of state is isothermal) and SNOOPY (Lesur & Longaretti 2009). As expected, the results of the spectral code SNOOPY are closer to the analytical results than the finite-volume code PLUTO; the difference is accounted for by compressibility.

Zoom In Zoom Out Reset image size

The early evolution of this sinusoidal, 2D flow leads to the development of sheets of cellular vortices of alternating polarities in both radius and height, as shown in Figure 2. The agreement between the analytical and numerical results in the linear regime is also excellent in terms of the structures formed. The extremum of the vorticity distribution is due to a combination of circulation and shear vorticity, so the peak of the color map does not coincide with the point around which the fluid seems to be circulating. Figure 3 shows the evolution of the instability into the nonlinear regime, where the vortices interact, grow, and merge.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The saturated state, depicted in the rightmost panel of Figure 3, consists of episodic fluctuations in the vertical component of the velocity, which leads to a decrease of v_z² to about 50% of its initial value. The radial component of the velocity also shows episodic fluctuations, which leads v_x² to increase to 80% of the initial kinetic energy (Figure 3). At later times, shock waves form, fine sheared structures appear, and kinetic energy is dissipated into thermal energy. Note that while the initial shear flow is aligned in the x-direction, the KH instability increases the flow speed in the transverse, z-direction, until it reaches the same order of magnitude, so that initially anisotropic flow becomes partially isotropized.

In the standard picture, parasitic instabilities are thought to feed off the exponentially growing background provided by the primary MRI mode, which corresponds to an exact solution of the MHD equations provided that compressibility is unimportant. In ideal MHD, the growth rate of these secondary instabilities is expected to be proportional to the amplitude of the MRI mode, and thus the growth of the parasites is superexponential. Therefore, the energy obtained by the secondary instabilities can become comparable to the energy in the MRI mode from which it feeds very quickly. Because the energy going into secondary instabilities is envisioned to be directly supplied by the primary mode, the fast growth of the secondary could potentially prevent the MRI from continuing growing, halting its exponential growth and leading thus to its saturation.

This picture assumes that the amplitude of the MRI is large enough so that several approximations can be invoked. Perhaps more importantly, it also assumes that the dynamics of the secondary perturbations is insensitive to the shearing background and that the MRI grows uninhibited. The first approximation implies that the parasitic perturbations can be described by a time-independent wavevector ${{{\boldsymbol{k}} }_{h}}=({{k}_{x}},{{k}_{y}})$, which characterizes the wavelength of the perturbations in horizontal planes, i.e., perpendicular to z-direction. However, because of advection by the background shear flow, the wavevectors of nonaxisymmetric perturbations will evolve in time according to ${{k}_{x}}(t)={{k}_{x}}-q{\Omega }t\;{{k}_{y}}$. The second approximation prevents any dynamic coupling between the primary mode and the secondary parasites. This approximation becomes increasingly worse as the amplitude of the parasite grows larger. In order to understand the role of secondary instabilities and how they mediate the transfer of power from axisymmetric MRI modes into nonaxisymmetric structures, we set up a tailored numerical experiment that is not limited by the assumptions usually invoked in the semi-analytical studies described above. For the sake of simplicity, we focus on following the dynamical evolution of the fastest-growing MRI mode.

The incompressible simulations are carried out using the SNOOPY spectral code (Lesur & Longaretti 2009), which solves the MHD equations in the constant-density regime and incorporates a finite physical resistivity and viscosity. The resistivity is such that ${{{\Lambda }}_{\eta }}\equiv \bar{v}_{Az}^{2}/\eta {\Omega }={{10}^{4}}$, and the magnetic Prandtl number is unity. For practical purposes, this is rather close to ideal MHD and corresponds to the regime where the secondary KH mode is expected to dominate.

We set up the numerical experiment as follows. In order to maximize the numerical resolution, we set the vertical extent of the simulation domain to be equal to the most unstable MRI wavelength, so that ${{L}_{z}}=2\pi {{K}_{{\rm max} }}=\sqrt{15/16}\;{{\bar{v}}_{A,z}}/{\Omega }$. This motivates choosing L_z and ${{{\Omega }}^{-1}}$ as the units of length and time. All the results in the paper correspond to a simulation with a resolution of 256 × 256 × 128 cells, which was carried out for $260\;{{{\Omega }}^{-1}}$. The only exception are the results where we used the lower resolution of 128 × 128 × 64 cells in order to run the simulation for $400\;{{{\Omega }}^{-1}}$.

We initialize the simulation by exciting the fastest-growing MRI mode at t = 0, i.e.,

$$\begin{eqnarray}&&{\boldsymbol{v}} (x,y,z,t=0)=\delta v\left( \begin{array}{ccccccccccccccc} {\rm sin} (Kz) \\ {\rm sin} (Kz) \\ 0 \\ \end{array} \right),\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&{\boldsymbol{B}} (x,y,z,t=0)=\delta B\left( \begin{array}{ccccccccccccccc} {\rm cos} (Kz) \\ -{\rm cos} (Kz) \\ 0 \\ \end{array} \right),\end{eqnarray} \tag{ 8 }$$

where $\delta v={{10}^{-3}}{{L}_{z}}{\Omega }$ and $\delta B/\sqrt{4\pi \rho }=\sqrt{5/3}\;\delta v$. As expected, this initial perturbation grows exponentially at a rate ${{{\Gamma }}_{{\rm max} }}=3/4{\Omega }$ (Figure 4) with the mode structure predicted by linear theory. In incompressible dynamics this exponential growth can continue for several orbits unless the fluid is perturbed.

Zoom In Zoom Out Reset image size

In order to break down the exponential growth of the MRI, we seed a perturbation with a horizontal wavevector given by ${{{\boldsymbol{k}} }_{h}}=(1,1,0)$.¹ Note that this wavevector is chosen in the direction where the shear induced by the MRI is fastest, i.e., at 45° with respect to the radial direction. Because this mode is nonaxisymmetric, its excitation requires some care for the following reason. If a nonaxisymmetric perturbation is excited at time t = 0, since the radial wavenumber has a linear dependence on time, it reaches the Nyquist wavenumber in ${{t}_{{\rm alias}}}={{{\rm log} }_{2}}({{N}_{x,g}})\;{{{\Omega }}^{-1}}$ (where ${{N}_{x,g}}$ is the number of grid cells in the x-direction), where it is aliased to k_x = 1 before it can grow significantly (Lyra & Klahr 2011). To circumvent this, we excite the perturbation

$$\begin{eqnarray}&&\delta {\boldsymbol{v}} (x,y,z,t=6)=\epsilon \;\delta v{\rm sin} ({{{\boldsymbol{k}} }_{h}}\cdot {\boldsymbol{x}} )\;\check{{\boldsymbol{z}} },\end{eqnarray} \tag{ 9 }$$

with $\epsilon ={{10}^{-4}}$ at time $t=6{{{\Omega }}^{-1}}$. This allows the perturbation to seed the parasitic mode, avoiding the effects of aliasing due to finite grid resolution. This secondary instability grows superexponentially (Figure 4), exhibiting the characteristics expected of a KH instability feeding off the velocity shear induced by the primary MRI mode (see Figure 5 and Pessah 2010). The subsequent evolution of this 3D, nonaxisymmetric perturbation, in the plane parallel to ${{{\boldsymbol{k}} }_{h}}$, i.e., where the MRI-induced velocity shear is largest, is shown in Figure 6. Note that the rightmost panel there shows a close-up of Figure 4 at the time where the amplitude of the secondary instability approaches that of the MRI.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In order to analyze the dynamics in detail, it is useful to divide the evolution of the flow into three stages: (i) an exponential growth phase driven by the MRI in which the amplitude of the fast-growing parasite is too small to modify the primary mode; (ii) a transitionstage in which the energy of the parasite is comparable to the energy of the primary MRI mode, which starts to depart significantly from the initial exact solution; and (iii) a turbulent state in which all flow variables fluctuate in time. In what follows, we describe briefly each of these stages.

The initial conditions in Equations (7) and (8) lead to exponential growth of the exact solution that grows for several ${{{\Omega }}^{-1}}$ timescales. The amplitude of this perturbation grows by more than two orders of magnitude, retaining the initial MRI mode structure in which the sinusoidal velocity field, with $\delta {{v}_{x}}=\delta {{v}_{y}}$, is orthogonal to the co-sinusoidal magnetic field perturbations, with $\delta {{B}_{x}}=\delta {{B}_{y}}$. This is illustrated in the leftmost panel of Figure 7, which shows the magnitude of the radial and azimuthal components of the velocity and magnetic field (Pessah & Chan 2008) as a scatter plot (see, e.g., Bodo et al. 2008).

Zoom In Zoom Out Reset image size

In the incompressible regime under consideration, the exact MRI mode does not lead to velocity or magnetic field perturbations in the vertical direction. Therefore, keeping track of $\delta {{v}_{z}}$ and $\delta {{B}_{z}}$ allows us to quantify how the evolution of the flow departs from the exact MRI mode excited at t = 0. Figure 4 illustrates that the vertical component of the magnetic field, which gives rise to the MRI, remains constant until the secondary perturbation (9) is excited at $t=6\;{{{\Omega }}^{-1}}$. The initial perturbation seeding $\delta {{v}_{z}}$ grows superexponentially together with $\delta {{B}_{z}}$, affecting all other field components, which start to depart significantly from the MRI behavior predicted by linear theory soon after. This is also evident from the middle and rightmost panels in Figure 7.

The shear resulting from the exponentially growing sinusoidal velocity field that results from the MRI is maximum along the direction x = y, which is the direction in which we have chosen ${{{\boldsymbol{k}} }_{h}}$. This plays an important role in the 3D, nonaxisymmetric nature of the parasitic modes. The evolution of the parasitic mode is illustrated in Figure 6. The first four panels, from left to right, show the projected component of the vorticity into the $({{x}_{h}},z)$ plane, with ${{x}_{h}}={\boldsymbol{x}} \cdot {{{\boldsymbol{k}} }_{h}}$. The initial vorticity perturbation, which is of the form $\nabla \times {{{\boldsymbol{v}} }_{h}}\propto {\rm cos} ({{{\boldsymbol{k}} }_{h}}\cdot {\boldsymbol{x}} )\;{{{\boldsymbol{\hat{k}}} }_{p}}$, with ${{{\boldsymbol{\hat{k}}} }_{p}}={{{\boldsymbol{\hat{k}}} }_{z}}\times {{{\boldsymbol{\hat{k}}} }_{h}}$, gives rise to a sheet of cellular vortices, as predicted in Pessah (2010). The first and second panels, from left to right, show that the cellular vortices begin to interact in the nonlinear stage and vorticity concentrates near the boundary of the vortex sheets, at $Kz=0,\pm 0.5$. The sheets of vortices can be seen in the 3D slices of vorticity and velocity field lines projected into the $({{x}_{h}},z)$ plane (see Figure 5), where the familiar parasite structure exhibits extrema visible at $z=-0.25,0.25$.

The vorticity plots in the third and fourth panels of Figure 6 illustrate how the interface undulations roll over into the classical cat's eye associated with KH vortices, which appears between the two shearing layers, similar to those seen in Frank et al. (1996), Jones et al. (1997), and Malagoli et al. (1996). At this time (6.3–6.4 ${{{\Omega }}^{-1}}$) the peak value of the magnetic field B_z exceeds the peak value of the vertical velocity component, v_z. The radial component of the magnetic field has also started to depart from its linear growth and exceeds the MRI growth rate, as can be seen from Figure 4. As a result of the superexponential growth of the KH instability, the flow becomes completely disrupted and transitions rapidly into a turbulent state.

It is instructive to compare the evolution of the secondary KH instabilities that feed off the exponentially growing sinusoidal shear provided by the MRI and the constant, imposed sinusoidal shear addressed in Section 3.2. Because the growth of the MRI mode induces magnetic field perturbations that are orthogonal to the velocity fluctuations, the nature of the motions in the plane $({{x}_{h}},z)$ is essentially hydrodynamic. Therefore, it is not surprising that the perturbations (9) lead to the excitation of KH instabilities similar to the ones described in Section 3.2. However, because the amplitude of the velocity shear off of which the secondary instability feeds is growing exponentially in time, their late temporal evolution exhibits significant differences with respect to the KH instabilities triggered by the constant-amplitude sinusoidal profile (compare Figures 3 and 6).

In the linear stage of the MRI, the flow is strongly anisotropic, as the vertical components of the magnetic field and velocity are greatly exceeded by their radial and axial counterparts. The departure from the linear state can be seen in the scatter plots of ${{B}_{x,y}}$ and ${{v}_{x,y}}$ (see Figure 7); the second panel shows the departure from the linear mode at $t=6.5\;{{{\Omega }}^{-1}}$, shortly after the perturbation was excited. The magnetic field components B_x and B_y no longer increase exponentially, but instead a parasite mode starts to take over. The ${{B}_{x}},{{B}_{y}}$ structures, which grow at right angles to the initial MRI flow, at $({{B}_{x}},{{B}_{y}})$ = (−2, 1) and (2, −1) in the second panel, then correspond to the growth of parasitic modes. At saturation time, all three components of the velocity and magnetic fields are approaching the same magnitude. In essence, the growth of parasite modes has a strong effect on the anisotropy of the flow, reducing it and attempting to return to isotropy, as can also be seen in Figure 7.

The MRI mode starts to saturate as energy is drained from the radial and azimuthal components of velocity and magnetic fields and fed into the vertical components of velocity and magnetic fields. Thus, a reduction in the strong anisotropy, from ${{v}_{x}}\sim {{v}_{y}}\gg {{v}_{z}}$ to ${{v}_{x}}\sim {{v}_{y}}\sim {{v}_{z}}$, appears to accompany the parasite modal growth.

After the parasite mode has grown to its maximum value at $t=6.6\;{{{\Omega }}^{-1}}$, and several $\;{{{\Omega }}^{-1}}$ have elapsed, a state of fully developed, anisotropic turbulence is reached. The analysis of the MRI-driven turbulent state has been the subject of a large body of literature. Here we point out that the turbulent regime observed in our numerical experiment is characterized by sequences of peaks in the v_x component followed by peaks in v_z. These quasi-periodic oscillations are evident in Figure 8. The radial component v_x always increases more slowly, followed by a rapid increase in the vertical component v_z. As both components attain their respective local maxima, they decrease in tandem and at the same rate. Interestingly, this echoes the behavior seen in the linear stage, where the growth in v_x and v_y (driven by the MRI mode) triggers the growth in v_z (resulting from the excitation of parasite modes), which then prevents the amplitude of the MRI from growing further.

In order to quantify this periodicity, we analyze the time series provided by v_x(t). We consider a total of $350\;{{{\Omega }}^{-1}}$, discarding the initial 50$\;{{{\Omega }}^{-1}}$ in which the simulation is growing exponentially and then transitioning toward fully turbulent flow. We use Bartlett's method of estimating power spectra via averaging periodograms (Bartlett 1948) in order to reduce their variance. This is achieved by splitting the time series for ${{v}_{x}}(50\lt t\lt 300)$ into four independent samples, calculating the Fourier amplitudes of these time series, and averaging their square for each frequency (see Figure 8). The bulk of the power is concentrated between 9 and $13\;{{{\Omega }}^{-1}}$. A spike is also seen at two-thirds of a period, which corresponds to one shearing time.

In the saturated state, the flow is characterized by highly anisotropic turbulence, which exhibits quasi-episodic bursts. Bodo et al. (2008) attributed the presence (absence) of these bursts to the absence (presence) of parasite modes. In the linear stage, it has been possible to directly subtract the background MRI exact solution, in order to identify signatures and compare with those predicted in semi-analytical calculations (Pessah 2010). However, in the turbulent state, the same approach becomes more difficult. Exact solutions are not easy to identify or subtract in a consistent way. Therefore, we pursue an alternative approach, in order to engage directly with the fluctuations in the Reynolds and Maxwell stresses that are associated with angular momentum transport in the turbulent state.

In order to obtain a better understanding of the anisotropy of the flow and its temporal evolution, we employ the turbulent stress invariant analysis introduced by Lumley & Newman (1977) in the context of hydrodynamics and apply it to study MRI-driven turbulence.

The Reynolds and Maxwell stress tensors are defined according to

$$\begin{eqnarray}&&{{R}_{ij}}=\rho \;\langle \delta {{v}_{i}}\;\delta {{v}_{j}}\rangle ,\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}&&{{M}_{ij}}=\frac{\langle \delta {{B}_{i}}\;\delta {{B}_{j}}\rangle }{4\pi },\end{eqnarray} \tag{ 11 }$$

where the brackets represent an appropriate spatial average. Each of these real and symmetric tensors can be divided into a deviator and an isotropic part. By normalizing the tensors using their trace, associated with the corresponding energy densities, and then subtracting the isotropic component, we arrive at a single tensor quantity that parametrizes the deviation from isotropy. The Reynolds stress anisotropy tensor (Lumley & Newman 1977) is defined as

$$\begin{eqnarray}&&{{\mathcal{R}}_{ij}}=\frac{{{R}_{ij}}}{{\rm Tr}({{R}_{ij}})}-\frac{1}{3}{{\delta }_{ij}},\end{eqnarray} \tag{ 12 }$$

where Tr stands for the trace and ${{\delta }_{ij}}$ is the Kronecker delta. The Reynolds stress anisotropy tensor provides an important diagnostic and has been used extensively in numerical and experimental studies of hydrodynamic turbulence (Biferale & Procaccia 2005). Here we extend this idea in order to study MRI-driven turbulence and define the Maxwell stress anisotropy tensor as

$$\begin{eqnarray}&&{{\mathcal{M}}_{ij}}=\frac{{{M}_{ij}}}{{\rm Tr}({{M}_{ij}})}-\frac{1}{3}{{\delta }_{ij}}.\end{eqnarray} \tag{ 13 }$$

The evolution of the global distribution of anisotropy when parasitic modes grow by feeding off the MRI can be quantified via ${{\mathcal{R}}_{ij}}$ in Equation (10) by averaging over cubic volumes composed of eight cells. At a given time, this procedure leads to a set of ellipsoidal glyphs, as shown in Figure 9 for three different time snapshots. The orientation and shape of each glyph are controlled by the eigenvectors of $\boldsymbol{\mathcal{R}}$ and their corresponding eigenvalues. The latter, defined in descending order as ${{\lambda }_{1}},{{\lambda }_{2}}$, and ${{\lambda }_{3}}$, correspond to the principal, medium, and minor axes of the ellipsoidal glyph. The direction of the glyph is aligned along the principal stress eigenvector, i.e., the eigenvector with the largest associated eigenvalue. The leftmost panel shows that all the tensors are aligned in the direction along which the MRI-induced shear is strongest, i.e., x = y, indicating that in the late linear stage of the MRI we have highly organized and highly anisotropic flow. The middle panel shows $\boldsymbol{\mathcal{R}}$, after the secondary perturbation has been excited and the parasites have started to grow. The initially highly anisotropic flow has now been partially isotropized, and the principal stress eigenvectors of $\boldsymbol{\mathcal{R}}$ are no longer parallel over the domain. The rightmost panel shows the partially isotropized field in saturation. In this last panel, the system has not yet reached fully developed turbulence, but it has reached the closest approach to isotropy.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In order to probe how the anisotropy on larger scales varies as a function of time in a quantitive way, it is useful to employ invariant quantities. For a traceless tensor, ${{\lambda }_{1}}+{{\lambda }_{2}}+{{\lambda }_{3}}=0$, the degree of anisotropy can be characterized by calculating the two nonzero invariants

$$\begin{eqnarray}&&\chi _{2}^{\mathcal{R}}=-(\lambda _{1}^{2}+\lambda _{2}^{2}+\lambda _{3}^{2}),\end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray}&&\chi _{3}^{\mathcal{R}}={{\lambda }_{1}}\;{{\lambda }_{2}}\;{{\lambda }_{3}}.\end{eqnarray} \tag{ 15 }$$

By plotting all the possible values of these two independent invariants of the Reynolds stress in the plane spanned by $({{\chi }_{2}},{{\chi }_{3}})$, Lumley & Newman (1977) found that every realizable Reynolds stress is constrained to lie within three loci that intersect at three vertices. The enclosed deltoid, known as the Lumley triangle, represents all possible values for $\boldsymbol{\mathcal{R}}$. The upper right vertex corresponds to the case of ${{\lambda }_{1}}\gt {{\lambda }_{2}},{{\lambda }_{3}}$, and the origin corresponds to isotropy, ${{\lambda }_{1}}={{\lambda }_{2}}={{\lambda }_{3}}$. The upper left vertex corresponds to two-component turbulence, with ${{\lambda }_{1}},{{\lambda }_{2}}\gt {{\lambda }_{3}}$. The vertical axis in these plots corresponds to the deviation from isotropy, from zero (fully isotropic) through 0.083 (two-component isotropic) to 0.333 (one-component). The horizontal axis in the plots correspond to the determinant, which distinguishes between those cases that have the same magnitude of deviation but a different sign.

We provide a quantitative representation of the time evolution of anisotropy on larger scales, by showing in Figure 10 how the Lumley triangles associated with the Reynolds and Maxwell stresses are populated as a function of time. The plots are constructed by taking contours of a 2D histogram of the distribution of vertically averaged tensors. The left column shows the Reynolds stress tensor invariant map at five successive times. As time progresses, the distribution of anisotropy fluctuates, initially increasing between $109\lt t\;{{{\Omega }}^{-1}}\lt 113$ and later decreasing between $113\lt t\;{{{\Omega }}^{-1}}\lt 117$. The left column of Figure 10 shows the normalized Reynolds stress tensor invariant map, with the second invariant in the vertical axis and the third invariant (determinant) in the horizontal axis. The color map shows the density of points that have associated the pair of values $({{\chi }_{2}},{{\chi }_{3}})$ when each stress tensor is averaged in the vertical direction. At $t=119\;{{{\Omega }}^{-1}}$, the stress tensors have a peak that is far from the origin. The highest density (shown in red) is offset from the origin, which is consistent with anisotropic turbulence. The majority of points are distributed close to the 1D–3D locus, with almost no points on the 1D–2D locus. As time progress, at $t=111\;{{{\Omega }}^{-1}}$ and $t=113\;{{{\Omega }}^{-1}}$ the distribution of anisotropy shifts along the 1D–3D locus, further from isotropy (i.e., the origin). At $t=113\;{{{\Omega }}^{-1}}$, the peak of anisotropy is reached and almost no power is present near the origin. At $t=115\;{{{\Omega }}^{-1}}$, the distribution of invariants shifts back toward the origin, indicating a return toward isotropy. The peak of the distribution remains far from the origin, however. At $t=117\;{{{\Omega }}^{-1}}$, the distribution shifts even further toward the origin.

Zoom In Zoom Out Reset image size

This effect is even more pronounced in the Maxwell stress tensor (right column of Figure 10), which shows a strong increase in anisotropy between $109\lt t\;{{{\Omega }}^{-1}}\lt 113$, which then reverses between $113\lt t\;{{{\Omega }}^{-1}}\lt 117$. The examination of other time frames (not reproduced here) reveals that this pattern of increase and decrease in anisotropy is quasi-periodic, in agreement with Figure 8.

A number of works have studied the Fourier spectrum of MRI-driven turbulence (for early studies see, e.g., Hawley et al. 1995; for more recent work see Fromang 2010; Lesur & Longaretti 2011; Nauman & Blackman 2014). In order to visualize the distribution of power as a function of scale, it is customary to plot either spherical shell-averaged spectra or 2D slices in the planes perpendicular to the directions given by ${{{\boldsymbol{\hat{k}}} }_{x}},{{{\boldsymbol{\hat{k}}} }_{y}}$, and ${{{\boldsymbol{\hat{k}}} }_{z}}$. The second approach reveals that the turbulence is highly anisotropic, rendering the first approach suspect. In either case, temporal averages over many orbital timescales are usually invoked in order to obtain smooth spectra. However, as we have shown in the previous section, the degree of anisotropy of the turbulence varies with time significantly.

In order to gain a deeper understanding of the temporal evolution of the anisotropy of the turbulence in Fourier space, we pursue a different approach, as follows. We demonstrate the procedure by computing the Fourier spectrum² associated with the magnetic energy density in the early stages after the saturation of the linear instability and in the turbulent state. It is useful to visualize the results by selecting 2D cuts in the 3D Fourier space. The anisotropic nature of the turbulence in the $({{k}_{x}},{{k}_{y}})$ plane presents a natural direction to define two perpendicular directions in this plane. We define ${{{\boldsymbol{\hat{k}}} }_{1}}$ as the unit wavevector corresponding to the direction of maximum anisotropy in the plane spanned by $({{k}_{x}},{{k}_{y}})$. The unit vectors ${{{\boldsymbol{\hat{k}}} }_{3}}={{{\boldsymbol{\hat{k}}} }_{z}}$ and ${{\mathop{{}}\limits^{{}}}_{{}}}{{{\boldsymbol{\hat{k}}} }_{3}}\times {{{\boldsymbol{\hat{k}}} }_{1}}$ complete the orthonormal triad. Note that the directions provided by ${{{\boldsymbol{\hat{k}}} }_{1}}$ and ${{{\boldsymbol{\hat{k}}} }_{2}}$ depend in general on time!

Figure 11 shows the 2D cuts along the three different planes perpendicular to ${{{\boldsymbol{\hat{k}}} }_{1}}$, ${{{\boldsymbol{\hat{k}}} }_{2}}$, and ${{{\boldsymbol{\hat{k}}} }_{z}}$ for two different times, corresponding to the early times after the saturation of the linear phase of the MRI ($t=10\;{{{\Omega }}^{-1}}$, upper row) and the fully developed turbulent state ($t=260\;{{{\Omega }}^{-1}}$, lower row). In order to quantify the degree of anisotropy, at a given time, we performed a 2D Gaussian fit, which is overplotted in each panel. The Gaussian fit is carried out by identifying the peak of the power spectral density and then fitting a 1D Gaussian along the direction of maximum anisotropy and another 1D Gaussian in the direction perpendicular to it. The width of these two Gaussians determines the semimajor axis of the 2D ellipses that are plotted over the data.

Zoom In Zoom Out Reset image size

At time $t=10\;{{{\Omega }}^{-1}}$, the parasite mode has reached the peak of its super-exponential growth, as attested by the presence of significant power in non-axisymmetric modes. The distribution of power in the $({{k}_{2}},{{k}_{z}})$ plane is mostly isotropic, the 2D Gaussian fit has an aspect ratio of 1.12. On the other hand, the Gaussian fits performed in the $({{k}_{1}},{{k}_{z}})$ and $({{k}_{x}},{{k}_{y}})$ planes show aspect ratios of ∼1.88, 3.12 respectively. In the fully developed turbulent state, at $t=260\;{{{\Omega }}^{-1}}$, the flow is highly anisotropic in the $({{k}_{1}},{{k}_{z}})$ plane. The aspect ratios of the Gaussian fits in the $({{k}_{1}},{{k}_{z}})$, $({{k}_{2}},{{k}_{z}})$, and $({{k}_{x}},{{k}_{y}})$ planes are, respectively, 2.63, 1.16, and 3.59. The anisotropy of the spectrum has therefore increased significantly from the late linear evolution of the MRI to the ensuing turbulent state.

The contrast between the energy spectral density obtained using spherical shell-averaged values versus 1D profiles in each of the three planes described above is depicted in the rightmost panel in Figure 11. At early times, i.e., soon after the parasitic modes have disrupted the highly organized flow driven by the exponential growth of the MRI, the 1D spectra along the directions ${{{\boldsymbol{\hat{k}}} }_{1}}$ and ${{{\boldsymbol{\hat{k}}} }_{2}}$ is not exceedingly different, although it is clear that there is more power along ${{{\boldsymbol{\hat{k}}} }_{1}}$. However, at later times, in the fully developed turbulent regime, the turbulent power along ${{{\boldsymbol{\hat{k}}} }_{1}}$ is almost two orders of magnitude higher than the power in the direction perpendicular to it throughout most of the scales. In this case, the shell average is unable to provide a good representation of the anisotropic power spectrum of the turbulence. In both cases, averaging over spherical shells masks the spread of almost two orders of magnitude in power along the directions given by ${{{\boldsymbol{\hat{k}}} }_{1}}$ and ${{{\boldsymbol{\hat{k}}} }_{2}}$.

The distribution of power characterizing the Maxwell stress in Fourier space can be analyzed in a similar way that we have analyzed the magnetic energy density. The results are shown in Figure 12. As observed for the magnetic energy density, the Maxwell stress is also highly anisotropic and the degree of anisotropy varies significantly with time. Note that the fast decrease in power along the k_x and k_y directions, i.e., the usual directions along which 1D PSD cuts are plotted, is even steeper for the Maxwell stress than for the magnetic energy. This implies that if we are interested in characterizing the degree of anisotropy of the stress is even more critical to analyze the spectra along ${{{\boldsymbol{\hat{k}}} }_{1}}$ and ${{{\boldsymbol{\hat{k}}} }_{2}}$, which provide a better characterization of the 3D distribution of power in Fourier space.

Zoom In Zoom Out Reset image size

It is evident that the distribution of power in Fourier space is highly anisotropic and time dependent. Therefore, in order to provide a robust characterization of the turbulent spectrum we proceed as follows. At a fixed time, we obtain the 3D PSD in Fourier space and define the direction of maximum anisotropy in horizontal planes. This is accomplished by fitting a 2D Gaussian ellipse to the distribution in the $({{k}_{x}},{{k}_{y}})$ plane and defining the direction of maximum anisotropy as ${{{\boldsymbol{k}} }_{1}}$ and the two orthogonal directions ${{{\boldsymbol{k}} }_{2}}$ and ${{{\boldsymbol{k}} }_{3}}$. The 1D power spectra along each of these directions are readily obtained. This procedure is repeated for different times, and the resulting 1D spectra can be properly averaged along each independent direction.

The average power spectra of the magnetic energy density obtained in this way are shown in Figure 13. Most of the power is concentrated in large scales, and a power-law behavior is observed beyond a characteristic wavenumber, which is different along each direction. In order to quantify the steepness of these profiles, we fitted power-law ${{k}^{-{{\alpha }_{i}}}}$ distributions to $k|{\rm DFT}(B){{|}^{2}}$ along each independent direction i = 1, 2, 3. We obtained ${{\alpha }_{1}}\sim 2.8,{{\alpha }_{2}}\sim 2.9,{{\alpha }_{z}}\sim 3.2$. For comparison, we also plotted the power spectra that result from a standard spherical shell average. This corresponds to ${{\alpha }_{r}}=3$. The orientation of the ellipsoid in Fourier space can be characterized by the angle subtended between ${{{\boldsymbol{k}} }_{1}}$ and ${{{\boldsymbol{k}} }_{x}}$. The average value of this angle is about 35°. The distribution of the values of this angle is shown in the inset.

Zoom In Zoom Out Reset image size

The fact that the distribution of power in MRI-driven turbulence is anisotropic in Fourier space was already indicated by Hawley et al. (1995). Their Figure 4 shows that most of the power in the magnetic field resides in modes along k_z, in qualitative agreement with our results. Hawley et al. (1995) also note a steepening at higher values of k due to numerical diffusion, which we do not find in our simulations, possibly owing to the incompressible pseudo-spectral method being less diffusive than the second-order scheme (see, e.g., Figure 1 of Brandenburg 2001). Our results can also be compared with those of Lesur & Longaretti (2011), who defined a spherical average and found a 3/2 slope for both the magnetic and kinetic energy spectra in the $1\lt k\lt 10$ range, with a subsequent steeper falloff in the range $10\lt k\lt 100$. In those simulations, there does not seem to be a clear sign of an inertial range. The simulations in Lesur & Longaretti (2011) used a higher value of $\beta =q{\Omega }L/{{B}^{2}}={{10}^{3}}$ than our β ∼ 10, a slightly different aspect ratio of 4 × 4 × 1, and a broadband noise excitation. The fact that we have chosen the value of β such that a single wavelength of the fastest-growing MRI mode fits in the domain may have helped to set an effectively larger injection scale, possibly resulting in the presence of a larger inertial range, as seen in Figure 13.

It is worth pointing out that most of the turbulent power, in both the magnetic energy density and the Maxwell stress, resides in modes with wavevectors along ${{{\boldsymbol{k}} }_{z}}$, i.e., the direction perpendicular to the shear, with the contribution from nonaxisymmetric modes being subdominant. This could be due to the fact that we have considered a rather strong vertical magnetic field (see, e.g., Pessah et al. 2007). Whether this ordering of power in Fourier space is a generic feature of MRI-driven turbulence in more general settings, e.g., in the presence of weaker fields and/or more general magnetic field geometries, or stratification will be addressed in subsequent work.