# The Weakly Nonlinear Magnetorotational Instability in a local geometry

###### Abstract

The magnetorotational instability (MRI) is a fundamental process of accretion disk physics, but its saturation mechanism remains poorly understood despite considerable theoretical and computational effort. We present a multiple scales analysis of the non-ideal MRI in the weakly nonlinear regime – that is, when the most unstable MRI mode has a growth rate asymptotically approaching zero from above. Here, we develop our theory in a local, Cartesian channel. Our results confirm the finding by Umurhan et al. (2007) that the perturbation amplitude follows a Ginzburg-Landau equation. We further find that the Ginzburg-Landau equation will arise for the local MRI system with shear-periodic boundary conditions when the effects of ambipolar diffusion are considered. A detailed force balance for the saturated azimuthal velocity and vertical magnetic field demonstrates that even when diffusive effects are important, the bulk flow saturates via the combined processes of reducing the background shear and rearranging and strengthening the background vertical magnetic field. We directly simulate the Ginzburg-Landau amplitude evolution for our system and demonstrate the pattern formation our model predicts on long length and time scales. We compare the weakly nonlinear theory results to a direct numerical simulation of the MRI in a thin-gap Taylor Couette flow.

^{†}

^{†}journal: ApJ\correspondingauthor

S. E. Clark

0000-0002-7633-3376]S. E. Clark

0000-0001-8531-6570]Jeffrey S. Oishi \nocollaboration

## 1 Introduction

For matter to accrete from a disk onto a central object, angular momentum must be transported radially outward in the disk. The transport mechanism is likely turbulent, as molecular viscosity alone cannot account for the needed angular momentum transfer, and likely magnetic, as this turbulence is excited even in hydrodynamically stable disks (Shakura & Sunyaev, 1973). Discovered by Chandrasekhar (1960) and Velikhov (1959) in a global geometry, the magnetorotational instability (MRI) was subsequently rediscovered and applied to accretion disks by Balbus & Hawley (1991). Since then, the MRI remains the leading explanation for rapid angular momentum transport in astrophysical disks. The instability in its simplest geometry arises when a differentially rotating disk is threaded by a vertical magnetic field. The presence of the magnetic field linearly destabilizes the disk gas, driving turbulence and angular momentum transport (e.g. Hawley et al., 2011; Parkin, 2014; Parkin & Bicknell, 2013). The MRI likely plays a role in a diverse host of astrophysical systems, including protoplanetary disks (e.g. Bai, 2015) and black hole accretion disks (e.g. Schnittman et al., 2013), as well as stellar interiors (e.g. Wheeler et al., 2015). Despite its importance, many aspects of the MRI remain poorly understood. In particular, the nonlinear saturation mechanism for the MRI is an open question, and a formidable challenge. MRI saturation has been tackled almost exclusively with simulation, with a few notable exceptions detailed below. In this work we analytically investigate the weakly nonlinear saturation of the MRI.

Weakly nonlinear analysis is a perturbative method used to examine the asymptotic behavior of a system near threshold – that is, when the system is just barely unstable to its most unstable mode. The analytical technique follows the multiscale evolution of fluid variables in a perturbation expansion, allowing the controlled interaction of modes between orders in a perturbation series (Bender & Orszag, 1978). Weakly nonlinear analysis can be a powerful technique for analytically examining systems which in their full generality exhibit such complicated nonlinear behavior that their study is relegated primarily to the simulation domain. The MRI is one such phenomenon: while there is a rich literature analytically examining the linear MRI, analytical treatments of the nonlinear system are relatively few. The weakly nonlinear treatment of the MRI was pioneered by Knobloch & Julien (2005) and Umurhan et al. (2007b, hereafter URM07; see also ). The latter authors undertook the first weakly nonlinear analysis of the MRI in a thin-gap Taylor-Couette (TC) flow with strong dissipation (as is appropriate to liquid metal experiments), and found that the marginal MRI system approaches saturation in a manner analogous to that of Rayleigh-Bénard convection. Weakly nonlinear analysis was instrumental in our understanding of Rayleigh-Bénard convection saturation (Newell & Whitehead, 1969), and the similarities between convection and the local MRI are the result of important shared symmetries between the systems. The success of URM07 in modeling the MRI system near threshold merits further consideration, but we are unaware of any other attempts to expand upon their theoretical framework. In this work we rederive the theory of URM07, and expand upon their findings. Our focus here is on fully characterizing the local MRI system, both for independent theoretical interest and to have a robust comparison point for extensions of this theory into more complicated geometries. In a companion paper we derive for the first time the weakly nonlinear theory of the standard and helical MRI in a global, cylindrical TC flow (Clark & Oishi, 2017b, hereafter Paper II). The thin- and wide-gap treatments complement one another theoretically, and both are important regimes for comparison with simulation.

This work examines TC flow in the thin-gap regime, an idealization in which the radial extent of the channel is very small compared to its distance from the center of rotation, i.e. (r_{2}-r_{1})\ll\frac{1}{2}(r_{1}+r_{2}) where r_{1} and r_{2} are the radii of the inner and outer flow boundaries, respectively. The thin-gap approximation eliminates curvature terms, so the domain geometry is Cartesian rather than cylindrical. The excluded curvature terms have an explicit dependence on r, so they make the problem more challenging both analytically and numerically. In particular, in the wide-gap geometry (i.e. true Taylor-Couette flow) the base angular velocity is a function of r, where in the thin-gap approximation the shear flow reduces to a linear profile. The equations of motion in thin-gap TC flow are thus identical to the MRI in a local shearing box, which differs from our fiducial setup only in the application of periodic boundary conditions.

We note several other important analytical studies of MRI saturation. Knobloch & Julien (2005) analyze the MRI in the strongly nonlinear regime, by following the already-developed MRI modes into asymptotic saturation. They consider a thin-gap regime as well, and so their theory may be considered the strongly nonlinear analogue to the one developed here. Vasil (2015) examines the weakly nonlinear MRI in a thin-gap regime in a minimal model, finding deep mathematical similarities between the MRI system and the elastodynamic instability of a buckling beam. We discuss these results and their relation to ours in Section 6. Several authors have investigated the behavior of the MRI when the boundary conditions are shear periodic, and so the MRI has no mechanism by which to modify the background shear flow profile. In this approximation linear MRI growth is dominated by channel modes, a type of MRI mode that, for periodic boundary conditions, are exact solutions of both the linear and nonlinear MRI equations (Goodman & Xu, 1994). In this regime the MRI saturates via parasitic instabilities, which feed off and destroy the primary MRI modes. Analytical investigation of this case reveals that MRI saturation can be caused by parasitic Kelvin-Helmholtz and tearing mode instabilities, depending on parameter regime (Pessah, 2010). The theory of MRI channel mode parasites is robust (see also Pessah & Goodman, 2009; Latter et al., 2010; Rembiasz et al., 2016), but their importance may be overestimated by the local approximation (Latter et al., 2015), and not germane to global analyses like the one presented here. Latter et al. (2015) gives a detailed analysis of the relation between local and global linear MRI modes. In this work we describe the applicability of our weakly nonlinear theory to shear-periodic boxes. We find that under certain conditions the weakly nonlinear mode interaction described here may provide an alternative MRI saturation mechanism in the shearing box that does not rely on parasitic instabilities.

We begin with an overview of our basic model equations for the local MRI in Section 2 and then describe our weakly nonlinear analysis and give results for the thin-gap TC flow in Section 3. In Section 4 we detail the conditions under which our theory applies to the case where the boundary conditions are shear periodic, namely the consideration of ambipolar diffusion. We compare our results to a direct numerical simulation in Section 5. We then place our results in the context of previous results from both analytic and computational studies in Section 6 and draw conclusions in Section 7.

## 2 Equations

The evolution of a conducting fluid is governed by the momentum and induction equations,

\displaystyle\begin{split}&\displaystyle\partial_{t}\mathbf{u}\,+\,\mathbf{u}% \cdot\nabla\mathbf{u}\,=\,-\frac{1}{\rho}\nabla P\,-\,\nabla\Phi\,+\,\frac{1}{% \rho}\left(\mathbf{J}\times\mathbf{B}\right)\\ &\displaystyle+\,\nu\nabla^{2}\mathbf{u}\,-\,2\mathbf{\Omega}\times\mathbf{u}% \,-\,\mathbf{\Omega}\times\left(\mathbf{\Omega}\times\mathbf{r}\right),\\ \end{split} | (1) |

\partial_{t}\mathbf{B}=\nabla\times\left(\mathbf{u}\times\mathbf{B}\right)+% \eta\nabla^{2}\mathbf{B}, | (2) |

