MAGNETOROTATIONAL INSTABILITY: NONMODAL GROWTH AND THE RELATIONSHIP OF GLOBAL MODES TO THE SHEARING BOX

To introduce ideas used in the remainder of this work, we give here a very simple example showing the physical origins of nonmodal growth of the simplest axisymmetric MRI. While the general idea (which is nothing but the standard Ω effect) has been discussed in previous works in a somewhat different context (Rincon et al. 2007, 2008), we feel that its presentation as a linear instability is a useful starting point for our examination of more complicated non-axisymmetric situations later in the text.

Consider a MHD system with a background linear shear flow and impose an initial perturbation to the magnetic field (the lack of a velocity perturbation renders the presence of a background magnetic field irrelevant). For perturbations that depend only on the vertical coordinate as $\boldsymbol {B}(z,t)=\boldsymbol {B}\exp (i k_z z)$, the induction equation,

$$\begin{equation} \frac{\partial \boldsymbol {B}}{\partial t}=\nabla \times \left(\boldsymbol {U}\times \boldsymbol {B}\right) + \bar{\eta }\nabla ^2 \boldsymbol {B}, \end{equation} \tag{ 1 }$$

with $\boldsymbol {U}=\left(0,-q x,0\right)$ becomes

$$\begin{equation} \frac{\partial }{\partial t}\left(\begin{array}{c}B_{x}\\ B_{y} \end{array}\right)=\left(\begin{array}{cc}-\bar{\eta }k_{z}^{2} & 0\\ -q & -\bar{\eta }k_{z}^{2} \end{array}\right)\left(\begin{array}{c}B_{x}\\ B_{y} \end{array}\right). \end{equation} \tag{ 2 }$$

This system is perhaps the simplest paradigm of nonmodal stability theory, appearing in many introductory treatments due to its tendency to exhibit strong transient growth at small $\bar{\eta } k_z^2$ (Trefethen & Embree 2005). More precisely, although the eigenvalues of the system ($-\bar{\eta } k_z^2$ repeated) indicate it is stable, in the limit $\bar{\eta } k_z^2\rightarrow 0$ the system can grow many orders of magnitude before eventually decaying exponentially, as illustrated in Figure 1. Indeed, with $\bar{\eta } k_z^2= 0$ the solution, B_x(t) = B_x(0), B_y(t) = B_y(0) − qtB_x(0), can grow indefinitely, the physical mechanism being simple advection of the initial perturbation by the shear (Ω effect). Of course, over long timescales this algebraic growth is dwarfed by the standard MRI (if there is a vertical background field), which can grow as |B(t)|² ∼ exp (qt/2) in these units. Nonetheless, it is interesting to note that with the initial conditions B_x(0) = −B_y(0) the magnetic energy growth $\partial _t \ln \left(B_x^2+B_y^2\right)$ for Equation (2) at t = 0 is q, the same as for the standard MRI. This result—over short timescales the MRI energy growth rate is q—holds for all axisymmetric and non-axisymmetric MRI modes given an appropriate choice of initial conditions (SB14).

Zoom In Zoom Out Reset image size

What is the significance of this transient growth if such perturbations eventually decay? In the case of perturbations on top of a quiescent (non-turbulent) disk, one could imagine that some nonlinear effect might become important as the azimuthal magnetic field increased in magnitude, particularly since such a field would have the radial and vertical structure of the initial perturbation if dissipation were sufficiently low. However, the question of what nonlinear mechanisms might be at work to allow for the sustenance of a turbulent state are subtle and very difficult to answer conclusively (Rincon et al. 2008; Riols et al. 2013) and we do not speculate on such ideas in this work. Another, less obvious, situation where the shorter time behavior of such perturbations may be of interest is if the plasma is already turbulent. Neglecting the issue of why the turbulence is sustaining in the first place, it is of note that this very simple nonmodal growth mechanism could be injecting energy into the turbulence at a similar rate to the most unstable MRI modes, since both have the same growth rate over short times. To explore these ideas more rigorously, we have recently been studying (Squire & Bhattacharjee 2014b) the interaction of the linear MRI system with space-time-dependent mean-fields in a self-consistent quasi-linear theory (Farrell & Ioannou 2003). The nonmodal growth of MRI is critical in this endeavor, and a more traditional quasi-linear approach (e.g., Sagdeev & Galeev 1969) based on eigenmode growth rates would be quite unsuitable for this problem. Finally, we note that simple applications of linear nonmodal ideas to wall bounded hydrodynamic flows (using parameterized versions of the turbulence induced viscosity) have shown that the fastest growing structures are remarkably similar to many of the large-scale structures observed in full nonlinear turbulence simulations (Butler & Farrell 1993; Cossu et al. 2008).

