Active Modes and Dynamical Balances in MRI Turbulence of Keplerian Disks with a Net Vertical Magnetic Field

To analyze spectral dynamics, we decompose the perturbations $f\equiv ({\boldsymbol{u}},p,\theta ,{\boldsymbol{b}})$ into Fourier harmonics/modes,

$$\begin{eqnarray}&&f({\boldsymbol{r}},t)=\int \bar{f}({\boldsymbol{k}},t)\exp (i{\boldsymbol{k}}\cdot {\boldsymbol{r}}){d}^{3}{\boldsymbol{k}},\end{eqnarray} \tag{ 6 }$$

where $\bar{f}\equiv (\bar{{\boldsymbol{u}}},\bar{p},\bar{\theta },\bar{{\boldsymbol{b}}})$ are the corresponding Fourier transforms. The derivation of spectral equations of perturbations is a technical task and presented in the Appendix: substituting decomposition (6) into perturbation Equations (34)–(42), we obtain the equations for the spectral velocity (Equations (58)–(60)), logarithmic density (Equation (46)), and magnetic field (Equations (47)–(49)). Below, we give the final set of equations for the quadratic forms of these quantities in Fourier space.

Multiplying Equations (58)–(60), respectively, by ${\bar{u}}_{x}^{* }$, ${\bar{u}}_{y}^{* }$, and ${\bar{u}}_{z}^{* }$ and adding up their complex conjugates, we obtain

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{\partial }{\partial t}\displaystyle \frac{| {\bar{u}}_{x}{| }^{2}}{2} & = & -{{qk}}_{y}\displaystyle \frac{\partial }{\partial {k}_{x}}\displaystyle \frac{| {\bar{u}}_{x}{| }^{2}}{2}\\ & & +\,{{ \mathcal H }}_{x}+{{ \mathcal I }}_{x}^{(u\theta )}+{{ \mathcal I }}_{x}^{({ub})}+{{ \mathcal D }}_{x}^{(u)}+{{ \mathcal N }}_{x}^{(u)},\end{array}\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{\partial }{\partial t}\displaystyle \frac{| {\bar{u}}_{y}{| }^{2}}{2} & = & -{{qk}}_{y}\displaystyle \frac{\partial }{\partial {k}_{x}}\displaystyle \frac{| {\bar{u}}_{y}{| }^{2}}{2}\\ & & +\,{{ \mathcal H }}_{y}+{{ \mathcal I }}_{y}^{(u\theta )}+{{ \mathcal I }}_{y}^{({ub})}+{{ \mathcal D }}_{y}^{(u)}+{{ \mathcal N }}_{y}^{(u)},\end{array}\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{\partial }{\partial t}\displaystyle \frac{| {\bar{u}}_{z}{| }^{2}}{2} & = & -{{qk}}_{y}\displaystyle \frac{\partial }{\partial {k}_{x}}\displaystyle \frac{| {\bar{u}}_{z}{| }^{2}}{2}\\ & & +\,{{ \mathcal H }}_{z}+{{ \mathcal I }}_{z}^{(u\theta )}+{{ \mathcal I }}_{z}^{({ub})}+{{ \mathcal D }}_{z}^{(u)}+{{ \mathcal N }}_{z}^{(u)},\end{array}\end{eqnarray} \tag{ 9 }$$

where the terms of linear origin are

$$\begin{eqnarray}&&{{ \mathcal H }}_{x}=\left(1-\displaystyle \frac{{k}_{x}^{2}}{{k}^{2}}\right)({\bar{u}}_{x}{\bar{u}}_{y}^{* }+{\bar{u}}_{x}^{* }{\bar{u}}_{y})+2(1-q)\displaystyle \frac{{k}_{x}{k}_{y}}{{k}^{2}}| {\bar{u}}_{x}{| }^{2},\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}\begin{array}{rcl}{{ \mathcal H }}_{y} & = & \displaystyle \frac{1}{2}\left[q-2-2(q-1)\displaystyle \frac{{k}_{y}^{2}}{{k}^{2}}\right]({\bar{u}}_{x}{\bar{u}}_{y}^{* }+{\bar{u}}_{x}^{* }{\bar{u}}_{y})\\ & & -\,2\displaystyle \frac{{k}_{x}{k}_{y}}{{k}^{2}}| {\bar{u}}_{y}{| }^{2},\end{array}\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}&&{{ \mathcal H }}_{z}=(1-q)\displaystyle \frac{{k}_{y}{k}_{z}}{{k}^{2}}({\bar{u}}_{x}{\bar{u}}_{z}^{* }+{\bar{u}}_{x}^{* }{\bar{u}}_{z})-\displaystyle \frac{{k}_{x}{k}_{z}}{{k}^{2}}({\bar{u}}_{y}{\bar{u}}_{z}^{* }+{\bar{u}}_{y}^{* }{\bar{u}}_{z}),\end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray}&&{{ \mathcal I }}_{i}^{(u\theta )}=\left(\displaystyle \frac{{k}_{i}{k}_{z}}{{k}^{2}}-{\delta }_{{iz}}\right)\displaystyle \frac{\bar{\theta }{\bar{u}}_{i}^{* }+{\bar{\theta }}^{* }{\bar{u}}_{i}}{2},\end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray}&&{{ \mathcal I }}_{i}^{({ub})}=\displaystyle \frac{i}{2}{k}_{z}{B}_{0z}({\bar{u}}_{i}^{* }{\bar{b}}_{i}-{\bar{u}}_{i}{\bar{b}}_{i}^{* }),\end{eqnarray} \tag{ 14 }$$

and

$$\begin{eqnarray}&&{{ \mathcal D }}_{i}^{(u)}=-\displaystyle \frac{{k}^{2}}{\mathrm{Re}}| {\bar{u}}_{i}{| }^{2},\end{eqnarray} \tag{ 15 }$$

and the modified nonlinear transfer functions for the quadratic forms of the velocity components are

$$\begin{eqnarray}&&{{ \mathcal N }}_{i}^{(u)}=\displaystyle \frac{1}{2}({\bar{u}}_{i}{Q}_{i}^{* }+{\bar{u}}_{i}^{* }{Q}_{i}).\end{eqnarray} \tag{ 16 }$$

Here the index i = x, y, z henceforth, ${\delta }_{{iz}}$ is the Kronecker delta, and Q_i (given by Equation (61)) describes the nonlinear transfers via triad interactions for the spectral velocities ${\bar{u}}_{i}$ in Equations (58)–(60). It is readily shown that the sum of ${{ \mathcal H }}_{i}$ is equal to the spectrum of the Reynolds stress multiplied by the shear parameter q, ${ \mathcal H }={{ \mathcal H }}_{x}+{{ \mathcal H }}_{y}+{{ \mathcal H }}_{z}=q({\bar{u}}_{x}{\bar{u}}_{y}^{* }+{\bar{u}}_{x}^{* }{\bar{u}}_{y})/2$.

Next, multiplying Equation (46) by ${\bar{\theta }}^{* }$ and adding its complex conjugate, we get

$$\begin{eqnarray}&&\displaystyle \frac{\partial }{\partial t}\displaystyle \frac{| \bar{\theta }{| }^{2}}{2}=-{{qk}}_{y}\displaystyle \frac{\partial }{\partial {k}_{x}}\displaystyle \frac{| \bar{\theta }{| }^{2}}{2}+{{ \mathcal I }}^{(\theta u)}+{{ \mathcal D }}^{(\theta )}+{{ \mathcal N }}^{(\theta )},\end{eqnarray} \tag{ 17 }$$

where the terms of linear origin are

$$\begin{eqnarray}&&{{ \mathcal I }}^{(\theta u)}=\displaystyle \frac{{N}^{2}}{2}({\bar{u}}_{z}{\bar{\theta }}^{* }+{\bar{u}}_{z}^{* }\bar{\theta }),\end{eqnarray} \tag{ 18 }$$

$$\begin{eqnarray}&&{{ \mathcal D }}^{(\theta )}=-\displaystyle \frac{{k}^{2}}{\mathrm{Pe}}| \bar{\theta }{| }^{2},\end{eqnarray} \tag{ 19 }$$

and the modified nonlinear transfer function for the quadratic form of the entropy is

