The magnetorotational instability (MRI) occurs when a weak magnetic field destabilizes a rotating, electrically conducting fluid with inwardly increasing angular velocity. The MRI is essential to astrophysical disc theory where the shear is typically Keplerian. Internal shear layers in stars may also be MRI-unstable, and they take a wide range of profiles, including near-critical. We show that the fastest growing modes of an ideal magnetofluid are three-dimensional provided the shear rate, S, is near the two-dimensional onset value, S_c. For a Keplerian shear, three-dimensional modes are unstable above S ≈ 0.10S_c, and dominate the two-dimensional modes until S ≈ 2.05S_c. These three-dimensional modes dominate for shear profiles relevant to stars and at magnetic Prandtl numbers relevant to liquid-metal laboratory experiments. Significant numbers of rapidly growing three-dimensional modes remainy well past 2.05S_c. These finding are significant in three ways. First, weakly nonlinear theory suggests that the MRI saturates by pushing the shear rate to its critical value. This can happen for systems, such as stars and laboratory experiments, that can rearrange their angular velocity profiles. Second, the non-normal character and large transient growth of MRI modes should be important whenever three-dimensionality exists. Finally, three-dimensional growth suggests direct dynamo action driven from the linear instability.

1. Introduction

The magnetorotational instability (MRI) is extremely important in astrophysical fluid dynamics. A weak magnetic field catalyses turbulence in a Keplerian shear by changing the stability criterion for differentially rotating flows from a negative angular momentum gradient to a negative angular velocity gradient [1,2]. This discovery explained the ubiquitous accretion onto compact objects at rates compatible with observations, and may also influence the formation of planets [3]. In discs, the gravitational field dominates the local plasma dynamics, and thus the MRI cannot significantly affect the background shear; it must saturate by other means [4]. However, stars and liquid-metal Taylor–Couette experiments have differential rotation profiles driven by much weaker stresses. Where the MRI is active in these flows, it saturates by pushing the background shear close to critical [5–7], analogous to convection mixing entropy. Stellar interiors have extremely high fluid and magnetic Reynolds numbers, but can operate at or near the critical shear rate for the MRI. This finite critical shear results from a finite-channel cut-off in the radial direction. Despite the extensive literature on accretion discs (strong shear, low dissipation) and liquid-metal experiments (weak shear, large dissipation), even the linear MRI is not well understood in the weak-shear (e.g. $S ≳ S_{c}$ ), low-dissipation regime.

Here, we investigate the stability of three-dimensional perturbations near the two-dimensional critical shear rate S_c for a nearly inviscid, ideal magnetohydrodynamic (MHD) flow. Throughout this paper, we use ‘three-dimensional’ to refer to non-axisymmetric perturbations (and their local Cartesian equivalents), while ‘two-dimensional’ refers to axisymmetric perturbations. In both cases, we retain all three components of velocity and magnetic fields. We find that the first destabilized modes are three dimensions, and thus could act as a dynamo even in the absence of secondary instability. These results also suggest that the non-normal growth of the MRI is always important, even when axisymmetric modes dominate.

2. Methods

We numerically solve the linearized MHD equations in rotating plane Couette geometry. This corresponds to a Cartesian frame rotating with angular frequency, Ω, and a linear background shear, V_y = Sx [5]; it is also the narrow-gap limit of the Taylor–Couette geometry. We cast the Navier–Stokes equation in the form

\frac{D v}{D t} + f \hat{z} \times v + S v_{x} \hat{y} + \nabla p + ν \nabla \times ω = B_{0} \partial_{z} b,

2.1

where

ω = \nabla \times v and \frac{D}{D t} = \partial_{t} + S x \partial_{y} .

2.2

We write the induction equation in terms of the x-component of the magnetic field

\frac{D b_{x}}{D t} + η (\partial_{y} j_{z} - \partial_{z} j_{y}) = B_{0} \partial_{z} v_{x},

2.3

and the x-component of the current density (j_x = ∂_yb_z − ∂_zb_y)

\frac{D j_{x}}{D t} - η \nabla^{2} j_{x} = B_{0} \partial_{z} ω_{x} - S \partial_{z} b_{x} .

2.4

We explicitly enforce the divergence-free condition on the velocity and magnetic fields

\nabla \cdot v = \nabla \cdot b = 0.

2.5

The spatial domain is a doubly periodic channel in y, z with width −d/2 ≤ x ≤ d/2. The boundary conditions are impenetrable stress-free and perfectly conducting; v_x = ω_y = ω_z = b_x = ∂_xj_x = 0 at x = ±d/2.

The main input parameters are the Coriolis parameter, f = 2Ω; the background shear rate, S = dV_y/dx < 0; and the vertical magnetic field, B₀ (in Alfvén units μ₀ρ₀ = 1). Accretion-disc modelling usually considers the Rossby number, Ro = −S/f < 1. The regime Ro ≥ 1 corresponds to purely hydrodynamical Rayleigh unstable shear. Unless otherwise stated Ro = 3/4 (Keplerian). The solution also depends on the viscosity ν and resistivity η; in our non-dimensionalization, these are equivalent to the inverse Reynolds and magnetic Reynolds numbers, respectively. That is, Re = 1/ν and Rm = 1/η.