Rather than examining a particular case in detail, we have structured this paper to survey several different ways that nonmodal methods can be useful for the analysis of MRI. This choice was made because the techniques are useful in understanding both global (i.e., r dependent) and simplified local versions of MRI, as well as the connection between them. After describing our models and the fundamentals of nonmodal stability theory (Sections 2 and 3, respectively), we give a basic explanation of the relationship between MRI eigenmodes and nonmodal structures in Section 4. This is done for non-axisymmetric modes in both local and global models with hard-wall boundary conditions, to illustrate the origins and fundamental importance of structures that shear with the background flow (shear waves). Having seen the importance of such structures, we then illustrate the utility of the local model in Section 5 by directly comparing global nonmodal structures to the shearing wave equations (see Section 2.2). This relationship between the global and local pictures implies that the shearing wave equations should themselves be interpreted from the nonmodal standpoint and this is the purpose of Section 6. We illustrate how such an interpretation of the equations can be fruitful, perhaps allowing simple quasi-linear interpretations of MRI turbulence.

Our starting point is the incompressible, resistive MHD model,

$$\begin{eqnarray} \frac{\partial \boldsymbol {u}}{\partial t}+\left(\boldsymbol {u}\cdot \nabla \right)\boldsymbol {u} &=&-\nabla p+\nabla \times \boldsymbol {B}\times \boldsymbol {B} - \nabla \Phi +\bar{\nu }\nabla ^{2}\boldsymbol {u},\nonumber \\ \frac{\partial \boldsymbol {B}}{\partial t}+\left(\boldsymbol {u}\cdot \nabla \right)\boldsymbol {B} &=& \left(\boldsymbol {B}\cdot \nabla \right)\boldsymbol {u} + \bar{\eta }\nabla ^2\boldsymbol {B},\nonumber \\ [2ex] \nabla \cdot \boldsymbol {u}&=& 0,\;\;\;\nabla \cdot \boldsymbol {B}=0. \end{eqnarray} \tag{ 3 }$$

In Equation (3), as for the remainder of the article, we use dimensionless variables: $\boldsymbol {u}=\boldsymbol {u}_{si}/u_0,\: \boldsymbol {B}=\boldsymbol {B}_{si}/(u_0 \sqrt{\mu _0 \rho _0} ),\: p=p_{si}/(u_0^2 \rho _0 ),\, {\rm and}\, \Phi =\Phi _{si}/(u_0^2 \rho _0)$, where $\boldsymbol {u}_{si}$, $\boldsymbol {B}_{si}$, p_si, and Φ_si are, respectively, the fluid velocity, magnetic field, pressure, and gravitational potential in SI units, and u₀, ρ₀, and μ₀ are a characteristic velocity, the density (considered constant for simplicity), and the vacuum permeability. Lengths have been scaled by characteristic scale L₀ in Equation (3), and time is scaled by L₀/u₀. The fluid and magnetic diffusivities, $\bar{\nu }$ and $\bar{\eta }$, are defined as $\bar{\nu }=\nu /\left(u_0 L_0\right)$ and $\bar{\eta }=\eta /\left(u_0 L_0\right)$, where ν and η are the kinematic viscosity and resistivity of the plasma. Since most parameters in our problem are of the order of one, $\bar{\nu }$ and $\bar{\eta }$ are approximately the inverses of the fluid and magnetic Reynolds numbers, respectively.

For all global calculations we use a simplified version of the equilibrium in cylindrical coordinates proposed by Kersalé et al. (2004). This model includes a very small radial inflow velocity,

$$\begin{equation*} U_r=\alpha /r, \end{equation*}$$

driven by the viscosity acting on the azimuthal component of the velocity,

$$\begin{equation*} U_\theta =U_0 r^{1+\alpha /\bar{\nu }}. \end{equation*}$$

We take α to be $-3/2 \, \bar{\nu }$ to give a Keplerian rotation profile and set U₀ = 1 in keeping with our normalization. For simplicity, we use the magnetic field

$$\begin{equation*} \boldsymbol {B}_0=\left(0,r B_{0\theta },B_{0 z}\right), \end{equation*}$$

with B_0θ and B_0z constant. The pressure is determined through the equilibrium equation, and Φ = −1/r. Note that the equilibrium is determined by only four free parameters $B_{0\theta },\:B_{0 z},\:\bar{\nu },\: \hbox{and } \bar{\eta }$. For all calculations presented here, we use the domain (0.25, 2.25) in r. While not given here, we have also carried out calculations with more general profiles and results seem to be quite similar.

The global linear equations are obtained by linearizing Equation (3) about the background profile, i.e.,

$$\begin{eqnarray} \boldsymbol {u}&=&\boldsymbol {u}_0+\boldsymbol {u}^{\prime }, \;\; \boldsymbol {B}=\boldsymbol {B}_0+\boldsymbol {B}^{\prime },\nonumber \\ p&=&p_0+p^{\prime }, \end{eqnarray} \tag{ 4 }$$

and inserting the ansatz