$$\begin{eqnarray}&&{{ \mathcal N }}^{(\theta )}=\displaystyle \frac{i}{2}{\bar{\theta }}^{* }({k}_{x}{N}_{x}^{(\theta )}+{k}_{y}{N}_{y}^{(\theta )}+{k}_{z}{N}_{z}^{(\theta )})+{\rm{c}}.{\rm{c}}.,\end{eqnarray} \tag{ 20 }$$

where ${N}_{i}^{(\theta )}$ (given by Equation (53)) describes the nonlinear transfers for the entropy $\bar{\theta }$ in Equation (46).

Finally, multiplying Equations (47)–(49), respectively, by ${\bar{b}}_{x}^{* }$, ${\bar{b}}_{y}^{* }$, and ${\bar{b}}_{z}^{* }$ and adding up their complex conjugates, we obtain

$$\begin{eqnarray}&&\displaystyle \frac{\partial }{\partial t}\displaystyle \frac{| {\bar{b}}_{x}{| }^{2}}{2}=-{{qk}}_{y}\displaystyle \frac{\partial }{\partial {k}_{x}}\displaystyle \frac{| {\bar{b}}_{x}{| }^{2}}{2}+{{ \mathcal I }}_{x}^{({bu})}+{{ \mathcal D }}_{x}^{(b)}+{{ \mathcal N }}_{x}^{(b)},\end{eqnarray} \tag{ 21 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial }{\partial t}\displaystyle \frac{| {\bar{b}}_{y}{| }^{2}}{2}=-{{qk}}_{y}\displaystyle \frac{\partial }{\partial {k}_{x}}\displaystyle \frac{| {\bar{b}}_{y}{| }^{2}}{2}+{ \mathcal M }+{{ \mathcal I }}_{y}^{({bu})}+{{ \mathcal D }}_{y}^{(b)}+{{ \mathcal N }}_{y}^{(b)},\end{eqnarray} \tag{ 22 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial }{\partial t}\displaystyle \frac{| {\bar{b}}_{z}{| }^{2}}{2}=-{{qk}}_{y}\displaystyle \frac{\partial }{\partial {k}_{x}}\displaystyle \frac{| {\bar{b}}_{z}{| }^{2}}{2}+{{ \mathcal I }}_{z}^{({bu})}+{{ \mathcal D }}_{z}^{(b)}+{{ \mathcal N }}_{z}^{(b)},\end{eqnarray} \tag{ 23 }$$

where ${ \mathcal M }$ is the spectrum of Maxwell stress multiplied by q,

$$\begin{eqnarray}&&{ \mathcal M }=-\displaystyle \frac{q}{2}({\bar{b}}_{x}{\bar{b}}_{y}^{* }+{\bar{b}}_{x}^{* }{\bar{b}}_{y}),\end{eqnarray} \tag{ 24 }$$

$$\begin{eqnarray}&&{{ \mathcal I }}_{i}^{({bu})}=-{{ \mathcal I }}_{i}^{({ub})}=\displaystyle \frac{i}{2}{k}_{z}{B}_{0z}({\bar{u}}_{i}{\bar{b}}_{i}^{* }-{\bar{u}}_{i}^{* }{\bar{b}}_{i}),\end{eqnarray} \tag{ 25 }$$

$$\begin{eqnarray}&&{{ \mathcal D }}_{i}^{(b)}=-\displaystyle \frac{{k}^{2}}{\mathrm{Rm}}| {\bar{b}}_{i}{| }^{2},\end{eqnarray} \tag{ 26 }$$

and the modified nonlinear terms for the quadratic forms of the magnetic field components are

$$\begin{eqnarray}&&{{ \mathcal N }}_{x}^{(b)}=\displaystyle \frac{i}{2}{\bar{b}}_{x}^{* }[{k}_{y}{\bar{F}}_{z}-{k}_{z}{\bar{F}}_{y}]+{\rm{c}}.{\rm{c}}.,\end{eqnarray} \tag{ 27 }$$

$$\begin{eqnarray}&&{{ \mathcal N }}_{y}^{(b)}=\displaystyle \frac{i}{2}{\bar{b}}_{y}^{* }[{k}_{z}{\bar{F}}_{x}-{k}_{x}{\bar{F}}_{z}]+{\rm{c}}.{\rm{c}}.,\end{eqnarray} \tag{ 28 }$$

$$\begin{eqnarray}&&{{ \mathcal N }}_{z}^{(b)}=\displaystyle \frac{i}{2}{\bar{b}}_{z}^{* }[{k}_{x}{\bar{F}}_{y}-{k}_{y}{\bar{F}}_{x}]+{\rm{c}}.{\rm{c}}.,\end{eqnarray} \tag{ 29 }$$

where ${\bar{F}}_{x},{\bar{F}}_{y},{\bar{F}}_{z}$ (given by Equations (54)–(56)) are the Fourier transforms of the respective components of the perturbed electromotive force and describe nonlinear transfers for the spectral magnetic field components in Equations (47)–(49).

Our analysis is based on Equations (7)–(9), (17), and (21)–(23), which describe the processes of linear $[{{ \mathcal H }}_{i}({\boldsymbol{k}},t)$, ${{ \mathcal I }}_{i}^{(u\theta )}({\boldsymbol{k}},t)$, ${{ \mathcal I }}^{(\theta u)}({\boldsymbol{k}},t)$, ${{ \mathcal I }}_{i}^{({ub})}({\boldsymbol{k}},t)$, ${{ \mathcal I }}_{i}^{({bu})}({\boldsymbol{k}},t)$, and ${ \mathcal M }({\boldsymbol{k}},t)]$ and nonlinear $[{{ \mathcal N }}_{i}^{(u)}({\boldsymbol{k}},t)$, ${{ \mathcal N }}^{(\theta )}({\boldsymbol{k}},t)$, and ${{ \mathcal N }}_{i}^{(b)}({\boldsymbol{k}},t)$] origin. The terms ${{ \mathcal D }}_{i}^{(u)}({\boldsymbol{k}},t)$, ${{ \mathcal D }}^{(\theta )}({\boldsymbol{k}},t)$, and ${{ \mathcal D }}_{i}^{(b)}({\boldsymbol{k}},t)$ describe, respectively, the effects of viscous, thermal, and resistive dissipation as a function of wavenumber and are negative definite. Except for the orientation of the background field, these basic dynamical equations in Fourier space, and hence the above terms, are formally analogous to the corresponding ones in the case of a net azimuthal magnetic field derived in Paper I (Equations (11)–(17) there), which also gives the description/meaning of each dynamical term in these spectral equations. The only, though principal, difference lies in the kinetic–magnetic cross terms ${{ \mathcal I }}_{i}^{({ub})}({\boldsymbol{k}},t)=-{{ \mathcal I }}_{i}^{({bu})}({\boldsymbol{k}},t)$; however, this difference drastically changes the dynamical processes, resulting in entirely different physics of the turbulence for azimuthal and vertical orientations of the background field. These cross terms describe, respectively, the influence of the i-component of the magnetic field (${\bar{b}}_{i}$) on the same component of the velocity (${\bar{u}}_{i}$) and vice versa for each mode. They are of linear origin and, in the present case with vertical magnetic field, ${{ \mathcal I }}_{i}^{({ub})}({\boldsymbol{k}},t)\propto {k}_{z}{B}_{0z}$, while in the case with azimuthal magnetic field, ${{ \mathcal I }}_{i}^{({ub})}({\boldsymbol{k}},t)\propto {k}_{y}{B}_{0y}$. Consequently, in the presence of net vertical field, this mutual influence is nonzero for axisymmetric (channel) modes with k_y = 0, which are not sheared by the flow, leading to their continual exponential growth, i.e., to the classical MRI in the linear regime. As discussed in the Introduction, in the given magnetized disk shear flow, in addition to this exponential growth, which dominates in fact at larger times, there exists shear-induced linear nonmodal growth of nonaxisymmetric (with ${k}_{y}\ne 0$) as well as axisymmetric modes themselves, which dominates instead at finite (dynamical) times (Mamatsashvili et al. 2013; Squire & Bhattacharjee 2014). Assessment of the relative roles of the axisymmetric and nonaxisymmetric modes in the vertical field MRI problem is one of the key questions of our study. We note that despite the transient nature of the nonmodal growth, it is not of secondary importance in forming spectral characteristics of the turbulence.