where P is the gas pressure, \nu is the kinematic viscosity, \eta is the microscopic diffusivity, \nabla\Phi is the gravitational force per unit mass, and the current density is \mathbf{J}=\nabla\times\mathbf{B}. Equations 1 and 2 are subject to the incompressibility and magnetic solenoid constraints,

\nabla\cdot\mathbf{u}=0 | (3) |

\nabla\cdot\mathbf{B}=0. | (4) |

We axisymmetrically perturb all three vector components of each of the fluid quantities. We nondimensionalize the equations, with lengths nondimensionalized by L, time by \Omega_{0}, velocities by \Omega_{0}L, magnetic fields by B_{0}, and pressure by \Omega_{0}^{2}L^{2}\rho_{0}, where L is the channel width, \Omega_{0} is the rotation rate at the center of the channel, and \rho_{0} is the constant pressure in the base state (see Figure 1). We define the Reynolds number, \mathrm{Re}\equiv{\Omega_{0}L^{2}}/{\nu}, magnetic Reynolds number, \mathrm{Rm}\equiv{\Omega_{0}L^{2}}/{\eta}, and Cowling number, \mathrm{Co}\equiv{2v_{A}^{2}}/{\Omega_{0}^{2}r_{0}^{2}}, where the Alfvén speed v_{A} is v_{A}^{2}={B_{0}^{2}}/{\rho_{0}}. The fluid symbols \mathbf{u}, \mathbf{B}, etc. will henceforth be used to refer to the nondimensional, perturbed quantities.

We define the streamfunction \Psi and flux function A, where A is the familiar two-dimensional vector potential. \Psi and A are scalar fields. The curl of \Psi and the curl of A are defined as the velocity and magnetic field perturbation, respectively, and so \Psi and A automatically satisfy our constraints (Equations 3 and 4).

\Psi and A are thus related to the velocity and magnetic field perturbations, respectively, as

\mathbf{u}\,=\,\left[\begin{matrix}\partial_{z}\Psi\\ u_{y}\\ -\partial_{x}\Psi\\ \end{matrix}\right], | (5) |

\mathbf{B}\,=\,\left[\begin{matrix}\partial_{z}A\\ B_{y}\\ -\partial_{x}A\\ \end{matrix}\right]. | (6) |

Our final equation set is

\partial_{t}\nabla^{2}\Psi\,-\,2\partial_{z}u_{y}\,-\,\mathrm{Co}B_{0}\partial% _{z}\nabla^{2}A\,-\,\frac{1}{\mathrm{Re}}\nabla^{4}\Psi\,=\,\mathrm{Co}J\left(% A,\nabla^{2}A\right)\,-\,J\left(\Psi,\nabla^{2}\Psi\right)\, | (7) |

\partial_{t}u_{y}\,+\,\left(2-q\right)\Omega_{0}\partial_{z}\Psi\,-\,\mathrm{% Co}B_{0}\partial_{z}B_{y}\,-\,\frac{1}{\mathrm{Re}}\nabla^{2}u_{y}\,=\,\mathrm% {Co}J\left(A,B_{y}\right)\,-\,J\left(\Psi,u_{y}\right) | (8) |

\partial_{t}A\,-\,B_{0}\partial_{z}\Psi\,-\,\frac{1}{\mathrm{Rm}}\nabla^{2}A\,% =\,J\left(A,\Psi\right) | (9) |

\partial_{t}B_{y}\,+\,q\Omega_{0}\partial_{z}A\,-\,B_{0}\partial_{z}u_{y}\,-\,% \frac{1}{\mathrm{Rm}}\nabla^{2}B_{y}\,=\,J\left(A,u_{y}\right)\,-\,J\left(\Psi% ,B_{y}\right), | (10) |

where J is the Jacobian operator,

J\left(f,g\right)\equiv\partial_{z}f\partial_{x}g-\partial_{x}f\partial_{z}g, | (11) |

and q\equiv-d\ln\Omega/\ln R=3/2 is the dimensionless shear parameter defining a rotation profile \Omega(r)=\Omega_{0}(r/r_{0})^{-q}, such that the background velocity profile is u_{0}=-q\Omega_{0}x.

The weakly nonlinear regime is where the MRI system is nonlinearly unstable to only the most unstable mode of the linear solution. We find the marginal state, where the most unstable linear MRI mode neither grows nor decays, for a set of dimensionless parameters, and then destabilize the system. We examine the system for fiducial parameters comparable to URM07, namely \mathrm{Pm}=1.0\times 10^{-3}, \mathrm{Co}=0.08, q=1.5. The system is marginal for a critical wavenumber k_{c}=0.75 and a critical magnetic Reynolds number \mathrm{Rm}_{c}=4.9.

Because we nondimensionalize B by the magnitude of the background field strength, B_{0}\equiv 1 in Equations 7 - 10. To excite the weakly nonlinear MRI, we tune the background magnetic field away from stability. We do so by substituting B=B_{0}\left(1+\epsilon^{2}\right). The degree of departure from the marginal state is measured by the small parameter \epsilon. An \mathcal{O}\left(\epsilon^{2}\right) strengthening of the background magnetic field destabilizes a finite band of wave modes with a width of \mathcal{O}\left(\epsilon\right), which interact nonlinearly. We note that this definition of \epsilon is opposite in sign to nearly all previous works (e.g. Umurhan et al., 2007a, b). Because in the ideal limit, the MRI can be tuned into instability by setting B_{0} to its critical value and then decreasing its value, it is natural to consider \epsilon^{2} as a weakening of the background field (as is done correctly in Vasil, 2015, for example). However, as we show in figure 2, for the dissipative case with \eta,\nu\neq 0, when all other parameters are critical, decreasing B_{0} leads to stability, while increasing it pushes the system into instability. Figure 2 is symmetric about B_{0}=0, as it must be, since the MRI is insensitive to the sign of the background field. There are several places at which the derivative of \gamma appears discontinuous; this is not physical but rather reflects the fact that we define \gamma as the growth rate of the most unstable mode. That is, it is the maximum real part of the eigenvalues of the linearized system (e.g. equation 13 with \mathbf{N}=0). Because there are four wave families in rotating incompressible MHD, each modified differently by changing B_{0}, when the growth rates of the individual modes cross, there appear piecewise continuous solutions. We highlight one such point in the inset in figure 2, where the MRI mode becomes more stable than another mode which is always stable. Since all of these piecewise discontinuities are below \gamma=0, they do not affect the analysis here.

The destabilizing substitution is made, and Equations 7 - 10 are rewritten such that the fluid variables are contained in a state vector

\mathbf{V}=\left[\Psi,u_{y},A,B_{y}\right]^{\mathrm{T}}. | (12) |

This yields the system of equations

\mathcal{D}\partial_{t}\mathbf{V}+\mathcal{L}\mathbf{V}+\epsilon^{2}\widetilde% {\mathcal{G}}=\mathbf{N}, | (13) |

where we leave the definition of the matrices \mathcal{D}, \mathcal{L}, and \widetilde{\mathcal{G}} to Appendix A, and the detailed form of the nonlinear vector \mathbf{N} to Appendix B. We solve this system subject to no-slip, perfectly conducting radial boundary conditions, defined as

\Psi=\partial_{x}\Psi=u_{y}=A=\partial_{x}B_{y}=0. | (14) |

## 3 Weakly Nonlinear Analysis

We conduct a formal multiple scales analysis of this system. Our perturbations are characterized in terms of fast- and slow-moving variables, that we treat as independent in order to simultaneously track the evolution of the system on two scales. The relative scalings of the fast and slow variables are chosen such that each of the temporal and spatial eigenvalues appear at the same lowest order in the linear dispersion relation (Appendix C). The scalings are

X\equiv\epsilon x,\,\,Y\equiv\epsilon y,\,\,Z\equiv\epsilon z,\,\,T\equiv% \epsilon^{2}t. | (15) |

Note that these are the same scalings as apply to Rayleigh-Bénard convection and hydrodynamic TC flow. Our x dimension, the direction of angular momentum transport, is analogous to the direction of temperature transport in the convection problem. In analogy to these problems, we posit slow variation in both Z and T. Each operator in Equations 7 - 10 is expanded to reflect these scalings – for instance, \partial_{z} becomes \partial_{z}+\epsilon\partial_{Z}.

The multiple scale dependencies of our solution are encoded into an ansatz for the linear MRI solution at marginality,

\mathbf{V_{1}}=\alpha(T,Z)\mathbb{V}_{11}(x)e^{ik_{c}z}+c.c.+\beta(T,Z)\mathbb% {U}_{11}(x)\\ | (16) |