$$\begin{equation*} f\left(r,\theta,z,t\right)=f\left(r,t\right)e^{i m \theta + i k_z z}, \end{equation*}$$

for each of the variables $\boldsymbol {u}^{\prime }$, $\boldsymbol {B}^{\prime }$, and p'. Finally, we rewrite the equations in terms of the Orr–Sommerfeld-like variables

$$\begin{eqnarray} u_r &=& u^{\prime }_r,\;\;B_r=B^{\prime }_r,\nonumber \\ [2ex] \zeta &=& i k_z u^{\prime }_{\theta }-i\frac{m}{r}u^{\prime }_z,\;\; \eta =i k_z B^{\prime }_{\theta }-i\frac{m}{r}B^{\prime }_z, \end{eqnarray} \tag{ 5 }$$

and rearrange to eliminate as many derivatives as possible. This choice of variables eliminates the pressure and reduces the eight equations to four, at the cost of causing fourth-order derivatives of u_r to appear in the equations. Because of the length of the equations resulting from this variable choice, we present them in Appendix A.

We use the incompressible shearing box (SB) equations for the local studies presented in this work. These equations are derived from the global equations with a shearing background velocity profile by transforming into the rotating frame and considering a small patch of fluid, see Umurhan & Regev (2004). In dimensionless variables, they are

$$\begin{eqnarray} \frac{\partial \boldsymbol {u}}{\partial t}+\left(\boldsymbol {u}\cdot \nabla \right)\boldsymbol {u}+2\Omega {\hat{\boldsymbol z}}\times \boldsymbol {u}&=&-\nabla p+\nabla \times \boldsymbol {B}\times \boldsymbol {B}\nonumber \\ &&\;\;\;+2q \Omega ^2 x {\hat{\boldsymbol x}} +\bar{\nu }\nabla ^{2}\boldsymbol {u},\nonumber \\ \frac{\partial \boldsymbol {B}}{\partial t}+\left(\boldsymbol {u}\cdot \nabla \right)\boldsymbol {B}&=&\left(\boldsymbol {B}\cdot \nabla \right)\boldsymbol {u} + \bar{\eta }\nabla ^2\boldsymbol {B},\nonumber \\ [2ex] \nabla \cdot \boldsymbol {u}&=& 0,\;\;\;\nabla \cdot \boldsymbol {B}=0. \end{eqnarray} \tag{ 6 }$$

Here $\boldsymbol {u}=\boldsymbol {u}_{si}/u_0,\: \boldsymbol {B}=\boldsymbol {B}_{si}/(u_0 \sqrt{\mu _0 \rho _0} ),\; {\rm and}\; p=p_{si}/(u_0^2 \rho _0 )$, where $\boldsymbol {u}_{si}$, $\boldsymbol {B}_{si}$, and p_si are, respectively, the fluid velocity, magnetic field, and pressure in SI units, and ρ₀ and μ₀ are the density (considered constant) and the vacuum permeability. Lengths have been scaled by characteristic scale L₀ in Equation (6), time is scaled by 1/Ω (with Ω the local rotation frequency), and the velocity scale u₀ is L₀Ω. As such, Ω = 1 in Equation (6). We have kept it explicitly to show more clearly how the basic MHD equations have been altered by the rotation, but will not include it explicitly for the remainder of this work. The directions x, y, and z correspond, respectively, to the radial, azimuthal, and vertical directions from the global model. The parameter q = −dln Ω/dln r embodies the radial velocity shear, and for all examples in this work we set q = 3/2 as for Keplerian rotation. The fluid and magnetic diffusivities, $\bar{\nu }$ and $\bar{\eta }$, are defined as $\bar{\nu }=\nu /(\Omega L_0^2)$ and $\bar{\eta }=\eta /(\Omega L_0^2)$, where ν and η are the kinematic viscosity and resistivity of the plasma. Since most parameters in our problem are of the order of one, $\bar{\nu }$ and $\bar{\eta }$ are approximately the inverses of the fluid and magnetic Reynolds numbers, respectively. The background velocity is azimuthal with linear shear in the radial direction, $\boldsymbol {u}_0=-q\Omega x$, and the background magnetic field is taken to be constant, $\boldsymbol {B}_0=(0,B_{0y},B_{0z})$. As for the global case, we linearize the equations about the background, $\boldsymbol {u}=\boldsymbol {u}_0+\boldsymbol {u}^{\prime },\:\boldsymbol {B}=\boldsymbol {B}_0+\boldsymbol {B}^{\prime },\, {\rm and}\,p=p_0+p^{\prime }$, and Fourier analyze in y and z by inserting

$$\begin{equation} f\left(x,y,z,t\right)=f\left(x,t\right)e^{i k_y y +i k_z z} \end{equation} \tag{ 7 }$$

for each dependent variable. Changing into the variables