The course of events in the presence of a net azimuthal magnetic field, studied in Paper I, is essentially different: the cross terms for axisymmetric modes vanish, which leads to the absence of the classical MRI and therefore of the exponential growth of the channel mode—a key ingredient in the dynamics of a net vertical field MRI turbulence. In this case, the only source of energy for the turbulence remains linear transient growth of nonaxisymmetric modes (see also Balbus & Hawley 1992; Hawley et al. 1995; Papaloizou & Terquem 1997; Brandenburg & Dintrans 2006; Simon & Hawley 2009), and its sustenance, as demonstrated in Paper I, relies on the interplay between this linear nonmodal transient growth and nonlinear transverse cascade. As a result, as we will see below, the active modes and dynamical balances in the azimuthal and vertical field cases are also fundamentally different—it is the channel mode that mainly provides energy for the rest modes via the transverse cascade in the latter case.

Although our analysis is based on the above dynamical equations for quadratic forms of physical quantities in Fourier space, for completeness, we also present an equation for the normalized spectral kinetic energy of modes/harmonics, ${{ \mathcal E }}_{K}=(| {\bar{u}}_{x}{| }^{2}+| {\bar{u}}_{y}{| }^{2}+| {\bar{u}}_{z}{| }^{2})/2$:

$$\begin{eqnarray}&&\displaystyle \frac{\partial {{ \mathcal E }}_{K}}{\partial t}=-{{qk}}_{y}\displaystyle \frac{\partial {{ \mathcal E }}_{K}}{\partial {k}_{x}}+{{ \mathcal H }}_{E}+{{ \mathcal I }}_{E}^{(u\theta )}+{{ \mathcal I }}_{E}^{({ub})}+{{ \mathcal D }}_{E}^{(u)}+{{ \mathcal N }}_{E}^{(u)},\end{eqnarray} \tag{ 30 }$$

which is obtained by summing Equations (7)–(9). The terms on the right-hand side are the sum of the following terms:

$$\begin{eqnarray*}\begin{array}{rcl}{{ \mathcal H }}_{E} & = & \displaystyle \sum _{i}{{ \mathcal H }}_{i}=q({\bar{u}}_{x}{\bar{u}}_{y}^{* }+{\bar{u}}_{x}^{* }{\bar{u}}_{y})/2,\\ {{ \mathcal I }}_{E}^{(u\theta )} & = & \displaystyle \sum _{i}{{ \mathcal I }}_{i}^{(u\theta )},\,\,\,\,{{ \mathcal I }}_{E}^{({ub})}=\displaystyle \sum _{i}{{ \mathcal I }}_{i}^{({ub})},\\ {{ \mathcal D }}_{E}^{(u)} & = & \displaystyle \sum _{i}{{ \mathcal D }}_{i}^{(u)}=-\displaystyle \frac{2{k}^{2}}{\mathrm{Re}}{{ \mathcal E }}_{K},\,{{ \mathcal N }}_{E}^{(u)}=\displaystyle \sum _{i}{{ \mathcal N }}_{i}^{(u)}.\end{array}\end{eqnarray*}$$

Similarly, we can get the equation for the normalized spectral magnetic energy of modes, ${{ \mathcal E }}_{M}=(| {\bar{b}}_{x}{| }^{2}+| {\bar{b}}_{y}{| }^{2}+| {\bar{b}}_{z}{| }^{2})/2$, by summing Equations (21)–(23),

$$\begin{eqnarray}&&\displaystyle \frac{\partial {{ \mathcal E }}_{M}}{\partial t}=-{{qk}}_{y}\displaystyle \frac{\partial {{ \mathcal E }}_{M}}{\partial {k}_{x}}+{ \mathcal M }+{{ \mathcal I }}_{E}^{({bu})}+{{ \mathcal D }}_{E}^{(b)}+{{ \mathcal N }}_{E}^{(b)},\end{eqnarray} \tag{ 31 }$$

where

$$\begin{eqnarray*}\begin{array}{rcl}{ \mathcal M } & = & -q({\bar{b}}_{x}{\bar{b}}_{y}^{* }+{\bar{b}}_{x}^{* }{\bar{b}}_{y})/2,\,\,\,\,{{ \mathcal I }}_{E}^{({bu})}=\displaystyle \sum _{i}{{ \mathcal I }}_{i}^{({bu})}=-{{ \mathcal I }}_{E}^{({ub})},\\ {{ \mathcal D }}_{E}^{(b)} & = & \displaystyle \sum _{i}{{ \mathcal D }}_{i}^{(b)}=-\displaystyle \frac{2{k}^{2}}{\mathrm{Rm}}{{ \mathcal E }}_{M},\,\,\,{{ \mathcal N }}_{E}^{(b)}=\displaystyle \sum _{i}{{ \mathcal N }}_{i}^{(b)}.\end{array}\end{eqnarray*}$$

The equation of the normalized thermal energy of modes, ${{ \mathcal E }}_{\mathrm{th}}=| \theta {| }^{2}/2{N}^{2}$, is readily derived by dividing Equation (17) just by N². We do not give it here because the thermal energy is always small compared to the kinetic and magnetic energies (see below).

One can also easily write the equation for the total normalized spectral energy of modes, ${ \mathcal E }={{ \mathcal E }}_{K}+{{ \mathcal E }}_{\mathrm{th}}+{{ \mathcal E }}_{M}$,

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{\partial { \mathcal E }}{\partial t} & = & -{{qk}}_{y}\displaystyle \frac{\partial { \mathcal E }}{\partial {k}_{x}}+{{ \mathcal H }}_{E}+{ \mathcal M }+{{ \mathcal D }}_{E}^{(u)}+{N}^{2}{{ \mathcal D }}^{(\mathrm{th})}\\ & & +\,{{ \mathcal D }}_{E}^{(b)}+{{ \mathcal N }}_{E}^{(u)}+{N}^{2}{{ \mathcal N }}^{(\theta )}+{{ \mathcal N }}_{E}^{(b)}.\end{array}\end{eqnarray} \tag{ 32 }$$

In the simulation box, wavenumbers are discrete and determined by the box sizes L_i, k_i = 2πn_i/L_i, where n_i = 0, ±1, ±2 ... and the index i = x, y, z. For convenience, we normalize these wavenumbers by the grid cell sizes of Fourier space, ${\rm{\Delta }}{k}_{i}=2\pi /{L}_{i}$, that is, k_i/Δk_i → k_i, after which all the wavenumbers become integers.

We start the analysis of the 3D spectra of the energies by first examining their dependence on the vertical wavenumber. For this purpose, we integrate the spectra of the energies and stresses in the (k_x, k_y)-plane, ${\hat{{ \mathcal E }}}_{K,M,\mathrm{th}}({k}_{z})=\int {{ \mathcal E }}_{K,M,\mathrm{th}}{{dk}}_{x}{{dk}}_{y}$, $(\hat{{ \mathcal H }}({k}_{z}),\hat{{ \mathcal M }}({k}_{z}))=\int ({ \mathcal H },{ \mathcal M }){{dk}}_{x}{{dk}}_{y}$, then average over time from t = 100 to 650 and represent it as a function of k_z in Figure 4. The spectral magnetic energy is larger than the kinetic one except at k_z = 0, where they are comparable, and both dominate the spectral thermal energy; the spectral Maxwell stress, in turn, is higher than the Reynolds one. The maximum of the magnetic energy, as well as both stresses, comes at $| {k}_{z}| =1$, corresponding to the channel mode k_c = (0, 0, ±1). Thus, most of the energy supply by the stresses occurs in fact at $| {k}_{z}| =0,1,2$ rather than at $| {k}_{z}| =3$, which corresponds to the most unstable axisymmetric mode of MRI in the modal approach for our parameters (lower panel of Figure 2). This is in agreement with the above linear nonmodal analysis, which has demonstrated how the nonnormality eliminates the dominance of the $| {k}_{z}| =3$ mode in preference to lower k_z modes during finite times (Figures 1 and 2). This result again supports the statement made above that the nonmodal growth over short, dynamical/orbital timescales is more relevant in the turbulent state. A maximum of the kinetic energy spectrum is instead at k_z = 0 and is a bit higher than the magnetic energy at this k_z. It corresponds to the axisymmetric zonal flow mode with k_zf = (±1, 0, 0).

Zoom In Zoom Out Reset image size