where \alpha(T,Z) is a slowly-varying amplitude and c.c. denotes the complex conjugate. The x dependence is contained in \mathbb{V}_{11}=(\Psi_{11},u_{11},A_{11},B_{11})^{\mathbf{T}}, and must be solved subject to the radial boundary conditions. The periodic vertical boundary conditions allow us to posit the z dependence, where k_{c} is the value of the vertical wavenumber at marginality. As noted by URM07, there exists a spatially constant neutral mode solution to the B_{y} equation, with \mathbb{U}_{11}=(0,0,0,1)^{\mathbf{T}}. The amplitude \beta(T,Z) encodes the slow evolution of this mode. This spatially constant mode cannot contribute to the nonlinear saturation of the MRI because all of the nonlinearities involve derivatives. The long-term evolution of \beta(T,Z) is described by a simple diffusion equation that decouples from \alpha(T,Z), and so we neglect it in what follows.

The state vector is expanded in a perturbation series in orders of \epsilon,

\mathbf{V}=\epsilon\mathbf{V_{1}}+\epsilon^{2}\mathbf{V_{2}}+\epsilon^{3}% \mathbf{V_{3}}+h.o.t. | (17) |

Our perturbed system is then expressed order by order as

\displaystyle\mathcal{O}(\epsilon): | \displaystyle\mathcal{L}\mathbf{V_{1}}+\mathcal{D}\partial_{t}\mathbf{V_{1}}=0 | (18) | |||

\displaystyle\mathcal{O}(\epsilon^{2}): | \displaystyle\mathcal{L}\mathbf{V_{2}}+\widetilde{\mathcal{L}}_{1}\partial_{Z}% \mathbf{V_{1}}=\mathbf{N_{2}} | (19) | |||

\displaystyle\mathcal{O}(\epsilon^{3}): | \displaystyle\mathcal{L}\mathbf{V_{3}}+\mathcal{D}\partial_{T}\mathbf{V_{1}}+% \widetilde{\mathcal{L}}_{1}\partial_{Z}\mathbf{V_{2}} | (20) | |||

\displaystyle+\widetilde{\mathcal{L}}_{2}\partial_{Z}^{2}\mathbf{V_{1}}+% \widetilde{\mathcal{G}}\mathbf{V_{1}}=\mathbf{N_{3}}. | (21) |

The partial differential equations that comprise Equations 18 to 21 are solved in succession. The practical advantage of our ansatz construction (Equation 16) is clear: the separable x-dependence means that the radial boundary conditions are solved in only one dimension. Thus our analytical framework is able to side-step many of the resolution issues faced by multidimensional simulations. We are able to resolve even small-scale structure in the boundary layers of our domain, because we need only resolve it in one dimension. We solve the radial component of each equation using the open source pseudospectral code Dedalus. We compute the radial components on a grid of Chebyshev polynomials, as is appropriate for bounded one-dimensional domains (Boyd, 2001, e.g.). The nonuniform spacing of the Chebyshev grid allows us to resolve the boundary layers well on a 128-point grid.

To close the perturbation series we enforce a solvability criterion on Equation 21 (see Appendix A). This leads to an amplitude equation for \alpha(T,Z) that governs the slow length- and timescale evolution of the system. This amplitude equation is

\partial_{T}\alpha=b\alpha+h\partial_{Z}^{2}\alpha-c\alpha\left|\alpha^{2}% \right|, | (22) |

a real Ginzburg-Landau equation. The saturated solution to Equation 22 is evidently \alpha_{saturation}=\pm\sqrt{b/c}. We plot the first order, second order, and total perturbation structure of the fluid variables in Figures 4 and 4 with a constant \alpha_{saturation}. This is the Ginzburg-Landau equation that was previously found by URM07. Those authors investigated the behavior of this MRI system as a function of \mathrm{Pm}. By analyzing the system over several orders of magnitude in \mathrm{Pm}, we reproduce the URM07 result that the analytic saturation amplitude scales as \alpha_{saturation}^{2}\propto\mathrm{Pm}^{4/3} in a thin-gap geometry when \mathrm{Pm}\ll 1.

## 4 Shearing box and ambipolar diffusion

Many studies of the MRI consider the instability in a shearing box, i.e. a wall-less local approximation that is meant to represent a small section of a disk. The shearing box is the limit in which Equations 7 - 10 are subjected to shear periodic radial boundary conditions rather than Equation 14 (e.g. Regev & Umurhan, 2008). The periodic nature of the shearing box allows us to decompose the fluid perturbations into Fourier modes proportional to e^{ik_{x}x+ik_{z}z}. This makes the shearing box MRI straightforward to treat analytically. However, as noted above, the fastest-growing linear MRI modes in the shearing box are also exact solutions of the nonlinear MRI equations – that is, J(\hat{\psi_{0}},\hat{\psi_{0}})=J(\hat{\psi_{0}},\nabla^{2}\hat{\psi_{0}})=0 for \hat{\psi_{0}}\propto e^{ik_{x}x+ik_{z}z}. While this may be an appealing trait for analytic simplicity, it leads to the unphysical conclusion that the fastest growing modes will never nonlinearly interact (Goodman & Xu, 1994). This ‘nonlinear property’ will not be satisfied for two MRI modes with nonparallel wavenumbers, but with vertically periodic boundary conditions and a vertical background magnetic field the most unstable mode has a strictly axial vertical wavenumber. Thus a formal weakly nonlinear analysis cannot be conducted, as the most unstable mode will never nonlinearly interact with itself or its complex conjugate. Similarly, we cannot analytically examine interactions between MRI channel modes and damped eigenmodes belonging to other wave families. This is analytically examined for other plasma instabilities by tracking the amplitudes of growing, marginal, and damped eigenmodes simultaneously (e.g. Makwana et al., 2011). While the shearing box approximation allows the projection of the perturbed MRI equations into the basis set of linear eigenmodes, nonlinear coupling between modes will remain zero.

The nonlinear property of primary MRI modes in the shearing box motivates the addition of radial boundaries, such that the nonlinear evolution of the weakly nonlinear MRI can be properly considered. It also raises the question of whether some additional nonlinear mechanism can be introduced such that the fastest-growing modes are no longer nonlinear solutions to the shearing box equations. It has already been shown that the Hall effect does not negate the nonlinear property of primary MRI modes (Kunz & Lesur, 2013). However, it seems to have been overlooked in the literature that these linear modes are not solutions of the nonlinear ambipolar diffusion term, which is proportional to

\nabla\times((\mathbf{J}\times\mathbf{B})\times\mathbf{B}). | (23) |

Furthermore, the radial wavenumber of the fastest-growing linear MRI mode in a shearing box with ambipolar diffusion is nonzero when a constant azimuthal background field is considered in addition to an axial one (Kunz & Balbus, 2004). This means that, in the presence of ambipolar diffusion, we can derive the weakly nonlinear envelope equation for the MRI in the shearing box. Ambipolar diffusion adds both linear and nonlinear terms to Equations 18 to 21, but does not change their Z or T dependence. The constant azimuthal background field component does not contribute to any other terms in the local MRI equations. Thus, the slow-scale evolution of the MRI in a shearing box with ambipolar diffusion is also governed by a Ginzburg-Landau equation.

The Ginzburg-Landau form of the amplitude equation can be found in any system with Euclidean symmetry and a quadratic maximum in growth rate with respect to the wavenumber (Hoyle, 2006). In this case, the Euclidean symmetry comes from axisymmetry in the x-z plane, and the quadratic maximum is a consequence of the linear dispersion relation given in Appendix C. In Paper II, we show that the same symmetry occurs in the axisymmetric global geometry as well. The Ginzburg-Landau equation arises due to symmetries in the local MRI equations, irrespective of the boundary conditions to which they are subjected. This means that the local MRI is able to saturate via nonlinear mode interaction so long as the primary MRI modes are not exact solutions of the nonlinear terms. This can be achieved by considering the effects of ambipolar diffusion when the boundary conditions are shear periodic, or by enforcing wall-like radial boundary conditions. Both constructions require the most unstable mode to have nonconstant radial structure. Physically, this radial variation impedes the free exchange of angular momentum facilitated by the uniform stretching of channel modes.

## 5 Direct Numerical Simulation

Here, we make a preliminary test of our weakly nonlinear theory by comparing it to direct numerical simulation. Using Dedalus, we solved the full, nonlinear equations 7 - 10 with all parameters (\mathrm{Rm}, Q) equal to their critical values except the background magnetic field, which we set to B_{0}=1+\epsilon^{2}. We thus drive the system MRI unstable in the same way as in our theory. The computational requirements of low \mathrm{Pm} simulations are quite intense in both time and space. Despite being virtually smooth, the solutions require a resolution of 192\times 1536 grid points at \mathrm{Pm}=10^{-2}. Because the system has such a small growth rate, it takes hundreds of orbits for the system to reach saturation, as compared to the few orbits typical of high \mathrm{Rm} simulations (e.g. Lesur & Longaretti, 2007). As a result, we make our comparison at \mathrm{Pm}=10^{-2}, which provides a good tradeoff between probing relatively low \mathrm{Pm} while keeping the computational time for these exploratory simulations modest.