$$\begin{eqnarray} u&=&u^{\prime }_x,\;\;B=B^{\prime }_x,\nonumber \\ \zeta &=&i k_z u^{\prime }_{y}-i k_y u^{\prime }_z,\;\; \eta =i k_z B^{\prime }_{y}-i k_y B^{\prime }_z, \end{eqnarray} \tag{ 8 }$$

and simplifying, we obtain the four linear partial differential equations in x and t,

$$\begin{eqnarray} && {\displaystyle\frac{\partial }{\partial t}\left(\begin{array}{@{}c@{}}\nabla ^2 u \\ \zeta \\ B \\ \eta \\ \end{array}\right) = } \nonumber\\ && {\left(\begin{array}{@{}cccc@{}}\bar{\nu } \nabla ^4+i q x k_y \nabla ^2 & 2 i k_z & i F \nabla ^2 & 0 \\ i (q-2) k_z & \bar{\nu } \nabla ^2+i q x k_y & 0 & i F \\ i F & 0 & \bar{\eta } \nabla ^2+i q x k_y & 0 \\ 0 & i F & i q k_z & \bar{\eta } \nabla ^2+i q x k_y \\ \end{array}\right) \!\cdot\! \left(\begin{array}{@{}c@{}}u \\ \zeta \\ B \\ \eta \\ \end{array}\right), } \nonumber\\ \end{eqnarray} \tag{ 9 }$$

where F ≡ k_yB_0y + k_zB_0z and $\nabla ^2\equiv -k_y^2-k_z^2+\partial ^2/\partial x^2$. Since these equations have no time dependence they can be Fourier analyzed in time using ∂/∂t → −iω to obtain a set of linear eigenvalue ordinary differential equations (ODEs); however, since much of this work focuses on nonmodal stability methods, we prefer to keep the time-dependence general even though they have been solved computationally from the eigenvalue standpoint (see Section 3).

Shearing wave equations. A common way to study the local non-axisymmetric linear MRI has been using a decomposition in terms of shearing waves. Shearing waves are simply sinusoidal waves that are static in the frame of the background flow; they have also been termed spatial Fourier harmonics or Kelvin waves by various authors. As part of this work, we compare the solutions obtained from assuming such a decomposition with global nonmodal stability calculations, showing excellent agreement.

The shearing wave equations are straightforwardly derived by inserting the ansatz

$$\begin{equation} f(x,t)=f(t)e^{i q k_y (t-t_{SW}) x}, \end{equation} \tag{ 10 }$$

for each dependent variable in Equation (9), where the initial orientation of the wave is determined by the parameter t_SW through k_x(0) = −qk_y(0 − t_SW). This yields the set of ODEs in time

$$\begin{equation} {\displaystyle\frac{\partial }{\partial t}U\!\left(t\right)= \left(\begin{array}{@{}cccc@{}}-\bar{\nu } k^2-2 q k_x k_y/k^2 & -2 i k_z/k^2 & i F & 0 \\ i (q-2)k_z & - \bar{\nu } k^2 & 0 & i F \\ i F & 0 & -\bar{\eta } k^2 & 0 \\ 0 & i F & -i q k_z & -\bar{\eta }k^2 \\ \end{array} \right) \cdot U\!\left(t\right), } \end{equation} \tag{ 11 }$$

with U(t) = (u, ζ, B, η) and $k^2=\sqrt{{{k_x^2+k_y^2+k_z^2}}}$. Due to the time dependence of k_x and k, Equation (11) cannot be usefully Fourier analyzed in time and must be solved numerically in general, although various analytic results have been obtained in previous works (Balbus & Hawley 1992; Terquem & Papaloizou 1996; Johnson 2007; Mamatsashvili et al. 2013; Squire & Bhattacharjee 2014a). It so happens that Equation (11) is actually nonlinearly valid (Goodman & Xu 1994; Balbus & Hawley 2006) due to rather fortuitous cancellations of nonlinear terms. As such, they can be derived by simply inserting the shearing wave ansatz directly into the nonlinear equations (Equation (6)) and changing variables (Equation (8)), skipping the linearization step entirely.

We start by considering the simplest possible background field configuration in the local box, a constant magnetic field in the z direction. However, in contrast to standard local stability approaches, we solve the full local differential equations (Equation (9)) with hard wall (perfectly conducting) boundary conditions. The reason for this choice is to illustrate the general irrelevance of the eigenmode at intermediate times; shearing wave structures are strongly apparent in the pseudo-mode, despite their incompatibility with the boundary conditions.

The transient and eigenmode growths for a weakly non-axisymmetric (k_y = 1) mode⁶ are illustrated in Figure 2 at B_0z = 1/10. To demonstrate the importance of the shearing wave, we also illustrate the time evolution of the pseudo-mode spatial structure (for t₀ = 10) in Figure 3. There are several important insights that can be gained from Figures 2 and 3.