The main energy-carrying dynamical processes concentrated mostly at small $| {k}_{z}| =0,1,2$ in the vital area are nearly unaffected by dissipation, as demonstrated in Figure 5 showing the 1D spectra of the Maxwell stress, $\hat{{ \mathcal M }}$, for the fiducial run and at twice higher Reynolds numbers. The difference between these two spectra is only appreciable outside the vital area, where the stress, as well as the energies, fall off by about three orders of magnitude; but these wavenumbers are not dynamically important anyway.

Zoom In Zoom Out Reset image size

Having examined the dependence of the spectra on the vertical wavenumber, we now present in Figure 6 slices of the 3D spectra of the kinetic, ${{ \mathcal E }}_{K}$, and magnetic, ${{ \mathcal E }}_{M}$, energies in the (k_x, k_y)-plane. Both spectra are fairly anisotropic due to the shear with a similar overall elliptical form. We checked that this spectral anisotropy in fact extends to larger wavenumbers, up to the dissipation scale. The maximum of the kinetic energy spectrum comes at the wavenumber k_zf = (±1, 0, 0) of the zonal flow mode, while the maximum of the magnetic energy spectrum comes at the channel mode wavenumber k_c = (0, 0, ±1). These spectra also indicate that there is a broad range of trailing nonaxisymmetric modes (red areas) with energies, on average, comparable to those of the channel mode and hence taking an active part in the turbulence dynamics (see also Longaretti & Lesur 2010). By contrast, in the case with azimuthal field in Paper I, although the maximum of the magnetic energy spectrum comes again at k_z = 1, in the (k_x, k_y)-plane, this maximum comes instead at nearby nonaxisymmetric modes with $| {k}_{y}| =1$, while the channel mode is not dynamically as significant. A more detailed analysis of the relative roles and dynamics of the channel, zonal flow, and other (nonaxisymmetric) modes is presented below. Although not shown here, the slices of the 3D spectrum of the Maxwell stress in the (k_x, k_y)-plane share a similar shape and type of anisotropy as the magnetic energy spectrum in Figure 6. Finally, we note that analogous anisotropic energy spectra were also observed in other local studies of MRI turbulence with a nonzero net vertical field (Hawley et al. 1995; Lesur & Longaretti 2011; Murphy & Pessah 2015).

Zoom In Zoom Out Reset image size

The dynamically important modes with the most effective amplification have been identified above based on the optimal nonmodal growth calculations (Figure 1), i.e., on the analysis of the linear dynamics. However, the overall dynamical picture of the turbulence is formed as a result of the interplay of linear and nonlinear processes. So, it is reasonable to introduce the notion of active modes in k-space—the energy-carrying modes playing a major role in the dynamics—separately for the kinetic and magnetic components. The active kinetic (magnetic) modes are labeled as those modes whose spectral kinetic (magnetic) energy is higher than 50% of the maximum spectral kinetic energy ${{ \mathcal E }}_{K,\max }$ (magnetic energy ${{ \mathcal E }}_{M,\max }$) at the same time. Both types of active modes in Fourier space at k_z = 0, 1 are shown in Figure 7 with colored dots. The color of each mode indicates the fraction of the total evolution time during which this mode retains the higher energy. These dynamically active modes are located anisotropically in the (k_x, k_y)-plane in the region of wavenumbers $| {k}_{x}| \leqslant 6,| {k}_{y}| \leqslant 2$, forming the so-called vital area of the turbulence in k-space. In this area, the most energetic dynamical processes that set the strength and statistical properties of the turbulence are concentrated. There are two—the channel k_c and the zonal flow k_zf—modes (bigger red dots, respectively, in panels (a) and (d) of Figure 7) that clearly stand out against other modes, as they retain, respectively, higher magnetic (channel mode) and kinetic (zonal flow mode) energies most of the time. As a result, these two modes play a key role in shaping the turbulence dynamics and deserve a more detailed analysis, which is given in the following subsections.

Zoom In Zoom Out Reset image size

Apart from the channel and zonal flow modes, there are a sufficient number of other active modes (blue dots in Figure 7) that, individually, are not so significant compared to these two modes, since they retain the higher kinetic or magnetic energies for much less time. However, the combined action of the multitude of these modes is competitive to the channel and zonal flow modes. We label them as the rest modes and describe them in Section 4.5 below. By contrast, in the case of a purely azimuthal field (Paper I), these rest nonaxisymmetric modes (mainly those with $| {k}_{y}| =1$ and k_z = 0, 1) prevail not only collectively but also individually—they carry energies higher than or comparable to those of the channel mode for longer times than those in the present case with net vertical field. This results in different properties of the turbulence defined instead by dominant individual nonaxisymmetric rest modes. The other modes (including those with $| {k}_{z}| \geqslant 2$ not shown in Figure 7) lie outside the vital area, playing only a minor role in the dynamics.

We have seen above that the channel mode—the harmonic with wavenumber k_c = (0, 0, ±1) that is uniform in the horizontal (x, y)-plane and has the largest vertical wavelength (correspondingly, the smallest wavenumber, k_z = ±1) in the domain—is a key participant in the turbulence dynamics. It carries higher energy most of the time in the turbulent state among other active modes in the vital area and corresponds to the maximum of the spectral magnetic energy and stresses in Fourier space at any moment (see Figure 10 below). So, we first analyze the dynamics of this mode and then move to the zonal flow and the rest modes.

Figure 8 shows the evolution of the main spectral amplitudes of the magnetic field, $| {\bar{b}}_{x}({{\boldsymbol{k}}}_{c})| ,| {\bar{b}}_{y}({{\boldsymbol{k}}}_{c})|$, and velocity, $| {\bar{u}}_{x}({{\boldsymbol{k}}}_{c})| ,\,| {\bar{u}}_{y}({{\boldsymbol{k}}}_{c})|$, of the channel mode. Also plotted are the time-histories of the corresponding linear—stresses, ${ \mathcal M }({{\boldsymbol{k}}}_{c})$, ${ \mathcal H }({{\boldsymbol{k}}}_{c})$, and exchange, ${{ \mathcal I }}^{({ub})}({{\boldsymbol{k}}}_{c})$, ${{ \mathcal I }}^{({bu})}({{\boldsymbol{k}}}_{c})$—and nonlinear, ${{ \mathcal N }}^{(u)}({{\boldsymbol{k}}}_{c})$, ${{ \mathcal N }}^{(b)}({{\boldsymbol{k}}}_{c})$, terms in the above spectral Equations (7)–(8) and (21)–(22), which govern the dynamics of this mode. The action of the stresses and exchange terms together describes the (nonmodal) effect of MRI on the channel mode and is mainly responsible for its amplification. The action of the nonlinear terms describes, in turn, interaction of the channel mode with other active modes with different wavenumbers and components (this process is covered in detail in Section 4.6; see also Figures 12 and 14 below). These linear and nonlinear terms jointly determine the temporal evolution of the channel mode's velocity and magnetic field with characteristic recurrent bursts seen in Figure 8. It is clear from this figure that the typical variation, or dynamical time, of the nonlinear terms is indeed of the order of the orbital time, indicating again that the nonmodal physics of MRI is more relevant in the state of developed turbulence. Analysis of these terms allows us to understand nuances of the channel mode dynamics and ultimately the behavior of the total stress and the energy (Figure 3), because, as we will see below (Section 4.5), the channel mode bursts manifest themselves in the time-development of the volume-averaged magnetic energy and the Maxwell stress. These terms operate differently for velocity and magnetic fields.

Zoom In Zoom Out Reset image size