We initialize the runs with the linear eigenvectors of the MRI unstable mode (also computed by Dedalus; see section 3) multiplied by an initial amplitude A_{0}=10^{-3}. Doing so requires considerably less run time, as the MRI unstable mode starts growing immediately from A_{0}. By contrast, initializing random noise in \psi with amplitude A_{0} would give the unstable mode a much smaller amplitude. Nevertheless, we have confirmed that simulations with eigenvector initial conditions have similar evolutions to those with noise initial conditions once each enter linear growth.

We analyze the average energy and angular momentum transport in the simulation domain (Figure 6). The saturation amplitude predicted by the weakly nonlinear theory depends on the choice of normalization of the linear eigenvectors. The eigenvectors of the linear problem are only determined up to an arbitrary normalization, and the nonlinear coefficient of the Ginzburg-Landau equation is sensitive to this normalization. The undetermined factor is typically assigned by comparison with direct numerical simulation or laboratory experiment (e.g. Deyirmenjian, 1997; Recktenwald et al., 1993). Here we determine the constant by requiring that the maximum amplitude of B_{y} be equal in both theory and simulation. With this normalization choice all plotted quantities agree to within \sim 25\%. The theory and simulation are thus in reasonably good agreement considering that the weakly nonlinear theory applies rigorously to a channel of infinite height, while the simulation was carried out in a box with a vertical extent of only two critical wavelengths. We defer further comparisons between simulation and theory, including an analysis of the effect of the box height on the simulated flow, to future work.

## 6 Discussion

Here our focus is on a physical description of the saturation mechanism. Figure 7 shows saturated radial profiles of u_{0}-u_{y}=-q\Omega_{0}x-u_{y} and each term in the steady state force balance (i.e. Equation 8 with \partial_{t}u_{y}=0). In the bulk of the fluid away from the boundary layers, the saturated state shows reduced shear, with little diffusive contribution. This demonstrates that even in a case where diffusive effects are important, the bulk of the fluid saturates by balancing shear and magnetic tension. As discussed at length in Vasil (2015), when diffusive effects are not important, it is impossible to rearrange momentum without also rearranging the magnetic field. The Vasil (2015) model demonstrates saturation without diffusive effects; our results show that outside of the boundary layers, a simultaneous rearrangement of momentum and field occurs. In the boundary layers, the nonlinear advection balances viscous dissipation.

Figure 8 shows B_{z} and the terms corresponding to steady state inductive balance (\partial_{x} of Equation 10 with \partial_{t}B_{y}=0). Here, the instability acts to push the magnetic field toward the boundaries in both the bulk and the boundary layers. The radial average of the saturated B_{z} is B_{0}, i.e. B_{z} is marginally stable. Ebrahimi et al. (2009) considered the saturation of a single, strongly super-critical MRI mode allowed to interact nonlinearly only with itself and the mean. They considered two important cases, one in which the mean flow was forced to remain at its initial, quasi-Keplerian state for all time, and one in which the background flow was allowed to evolve. This is a crucial difference between the shearing box and our narrow-gap TC flow: perturbations in our simulation can adjust the background flow, whereas in a shearing box, the shear periodicity forbids perturbations from affecting the mean flow. In the case with a freely evolving background flow, Ebrahimi et al. (2009) found a saturated state quite similar to ours: field pushed to the boundaries, and a reduction in shear in the bulk of the flow. Their flows have less pronounced boundary layers, likely because of their much larger \mathrm{Pm}=0.1-1.

In the high \mathrm{Re} and \mathrm{Rm} limit, Vasil (2015) derives an amplitude equation considerably different than the one found here. By averaging in the z direction, the author computes a mean-field equation with striking similarity to the buckling of an elastic beam under load. The most salient feature of this equation is its non-local character. Unlike the present work, which focuses on Keplerian rotation profiles with q=3/2 with a critical background magnetic field strength, Vasil (2015) focuses on a fixed field strength and a weakly destabilized shear profile. These differences are minor, however: the destabilizing parameter \epsilon enters the analysis in the same quadratic proportion. Whether and how Vasil (2015)’s amplitude equation is equivalent to our own in the limit of dynamically important resistive and viscous effects is beyond the scope of this work. Nevertheless, the author identifies the nonlinear term responsible for saturation as consisting of flux and field transport and notes these are the only mechanisms able to produce saturation. Our results likewise demonstrate a combination of flux and field transport in the comparable region of our domain. This suggests that despite our formulation displaying different saturation dynamics (Ginzburg-Landau in our case; a network of coupled Duffing oscillators in Vasil 2015), there may indeed be an underlying unification.

The real Ginzburg-Landau equation describes the amplitude behavior of our system close to threshold. Although the form of the equation is generic to many systems, its coefficients depend on the specific physics of our system and govern its detailed evolution (see Appendix A). We simulate the evolution of the MRI amplitude equation by solving Equation 22 on a Fourier basis in Z using Dedalus. We initialize uniform random noise of amplitude -10^{-3} to +10^{-3} in Z, and timestep the system using a four-stage, third-order Runge-Kutta integrator. We evolve the system for 100\Omega_{0}^{-1} in timesteps of 0.02\Omega_{0}^{-1}. Results are shown in Figures 10 and 10, where the amplitude and phase structure over the vertical domain is plotted for every 20 timesteps. The system quickly organizes itself into rolls in Z bounded by the analytic saturation amplitude \alpha_{s}=\sqrt{b/c}. The specific geometry depends on the number of critical wavelengths \lambda_{crit}=2\pi/k_{c} that are initialized in Z. Figure 10 shows that a system with a height equal to two critical wavelengths will be modulated by simple rolls of sinusoidal amplitude. The saturation amplitude pattern becomes more complicated when more modes are allowed to interact. Figure 10 shows the evolution of a system of height 10\lambda_{crit}. While still bounded by \alpha_{s}, the saturation amplitude exhibits a nonlinear phase geometry due to the nonlinear interaction of modes in Z.

The weakly nonlinear theory predicts that the amplitude of the system is bounded by the saturation amplitude \alpha_{s}=\sqrt{b/c}, where b and c are coefficients corresponding to the linear growth term and nonlinear term of the Ginzburg-Landau equation, respectively. The coefficient b comes from the interaction between the background magnetic field and the linear MRI solution. The coefficient c describes the third-order nonlinear interaction between terms in the perturbation series. Physically, we see that the saturation amplitude is controlled by the strength of the mode interaction within our finite band of unstable modes. We stress that while the third-order nonlinear terms in the walled TC flow are strongly influenced by the boundary layers, this is not generically true of the MRI system. Indeed, in the shearing box MRI with ambipolar diffusion (the case sketched out in Section 4), boundary layers are impossible in the shear periodic flow. In this case the third-order nonlinear behavior of the system includes three-mode interactions from the cubic nonlinearity in the ambipolar diffusion term.

Figure 5 shows the total stress u_{x}u_{y}-\mathrm{Co}B_{x}B_{y} for the \mathrm{Pm}=10^{-2} model with \epsilon=0.5. The stress shows significantly more structure throughout the domain than the variables u_{x}, u_{y}, B_{x} and B_{y} that comprise it, demonstrating that a non-trivial correlation exists even in the weakly non-linear state. As in simulations at higher \mathrm{Rm}, figure 6 shows that the Maxwell stress dominates over the Reynolds stresses even though the kinetic energy significantly exceeds the magnetic energy.

## 7 Conclusion

In this paper we construct a weakly nonlinear analysis of the MRI using multiple scales analysis, leading to a real Ginzburg-Landau equation for the nonlinear amplitude, confirming the previous results of Umurhan et al. (2007b). We also confirm their results for the scaling of the analytic saturation amplitude with \mathrm{Pm}. We extend their results by constructing a detailed force and inductive balance for the saturated u_{y} and B_{z} components. In doing so, we find the saturated state is a complex balance in which reduction of shear and amplification and redistribution of B_{z} combine to saturate the instability. We perform numerical simulations of the amplitude equation and a direct numerical simulation of the MRI system. Using the former, we demonstrate that complex patterns can organize the flow on long length scales Z, though the maximum magnitude of the amplitude \alpha is well predicted by the steady state solution. The latter show that there is rough agreement for both total energy and average angular momentum transport between the weakly nonlinear theory and simulation for a representative case at \mathrm{Pm}=10^{-2}. We defer a full comparison between theory and simulation to later work. We describe the application of shear-periodic boundary conditions to the local MRI and find that with the inclusion of certain nonideal physical effects, namely ambipolar diffusion, our theory points to a new saturation avenue for the MRI in a shearing box. In Paper II, we make use of the techniques developed here to extend the weakly nonlinear analysis of the MRI to a full cylindrical geometry appropriate for a Taylor-Couette experiment.