• 1.  
  The maximum linear growth rate achievable (G(t) and G_Γ(t, 10) in Figure 2) is approximately twice as large as that of the eigenmode. In addition, this fast growth rate continues until the perturbation has been amplified by a factor of nearly 10⁵. In a system with relatively strong initial perturbations, this amplification is presumably sufficient for various nonlinear effects or secondary instabilities to become significant, while in the case of small initial perturbations the most unstable channel mode would be many orders of magnitude larger than non-axisymmetric modes by the time such nonlinear effects took over (e.g., by t = 30 the energy of most unstable channel mode would be 10¹³ times larger than that of the eigenmode in Figure 2). Thus, we surmise that the non-axisymmetric eigenvalue growth rate is largely irrelevant at these parameters. This important conclusion carries over to the global case (Figure 4).
• 2.  
  The pseudo-mode is a shearing wave, despite the presence of the hard wall boundary conditions. Considering that the most unstable eigenmodes are localized near the boundaries of the domain, it is perhaps initially surprising that the pseudo-mode is localized in the middle of the domain, at least until the transient growth subsides (around t = 20). Note that a very similar effect is seen in the hydrodynamic case, see, e.g., Mukhopadhyay et al. (2005). The general dominance of the shearing wave is nicely justified by our recent proof that shearing wave growth rates are always larger over short timescales than those of static structures (SB14).
• 3.  
  Unlike the (spectrally stable) hydrodynamic case, the time at which k_x ≈ 0 (i.e., the shearing wave is horizontal) does not correspond to any obvious change in the growth (in the hydrodynamic case k_x ≈ 0 when G(t) is maximum; Mukhopadhyay et al. 2005). In addition, we see little change in the initial shearing wave orientation (k_x(0)) with changes in dissipation, $\bar{\nu }$ and $\bar{\eta }$, in stark contrast to the hydrodynamic case. In fact, in all pseudo-mode calculations we have done for non-axisymmetric modes, the initial conditions satisfy k_x(0) ≈ k_y. This is partially explained by the calculations in SB14, where it is seen that the strongest growth over very short timescales (t = 0⁺) is for a shear wave with k_x(0) = ±k_y.
• 4.  
  At intermediate timescales the boundaries of the domain seem largely irrelevant. Indeed, it is a general feature of nonmodal stability that the transient growth is much less sensitive to modifications of the system than the eigenmode growth (Trefethen et al. 1993; Trefethen & Embree 2005). In this case, the modification is the change of boundary conditions from those that naturally accept shearing waves (e.g., SB boundary conditions) to those that do not (hard wall conditions).
• 5.  
  At late times the pseudo-mode starts to more closely resemble the most unstable eigenmodes (as might be expected) becoming more localized near the wall. As an interesting corollary of this, we note that the eventual decay of shearing waves due to the increasing k_x (Balbus & Hawley 1992; Brandenburg & Dintrans 2006) is not necessarily physically important, even discounting nonlinear effects. The reason is that the eigenmode growth can "take over" at large times, with the shearing structure transitioning into a non-shearing structure. Although we have not done the calculation, we conjecture that this could also be true when SB boundary conditions are utilized, with the very late time structure starting to resemble some type of time-periodic Floquet eigenmode (the SB system is periodic in time). Of course, such an (time-periodic) eigenmode could be stable and decay, though perhaps more slowly than a shearing wave.
• 6.  
  As the dissipation parameters ($\bar{\nu }$ and $\bar{\eta }$) are decreased the period over which the pseudo-mode resembles the shearing wave increases in time, thus leading to a larger total amplification of the disturbance. This is in spite of the fact that the non-axisymmetric spectral instability can disappear as the dissipation is decreased (Kitchatinov & Rudiger 2010). This is essentially implying that nonmodal effects become more important as $\bar{\nu },\bar{\eta }\rightarrow 0$. Note that the shearing wave growth will not continue forever even if $\bar{\nu }=\bar{\eta }=0$, as can be seen by solving the shearing wave equations (Equation (11)) in the dissipation-less limit (e.g., Brandenburg & Dintrans 2006).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The structure of the pseudo-mode in time does depend on the choice of when to maximize the growth, t₀. For instance, for large t₀ the structure is more localized near the boundaries at all times, but is still strongly shearing with the background flow. Note that the other variables (u, ζ, η) have very similar time-evolution (not shown). The illustrated transient growth is not simple stretching of the initial perturbation by the background flow (as in Section 1.1), which can be confirmed by noting that all three components of $\boldsymbol {B^{\prime }}$ grow at similar rates (not shown). We have also run calculations with different boundary conditions in x, including standard periodic conditions and the local equivalent of those advocated in Kersalé et al. (2004). We see that the structures observed in the pseudo-modes are always shearing waves in support of Item 4 above, so long as there are no strongly unphysical energy sources or sinks in the chosen boundary conditions.