Let us first consider the magnetic field components in Figure 8. Here $| {\bar{b}}_{x}|$ and $| {\bar{b}}_{y}|$ are amplified, respectively, by ${{ \mathcal I }}_{x}^{({bu})}$ and ${ \mathcal M }$, as they are always positive, and then drained, respectively, by the nonlinear terms ${{ \mathcal N }}_{x}^{(b)}$ and ${{ \mathcal N }}_{y}^{(b)}$, which are always negative, transferring energy to other wavenumbers and components. Consequently, the negative peaks of these nonlinear terms slightly lag the corresponding positive peaks of the linear terms. The azimuthal field, which is a dominant field component in the channel mode, is also drained to a lesser degree by the negative exchange term ${{ \mathcal I }}_{y}^{({bu})}$, transferring its energy to $| {\bar{u}}_{y}|$. As for the velocity components, $| {\bar{u}}_{x}|$ is amplified by positive ${{ \mathcal H }}_{x}$, then drained mostly by ${{ \mathcal N }}_{x}^{(u)}$, which is always negative and, to a lesser degree, by also negative ${{ \mathcal I }}_{x}^{({ub})}$; the peaks of the latter two sink terms lag the corresponding peaks of the linear source term ${{ \mathcal H }}_{x}$. Finally, $| {\bar{u}}_{y}|$ is mostly amplified by ${{ \mathcal I }}_{y}^{({ub})}$, taking energy from $| {\bar{b}}_{y}|$, and is drained by always negative ${{ \mathcal H }}_{y}$. Here ${{ \mathcal N }}_{y}^{(u)}$ acts differently from the above nonlinear terms—it alternates sign, providing for $| {\bar{u}}_{y}|$ either source, when positive, or sink, when negative. So, other modes via nonlinear interaction occasionally reinforce the growth of the channel mode's azimuthal velocity, which is primarily due to MRI, being represented by the positive exchange term ${{ \mathcal I }}_{y}^{({ub})}$.

In conclusion, the channel mode is supported by linear—nonmodal MRI growth—processes, while nonlinear processes mostly oppose this growth. The effect of the nonlinear terms is maximal during the peaks of $| {\bar{b}}_{x}({{\boldsymbol{k}}}_{c})|$ and $| {\bar{b}}_{y}({{\boldsymbol{k}}}_{c})|$; they halt the growth and lead to a fast drain of the channel mode energy, redistributing this energy to other modes. Thus, in a strict self-consistent approach of the net vertical field MRI turbulence, this redistribution of the channel mode energy to other modes is actually due to the nonlinear processes (in the form of the transverse cascade; Figures 12 and 14). However, in previous studies, the channel mode was considered, for simplicity, not as a variable/perturbation mode itself, but as a part of the basic/stationary dynamics (see, e.g., Goodman & Xu 1994; Latter et al. 2009; Pessah & Goodman 2009; Pessah 2010). In those analyses, the redistribution of the channel mode energy to other modes is classified as a linear process—called parasitic instability—and these modes (parasites; in our terms, "the rest modes"), which feed on the channel mode, are assumed to have smaller amplitudes than those of the latter. This simplified approach definitely gives a good feeling of the broadening of the perturbation spectrum at the expense of the channel mode. At the same time, this simplification might suffer somewhat from inaccuracy when used for the turbulent case, since the channel mode in fact is not stationary and consists of a chaotically repeated processes of bursts and subsequent drains—both the values of the peaks and the time intervals between them are chaotic.

The second mode that also plays an important role in the turbulence dynamics is the zonal flow mode—the harmonic with wavenumber k_zf = (±1, 0, 0) that is independent of the azimuthal and vertical coordinates (i.e., is axisymmetric and vertically uniform) and has the largest radial wavelength (correspondingly, the smallest wavenumber, k_x = ±1) in the domain. Zonal flow in MRI turbulence with zero and nonzero net vertical field was found in several studies (Johansen et al. 2009; Bai & Stone 2014; Simon & Armitage 2014); however, its dynamics was usually analyzed based on the simplified models. Here we trace the zonal flow evolution in Fourier space self-consistently (without invoking the simplified models) by examining its governing dynamical terms in the main spectral equations directly from the simulation data.

We have seen in Figure 7 that in the turbulent state, the zonal flow mode carries higher kinetic energy among other active modes and corresponds to the maximum of the spectral kinetic energy in Fourier space most of the evolution time (see Figure 10 below). However, this mode does not contribute to the spectral Reynolds and Maxwell stresses, because its radial velocity and magnetic field components are identically zero, ${\bar{u}}_{x}({{\boldsymbol{k}}}_{{zf}})={\bar{b}}_{x}({{\boldsymbol{k}}}_{{zf}})=0$, due to incompressibility (Equation (50)) and divergence-free (Equation (51)) conditions. Besides, we checked that the vertical velocity, ${\bar{u}}_{z}({{\boldsymbol{k}}}_{{zf}})$, and magnetic field, ${\bar{b}}_{z}({{\boldsymbol{k}}}_{{zf}})$, components are also much smaller than the respective azimuthal ones. So, in Figure 9, we present only the evolution of the dominant spectral amplitudes of the azimuthal magnetic field, $| {\bar{b}}_{y}({{\boldsymbol{k}}}_{{zf}})|$, and velocity, $| {\bar{u}}_{y}({{\boldsymbol{k}}}_{{zf}})|$, in this mode. Here $| {u}_{y}({{\boldsymbol{k}}}_{{zf}})|$ changes on a longer timescale corresponding to axisymmetric zonal flow in physical space, while $| {b}_{y}({{\boldsymbol{k}}}_{{zf}})|$ exhibits both short timescale and long timescale variations. It is readily seen that the linear—stresses and exchange—terms are identically zero for the zonal flow mode in Equations (8)and (22), ${{ \mathcal H }}_{y}({{\boldsymbol{k}}}_{{zf}})={ \mathcal M }({{\boldsymbol{k}}}_{{zf}})=0$, ${{ \mathcal I }}_{y}^{({ub})}({{\boldsymbol{k}}}_{{zf}})={{ \mathcal I }}_{y}^{(u\theta )}({{\boldsymbol{k}}}_{{zf}})=0$. Thus, the linear processes (i.e., MRI) do not affect this mode. Therefore, it can be supported only by the nonlinear terms ${{ \mathcal N }}_{y}^{(u)},{{ \mathcal N }}_{y}^{(b)}$. Since ${\bar{u}}_{y}$ is the dominant component in this mode, we look more into its nonlinear term in order to pin down a mechanism of zonal flow generation. This term consists of the magnetic, ${{ \mathcal N }}_{y}^{(u,\mathrm{mag})}$, and hydrodynamic, ${{ \mathcal N }}_{y}^{(u,\mathrm{kin})}$, contributions,

$$\begin{eqnarray}&&{{ \mathcal N }}_{y}^{(u)}={{ \mathcal N }}_{y}^{(u,\mathrm{mag})}+{{ \mathcal N }}_{y}^{(u,\mathrm{kin})},\end{eqnarray} \tag{ 33 }$$

which for k_zf = (±1, 0, 0) have the following forms:

$$\begin{eqnarray*}&&\begin{array}{l}{{ \mathcal N }}_{y}^{(u,\mathrm{mag})}({{\boldsymbol{k}}}_{{zf}})\\ \quad =\,\displaystyle \frac{i}{2}{\bar{u}}_{y}^{* }({{\boldsymbol{k}}}_{{zf}})\int {d}^{3}{\boldsymbol{k}}^{\prime} {\bar{b}}_{y}({\boldsymbol{k}}^{\prime} ){\bar{b}}_{x}({{\boldsymbol{k}}}_{{zf}}-{\boldsymbol{k}}^{\prime} )+{\rm{c}}.{\rm{c}}.,\end{array}\end{eqnarray*}$$

with the integrand composed of the turbulent magnetic stresses, and

$$\begin{eqnarray*}&&\begin{array}{l}{{ \mathcal N }}_{y}^{(u,\mathrm{kin})}({{\boldsymbol{k}}}_{{zf}})\\ =\,-\displaystyle \frac{i}{2}{\bar{u}}_{y}^{* }({{\boldsymbol{k}}}_{{zf}})\int {d}^{3}{\boldsymbol{k}}^{\prime} {\bar{u}}_{y}({\boldsymbol{k}}^{\prime} ){\bar{u}}_{x}({{\boldsymbol{k}}}_{{zf}}-{\boldsymbol{k}}^{\prime} )+{\rm{c}}.{\rm{c}}.,\end{array}\end{eqnarray*}$$