## 8 Acknowledgments

S.E.C. was supported by a National Science Foundation Graduate Research Fellowship under grant No. DGE-16-44869. J.S.O. acknowledges support from NASA grant NNX16AC92G. We thank the anonymous referee for thoughtful comments that greatly improved the manuscript. We also thank Mordecai-Mark Mac Low, Jeremy Goodman, John Krommes, Geoff Vasil, and Ellen Zweibel for useful discussion. \softwareDedalus: http://dedalus-project.org/

## Appendix A Detailed Equations

Here we detail the perturbation analysis described in Section 3. The perturbation series is described by Equations 18 - 21, where

\mathcal{L}=\mathcal{L}_{0}+\mathcal{L}_{1}\partial_{z}+\mathcal{L}_{2}% \partial_{z}^{2}+\mathcal{L}_{3}\partial_{z}^{3}+\mathcal{L}_{4}\partial_{z}^{% 4}, | (A1) |

\widetilde{\mathcal{L}}_{1}=\mathcal{L}_{1}+2\mathcal{L}_{2}\partial_{z}+3% \mathcal{L}_{3}\partial_{z}^{2}+4\mathcal{L}_{4}\partial_{z}^{3} | (A2) |

\widetilde{\mathcal{L}}_{2}=\mathcal{L}_{2}+3\mathcal{L}_{3}\partial_{z}+6% \mathcal{L}_{4}\partial_{z}^{2} | (A3) |

\widetilde{\mathcal{G}}=\mathcal{G}\partial_{z}+\mathcal{L}_{3}\partial_{z}^{3}, | (A4) |

and the constituent matrices are defined as

\mathcal{D}=\left[\begin{matrix}\nabla^{2}&0&0&0\\ 0&1&0&0\\ 0&0&1&0\\ 0&0&0&1\\ \end{matrix}\right] | (A5) |

\mathcal{L}_{0}=\left[\begin{matrix}-\frac{1}{\mathrm{Re}}\partial_{x}^{4}&0&0% &0\\ 0&-\frac{1}{\mathrm{Re}}\partial_{x}^{2}&0&0\\ 0&0&-\frac{1}{\mathrm{Rm}}\partial_{x}^{2}&0\\ 0&0&0&-\frac{1}{\mathrm{Rm}}\partial_{x}^{2}\\ \end{matrix}\right] | (A6) |

\mathcal{L}_{1}=\left[\begin{matrix}0&-2&-\mathrm{Co}\partial_{x}^{2}&0\\ (2-q)\Omega_{0}&0&0&-\mathrm{Co}\\ -1&0&0&0\\ 0&-1&q\Omega_{0}&0\\ \end{matrix}\right] | (A7) |

\mathcal{L}_{2}=\left[\begin{matrix}-2\frac{1}{\mathrm{Re}}\partial_{x}^{2}&0&% 0&0\\ 0&-\frac{1}{\mathrm{Re}}&0&0\\ 0&0&-\frac{1}{\mathrm{Rm}}&0\\ 0&0&0&-\frac{1}{\mathrm{Rm}}\\ \end{matrix}\right] | (A8) |

\mathcal{L}_{3}=\left[\begin{matrix}0&0&-\mathrm{Co}&0\\ 0&0&0&0\\ 0&0&0&0\\ 0&0&0&0\\ \end{matrix}\right] | (A9) |

\mathcal{L}_{4}=\left[\begin{matrix}-\frac{1}{\mathrm{Re}}&0&0&0\\ 0&0&0&0\\ 0&0&0&0\\ 0&0&0&0\\ \end{matrix}\right] | (A10) |

\mathcal{G}=\left[\begin{matrix}0&0&-\mathrm{Co}\partial_{x}^{2}&0\\ 0&0&0&-\mathrm{Co}\\ -1&0&0&0\\ 0&-1&0&0\\ \end{matrix}\right] | (A11) |

Once perturbed, the system is solved for successive orders of \epsilon (Equations 18 - 21). \mathcal{O}(\epsilon) is the linear system. At \mathcal{O}(\epsilon^{2}), first-order MRI modes nonlinearly interact with themselves and with their complex conjugates, and so the term \mathbf{N}_{2} in Equation 19 has the form

\mathbf{N}_{2}=|\alpha|^{2}\mathbf{N}_{20}+\alpha^{2}\mathbf{N}_{22}e^{2ik_{c}z} | (A12) |

(see Appendix B for the full form of \mathbf{N}_{20} and \mathbf{N}_{22}).

Note that, following the notation of Umurhan et al. (2007b), the subscripts refer to \epsilon order, z order, successively, such that \mathbf{N}_{22} is the second-order nonlinear term which corresponds to e^{2ik_{c}z} z-dependence.

Equation 19 is solved as three separate systems of equations, one for each possible z resonance:

\displaystyle\mathcal{L}\mathbf{V}_{20} | \displaystyle=\mathbf{N}_{20} | (A13) | ||

\displaystyle\mathcal{L}\mathbf{V}_{21} | \displaystyle=-\widetilde{\mathcal{L}}_{1}\partial_{Z}\mathbf{V}_{11} | (A14) | ||

\displaystyle\mathcal{L}\mathbf{V}_{22} | \displaystyle=\mathbf{N}_{22} | (A15) |

Finally, at \mathcal{O}(\epsilon^{3}) we eliminate secular terms to close the system. Secular terms are terms which are resonant with the solution to the homogenous linear equation (Equation 18), and which cause the higher-order solutions to grow without bound. The solvability criterion we enforce to eliminate these terms is the vanishing of the inner product of the solution to the adjoint linear homogenous equation \mathcal{L}^{\dagger}\mathbf{V}^{\dagger}=0 with the nonhomogenous terms in Equation 21, namely

\langle\mathbb{V}^{\dagger}|\mathcal{D}\mathbb{V}_{11}\rangle\partial_{T}% \alpha+\langle\mathbb{V}^{\dagger}|\widetilde{\mathcal{G}}\mathbb{V}_{11}% \rangle\alpha+\langle\mathbb{V}^{\dagger}|\widetilde{\mathcal{L}}_{1}\mathbb{V% }_{21}+\widetilde{\mathcal{L}}_{2}\mathbb{V}_{11}\rangle\partial_{Z}^{2}\alpha% =\langle\mathbb{V}^{\dagger}|\mathbf{N}_{31}\rangle\alpha|\alpha|^{2}. | (A16) |

This solvability criterion derives from a corollary to the Fredholm Alternative (see Paper II for a formal definition).

Equation A16 can be rewritten as Equation 22, the Ginzburg-Landau equation, where the coefficients are

b=\langle\mathbb{V}^{\dagger}|\widetilde{\mathcal{G}}\mathbb{V}_{11}\rangle/% \langle\mathbb{V}^{\dagger}|\mathcal{D}\mathbb{V}_{11}\rangle, | (A17) |

h=\langle\mathbb{V}^{\dagger}|\widetilde{\mathcal{L}}_{1}\mathbb{V}_{21}+% \widetilde{\mathcal{L}}_{2}\mathbb{V}_{11}\rangle/\langle\mathbb{V}^{\dagger}|% \mathcal{D}\mathbb{V}_{11}\rangle, | (A18) |

and

c=\langle\mathbb{V}^{\dagger}|\mathbf{N}_{31}\rangle/\langle\mathbb{V}^{% \dagger}|\mathcal{D}\mathbb{V}_{11}\rangle. | (A19) |

We define the adjoint operator \mathcal{L}^{\dagger} and solution \mathbf{V}^{\dagger} as

\langle\mathbf{V^{\dagger}}|\mathcal{L}\mathbf{V}\rangle=\langle\mathcal{L}^{% \dagger}\mathbf{V}^{\dagger}|\mathbf{V}\rangle, | (A20) |

where the inner product is defined as

\langle\mathbf{V^{\dagger}}|\mathcal{L}\mathbf{V}\rangle=\frac{k_{c}}{2\pi}% \int_{-\pi/k_{c}}^{\pi/k_{c}}\int_{x_{1}}^{x_{2}}\mathbf{V}^{\dagger*}\cdot% \mathcal{L}\mathbf{V}\,\mathrm{d}x\mathrm{d}z. | (A21) |

The solution to the adjoint homogenous equation has the form

\mathbf{V^{\dagger}}=\mathbb{V}^{\dagger}(x)e^{ik_{c}z}+c.c. | (A22) |

As noted by URM07, a second amplitude equation for a spatially constant azimuthal magnetic field mode arises from the terms in the \mathcal{O}(\epsilon^{3}) equation which contain no z dependence. This is a diffusion equation, so the neutral mode simply decays away.

## Appendix B Expansion of Nonlinear Terms