Finally, we note that transient growth is not limited to non-axisymmetric modes, but can also be significant for the axisymmetric channel mode (k_y = k_x = 0) in the chosen energy norm. To be precise, some transient growth is possible even with periodic boundary conditions, whenever the vertical wavenumber is different from the wavenumber that gives maximum eigenmode growth, $k_z=1/B_{0z} \sqrt{15/16}$. In the local case, there is no substantial difference in spatial structure between the eigenmodes and pseudo-modes with hard-wall boundary conditions, but the ratio of components (u, ζ, B, η) is different. Note that one can straightforwardly choose a simple energy-like norm that removes the transient growth of axisymmetric modes, at least in the two-dimensional (2D) hydrodynamic case⁷ (Zhuravlev & Razdoburdin 2014). However, as illustrated by the introductory example (Section 1.1), transient growth of the axisymmetric instability is a very real physical effect. We give an example of global axisymmetric pseudo-mode growth in Section 5.

To illustrate that the prevalence of shearing wave structures is by no means unique to the local model, in Figures 4 and 5 we display the pseudo-mode growth and structure for a weakly non-axisymmetric mode in a weak purely vertical field. Note that the chosen k_z is around the lower limit of what might be physically relevant in a thin accretion disk (Kersalé et al. 2004). We see that all of the same conclusions that held in the local computation carry over to the global case. In fact, generally we have observed a greater prominence of transient effects in the global equations than the local equations, probably due to a greater propensity for pseudo-modes and eigenmodes to be localized in very different regions. This is certainly the case here, as evidenced by comparison of Figures 5 and 6 (the most unstable eigenmode); while the eigenmode is strongly localized near the outer boundary, the pseudo-mode is far removed from this. Of course, at very large times (not shown), the pseudo-mode moves out in radius and starts to more closely resemble the eigenmode. Terquem & Papaloizou (1996) noted a similar difference between the localization of eigenmodes and that of structures emerging from random noise (in a toroidal field with no nonlinear effects). They explain these findings in terms of the local growth rates, a connection that we make in Section 5. Finally, we note the extreme difference in growth rate between the nonmodal structures and eigenmodes (Figure 4). The pseudo-mode grows approximately six times faster than the least stable eigenmodes and reaches an amplification of 10⁵ before this fast growth shows any sign of slowing.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

As a final example to illustrate the greater relevance of pseudo-modes over eigenmodes, we initialize with random realizations of Gaussian noise and examine growth rates and prominent structures. This calculation can mitigate fears that the pseudo-mode structures might be less likely to be excited for some reason, and that total growth may not always be a good indicator of dynamical importance in a physical situation. We present an example of this calculation in Figure 7, for local parameters very similar to those of Figure 2. After an initial dip due to damped modes in the initial conditions,⁸ we see the growth curve follow that of G(t) very closely. In fact, for these parameters we see that even the minimum growth seen out of 100 realizations has overtaken that of the most unstable eigenmode by late times, i.e., the most unstable eigenmode is statistically a particularly bad choice of initial condition for the total amplification of the disturbance. Observing the structure of random realizations (not shown) we see a strong dominance of shearing waves at later times.

Zoom In Zoom Out Reset image size

In Figure 8, we illustrate the comparison of shearing waves with global pseudo-mode energy growth using the procedure outlined above. The parameters chosen are those for a weakly non-axisymmetric mode in a vertical field (the same as Figure 5), with two values of r₀ chosen for comparison. We see excellent agreement, although unsurprisingly the growth is most similar where it is strongest, around the maximum of the pseudo-mode (r₀ = 1). Moving very far from the maximum (e.g., r₀ = 2, not shown), we see rather poor agreement, presumably due to noise and errors in the numerical result. We have run many other similar computations and see excellent agreement across a wide range of parameters.

Zoom In Zoom Out Reset image size

As an interesting corollary of such results, one can approximately predict the structure of the global pseudo-mode using the shearing wave equations. The basic idea is to solve the shearing wave equations at each point in the global domain, maximizing the growth at a chosen t₀ using the nonmodal technique. Examining the amplification as a function of radius gives an approximation of the structure of the global pseudo-mode. While an exact comparison is tricky due to the choice of t₀ in the shearing wave equations, we have considered a range of parameters (not shown here) and the agreement generally appears rather good. In particular, the prediction of the spatial location of the pseudo-mode maximum is quite accurate. Such computations present further evidence that the local shear wave approximation is accurate in many cases (Papaloizou & Terquem 1997), and will be more meaningful than global eigenmodes over moderate timescales.