with the integrand composed of the turbulent hydrodynamic stresses. The time-histories of these hydrodynamic and magnetic nonlinear parts, as well as ${{ \mathcal N }}_{y}^{(b)}$ for the zonal flow mode, are also shown in Figure 9. It is seen that ${{ \mathcal N }}_{y}^{(b)}({{\boldsymbol{k}}}_{{zf}})$ changes sign, acting as a source for $| {\bar{b}}_{y}({{\boldsymbol{k}}}_{{zf}})|$, when positive, and as a sink, when negative. On the other hand, the magnetic part is positive most of the time, producing and amplifying $| {\bar{u}}_{y}({{\boldsymbol{k}}}_{{zf}})|$. Thus, the zonal flow mode is supported by ${{ \mathcal N }}_{y}^{(u,\mathrm{mag})}$, which describes the action of the total azimuthal magnetic tension exerted by all the rest modes on this mode. By contrast, the hydrodynamic part ${{ \mathcal N }}_{y}^{(u,\mathrm{kin})}({{\boldsymbol{k}}}_{{zf}})$, which describes the total azimuthal hydrodynamic tension exerted by all other modes, is always negative and drains $| {\bar{u}}_{y}({{\boldsymbol{k}}}_{{zf}})|$. As a result, the peaks of the hydrodynamic part ${{ \mathcal N }}_{y}^{(u,\mathrm{kin})}$ slightly lag behind the peaks of the magnetic part ${{ \mathcal N }}_{y}^{(u,\mathrm{mag})}$, which initiates the growth of the azimuthal velocity.

Zoom In Zoom Out Reset image size

So far we have described two—the channel and zonal flow—modes that are the main participants of the dynamical processes. In addition to these modes, as classified in Section 4.2, there are a sufficient number of active modes in the system (blue dots in Figure 7) that are not individually significant but collectively have sufficiently large kinetic and magnetic energies—more than the channel and/or zonal flow modes. Consequently, the combined influence of these modes on the overall dynamics of the turbulence becomes important. This set of active modes has been labeled as the rest modes, which, as mentioned above, are also called parasitic modes in the studies of net vertical field MRI turbulence. As we checked, the contribution of other larger wavenumber modes, whose kinetic (magnetic) energy always remains less than 50% of the maximum spectral kinetic (magnetic) energy, in the energy balances is small compared to that of these three main types of modes and therefore is neglected below.¹⁰

Figure 10 allows us to compare the total magnetic, ${{ \mathcal E }}_{M,r}$ (left panel), and kinetic, ${{ \mathcal E }}_{K,r}$ (right panel), energies of the rest modes with the corresponding energies of the channel and zonal flow modes at different stages of the evolution. The magnetic energy of the zonal flow mode, ${{ \mathcal E }}_{M,{zf}}={{ \mathcal E }}_{M}(-{{\boldsymbol{k}}}_{{zf}})+{{ \mathcal E }}_{M}({{\boldsymbol{k}}}_{{zf}})\approx 0$, is small compared to that of the channel and rest modes and is not shown in the left plot. In this case, the total volume-averaged magnetic energy is approximately the sum of the magnetic energies of the channel mode, ${{ \mathcal E }}_{M,c}={{ \mathcal E }}_{M}(-{{\boldsymbol{k}}}_{c})+{{ \mathcal E }}_{M}({{\boldsymbol{k}}}_{c})$, and the rest modes, $\langle {E}_{M}\rangle \approx {{ \mathcal E }}_{M,c}+{{ \mathcal E }}_{M,r}$ (neglecting the very small contribution from the zonal flow mode), whereas in the total kinetic energy, all three types of modes contribute: $\langle {E}_{K}\rangle \approx {{ \mathcal E }}_{K,c}+{{ \mathcal E }}_{K,{zf}}+{{ \mathcal E }}_{K,r}$ with ${{ \mathcal E }}_{K,c}={{ \mathcal E }}_{K}(-{{\boldsymbol{k}}}_{c})+{{ \mathcal E }}_{K}({{\boldsymbol{k}}}_{c})$ and ${{ \mathcal E }}_{K,{zf}}={{ \mathcal E }}_{K}(-{{\boldsymbol{k}}}_{{zf}})+{{ \mathcal E }}_{K}({{\boldsymbol{k}}}_{{zf}})$ being the kinetic energies of the channel and zonal flow modes, respectively. In Figure 10, ${{ \mathcal E }}_{M,\max }$ and ${{ \mathcal E }}_{K,\max }$ are, respectively, the maximum values of the spectral magnetic and kinetic energies.¹¹ The left panel shows that this maximum value of the magnetic energy falls on the channel mode, i.e., ${{ \mathcal E }}_{M,\max }={{ \mathcal E }}_{M,c}$, during an entire course of the evolution. In other words, the channel mode is always magnetically the strongest among all the individual active rest modes (see also Longaretti & Lesur 2010), although the total magnetic energy of the rest modes dominates the magnetic energy of the channel mode, ${{ \mathcal E }}_{M,r}\gt {{ \mathcal E }}_{M,c}$, especially in the quiescent intervals between the bursts, when the magnetic energy of the channel mode is relatively low. This trend is also seen in the transport due to these modes (see Figure 11 below) and points to the importance of the rest modes in the turbulence dynamics. Despite this dominance, however, it is seen in the left panel of Figure 10 that the peaks of the rest modes' magnetic energy always occur shortly after the corresponding peaks of the channel mode's magnetic energy, indicating that the former increases, or is driven by, the latter (see Section 4.6). Thus, the typical burst-like behavior of the nonzero net vertical field MRI turbulence is closely related to the manifestation of the channel mode dynamics.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The role of the rest modes is also seen in the evolution of the kinetic energies in the right panel of Figure 10. As in the case of the magnetic energy, the contribution of the total kinetic energy of the rest modes, ${{ \mathcal E }}_{K,r}$, in the volume-averaged total kinetic energy is always dominant. The kinetic energy of the zonal flow mode, ${{ \mathcal E }}_{K,{zf}}$, is comparable to the kinetic energy of the rest modes only occasionally, while the contribution of the kinetic energy of the channel mode, ${{ \mathcal E }}_{K,c}$, in the overall balance of the total kinetic energy is even less. However, also in this case, the maximum value of the spectral kinetic energy, ${{ \mathcal E }}_{K,\max }$, always coincides with the kinetic energy of either the channel mode or the zonal flow mode, whichever is larger at a given moment. Moreover, on closer inspection of the kinetic energy curves, it appears that the peaks of the total kinetic energy of the rest modes in fact tend to mainly follow the corresponding peaks of the channel mode. Note also that often, the zonal flow energy increases and has a peak when the channel mode energy decreases and has a minimum; that is, the zonal flow tends to grow when the channel mode is at its minimum, and vice versa. This anticorrelation between the zonal flow and magnetic activity, which was also reported in Johansen et al. (2009), is further explored in Fourier space in the next subsection.

Finally, we characterize the angular momentum transport due to the rest modes and compare it to that of the channel mode. The transport is measured by the normalized sum of the Maxwell and Reynolds stresses, i.e., by the usual α parameter, which, in our nondimensional units, is defined for each harmonic as

$$\begin{eqnarray*}&&\bar{\alpha }=\displaystyle \frac{1}{2}({\bar{u}}_{x}{\bar{u}}_{y}^{* }-{\bar{b}}_{x}{\bar{b}}_{y}^{* })+{\rm{c}}.{\rm{c}}.\end{eqnarray*}$$

Figure 11 shows the time-histories of the transport parameter for the rest modes, separately for nonaxisymmetric ones, ${\alpha }_{r,\mathrm{non} \mbox{-} \mathrm{axisym}}={\sum }_{{k}_{y}\ne 0}\bar{\alpha }$, and axisymmetric ones, ${\alpha }_{r,\mathrm{axisym}}={\sum }_{{k}_{y}=0}\bar{\alpha }$, as well as for the channel mode, ${\alpha }_{c}=\bar{\alpha }({{\boldsymbol{k}}}_{c})$, and the total transport due to all the active modes in the box, $\alpha ={\sum }_{{\boldsymbol{k}}}\bar{\alpha }\approx {\alpha }_{r,\mathrm{non} \mbox{-} \mathrm{axisym}}+{\alpha }_{r,\mathrm{axisym}}+{\alpha }_{c}$. (The zonal flow mode does not contribute to the transport, since ${\bar{u}}_{x}({{\boldsymbol{k}}}_{{zf}})={\bar{b}}_{x}({{\boldsymbol{k}}}_{{zf}})=0$.) Like that for the energies above, the transport due to the nonaxisymmetric rest modes always dominates the transport due to the channel and axisymmetric rest modes, ${\alpha }_{r,\mathrm{non} \mbox{-} \mathrm{axisym}}\gt {\alpha }_{c},{\alpha }_{r,\mathrm{axisym}}$. This is consistent with the results of Longaretti & Lesur (2010), who also found the prevalence of nonaxisymmetric modes over the channel mode in the transport, although they arrived at the same conclusion by calculating the 1D Fourier transform of both the Maxwell and Reynolds stress as a function of k_y at the minimum and maximum of the transport, whereas we calculate the cumulative transport of the nonaxisymmetric rest modes. Note, however, that the peaks in α_{r, non-axisym} follow shortly after the peaks of the channel mode, α_c, despite the dominance of the former over the latter. At these moments, the channel mode typically exhibits more dramatic changes in the transport than the nonaxisymmetric modes. This again indicates that the channel mode, nonlinearly interacting with the rest modes, acts as a trigger for the bursts in the total transport, which, in turn, are primarily due to the collective activity of nonaxisymmetric modes.