At each order in our perturbation series, lower-order MRI modes nonlinearly interact. Thus there is a nonlinear term contribution at \mathcal{O}(\epsilon^{2}) and \mathcal{O}(\epsilon^{3}). Here we detail the form of these nonlinear terms.

The overall nonlinear contribution to our system, written as a vector \mathbf{N} in Equation 13, is

\mathbf{N}=\epsilon^{2}\mathbf{N_{2}}\,+\,\epsilon^{3}\mathbf{N_{3}}\,+\,% \mathcal{O}(\epsilon^{4}) | (B1) |

where

N_{2}^{(\Psi)}=J(\Psi_{1},\nabla^{2}\Psi_{1})\,-\,\mathrm{Co}J(A_{1},\nabla^{2% }A_{1}) | (B2) |

N_{2}^{(u)}=J(\Psi_{1},u_{1})\,-\,\mathrm{Co}J(A_{1},B_{1}) | (B3) |

N_{2}^{(A)}=-J(A_{1},\Psi_{1}) | (B4) |

N_{2}^{(B)}=J(\Psi_{1},B_{1})\,-\,J(A_{1},u_{1}) | (B5) |

and

\begin{split}\displaystyle N_{3}^{(\Psi)}&\displaystyle=J(\Psi_{1},\nabla^{2}% \Psi_{2})\,-\,\mathrm{Co}J(A_{1},\nabla^{2}A_{2})\,+\,J(\Psi_{2},\nabla^{2}% \Psi_{1})-\,\mathrm{Co}J(A_{2},\nabla^{2}A_{1})\,+\,2J(\Psi_{1},\partial_{z}% \partial_{Z}\Psi_{1})\\ &\displaystyle-\,2\mathrm{Co}J(A_{1},\partial_{z}\partial_{Z}A_{1})\,+\,% \widetilde{J}(\Psi_{1},\nabla^{2}\Psi_{1})\,-\,\mathrm{Co}\widetilde{J}(A_{1},% \nabla^{2}A_{1})\\ \end{split} | (B6) |

N_{3}^{(u)}=J(\Psi_{1},u_{2})\,+\,J(\Psi_{2},u_{1})\,+\,\widetilde{J}(\Psi_{1}% ,u_{1})\,-\,\mathrm{Co}J(A_{1},B_{2})\,-\,\mathrm{Co}J(A_{2},B_{1})\,-\,% \mathrm{Co}\widetilde{J}(A_{1},B_{1}) | (B7) |

N_{3}^{(A)}=-J(A_{1},\Psi_{2})\,-\,J(A_{2},\Psi_{1})\,-\,\widetilde{J}(A_{1},% \Psi_{1}) | (B8) |

N_{3}^{(B)}=\,J(\Psi_{1},B_{2})\,+\,J(\Psi_{2},B_{1})\,+\,\widetilde{J}(\Psi_{% 1},B_{1})\,-\,J(A_{1},u_{2})\,-\,J(A_{2},u_{1})\,-\,\widetilde{J}(A_{1},u_{1})% .\\ | (B9) |

\mathbf{N_{2}} and \mathbf{N_{3}} expand to become

\mathbf{N_{2}}=\alpha^{2}\mathbb{N}_{22}e^{i2k_{c}z}+\left|\alpha\right|^{2}% \mathbb{N}_{20}+c.c. | (B10) |

and

\mathbf{N_{3}}=\alpha^{3}\mathbb{N}_{33}e^{i3k_{c}z}+\alpha\partial_{Z}\alpha% \mathbb{N}_{32}e^{i2k_{c}z}+\alpha\left|\alpha\right|^{2}\mathbb{N}_{31}e^{ik_% {c}z}+\alpha\partial_{Z}\beta\mathbb{\widetilde{N}}_{31}e^{ik_{c}z}+\alpha^{*}% \partial_{Z}\alpha\mathbb{N}_{30}+c.c. | (B11) |

The second order nonlinear terms are

\begin{split}\displaystyle N_{22}^{(\Psi)}=&\displaystyle\,ik_{c}\Psi_{11}% \cdot\left(\partial_{x}^{3}\Psi_{11}-k_{c}^{2}\partial_{x}\Psi_{11}\right)-% \partial_{x}\Psi_{11}\cdot\left(ik_{c}\partial_{x}^{2}\Psi_{11}-ik_{c}^{3}\Psi% _{11}\right)\\ &\displaystyle+\mathrm{Co}\partial_{x}A_{11}\cdot\left(ik_{c}\partial_{x}^{2}A% _{11}-ik_{c}^{3}A_{11}\right)-\mathrm{Co}ik_{c}A_{11}\cdot\left(\partial_{x}^{% 3}A_{11}-k_{c}^{2}\partial_{x}A_{11}\right)\end{split} | (B12) |

\begin{split}\displaystyle N_{22}^{(u)}=&\displaystyle\,ik_{c}\Psi_{11}\cdot% \partial_{x}u_{11}-\partial_{x}\Psi_{11}\cdot ik_{c}u_{11}-\mathrm{Co}ik_{c}A_% {11}\cdot\partial_{x}B_{11}+\mathrm{Co}\partial_{x}A_{11}\cdot ik_{c}B_{11}% \end{split} | (B13) |

N_{22}^{(A)}=-ik_{c}A_{11}\cdot\partial_{x}\Psi_{11}+\partial_{x}A_{11}\cdot ik% _{c}\Psi_{11} | (B14) |

N_{22}^{(B)}=ik_{c}\Psi_{11}\cdot\partial_{x}B_{11}-\partial_{x}\Psi_{11}\cdot ik% _{c}B_{11}-ik_{c}A_{11}\cdot\partial_{x}u_{11}+\partial_{x}A_{11}\cdot ik_{c}u% _{11} | (B15) |

\begin{split}\displaystyle N_{20}^{(\Psi)}=&\displaystyle\,ik_{c}\Psi_{11}% \cdot\left(\partial_{x}^{3}\Psi_{11}^{*}-k_{c}^{2}\partial_{x}\Psi_{11}^{*}% \right)-\partial_{x}\Psi_{11}\cdot\left(ik_{c}^{3}\Psi_{11}^{*}-ik_{c}\partial% _{x}^{2}\Psi_{11}^{*}\right)\\ &\displaystyle+\mathrm{Co}\partial_{x}A_{11}\cdot\left(ik_{c}^{3}A_{11}^{*}-ik% _{c}\partial_{x}^{2}A_{11}^{*}\right)-\mathrm{Co}ik_{c}A_{11}\cdot\left(% \partial_{x}^{3}A_{11}^{*}-k_{c}^{2}\partial_{x}A_{11}^{*}\right)\end{split} | (B16) |

N_{20}^{(u)}=ik_{c}\Psi_{11}\cdot\partial_{x}u_{11}^{*}+\partial_{x}\Psi_{11}% \cdot ik_{c}u_{11}^{*}-\mathrm{Co}ik_{c}A_{11}\cdot\partial_{x}B_{11}^{*}-% \mathrm{Co}\partial_{x}A_{11}\cdot ik_{c}B_{11}^{*} | (B17) |

N_{20}^{(A)}=-ik_{c}A_{11}\cdot\partial_{x}\Psi_{11}^{*}-\partial_{x}A_{11}% \cdot ik_{c}\Psi_{11}^{*} | (B18) |

N_{20}^{(B)}=ik_{c}\Psi_{11}\cdot\partial_{x}B_{11}^{*}+\partial_{x}\Psi_{11}% \cdot ik_{c}B_{11}^{*}-ik_{c}A_{11}\cdot\partial_{x}u_{11}^{*}-\partial_{x}A_{% 11}\cdot ik_{c}u_{11}^{*} | (B19) |

and the third order nonlinear terms become