The MRI is a weak-field instability; in the inviscid, ideal case the critical shear rate for axisymmetric instability (i.e. in two dimensions) is [5]

S_{c} = - \frac{π^{2} B_{0}^{2}}{f d^{2}} .

2.6

We use S/S_c as our instability control parameter. Because S/S_c serves as a ratio of the dimensionless magnetic field strength $B_{0}^{2} / f d$ to the length scale d, this serves to specify the background field strength.

We assume harmonic perturbations in y and z (e.g. pressure) $p = \hat{p} (x) e^{i (k_{y} y + k_{z} z) + σ t}$ . We use a complex-valued growth rate σ = γ + iω with γ, ω both real. The system reduces to a σ-eigenvalue problem of 10 first-order ODEs in x with Dirichlet boundary conditions. We pose and solve equations (2.1)–(2.5) using the Dedalus framework [8].

Our main interest is in ideal (η = 0), inviscid (ν = 0) conditions. However, we set η = ν = 10⁻⁵ to avoid critical layers in the stable solution branch. We confirmed that our results for unstable solutions are insensitive to small diffusion. For each (k_y, k_z) pair, we solve the eigenvalue problem using n_x = 128 modes; all our results are identical at double the resolution.¹

For both the ideal and non-ideal MHD equations, we solve the eigenvalue problem in the x-direction using the EigenValueProblem solver in Dedalus for a grid of N_y × N_z modes in the y- and z-directions. For most of our runs, we use a targeted, sparse eigenvalue solver to find the 15 modes closest to a guess for the maximum growth rate. In run 3, we have confirmed that dense solvers retrieve identical results. For nearly ideal runs, we use the ideal two-dimensional growth rate as input for the smallest k_y > 0 mode at each k_z and then use the output from each previous k_y as an input guess for the next mode. For those runs that have significant η, we instead use a dense solver for the k_y = 0 modes and then step forward in k_y at each k_z as before. Our solver is embarrassingly parallelized over the k_z modes. For the spectrum in figure 1, we used a dense eigenvalue solver at $(k_{y}^{max}, k_{z}^{max})$ for S/S_c = 1.002.