So far, we have separately described the dynamics of the main—the channel, zonal flow, and rest modes. Now we move to the analysis of their interlaced dynamics in Fourier space, which determines the properties and balances of MRI turbulence with net vertical magnetic field. First of all, it is a characteristic burst-like behavior of the volume-averaged energies and stresses (Figures 3), which is related to the manifestations of these modes' activity and nonlinear interaction. As was discussed in Section 4.3, the short bursts of the channel mode arise as a result of the competition between the linear nonmodal MRI amplification and nonlinear redistribution of the channel mode's energy in k-space. Each such burst event starts with the linear amplification, which initially overwhelms the nonlinear redistribution to other modes, and therefore the energy of the channel mode increases. At the same time, the effect of the nonlinear terms on the channel mode, described by ${{ \mathcal N }}^{(b)}({{\boldsymbol{k}}}_{c}),{{ \mathcal N }}^{(u)}({{\boldsymbol{k}}}_{c})$, increases as well, as it gradually loses its energy to the rest modes due to nonlinear transfers. At a certain amplitude of the channel mode, these nonlinear terms become comparable to the linear ones that are responsible for MRI, halting the further growth of the channel mode and starting its fast drain (Figure 8). This moment, corresponding to the peak of the channel mode energy and stress, is a critical stage in its evolution, and therefore it is important to examine in more detail how nonlinear processes redistribute/transfer the energy of the channel mode in Fourier space at this time. For this purpose, we present the map of the nonlinear transfer terms in the (k_x, k_y)-plane at k_z = 0, 1, 2 at two key stages of the burst dynamics: when the channel mode energy has a peak (Figure 12) and when its energy is at a minimum (Figure 13). The nonlinear redistributions quite differ from each other at these "top" and "bottom" points of the channel mode bursts. In both moments, they are notably anisotropic due to the shear, i.e., depend on a wavenumber polar angle, being preferably concentrated in the first and third quadrants, and thus realize the transverse cascade of the modes in Fourier space. Figure 12, corresponding to the "top" in the vicinity of the moment t = 294 of the highest magnetic energy peak of the channel mode, shows strong negative peaks of the nonlinear transfer terms ${{ \mathcal N }}_{x}^{(b)}$, ${{ \mathcal N }}_{y}^{(b)}$, and ${{ \mathcal N }}_{x}^{(u)}$ for the channel mode (i.e., at k_x = k_y = 0, k_z = 1), as the nonlinearity leads its fast drain at this moment. The red and yellow areas in this figure represent the subset of the rest modes, mostly with k_z = 0, 2, which receive energy via nonlinear transfers at this moment of the channel mode's peak activity. Note that at this time, the nonlinearity transfers the channel mode energy not to the zonal flow mode (since ${{ \mathcal N }}_{y}^{(u)}({{\boldsymbol{k}}}_{{zf}})\lt 0$ at this moment) but to the rest modes. In other words, the rest modes can be considered the main "culprits/parasites" in the decline of the channel mode in the developed turbulent state. The action of the above nonlinear terms in Fourier space during other peaks of the channel mode are qualitatively similar to that given in Figure 12 at around t = 294; however, which modes, of the total set of rest modes, concretely receive energy is, of course, different for each burst.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Thus, this unifying approach—the analysis of the dynamics and the nonlinear transfer functions in Fourier space—enables us to self-consistently describe the interaction between the channel and the rest (parasitic) modes when the standard assumptions in the usual treatment of these modes (Goodman & Xu 1994; Pessah & Goodman 2009) break down: the amplitudes of these two mode types become comparable, so that one can no longer separate the channel as a primary background and parasites as small perturbations on top of that. Overall, the rest (parasitic) modes play a major role in the channel mode decline. After each burst, the transverse cascade determines a specific spectrum of the parasitic modes in the turbulent state, which mostly gain energy as a result of draining the channel mode.

The nonlinearity at the "bottom" point acts quite differently when the channel mode is at its minimum. In this case, the main effect is the production of the zonal flow mode, so that in Figure 13, we show only the nonlinear transfer term for the azimuthal velocity, ${{ \mathcal N }}_{y}^{(u)}$, at the end ("bottom") of the channel burst at three different moments. It is seen that this term is positive at the wavenumber of the zonal flow mode, k_zf, at all three moments (red dots at k_x = ±1, k_y = 0 on the plots with k_z = 0), acting as the source for the zonal flow mode. As a result, the formation of the zonal flow mode starts after the end of the channel mode burst. For the channel mode wavenumber k_c, ${{ \mathcal N }}_{y}^{(u)}$ is negative at these moments, acting as a sink. As for the other nonlinear terms at the "bottom" points (not shown here), they look qualitatively similar to those in Figure 12, except that they no longer have a prominent negative peak at k_c, since the rest and zonal flow modes dominate instead in the dynamics at these times.

We have seen from Figures 12 and 13 that the shear-induced anisotropy of nonlinear transfers in Fourier space lies at the basis of the balances of dynamical processes in vertical field MRI turbulence. In particular, they determine the drain of the channel mode and the resulting anisotropic spectrum of the rest modes, which gain energy after each peak of the channel mode. However, the snapshots of the spectral dynamics of the modes shown in these figures, being at separate moments, are still irregular, somewhat obscuring this generic spectral anisotropy of the MRI turbulence due to the shear. They convey only a short-term behavior of the turbulence. The spectral anisotropy of MRI turbulence with a net vertical field was demonstrated by Hawley et al. (1995), Lesur & Longaretti (2011), and, recently, in more detail by Murphy & Pessah (2015), who, however, characterized only the anisotropic spectra of the magnetic energy and Maxwell stress in k-space, not the nonlinear processes. This inherent anisotropic character of the dynamical process of shear flows in Fourier space is best identified and described by considering a long-term evolution of the turbulence, i.e., by averaging the spectral quantities over time, as was done in the above studies. Such an averaging in the present case with a net vertical field indeed shows the dominance of the nonlinear transverse cascade, i.e., the angular (i.e., over wavevector angles) redistribution of the power in Fourier space in the overall dynamical balances (nonlinear interactions) between the channel and the rest modes, which is a main process underlying the turbulence. Previously, Lesur & Longaretti (2011) also investigated the nonlinear transfers in a net vertical field MRI turbulence; however, together with time-averaging, they also applied shell-averaging of the linear and nonlinear dynamical terms in k-space that naturally smears out the spectral anisotropy and hence the transverse cascade.

We demonstrate the anisotropic nature of the nonlinear transfers on longer timescales in the example of the radial component of the magnetic field, for which it is most apparent. Figure 14 shows the corresponding nonlinear term ${{ \mathcal N }}_{x}^{(b)}$ in the (k_x, k_y)-plane at different k_z = 0, 1, 2 averaged in time from t = 100 to 650. Like in a single instant shown in Figure 12, the action of this term is highly anisotropic, that is, depends strongly on the azimuthal angle in the (k_x, k_y)-plane due to the shear. However, this time-averaged distribution is much smoother than its instantaneous counterpart and clearly shows the transfer of modes over wavevector angles, i.e., the transverse cascade of $| {\bar{b}}_{x}{| }^{2}$ from the blue regions, where ${{ \mathcal N }}_{x}^{(b)}\lt 0$ and acts as a sink for it, to the red and yellow regions, where ${{ \mathcal N }}_{x}^{(b)}\gt 0$ and acts as a source/production. In the green regions (outside the vital area), these terms are small, although they preserve the same anisotropic shape. As in Figure 12, at k_z = 1, ${{ \mathcal N }}_{x}^{(b)}$ has a pronounced negative peak at k_x = k_y = 0 in Figure 14, implying that in this time-averaged sense, too, the dominant nonlinear process at k_z = 1 is draining of the channel mode energy. This energy is transferred transversely to the rest modes with other k_z (mainly k_z = 0, 1) and a range of k_x, k_y indicated by the red and yellow areas in Figure 14, where ${{ \mathcal N }}_{x}^{(b)}\gt 0$. Since the time-averaging procedure spans the whole duration of the turbulence in the simulations, these smooth regions in fact enclose all of the rest modes ever excited during the entire evolution.