\begin{split}\displaystyle N_{31}^{(\Psi)}=&\displaystyle\,ik_{c}\left(\Psi_{1% 1}\cdot\partial_{x}^{3}\Psi_{20}\right)+ik_{c}\left(\Psi_{11}\cdot\partial_{x}% ^{3}\Psi_{20}^{*}\right)-ik_{c}\left(\Psi_{11}^{*}\cdot\partial_{x}^{3}\Psi_{2% 2}\right)-i2k_{c}\left(\partial_{x}\Psi_{11}^{*}\cdot\partial_{x}^{2}\Psi_{22}% \right)\\ &\displaystyle+i8k_{c}^{3}\left(\partial_{x}\Psi_{11}^{*}\cdot\Psi_{22}\right)% +i4k_{c}^{3}\left(\Psi_{11}^{*}\cdot\partial_{x}\Psi_{22}\right)+\mathrm{Co}% \left[-ik_{c}\left(A_{11}\cdot\partial^{3}A_{20}\right)-ik_{c}\left(A_{11}% \cdot\partial_{x}^{3}A_{20}^{*}\right)\right]\\ &\displaystyle+\mathrm{Co}\left[ik_{c}\left(A_{11}^{*}\cdot\partial_{x}^{3}A_{% 22}\right)+i2k_{c}\left(\partial_{x}A_{11}^{*}\cdot\partial_{x}^{2}A_{22}% \right)-i8k_{c}^{3}\left(\partial_{x}A_{11}^{*}\cdot A_{22}\right)-i4k_{c}^{3}% \left(A_{11}^{*}\cdot\partial_{x}A_{22}\right)\right]\\ &\displaystyle+i2k_{c}\left(\Psi_{22}\cdot\partial_{x}^{3}\Psi_{11}^{*}\right)% -i2k_{c}^{3}\left(\Psi_{22}\cdot\partial_{x}\Psi_{11}^{*}\right)-ik_{c}\left(% \partial_{x}\Psi_{20}\cdot\partial_{x}^{2}\Psi_{11}\right)+ik_{c}\left(% \partial_{x}\Psi_{22}\cdot\partial_{x}^{2}\Psi_{11}^{*}\right)\\ &\displaystyle-ik_{c}\left(\partial_{x}\Psi_{20}^{*}\cdot\partial_{x}^{2}\Psi_% {11}\right)+ik_{c}^{3}\left(\partial_{x}\Psi_{20}\cdot\Psi_{11}\right)+ik_{c}^% {3}\left(\partial_{x}\Psi_{20}^{*}\cdot\Psi_{11}\right)-ik_{c}^{3}\left(% \partial_{x}\Psi_{22}\cdot\Psi_{11}^{*}\right)\\ &\displaystyle+\mathrm{Co}\left[-i2k_{c}\left(A_{22}\cdot\partial_{x}^{3}A_{11% }^{*}\right)+i2k_{c}^{3}\left(A_{22}\cdot\partial_{x}A_{11}^{*}\right)+ik_{c}% \left(\partial_{x}A_{20}\cdot\partial_{x}^{2}A_{11}\right)-ik_{c}\left(% \partial_{x}A_{22}\cdot\partial_{x}^{2}A_{11}^{*}\right)\right]\\ &\displaystyle+\mathrm{Co}\left[ik_{c}\left(\partial_{x}A_{20}^{*}\cdot% \partial_{x}^{2}A_{11}\right)-ik_{c}^{3}\left(\partial_{x}A_{20}\cdot A_{11}% \right)-ik_{c}^{3}\left(\partial_{x}A_{20}^{*}\cdot A_{11}\right)+ik_{c}^{3}% \left(\partial_{x}A_{22}\cdot A_{11}^{*}\right)\right]\end{split} | (B20) |

\begin{split}\displaystyle N_{31}^{(u)}=&\displaystyle ik_{c}\left(\Psi_{11}% \cdot\partial_{x}u_{20}\right)+ik_{c}\left(\Psi_{11}\cdot\partial_{x}u_{20}^{*% }\right)-ik_{c}\left(\Psi_{11}^{*}\cdot\partial_{x}u_{22}\right)-i2k_{c}\left(% \partial_{x}\Psi_{11}^{*}\cdot u_{22}\right)\\ &\displaystyle-ik_{c}\left(u_{11}\cdot\partial_{x}\Psi_{20}\right)-ik_{c}\left% (u_{11}\cdot\partial_{x}\Psi_{20}^{*}\right)+ik_{c}\left(u_{11}^{*}\cdot% \partial_{x}\Psi_{22}\right)+i2k_{c}\left(\partial_{x}u_{11}^{*}\cdot\Psi_{22}% \right)\\ &\displaystyle+\mathrm{Co}\left[-ik_{c}\left(A_{11}\cdot\partial_{x}B_{20}% \right)-ik_{c}\left(A_{11}\cdot\partial_{x}B_{20}^{*}\right)+ik_{c}\left(A_{11% }^{*}\cdot\partial_{x}B_{22}\right)+i2k_{c}\left(\partial_{x}A_{11}^{*}\cdot B% _{22}\right)\right]\\ &\displaystyle+\mathrm{Co}\left[ik_{c}\left(B_{11}\cdot\partial_{x}A_{20}% \right)+ik_{c}\left(B_{11}\cdot\partial_{x}A_{20}^{*}\right)-ik_{c}\left(B_{11% }^{*}\cdot\partial_{x}A_{20}\right)-i2k_{c}\left(\partial_{x}B_{11}^{*}\cdot A% _{22}\right)\right]\end{split} | (B21) |

\begin{split}\displaystyle N_{31}^{(A)}=&\displaystyle-ik_{c}\left(A_{11}\cdot% \partial_{x}\Psi_{20}\right)-ik_{c}\left(A_{11}\cdot\partial_{x}\Psi_{20}^{*}% \right)+ik_{c}\left(A_{11}^{*}\cdot\partial_{x}\Psi_{22}\right)+i2k_{c}\left(% \partial_{x}A_{11}^{*}\cdot\Psi_{22}\right)\\ &\displaystyle+ik_{c}\left(\Psi_{11}\cdot\partial_{x}A_{20}\right)+ik_{c}\left% (\Psi_{11}\cdot\partial_{x}A_{20}^{*}\right)-ik_{c}\left(\Psi_{11}^{*}\cdot% \partial_{x}A_{22}\right)-i2k_{c}\left(\partial_{x}\Psi_{11}^{*}\cdot A_{22}% \right)\end{split} | (B22) |

\begin{split}\displaystyle N_{31}^{(B)}=&\displaystyle ik_{c}\left(\Psi_{11}% \cdot\partial_{x}B_{20}\right)+ik_{c}\left(\Psi_{11}\cdot\partial_{x}B_{20}^{*% }\right)-ik_{c}\left(\Psi_{11}^{*}\cdot\partial_{x}B_{22}\right)-i2k_{c}\left(% \partial_{x}\Psi_{11}^{*}\cdot B_{22}\right)\\ &\displaystyle-ik_{c}\left(B_{11}\cdot\partial_{x}\Psi_{20}\right)-ik_{c}\left% (B_{11}\cdot\partial_{x}\Psi_{20}^{*}\right)+ik_{c}\left(B_{11}^{*}\cdot% \partial_{x}\Psi_{22}\right)+i2k_{c}\left(\partial_{x}B_{11}^{*}\cdot\Psi_{22}% \right)\\ &\displaystyle-ik_{c}\left(A_{11}\cdot\partial_{x}u_{20}\right)-ik_{c}\left(A_% {11}\cdot\partial_{x}u_{20}^{*}\right)+ik_{c}\left(A_{11}^{*}\cdot\partial_{x}% u_{22}\right)+i2k_{c}\left(\partial_{x}A_{11}^{*}\cdot u_{22}\right)\\ &\displaystyle ik_{c}\left(u_{11}\cdot\partial_{x}A_{20}\right)+ik_{c}\left(u_% {11}\cdot\partial_{x}A_{20}^{*}\right)-ik_{c}\left(u_{11}^{*}\cdot\partial_{x}% A_{22}\right)-i2k_{c}\left(\partial_{x}u_{11}^{*}\cdot A_{22}\right)\end{split} | (B23) |

## Appendix C Linear dispersion relation

The linear dispersion relation, which determines the variable scalings in the multiple scales analysis. This relation is found by perturbing the linear system (Equation 18) with a small perturbation of the form e^{\sigma t+ik_{x}x+ik_{z}z}. Note that the spatial eigenvalues appear as k_{z}^{2} and k_{x}^{2} at lowest order.