Figure 1. Growth rates for three-dimensional MRI modes. (a) Growth rates γ over a grid of k_y and k_z for four values of supercriticality S/S_c. The black contours are equally spaced between zero and the maximum growth rate. The grey contour highlights γ = 0; the red dots indicate the fastest growing mode. Two-dimensional modes occur at k_y = 0 on the bottom of each figure. At S/S_c = 0.640, there are no two-dimensional modes. (b) Growth rate versus k_y for k_z = 0.259 at S/S_c = 1.002 (highlighted by the blue line in (a; the orange line gives the asymptotic result from equation (3.30) and the blue line is the numerical results; the asymptotic form is valid as S/S_c → 1 and k_y ≪ 1). (c) The full discrete spectrum of the MRI for (k_y, k_z) ≃ (0.263, 0.447). The unstable mode is plotted in orange, its stable complex conjugate is blue and all other modes are grey. (d) The phase angle, $\arctan (k_{y} / k_{z})$ , as a function of S/S_c showing the exact transition to purely two-dimensional fastest modes (ϕ = 0) at $S / S_{c} ≳ 2.05$ . (Online version in colour.)

To ensure our results are converged, we have repeated runs at n_x = 256 Chebyshev modes as well as doubling the number of modes in y and z.

3. Results

(a) Growth rates and three-dimensionality

Our two main results are as follows. (i) When the first two-dimensional mode becomes unstable, there already exist three-dimensional modes with positive growth rate. (i) At sufficiently large criticality, the fastest growing mode becomes purely two dimensions. These results are pertinent in two ways. First, the MRI contains a substantial ‘Goldilocks regime’ with possible direct dynamo action, and this regime probably applies to stellar interiors and laboratory experiments. Second, our results accord with well-established results for accretion discs that expect two-dimensional primary modes. Figure 1a shows the growth rates for four values of S/S_c. At S/S_c = 1.002, the maximum growth rate (γ_max/|S| = 0.093) happens at (k_y, k_z) ≃ (0.263, 0.447). Two-dimensional instability occurs for modes with finite growth rates at k_y = 0 along the bottom of the figure. In all cases in figure 1a, the maximum growth rate does not occur along this line, indicating the dominance of three-dimensional modes. The first panel of figure 1a shows S/S_c = 0.640, which contains only three-dimensional instability: the zero growth rate contour (highlighted in grey) does not intercept the k_z-axis.

(b) Asymptotic calculation

The MRI’s preference for three-dimensional modes can be predicted analytically. We analyse equations (2.1)–(2.5) asymptotically close to S_c, and compute the leading-order correction to the growth rate from three-dimensional effects assuming S/S_c = 1 + ϵ²R, $k_{z} \sim O (ϵ)$ , $σ \sim k_{y} \sim O (ϵ^{2})$ . Here, we outline the asymptotic calculation for the leading-order correction to the two-dimensional growth rates when accounting for three-dimensional effects via finite k_y. Being in the rotating frame means the shear has no net

\int_{- d / 2}^{d / 2} V_{0, y} (x) d x = 0 ⟹ V_{0, y} (x) = S x .

3.1

The full linear ideal equations are

\begin{aligned} \partial_{t} v_{x} + S x \partial_{y} v_{x} - f v_{y} + \partial_{x} p - B_{0} \partial_{z} b_{x} = 0, \end{aligned}

3.2

\begin{aligned} \partial_{t} v_{y} + S x \partial_{y} v_{y} + (f + S) v_{x} + \partial_{y} p - B_{0} \partial_{z} b_{y} = 0, \end{aligned}

3.3

\begin{aligned} \partial_{t} v_{z} + S x \partial_{y} v_{z} + \partial_{z} p - B_{0} \partial_{z} b_{z} = 0, \end{aligned}

3.4

\begin{aligned} \partial_{x} v_{x} + \partial_{y} v_{y} + \partial_{z} v_{z} = 0, \end{aligned}

3.5

\begin{aligned} \partial_{t} b_{x} + S x \partial_{y} b_{x} - B_{0} \partial_{z} v_{x} = 0, \end{aligned}

3.6

\begin{aligned} \partial_{t} b_{y} + S x \partial_{y} b_{y} - B_{0} \partial_{z} v_{y} - S b_{x} = 0 \end{aligned}

3.7

\begin{aligned} and \partial_{t} b_{z} + S x \partial_{y} b_{z} - B_{0} \partial_{z} v_{z} = 0. \end{aligned}

3.8

The divergence-free condition for the velocity implies the same for the magnetic field. The system is second order in x because ∂_x only appears in two places. All variables take the form

p = \hat{P} (x) e^{i (ω t + k_{z} z + k_{y} y)},

3.9

where we note that, for this section, we use the complex frequency ω = −iσ for the time dependence.

To make analytical progress, define the frequency ‘parameters’

ω_{S} (x) \equiv ω + S x k_{y} and ω_{A} \equiv B_{0} k .

3.10

We can use equations (3.3)–(3.8) to find all amplitudes in terms of

{\hat{V}}_{x} (x)

and

{\hat{V}}_{x}^{'} (x)

\begin{aligned} {\hat{V}}_{y} (x) = \frac{i (k_{y} {\hat{V}}_{x}^{'} (x) + k_{z}^{2} ((f + S) ω_{S}^{2} - S ω_{A}^{2}) / (ω_{S} (ω_{S}^{2} - ω_{A}^{2})) {\hat{V}}_{x} (x))}{k_{z}^{2} + k_{y}^{2}}, \end{aligned}

3.11

\begin{aligned} {\hat{V}}_{z} (x) = \frac{i k_{z} ({\hat{V}}_{x}^{'} (x) - ((f + S) ω_{S}^{2} - S ω_{A}^{2}) / (ω_{S} (ω_{S}^{2} - ω_{A}^{2})) k_{y} {\hat{V}}_{x} (x))}{k_{z}^{2} + k_{y}^{2}}, \end{aligned}

3.12

\begin{aligned} \hat{P} (x) = \frac{i (ω_{A}^{2} - ω_{S}^{2}) ({\hat{V}}_{x}^{'} (x) - k_{y} ((f + S) ω_{S}^{2} - S ω_{A}^{2}) / (ω_{S} (ω_{S}^{2} - ω_{A}^{2})) {\hat{V}}_{x} (x))}{(k_{z}^{2} + k_{y}^{2}) ω_{S}}, \end{aligned}

3.13

\begin{aligned} {\hat{B}}_{x} (x) = \frac{ω_{A} {\hat{V}}_{x} (x)}{ω_{S}}, \end{aligned}

3.14

\begin{aligned} {\hat{B}}_{y} (x) = \frac{ω_{A} {\hat{V}}_{y} (x)}{ω_{S}} - \frac{i S ω_{A} {\hat{V}}_{x} (x)}{ω_{S}^{2}} \end{aligned}

3.15

\begin{aligned} and {\hat{B}}_{z} (x) = \frac{ω_{A} {\hat{V}}_{z} (x)}{ω_{S}} . \end{aligned}

3.16

One can substitute everything into equation (3.2) to get a second-order equation for

{\hat{V}}_{x} (x)

of the general form

{\hat{V}}_{x}^{″} (x) + 2 C_{1} (x) {\hat{V}}_{x}^{'} (x) + C_{0} (x) {\hat{V}}_{x} (x) = 0,

3.17

where

C_{1} (x) = \frac{S k_{y} ω_{A}^{2}}{ω_{S} (ω_{S}^{2} - ω_{A}^{2})} and C_{0} (x) = \frac{(f^{2} k_{z}^{2} ω_{S}^{2}) / (ω_{S}^{2} - ω_{A}^{2}) + f S k_{z}^{2} - (2 S^{2} k_{y}^{2} ω_{A}^{2}) / (ω_{S}^{2})}{ω_{S}^{2} - ω_{A}^{2}} - (k_{y}^{2} + k_{z}^{2}) .

3.18

We can eliminate the first-order term via

{\hat{V}}_{x} (x) = χ (x) ψ (x), where \frac{χ^{'} (x)}{χ (x)} = - C_{1} (x) .

3.19

Therefore

- ψ^{″} (x) + Φ (x) ψ (x) = 0, where Φ (x) = C_{1}^{'} (x) + C_{1} {(x)}^{2} - C_{0} (x) .

3.20

Finally, by putting everything together

Φ (x) = \frac{S (f k_{z}^{2} - S k_{y}^{2}) ω_{A}^{2} - f k_{z}^{2} (f + S) ω_{S}^{2}}{{(ω_{A}^{2} - ω_{S}^{2})}^{2}} + k_{y}^{2} + k_{z}^{2} .

3.21

The denominator of Φ(x) cannot vanish if ω_S = −iγ + Sxk_y, with Re(γ) ≠ 0. We want to solve the Schrödinger-type equation (3.20) with the boundary conditions

ψ (x = \pm d / 2) = 0.

3.22

The boundary conditions for ψ(x) are the same as V_x(x) because

χ (x) \propto \frac{ω_{S} (x)}{\sqrt{ω_{S} {(x)}^{2} - ω_{A}^{2}}} \neq 0 for γ \neq 0.

3.23

We solve perturbatively near the critical values for the two-dimensional instability. We rescale each parameter in terms of a bookkeeping parameter, ε,

S \to - \frac{π^{2} B_{0}^{2}}{d^{2} f} (1 + ε^{2} R), ω \to ε^{2} ω, k_{z} \to ε k_{z}, k_{y} \to ε^{2} k_{y},

3.24

keeping in mind

Ro \equiv \frac{π^{2} B_{0}^{2}}{d^{2} f^{2}} .

3.25

Expanding equation (3.20) to

O (ε^{2})

, we arrive at

Φ (x) = - \frac{π^{2}}{d^{2}} - \frac{π^{2} ε^{2}}{d^{2}} [R - \frac{d^{2}}{π^{2}} k_{z}^{2} + \frac{{Ro}^{2} B_{0}^{2} k_{y}^{2} + (1 + Ro) {(ω + f Ro x k_{y})}^{2}}{Ro B_{0}^{2} k_{z}^{2}}] .

3.26

The leading-order balance is

ψ_{0}^{″} (x) + \frac{π^{2}}{d^{2}} ψ_{0} (x) = 0 ⟹ ψ_{0} (x) = \cos (\frac{π x}{d}),

3.27

and next order is

ψ_{2}^{″} (x) + \frac{π^{2}}{d^{2}} ψ_{2} (x) = Φ_{2} (x) ψ_{0} (x) .

3.28

The solvability condition for ψ₂(x) is

\int_{- d / 2}^{d / 2} Φ_{2} (x) ψ_{0} {(x)}^{2} d x = 0.

3.29

This yields our final result

σ^{2} = - ω^{2} = \frac{Ro B_{0}^{2}}{1 + Ro} [k_{z}^{2} (R - \frac{d^{2}}{π^{2}} k_{z}^{2}) + Υ k_{y}^{2}] + \dots,

3.30

where

Υ = Ro + (1 + Ro) \frac{π^{2} - 6}{12} \approx 1.31 for Ro = 0.75,

3.31

and we have switched back to σ for time dependence. The first term in equation (3.30) results from the two-dimensional calculation. The second term is positive definite: it always leads to enhanced growth rates when k_y ≠ 0 and ultraviolet divergence in the absence of higher-order effects.

Figure 1b compares the numerical growth rates for S/S_c = 1.002 between 0 ≤ k_y ≤ 0.2 and the asymptotic approximation, showing good agreement where the latter is valid. Figure 1c shows the full spectrum for S/S_c = 1.002. The plot shows a purely growing/decaying complex-conjugate pair (orange/blue dots on the real axis) consistent with equation (3.30). The other stable modes (grey dots) are rotationally modified Alfvén waves found in left- and right-going pairs, consistent with the analytic predictions for two-dimensional stability calculations. We reject spurious eigenvalues by using eigentools² to solve equations (2.1)–(2.5) at two resolutions and retain only pairs equal to within one part in 10⁻⁶.

(c) Analysis of the transition to two dimensions

We can understand the transition from three- to two-dimensional instability semi-analytically. We do this by conducting a similar analysis to the asymptotic calculation in §3b near the fastest growing two-dimensional modes, looking for values for S/S_c where the analogous value of $Υ (S / S_{c})$ equals zero.

The growth rate can be written as a function of k_y, k_z and S/S_c near k_y ≈ 0. We therefore expand the whole problem in a power-series,

γ (\frac{S}{S_{c}}, k_{y}, k_{z}) = \sum_{n \geq 0} γ_{n} (\frac{S}{S_{c}}, k_{z}) k_{y}^{n} and ψ (x) = \sum_{n \geq 0} ψ_{n} (x) k_{y}^{n},

3.32

and use the same asymptotic techniques as before to obtain each progressive correction

\begin{aligned} ψ_{0}^{″} - Φ_{0} ψ_{0} = 0, \end{aligned}

3.33

\begin{aligned} ψ_{1}^{″} - Φ_{0} ψ_{1} = Φ_{1} ψ_{0} \end{aligned}

3.34

\begin{aligned} and ψ_{2}^{″} - Φ_{0} ψ_{2} = Φ_{2} ψ_{0} + Φ_{1} ψ_{1} . \end{aligned}

3.35

At leading order, we make use of the fact that all non-constant coefficients vanish for two-dimensional modes. The leading-order eigenfunction is ψ₀ = cos (πx/d), just as before. The following dispersion relations determine the maximum growth rate and the corresponding wavenumber

D_{0} (k_{z}) = \frac{f k_{z}^{2} (B_{0}^{2} S k_{z}^{2} + γ_{0}^{2} (f + S))}{{(B_{0}^{2} k_{z}^{2} + γ_{0}^{2})}^{2}} + \frac{π^{2}}{d^{2}} + k_{z}^{2} = 0

3.36

and

\frac{D_{0}^{'} (k_{z})}{2 k_{z}} = \frac{γ_{0}^{2} f (B_{0}^{2} k_{z}^{2} (S - f) + γ_{0}^{2} (f + S))}{{(B_{0}^{2} k_{z}^{2} + γ_{0}^{2})}^{3}} + 1 = 0.

3.37

We transform these relations into polynomials in terms of non-dimensional variables

y^{2} s^{2} {Ro}^{2} + y s (y s {Ro}^{2} + (2 - s) Ro + s) x + (2 y s Ro - s + 1) x^{2} + x^{3} = 0

3.38

and

y s (y s {Ro}^{2} + (2 - s) Ro + s) + 2 (2 y s Ro - s + 1) x + 3 x^{2} = 0,

3.39

where

x = \frac{d^{2} k_{z}^{2}}{π^{2}}, y = \frac{γ_{0}^{2}}{S^{2}} and s = \frac{S}{S_{c}} .

3.40

These can be computed from the parameters reported in figure 1 and table 1. The two polynomials are straightforward to solve numerically for x, y given s, Ro. However, we want to find critical values for s, where the fastest growing modes are two-dimensional.

Table 1. Eigenvalue calculations performed.
Collapse
run	S/S_c	ν	η	Ro	N_x	$N_{k_{y}}$	$N_{k_{z}}$	sparse/dense
1	1.02	10⁻⁵	10⁻⁵	0.75	128	200	200	sparse	resolution study
2	1.02	10⁻⁵	10⁻⁵	0.75	256	200	200	sparse
3	1.02	10⁻⁵	10⁻⁵	0.75	128	200	200	dense
4	1.02	10⁻⁵	10⁻⁵	0.75	128	512	512	sparse
5	1.02	10⁻⁵	10⁻⁵	0.75	128	100	100	sparse
6	0.2	10⁻⁵	10⁻⁵	0.75	128	200	200	sparse	S/S_c variation
7	0.3	10⁻⁵	10⁻⁵	0.75	128	200	200	sparse
8	0.4	10⁻⁵	10⁻⁵	0.75	128	200	200	sparse
9	0.5	10⁻⁵	10⁻⁵	0.75	128	200	200	sparse
10	0.64	10⁻⁵	10⁻⁵	0.75	128	200	200	sparse
11	1.002	10⁻⁵	10⁻⁵	0.75	128	200	200	sparse
12	1.002	10⁻⁵	10⁻⁵	0.75	128	200	200	dense
13	1.44	10⁻⁵	10⁻⁵	0.75	128	200	200	sparse
14	1.75	10⁻⁵	10⁻⁵	0.75	128	200	200	sparse
15	1.891	10⁻⁵	10⁻⁵	0.75	128	200	200	sparse
16	2.0	10⁻⁵	10⁻⁵	0.75	128	200	200	sparse
17	2.01	10⁻⁵	10⁻⁵	0.75	128	200	200	sparse
18	2.015	10⁻⁵	10⁻⁵	0.75	128	200	200	sparse
19	2.031	10⁻⁵	10⁻⁵	0.75	128	200	200	sparse
20	2.05	10⁻⁵	10⁻⁵	0.75	128	200	200	sparse
21	2.1	10⁻⁵	10⁻⁵	0.75	128	200	200	sparse
22	2.25	10⁻⁵	10⁻⁵	0.75	128	200	200	sparse
23	2.5	10⁻⁵	10⁻⁵	0.75	128	200	200	sparse
24	4	10⁻⁵	10⁻⁵	0.75	128	200	200	sparse
25	1.02	10⁻⁶	10⁻⁶	0.75	128	200	200	sparse	Reynolds number study
26	1.02	10⁻⁴	10⁻⁴	0.75	128	200	200	sparse
27	1.02	10⁻³	10⁻³	0.75	128	200	200	sparse
28	1.02	10⁻²	10⁻²	0.75	128	200	200	sparse
29	1.02	10⁻⁵	10⁻⁵	0.1	128	200	200	sparse	low Rossby
30	1.02	10⁻⁶	10⁻²	0.1	128	200	200	sparse	liquid-metal case

We find the higher-order γ_n(s, k_z) coefficients by applying order-by-order solvability conditions. We first find that γ₁(s, k_z) = 0 identically. This renders the first-order correction equation (3.34) solvable. The first-order eigenfunction takes the form

ψ_{1} (x) = A_{1} F (\frac{π x}{d}), where F (ϑ) = ϑ (ϑ \sin (ϑ) + \cos (ϑ)) - \frac{π^{2}}{4} \sin (ϑ),

3.41

and

A_{1} = \frac{i γ_{0} d^{3} f S k_{z}^{2} (B_{0}^{2} (f - S) k_{z}^{2} - γ_{0}^{2} (f + S))}{2 π^{3} {(B_{0}^{2} k_{z}^{2} + γ_{0}^{2})}^{3}} .

3.42

The structure of the growth rate is therefore

γ (s, k_{y}, k_{z}) = γ_{0} (s, x) + γ_{2} (s, x) k_{y}^{2} + \dots .

3.43

This is the same structure as the asymptotic calculation in §3b. Here, $Υ \propto γ_{2}$ effectively. The case with dominant two-dimensional modes corresponds to γ₂(s, x) = 0. The final result derives from the second-order solvability condition

\int_{- d / 2}^{d / 2} (Φ_{2} ψ_{0} + Φ_{1} ψ_{1}) ψ_{0} d x = 0,

3.44

which produces the degree-6 polynomial expression

\begin{aligned} y^{6} s^{6} {Ro}^{6} + y^{4} s^{5} {Ro}^{4} ((3 ξ s + 6 y - 1) Ro - 3 ξ s) x \\ + \frac{y^{3} s^{4} {Ro}^{2}}{4} ((60 y - 4 ξ (s - 8) s - 16 + 3 s^{2}) {Ro}^{2} + 2 s (4 ξ (s + 1) - 3 s) Ro - (4 ξ - 3) s^{2}) x^{2} \\ + \frac{y^{2} s^{3} Ro}{2} ((3 s^{2} - 4 ξ (s - 3) s + 40 y - 12) {Ro}^{2} + 24 ξ s Ro + (4 ξ - 3) s^{2}) x^{3} \\ + \frac{y s^{2}}{4} ((60 y - (4 ξ - 3) s^{2} - 16) {Ro}^{2} + 2 s (3 s - 4 ξ (s - 3)) Ro - (4 ξ - 3) s^{2}) x^{4} \\ + s ((6 y - ξ s - 1) Ro - ξ s) x^{5} + x^{6} = 0. \end{aligned}

3.45

The new coefficient is

ξ = \frac{π^{2} - 6}{12} = \frac{\int_{- d / 2}^{d / 2} x^{2} \cos^{2} ((π x / d)) d x}{d^{2} / π^{2} \int_{- d / 2}^{d / 2} \cos^{2} ((π x / d)) d x} \approx 0.322467.

3.46

This results from projecting the squared shear onto the leading-order eigenmodes, just as in the previous section.

We solve equations (3.38), (3.39) and (3.45) for x, y, s using Newton’s method for a variety of Ro, yielding

\begin{matrix} Ro & k_{z} d / π & γ_{max} / | S_{2 D} | & S_{2 D} / S_{c} \\ 1.00 & 0.689610 & 0.208727 & 2.09694 \\ 0.75 & 0.682819 & 0.214895 & 2.05136 \\ 0.50 & 0.676062 & 0.221309 & 2.00061 \\ 0.25 & 0.671528 & 0.228551 & 1.94963 \\ 0.10 & 0.673164 & 0.234376 & 1.92692 \\ 0.00 & 0.462324 & 0.240212 & 1.92465 \end{matrix}

The value of S/S_c = 2.05136 at the two- to three-dimensional transition when Ro = 0.75 matches the numerical results shown in figure 1 to all reported significant figures in the numerical calculation.

We conclude this section with some speculations concerning the higher-order corrections. Because of symmetry, γ₃(s, x) = 0 identically. This means that, to leading order,

γ (s, k_{y}, k_{z}) = γ_{0} (s, x) + γ_{2} (s, x) k_{y}^{2} + γ_{4} (s, x) k_{y}^{4} + \dots .

3.47

For S/S_c > S_2D/S_c, γ₂ > 0. This is sufficient to determine the most unstable three-dimensional wavenumber near the transition to dominant two-dimensional modes

k_{y}^{2} \approx - \frac{γ_{2}}{2 γ_{4}} \propto \frac{S_{2D}}{S_{c}} - \frac{S}{S_{c}} .

3.48

In general, γ₄ ≠ 0. We speculate that γ₄ < 0, though the calculation is outside the scope of this work. If there is a negative fourth-order feedback, then

k_{y} \sim \sqrt{\frac{S_{2 D}}{S_{c}}} - \frac{S}{S_{c}},

3.49

which is consistent with figure 1d near the critical point. However, a scenario where γ₄ > 0 would require behaviour that is inconsistent with our numerical calculations.

(d) Consistency with previous simulations

Numerical simulations of the MRI in accretion discs consistently show axisymmetric modes dominating the early evolution of the MRI before breaking down into three-dimensional MHD turbulence [9–11]. These ‘channel modes’ are exact non-linear solutions for x-shearing-periodic domains, and can only saturate via parasitic shear instabilities [12]. Impenetrable and stress-free boundary conditions are more applicable to stellar interiors. But even in discs, any kind of finite radial extent will cut off the unbounded growth of channel modes. We reconcile our results with earlier simulations by showing that the MRI indeed prefers two dimensions for the larger values of criticality found in disc simulations. Figure 1d shows the phase angle $ϕ = \arctan (k_{y} / k_{z})$ of the fastest growing mode as a function of S/S_c. The overall critical shear for three-dimensional modes is S_c,3D ≃ 0.102S_c. Above $S / S_{c} ≳ 2.05$ , ϕ becomes zero, indicating that axisymmetric modes have the fastest growth rates. For the fiducial run in [13], S/S_c ≃ 4.84 (most works use similar values). Thus, our theory predicts that axisymmetric modes should dominate the linear dynamics for parameters studied in prior numerical simulations. However, it may be that those simulations saw a predominance of axisymmetric modes because of the radial boundary conditions and not their large shear rates.

Even for S/S_c > 2.05 there are significant swathes of unstable three-dimensional modes with growth rates comparable to the maximum (figure 1a). Also, non-normality generically accompanies non-axisymmetry in shearing systems [14]. In non-normal systems, transient amplification can occur for stable modes, and can even cause turbulence. Therefore, significant three-dimensionality likely implies important non-normal behaviour near onset.

(e) Eigenvectors

Figure 2 shows the eigenvectors for the most unstable mode at S/S_c = 1.02 using ideal MHD. No critical layers can form for γ ≠ 0. We therefore solve equations (2.1)–(2.5) without dissipation as a second-order system in ∂_x and impose v_x = 0 at x = ±d/2. The eigenfunctions are indistinguishable from those at finite dissipation, lending additional confidence to our other results. The tilted structures of v_x and v_y as well as b_x and b_y imply non-trivial Reynolds and Maxwell stresses.

(f) Reynolds number and departure from ideal MHD

In order to ensure that our calculations are probing the ideal MHD regime that we are interested in, we have performed a series of runs (24–27) with η = ν ranging from 10⁻⁶ to 10⁻², equivalent to varying the Reynolds and magnetic Reynolds numbers Re = Rm = 10²–10⁶. Figure 3 shows the growth rates at S/S_c = 1.05 for five values of Re = Rm. For $Re ≳ 10^{4}$ , the maximum growth rate is located at the same position and the same growth rate; the unstable region expands, as expected. This suggests that our results at Re = 10⁵ are reasonably close to ideal MHD for the purposes of understanding the MRI’s preference for three-dimensional modes near shear onset.

Figure 3. Growth rates at S/S_c = 1.02 for varying diffusivities. Note that the top three panels are essentially identical, suggesting that above $Re ≳ 10^{4}$ we are reasonably close to ideal conditions. (Online version in colour.)

(g) Mean electromotive forces

We should expect the generation of a mean electromotive force (EMF) in the linear regime, i.e. 〈v × b〉 ∼ α〈B₀〉. This is significant because of the possibility of direct laminar dynamo action in regions of weak shear at large Reynolds number. Figure 4 shows the correlation (averaged over y, z) as a function of x over the domain for S/S_c = 1.02

EMF \propto \frac{{⟨ \hat{z} \cdot (v \times b) ⟩}_{y, z}}{v_{rms} b_{rms}},

3.50

where ‘rms’ denotes the average over the whole domain. The correlation is independent of the arbitrary normalization of the linear eigenvectors. The non-zero correlation shows the tendency for back reaction on the mean magnetic field.

Figure 4. Velocity–magnetic field correlations relevant to the mean EMF for the eigenvectors of the fastest growing mode at S/S_c = 1.02. (Online version in colour.)

MRI dynamos have been studied in a number of different contexts [15–17], but all except one focused on nonlinear, usually turbulent, dynamos. The only exception we are aware of [18] found exponential growth of mean magnetic fields during the linear growth phase of the MRI far from stability in numerical simulations. Our work explains this result in terms of purely linear dynamics: non-axisymmetric MRI unstable modes drive the dynamo growth of magnetic fields.

(h) Extensions to stellar and experimental conditions

Finally, we vary two important parameters: Ro and the magnetic Prandtl number Pm = ν/η. Figure 5 shows the growth rates for S/S_c = 1.02. Figure 5a shows Ro = 1/10 (addressing a possible range within stellar interiors), holding all other parameters equal to their fiducial values. The growth rates show similar behaviour to the Ro = 3/4 cases, i.e. dominated by three-dimensional modes with γ_max/|S| ≃ 0.1. Figure 5b shows ν = 10⁻⁶ and η = 0.1 with all other parameters equal to their fiducial values. This case is relevant to liquid metals, e.g. the Princeton MRI experiment [19] at Pm = 10⁻⁵. Our results are for the rotating plane Couette geometry (Taylor–Couette in the small-gap limit R₂ ≈ R₁) with highly idealized boundary conditions. It is nevertheless quite interesting that, near onset, the low-Pm MRI is only unstable to three-dimensional modes.

4. Conclusion

Our results show that three-dimensional modes grow faster than two-dimensional modes whenever the MRI is near its critical shear values. While we mainly focus on the ideal Keplerian case (Ro = 3/4); we also demonstrate robustness for low Rossby (Ro = 1/10) and magnetic Prandtl (Pm = 10⁻⁵) numbers. There are several important future directions this work suggests. First, the Sun possesses two internal shear layers with inwardly increasing shear, the high-latitude tachocline and the near-surface shear layer (NSSL). Past work has already pointed out the possibility of the small-scale MRI in the Sun using local analysis [20–22]. The NSSL, in particular, may also host slower MRI-driven dynamics, despite containing small-scale convection. The NSSL contains the strongest shear anywhere within the solar interior; Ω(r) ∝ 1/r, implying Ro = 1/2. It is therefore crucial to further elucidate the nonlinear saturation of the three-dimensional MRI, along with its robustness to convection. An important uncertainty with regard to stellar applications is the fact that S/S_c depends on the square of the magnetic field strength, a quantity that is rather uncertain. Second, our low-Pm results suggest that three-dimensional MRI modes may be the easiest to excite in liquid-metal experiments. Determining possible non-axisymmetric signatures in Taylor–Couette experiments requires follow-up work using more realistic boundary conditions and geometry. Our prior work on the axisymmetric MRI [6,7] shows only small differences between rotating plane and cylindrical Taylor–Couette geometries. We fully expect the general three-dimensional features of the MRI to persist in more complex applications.

Data accessibility

All code used for this project is available at https://github.com/jsoishi/mri_prefers_3D.

Authors' contributions

J.S.O. and G.M.V. led the project. J.S.O. prepared all figures. G.M.V. performed the asymptotic calculation. M.B. ran the bulk of the eigenvalue calculations and prepared several preliminary figures. A.S. performed preliminary eigenvalue calculations which led to the fiducial parameter choices used here. J.S.O., G.M.V., K.J.B., D.L. and B.P.B. developed Dedalus, which is used for all calculations in this paper. All authors contributed to the writing of the manuscript and physical interpretation of the results.

Competing interests

We declare we have no competing interest.

Funding

J.S.O., M.B. and B.P.B. acknowledge support from NASALWS grant no. NNX16AC92G. J.S.O. also acknowledges support from Research Corporation Scialog Collaborative Award(TDA) ID no. 24231.

Footnotes

1 See github.com/jsoishi/mri_prefers_3d for all code used in this paper.

2 https://github.com/dedalusproject/eigentools.

Acknowledgements

Computations were performed on the Leavitt cluster at the Bates College High Performance Computing Centre.

VIEW FULL TEXT DOWNLOAD PDF

This Issue

Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences cover image

January 2020

Volume 476Issue 2233

Article Information

DOI:https://doi.org/10.1098/rspa.2019.0622
PubMed:32082064
Published by:Royal Society
Online ISSN:1471-2946

History:

Manuscript received21/09/2019
Manuscript accepted18/11/2019
Published online08/01/2020
Published in print29/01/2020

License:

Citations and impact

Keywords

Subjects

Your password must have 8 characters or more and contain 3 of the following:

The magnetorotational instability prefers three dimensions

Abstract

1. Introduction

2. Methods

3. Results

(a) Growth rates and three-dimensionality

(b) Asymptotic calculation

(c) Analysis of the transition to two dimensions

(d) Consistency with previous simulations

(e) Eigenvectors

(f) Reynolds number and departure from ideal MHD

(g) Mean electromotive forces

(h) Extensions to stellar and experimental conditions

4. Conclusion

Data accessibility

Authors' contributions

Competing interests

Funding

Footnotes

Acknowledgements

This Issue