Zoom In Zoom Out Reset image size

It is well known from simulations of zero net flux MRI turbulence that vertical stratification favors notable dynamo action in magnetized disks, whereby a large-scale azimuthal magnetic field is generated and exhibits a "butterfly" diagram—regular quasi-periodic spatiotemporal oscillations (reversals; e.g., Brandenburg et al. 1995; Johansen et al. 2009; Davis et al. 2010; Gressel 2010; Shi et al. 2010; Simon et al. 2011; Bodo et al. 2012, 2014; Gressel & Pessah 2015). By contrast, although large-scale azimuthal field generation can still take place in unstratified zero net flux MRI turbulence in vertically sufficiently extended boxes, the spatiotemporal variation of this field is still more incoherent/irregular in the form of wandering patches (Lesur & Ogilvie 2008; Shi et al. 2016). A detailed physical mechanism of how stratification brings about such an organization of the magnetic field is not fully understood and requires further study. This is an important topical question, as it closely bears on the convergence problem in zero net flux MRI turbulence simulations in stratified and unstratified shearing boxes, which is currently under debate (see, e.g., Bodo et al. 2014; Shi et al. 2016; Ryan et al. 2017 and references therein). Despite the fact that the present problem, containing a nonzero net vertical flux, differs from the case of zero net flux usually considered in the context of the MRI dynamo, here we also observe a notable spatiotemporal organization of the large-scale azimuthal field in the turbulent state (see below).

An ordered cyclic behavior of characteristic quantities in MRI turbulence with a net vertical magnetic flux in the shearing box model, but with more complex isothermal vertical stratification and outflow boundary conditions in the vertical direction, was also found. In particular, Suzuki & Inutsuka (2009) analyzed quasi-periodic variations of the horizontally averaged vertical mass flux (wind) in MRI-turbulent disks, and Bai & Stone (2013), Fromang et al. (2013), Simon et al. (2013), and Salvesen et al. (2016) showed a similar oscillatory behavior for the mean azimuthal magnetic field, which is somewhat analogous to that found in our setup.

The disk flow configuration considered here is based on the Boussinesq approximation, where vertical stratification enters the governing perturbation Equations (58)–(60) in a simpler manner than in the above studies: through buoyancy-induced terms, proportional to the vertically uniform Brunt–Väisälä frequency squared, N². Still, this simplification has an advantage in that one could simply put only N² to zero in these dynamical equations without altering the other parameters and equilibrium configuration of the system. This allows us to directly compare to the unstratified case and more conveniently isolate the basic role of stratification in the dynamo action. In the linear regime, stratification has only a little influence on the dynamics of vertical field MRI (Mamatsashvili et al. 2013). Moreover, it is easily seen that the effect of stratification identically vanishes for the most effectively amplified channel mode with k_x = k_y = 0, whose azimuthal field, equal to the averaged in the horizontal (x, y)-plane total azimuthal field, in fact constitutes the large-scale dynamo field. So, the effect of stratification on the dynamo action is solely attributable to nonlinearity, which modifies the electromotive force.

Figure 15 shows the spacetime diagrams of the horizontally averaged azimuthal magnetic field, $\langle {b}_{y}\rangle$, in the stratified fiducial case with N² = 0.25 as used in the paper and the unstratified, N² = 0, case (Table 1). It is clearly seen that stratification causes a remarkable organization of the azimuthal field into a coherent wavelike pattern, in contrast to the unstratified case. In both cases, however, this horizontally averaged field is dominated by the k_z = ±1 harmonics;¹² i.e., it is associated with the dominant channel mode analyzed above, but its phase variation with t and z in the stratified case is much more regular than that in the unstratified one. The amplitude of this averaged azimuthal field is about an order of magnitude larger than that of the background vertical field B_0z. Figure 16 shows that the rms of the azimuthal field and the volume-averaged Maxwell stress in the stratified case are larger, with stronger/higher and more frequent bursts, than those in the unstratified case. Clearly, the spacetime diagram shown in Figure 15 in the stratified case with periodic vertical boundary conditions differs from a typical "butterfly" observed in the case of isothermal stratification in the abovementioned studies. This difference can be attributable to the different nature of the buoyancy term and vertical boundary conditions used, which jointly cause elevation of the azimuthal flux from the midplane. Nevertheless, this comparison clearly demonstrates the capability of stratification to order/regularize the spatiotemporal variation of the mean azimuthal field—the main field component in the accretion disk dynamo. Here we have presented only a primary manifestation of the effect of stratification in the net vertical field MRI turbulence. In its own right, it is the subject of a detailed and refined investigation.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Apart from the shear-induced anisotropy of the nonlinear processes analyzed in our paper, a net nonzero background magnetic field itself, by definition, gives rise to anisotropy of the nonlinear dynamics to its parallel and perpendicular directions. However, these two types of anisotropy, having different origins, essentially differ from each other.

Anisotropy of nonlinear cascade processes due to magnetic field (directed along the z-axis) in plasmas with static equilibrium is analyzed in Goldreich & Sridhar (1995). The source of energy is assumed to be isotropic in both physical and Fourier spaces, with some outer scale of turbulence, L, corresponding to the smallest wavenumber, k₀ ∼ 1/L. Nonlinear cascade processes take place in the inertial range, perpendicular and parallel to the field wavenumbers k_⊥, k_z ≳ k₀, which is in fact the region of activity of these processes only. Overall, the regions of the energy supply and nonlinear cascade in Fourier space are separated from each other. As a result, the nonlinear processes do not affect the turbulence energy supply—they only transfer energy, mainly to small perpendicular (to the magnetic field) wavelengths, i.e., to large k_⊥ rather than along k_z, thereby forming an anisotropic spectrum of the turbulence in the $({k}_{\perp },{k}_{z})$-plane. This nonlinear cascade mostly to large k_⊥, besides increasing the magnitude $k=\sqrt{{k}_{\perp }^{2}+{k}_{z}^{2}}$ of the total wavevector, k = (k_⊥,k_z) (direct cascade), by definition also implies a change of its orientation that results in some angular redistribution of harmonics (i.e., transverse cascade). At k_⊥/k_z ∼ 1, the direct and transverse cascades are comparable. As the cascade proceeds, the ratio k_⊥/k_z increases, leading k to change primarily in magnitude rather than in orientation, and hence the direct cascade becomes dominant. In any case, the angular redistribution (transverse cascade) in this case is of secondary importance; it does not affect the turbulence energy supply.

A completely different situation arises in the presence of the disk flow shear. The linear (nonmodal) energy supply processes of perturbations/turbulence in shear flows are strongly anisotropic in Fourier space and generally occur over a broad range of wavenumbers without leaving a free room (i.e., inertial range) for the action of nonlinear processes only. Due to the inherent anisotropy of the linear dynamics due to the shear, the nonlinear processes in this range of wavenumbers become strongly anisotropic; i.e., the transverse cascade is the dominant nonlinear process. Overall, a dynamical picture of the turbulence is formed as a result of the interplay of the linear and nonlinear (transverse cascade) processes. We call the area of the most intensive interplay of the linear and nonlinear processes the vital area of the turbulence. Of course, the nonlinear processes lead to energy exchange between different modes, redistributing perturbation energy in Fourier space while leaving the total energy unchanged. However, the nonlinear transverse cascade in shear flows repopulates the (transiently) growing harmonics and, in this way, indirectly contributes to the energy supply to turbulence. As we have shown above, in the present problem of net vertical field MRI turbulence in Keplerian shear flow, it significantly affects the dynamical "design" of the turbulence and determines the interplay of its "building blocks": the channel, zonal flow, and rest modes. This is clearly illustrated by Figures 12–14, showing the resulting anisotropic dynamics (transfers) in the (k_x, k_y)-plane perpendicular to the field. By contrast, in the absence of basic flow velocity shear, nonlinear cascades in this plane are isotropic in Goldreich & Sridhar theory.