\begin{split}&\displaystyle\frac{B_{0}^{4}k_{x}^{2}k_{z}^{4}}{16\pi^{2}}+\frac% {B_{0}^{4}k_{z}^{6}}{16\pi^{2}}-\frac{B_{0}^{2}\Omega_{0}k_{z}^{4}q}{2\pi}-2% \Omega_{0}k_{z}^{2}q\sigma^{2}-\frac{4\sigma}{\mathrm{Rm}}\Omega_{0}k_{x}^{2}k% _{z}^{2}q-\frac{4\sigma}{\mathrm{Rm}}\Omega_{0}k_{z}^{4}q-\frac{2\Omega_{0}}{% \mathrm{Rm}^{2}}k_{x}^{4}k_{z}^{2}q-\frac{4\Omega_{0}}{\mathrm{Rm}^{2}}k_{x}^{% 2}k_{z}^{4}q-\frac{2\Omega_{0}}{\mathrm{Rm}^{2}}k_{z}^{6}q\\ &\displaystyle-k_{x}^{2}\sigma^{4}-k_{z}^{2}\sigma^{4}+4k_{z}^{2}\sigma^{2}-% \frac{2\sigma^{3}}{\mathrm{Rm}}k_{x}^{4}-\frac{4\sigma^{3}}{\mathrm{Rm}}k_{x}^% {2}k_{z}^{2}+\frac{8\sigma}{\mathrm{Rm}}k_{x}^{2}k_{z}^{2}-\frac{2\sigma^{3}}{% \mathrm{Rm}}k_{z}^{4}+\frac{8\sigma}{\mathrm{Rm}}k_{z}^{4}-\frac{k_{x}^{6}% \sigma^{2}}{\mathrm{Rm}^{2}}-\frac{3\sigma^{2}}{\mathrm{Rm}^{2}}k_{x}^{4}k_{z}% ^{2}+\frac{4k_{x}^{4}}{\mathrm{Rm}^{2}}k_{z}^{2}\\ &\displaystyle-\frac{3\sigma^{2}}{\mathrm{Rm}^{2}}k_{x}^{2}k_{z}^{4}+\frac{8k_% {x}^{2}}{\mathrm{Rm}^{2}}k_{z}^{4}-\frac{k_{z}^{6}\sigma^{2}}{\mathrm{Rm}^{2}}% +\frac{4k_{z}^{6}}{\mathrm{Rm}^{2}}-\frac{2\sigma^{3}}{\mathrm{Re}}k_{x}^{4}-% \frac{4\sigma^{3}}{\mathrm{Re}}k_{x}^{2}k_{z}^{2}-\frac{2\sigma^{3}}{\mathrm{% Re}}k_{z}^{4}-\frac{4k_{x}^{6}\sigma^{2}}{\mathrm{Re}\mathrm{Rm}}-\frac{12k_{x% }^{4}k_{z}^{2}\sigma^{2}}{\mathrm{Re}\mathrm{Rm}}-\frac{12k_{x}^{2}k_{z}^{4}% \sigma^{2}}{\mathrm{Re}\mathrm{Rm}}\\ &\displaystyle-\frac{4k_{z}^{6}\sigma^{2}}{\mathrm{Re}\mathrm{Rm}}-\frac{2k_{x% }^{8}\sigma}{\mathrm{Re}\mathrm{Rm}^{2}}-\frac{8k_{x}^{6}k_{z}^{2}\sigma}{% \mathrm{Re}\mathrm{Rm}^{2}}-\frac{12k_{x}^{4}k_{z}^{4}\sigma}{\mathrm{Re}% \mathrm{Rm}^{2}}-\frac{8k_{x}^{2}k_{z}^{6}\sigma}{\mathrm{Re}\mathrm{Rm}^{2}}-% \frac{2k_{z}^{8}\sigma}{\mathrm{Re}\mathrm{Rm}^{2}}-\frac{k_{x}^{6}\sigma^{2}}% {\mathrm{Re}^{2}}-\frac{3\sigma^{2}}{\mathrm{Re}^{2}}k_{x}^{4}k_{z}^{2}-\frac{% 3\sigma^{2}}{\mathrm{Re}^{2}}k_{x}^{2}k_{z}^{4}-\frac{k_{z}^{6}\sigma^{2}}{% \mathrm{Re}^{2}}\\ &\displaystyle-\frac{2k_{x}^{8}\sigma}{\mathrm{Re}^{2}\mathrm{Rm}}-\frac{8k_{x% }^{6}k_{z}^{2}\sigma}{\mathrm{Re}^{2}\mathrm{Rm}}-\frac{12k_{x}^{4}k_{z}^{4}% \sigma}{\mathrm{Re}^{2}\mathrm{Rm}}-\frac{8k_{x}^{2}k_{z}^{6}\sigma}{\mathrm{% Re}^{2}\mathrm{Rm}}-\frac{2k_{z}^{8}\sigma}{\mathrm{Re}^{2}\mathrm{Rm}}-\frac{% k_{x}^{10}}{\mathrm{Re}^{2}\mathrm{Rm}^{2}}-\frac{5k_{x}^{8}k_{z}^{2}}{\mathrm% {Re}^{2}\mathrm{Rm}^{2}}-\frac{10k_{x}^{6}k_{z}^{4}}{\mathrm{Re}^{2}\mathrm{Rm% }^{2}}-\frac{10k_{x}^{4}k_{z}^{6}}{\mathrm{Re}^{2}\mathrm{Rm}^{2}}\\ &\displaystyle-\frac{5k_{x}^{2}k_{z}^{8}}{\mathrm{Re}^{2}\mathrm{Rm}^{2}}-% \frac{k_{z}^{10}}{\mathrm{Re}^{2}\mathrm{Rm}^{2}}=0\end{split} | (C1) |

## References

- Bai (2015) Bai, X.-N. 2015, ApJ, 798, 84
- Balbus & Hawley (1991) Balbus, S A and Hawley, J F, 1991, ApJ, 376, 214
- Bender & Orszag (1978) Bender, C. M., & Orszag, S. A. 1978, Advanced Mathematical Methods for Scientists and Engineers, New York: McGraw-Hill, 1978,
- Boyd (2001) Boyd, J P, 2001, Chebyshev and Fourier Spectral Methods, New York, Dover
- Chandrasekhar (1960) Chandrasekhar, S. 1960, Proceedings of the National Academy of Sciences of the United States of America, 46, 253
- Clark & Oishi (2017b) Clark, S.E. and Oishi, J.S., 2017, ApJ, accepted
- Deyirmenjian (1997) Deyirmenjian, Z., Daya, A., & Morris, S.W., 1997, Phys. Rev. E, 56, 1706
- Ebrahimi et al. (2009) Ebrahimi, F., Prager, S.C., Schnack, D.D, 2009, ApJ, 698, 233
- Goodman & Xu (1994) Goodman, J, Xu, G, 1994, ApJ, 432, 213
- Hawley et al. (2011) Hawley, J. F., Guan, X., & Krolik, J. H. 2011, ApJ, 738, 84
- Hoyle (2006) Hoyle, R. 2006, Pattern Formation, by Rebecca Hoyle, pp. 432. Cambridge University Press, March 2006. ISBN-10: 0521817501. ISBN-13: 9780521817509, 432
- Knobloch & Julien (2005) Knobloch, E, Julien K, 2005, Physics of Fluids, 17, 094106
- Kunz & Balbus (2004) Kunz, M. W. & Balbus, S. A. 2004, MNRAS, 348, 355
- Kunz & Lesur (2013) Kunz, M. W. & Lesur, G. 2013, MNRAS, 434, 2295
- Latter et al. (2015) Latter, H.N., Fromang, S., & Faure, J., 2015, MNRAS, 453, 3257
- Latter et al. (2010) Latter, H.N., Fromang, S., & Gressel, O., 2010, MNRAS, 406, 848
- Lesur & Longaretti (2007) Lesur, G., & Longaretti, P.-Y. 2007, MNRAS, 378, 1471
- Makwana et al. (2011) Makwana, K.D., Terry, P.W., Kim, J.H., & Hatch, D.R. 2011, Phys. of Plasmas, 18, 1
- Newell & Whitehead (1969) Newell, A.C. & Whitehead, J.A. 1969, JFM, 38, 279N
- Parkin (2014) Parkin, E. R. 2014, MNRAS, 438, 2513
- Parkin & Bicknell (2013) Parkin, E. R., & Bicknell, G. V. 2013, MNRAS, 435, 2281
- Pessah (2010) Pessah, M. E., 2010, ApJ, 716, 1012
- Pessah & Goodman (2009) Pessah, M. E., & Goodman, J., 2009, ApJ698, L72
- Recktenwald et al. (1993) Recktenwald, A., Lücke, M., Müller, H.W., 1993, Phys. Rev. E48, 4444
- Regev & Umurhan (2008) Regev, O., & Umurhan, O. M. 2008, A&A, 481, 21
- Rembiasz et al. (2016) Rembiasz, T., Obergaulinger, M., Cerdá-Durán, P., Müller, E., & Aloy, M.A., 2016, MNRAS456, 3782
- Schnittman et al. (2013) Schnittman, J. D., Krolik, J. H., & Noble, S. C. 2013, ApJ, 769, 156
- Shakura & Sunyaev (1973) Shakura, N. I. and Sunyaev, R. A., A&A, 1973, 24, 337
- Umurhan et al. (2007a) Umurhan, O. M., Regev, O., Menou, K., 2007, Phys. Rev. Lett., 98, 034501
- Umurhan et al. (2007b) Umurhan, O. M., Regev, O., Menou, K., 2007, Phys. Rev. E, 76, 036310
- Vasil (2015) Vasil, G. M. 2015, RSPSA, 471, 20140699
- Velikhov (1959) Velikhov, E. P., 1959, Sov. Phys. JETP, 36, 995.
- Wheeler et al. (2015) Wheeler, J. C., Kagan, D., & Chatzopoulos, E. 2015, ApJ, 799, 85