Figure 9 presents a similar comparison for the seldom studied case of an axisymmetric mode in a purely azimuthal field. While such a case could be argued to be somewhat pathological due to the importance of even a wisp of vertical field (at least without dissipation, we discuss this point more in Section 6.1, see also Balbus & Hawley 1998), it provides an interesting example. Despite the eigenmode being stable, there is rather strong growth, with the pseudo-mode amplified by ∼10³ by t = 7. The agreement with the shearing wave—in this case a simple channel mode with k_y = k_x = 0—is remarkably good. Of course, as discussed in Section 1.1, we are seeing simple advection of an initial field by the flow; nonetheless, it is comforting to see that the global pseudo-mode is locally behaving in the same way with very similar optimal initial conditions. For comparison, the inset to Figure 9 illustrates the radial structure of the pseudo-mode and least-stable eigenmode for the radial magnetic field. In this case, with a purely azimuthal field, the two are rather different; however, in the case of axisymmetric modes in a pure vertical field (not shown), the pseudo-mode generally closely resembles the eigenmode and the nonmodal growth is less significant due to the strong exponential growth of the standard MRI.

Zoom In Zoom Out Reset image size

As presented in Section 2.2, the shearing wave equations are derived through first applying a local expansion about the global equilibrium (Appendix B; also Umurhan & Regev 2004), then inserting the shearing wave ansatz. However, a more general way to obtain such equations is by directly inserting a shearing wave ansatz into the global equations, and only then applying the local expansion. For axisymmetric modes, the first step (insertion of a shear wave ansatz) is essentially a standard WKB expansion and has been used in many previous works. For example, in Blokland et al. (2005), the full WKB expansion (without a local approximation) is compared directly to r-dependent eigenmode solutions, showing excellent agreement. Noting that the standard WKB expansion has severe problems for non-axisymmetric modes (Knobloch 1992), we have shown how to extend such expansions to non-axisymmetry in SB14, by simply inserting the time-dependent ansatz

$$\begin{equation} f\!\left(r,t\right)\sim f\!\left(t\right)\,\exp \left(-i \frac{m}{r} U_0 r^{-q+1} (t-t_0)\right) \end{equation} \tag{ 20 }$$

for each variable and making the ordering assumptions (k_rr, k_zr, m) ∼ 1/, $\left(\bar{\nu },\:\bar{\eta }\right)\sim \epsilon ^2$ (see also Shtemler et al. 2012 for a different approach). Such an approximation is the natural extension of WKB to the case with global shear. Applying the local approximation to the resulting system of ODEs leads to the standard shearing wave equations (Equation (11)). Thought of in this way, we can consider the excellent agreement between global pseudo-modes and shearing waves to be a verification of the applicability of such WKB-like methods. To continue the analogy, in the same way that one might compare full eigenmodes to a WKB dispersion relations (i.e., WKB eigenmodes), the correct way to analyze such shearing wave equations is using nonmodal stability techniques (as carried out above).

This more global way of considering the problem may have several advantages. First, it is straightforward to extend the shearing wave equations to much more complicated domains and physical models. For example, strong magnetic fields, compressibility, stratification, or more complex diffusion operators (e.g., Pessah & Psaltis 2005; Heinemann & Papaloizou 2009; Salhi et al. 2012; Rosin & Mestel 2012) could easily be accounted for in the shearing wave equations.⁹ Second, the approach elucidates the connection between previous results that illustrate the quality of WKB methods for axisymmetric modes, and our results, which show the accuracy of the shearing wave equations over moderate timescales. Of course, we have primarily explored global models in which the local approximation (Appendix B) is accurate. In future work it would be interesting to study global models that include more complex physical effects, comparing global pseudo-modes to nonmodal solutions of the extended shearing wave equations derived directly from the chosen global equations.

Finally we clarify here that shearing wave equations can only ever give a good approximation to the global pseudo-mode behavior over moderate timescales. The reason is that eventually the eigenmode will take over, since the structures in a shearing wave necessarily move to smaller scales in time. This causes dissipative effects to dominate and the shearing wave to damp, even when the global system has one or more unstable eigenmodes. This effect is clearly seen in the last pane of Figure 3, where the pseudo-mode eventually starts to resemble the least stable eigenmodes. From a practical standpoint, we have noted that the shearing wave equations accurately represent the global pseudo-mode up until their solution starts to decrease in time.

Our study here focuses on how the local MRI changes with imposed vertical field while in a strong background azimuthal field. There are two primary motivations behind this choice of problem.

• 1.  
  Using analyses based on eigenmodes (or similar ideas for time-dependent shear waves, e.g., Balbus & Hawley 1992; Johnson 2007) MRI behaves a little unusually in an azimuthal field in the limit B_0z → 0 (Balbus & Hawley 1998). In particular, the growth rate is very sensitive to even a minute vertical field and enormous changes in the mode structure are seen for tiny changes in vertical field. Here we show that this problem is, unsurprisingly, very strongly dependent on the timescale considered: over shorter timescales the behavior is quite smooth as B_0z → 0.
• 2.  
  This system is really the simplest one could study that may have some relevance to unstratified SB turbulence simulations. In particular, the strong azimuthal field could be generated by an MRI dynamo (e.g., Käpylä & Korpi 2011; Lesur & Ogilvie 2008a), while the vertical field comes from a net-flux threading the domain.¹⁰ Of particular relevance may be the work of Longaretti & Lesur (2010), where the authors study how various characteristics SB turbulence with net magnetic flux (i.e., mean B_0z) change with parameters. Here we illustrate that trends in their turbulent simulations seem to be well matched by the linear physics, so long as nonmodal analysis techniques are used.

We illustrate these ideas in Figures 10 (short time growth) and 11 (long time growth). These each show the maximum amplification of a disturbance as a function of (k_y, k_z), at fixed azimuthal magnetic field, as the vertical field is decreased from left to right. At each point (k_y, k_z), we additionally maximize the growth over the initial orientation of the shearing wave, k_x(0); thus, the contours represent the maximum growth possible at the chosen (k_y, k_z). This is really for ease of presentation and there could certainly be interesting information in the k_x(0) structure that could be studied in future work. (Such plots are similar in spirit to hydrodynamic results given in Yecko 2004; Mukhopadhyay et al. 2005). We use a rather large dissipation ($\bar{\nu }=\bar{\eta }=1/3000$) to have some relevance to nonlinear simulations. Of course, in a SB the (k_x(0), k_y, k_z) is necessarily discretized based on the box size; nonetheless, the continuous k results presented here can either be considered as pertaining to a continuous range of box sizes or, more usefully, to different dissipation values and magnetic fields through a rescaling of the SB equations as outlined in Appendix B.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The enormous difference between Figures 10 and 11 is a stark illustration of the importance of correctly choosing the relevant timescale for a given situation. Over the long timescales illustrated in Figure 11, we are essentially seeing eigenmode behavior, with very little contribution from non-axisymmetric modes (this is more severe than it would be at lower dissipation). In addition, the change in behavior with B_0z is extreme; a change in amplification by 14 orders of magnitude with a one order of magnitude change in B_0z. In contrast, over moderate timescales t = 0 → π (Figure 10) the change with B_0z is rather smooth, even as it vanishes completely (Figure 10(d)). The growth in the case of B_0z = 0 is still substantial, with non-axisymmetric modes being amplified by a factor of 40, around a third of the amplification of the fastest growing channel mode. Note that the general trend of increasing non-axisymmetry with decreasing vertical field¹¹ matches the characteristics of nonlinear turbulence (e.g., Figure 9 from Longaretti & Lesur 2010) rather well. We have also considered the change in the mode structure with dissipation parameters (not shown) and do not see the contradictions between linear and nonlinear results that are discussed in Longaretti & Lesur (2010). Also of interest are the results of Lesur & Longaretti (2011), where it is shown numerically that the energy injection spectrum in net-flux MRI turbulence is broadly distributed across a wide range of wave-numbers. With these results in mind, it seems likely that nonmodal analyses could be useful in studying aspects of MRI turbulence from a linear standpoint, since growth over short timescales is almost certainly more relevant to turbulent situations than the t → ∞ limit explored by eigenmode analyses (Friedman & Carter 2014). Taking this to the extreme, we have proved in SB14 that the energy growth in t = 0⁺ limit is as fast as the fastest growing channel mode for all (k_y, k_z); that is, as t → 0, amplification plots such as Figures 10 and 11 become completely homogenous, with no preference for any wave-number over any other.

Another interesting case that is simple to analyze using nonmodal techniques is MRI with no background magnetic field at all. In this case, the system is spectrally stable; nevertheless, there can be significant growth over a wide range of wave-numbers, which can be sufficient to cause a transition to turbulence given large enough initial conditions (Rempel et al. 2010; Riols et al. 2013). In Figure 12, we illustrate the maximum amplification of perturbations with no background field in (a) the MHD case and (b) the well-studied hydrodynamic (HD) case with Keplerian shear. It is interesting to note the enormous change afforded by adding in magnetic perturbations, not in the magnitude of the maximum amplification, but in the range of wave numbers that can be strongly amplified. In particular, while vertical perturbations (non-zero k_z) are strongly suppressed in the HD case, these can grow rather strongly in the MHD system. The MHD growth mechanism is simple advection of the initial perturbation, as in the introductory example Section 1.1 (with some added dissipation for non-axisymmetric modes due to shearing), so the time-dependence of the perturbation growth strongly resembles Figure 1.

Zoom In Zoom Out Reset image size

Of course, one can never hope to understand the transition to turbulence with purely linear physics, and the addition of nonlinear interactions increases the complexity of the problem enormously (e.g., Shen et al. 2006; Lithwick 2007). Nonetheless, knowing a priori that the zero-net flux MHD system readily transitions to turbulence, while the hydrodynamic system appears stay laminar even at high Reynolds numbers (Lesur & Longaretti 2005; Balbus & Hawley 2006), it is interesting to note the difference from the purely linear perspective.