## Abstract

The magnetorotational instability is believed to play an important role in accretion disc physics in extracting angular momentum from the disc and allowing accretion to take place. For this reason the instability has been the subject of numerous numerical simulations and, increasingly, laboratory experiments. In this review, recent developments in both areas are surveyed, and a new theoretical approach to understanding the nonlinear processes involved in the saturation of the instability is outlined.

## 1. Introduction

Accretion is a process of fundamental importance in astrophysics. However, accretion can only occur in the presence of an efficient mechanism for angular momentum extraction. Many processes have been suggested that might lead to efficient accretion. These include barotropic (Arlt & Urpin 2004; Dubrulle *et al.* 2005*b*) and baroclinic instabilities (Knobloch & Spruit 1986), vortices (Barranco & Marcus 2005; Petersen *et al.* 2007*a*,*b*; Lithwick 2009), convection (Stone & Balbus 1996), spiral shock waves (Kaisig 1989; Roźycka & Spruit 1993) and finite-amplitude shear instabilities (Lerner & Knobloch 1988; Longaretti 2002; Dubrulle *et al.* 2005*a*; Lesur & Longaretti 2005), as well as different types of global instabilities associated with over-reflection. Significant efforts to identify finite-amplitude instabilities in a Couette–Taylor experiment have also been undertaken (Wendt 1933; Taylor 1936; Richard & Zahn 1999; Dubrulle *et al.* 2005*a*), but the general consensus is that none of these mechanisms leads to substantial outward angular momentum transport (Ji *et al.* 2006).

In this review we focus on another possibility, namely the magnetorotational instability (MRI). This instability, also originally identified in the Couette–Taylor geometry (Velikhov 1959; Chandrasekhar 1960), occurs in the Rayleigh-stable regime whenever a weak poloidal magnetic field is present, provided the angular velocity *Ω** in the disc decreases outwards (*dΩ**/*dr**<0).^{1} This is the case in thin Keplerian discs. Particularly noteworthy are the facts that the instability is linear, occurs on a dynamic time scale, *τ**_{MRI}∼(*Ω**)^{−1}, and is fundamentally axisymmetric. The latter is helpful for both theory and numerical simulations. When the imposed magnetic field is weak, the instability has a small wavelength in the direction parallel to the rotation and magnetic field, and takes the form of thin sheets of matter moving alternately radially inwards and outwards and redistributing angular momentum (Balbus & Hawley 1991, 1998).

## 2. Ideal magnetorotational instability

The MRI is easiest to understand in the context of ideal incompressible flow in an imposed magnetic field. The basic magnetohydrodynamic (MHD) equations are
These equations have a basic axisymmetric solution of the form
in (*r*,*ϕ*,*z*) coordinates. In a thin accretion disc the azimuthal flow *V* (*r*) is determined by the gravitational field of the central object; in the present formulation, this field is subsumed in the generalized pressure gradient ∇*p* and the flow *V* (*r*) is taken to be an arbitrary function of the radial coordinate *r*.

We write infinitesimal axisymmetric perturbations of this state in the form (*u*,*v*,*w*), (*a*,*b*,*c*), where , etc., *n* is the vertical wavenumber and *ω* is a (complex) frequency. The frequency *ω* solves the eigenvalue problem (Acheson 1973)
2.1
where
are the (squares of the) Alfvén speeds associated with the imposed azimuthal and vertical magnetic fields. This eigenvalue problem is to be solved subject to imposed boundary conditions in the radial direction. In general, the associated eigenvalue depends on the choice of these boundary conditions and there are then two types of modes, body modes for which the eigenfunction extends across the disc, and wall modes that are localized near the boundaries. The former are, in general, insensitive to the boundary conditions adopted, whereas the latter are not. The MRI is a body mode.

### (a) Standard magnetorotational instability

In the simplest case we may take *u*(*a*)=*u*(*b*)=0 and suppose that *V*_{z}=*const*.≠0, *V*_{ϕ}=0. Multiplying equation (2.1) by *u**(*r*), the complex conjugate of *u*(*r*), and integrating by parts we obtain the *identity*
where
From this identity it follows that *ω*^{2} must be real (Chandrasekhar 1960); moreover, we see that if the disc is hydrodynamically stable, i.e. *d*(*r*^{2}*V* ^{2})/*dr*>0, a sufficient condition for stability is that *d*(*V* ^{2}/*r*^{2})/*dr*>0 in *a*<*r*<*b*, while for instability it is necessary that *d*(*V* ^{2}/*r*^{2})/*dr*<0 somewhere in *a*<*r*<*b*. This is the classical MRI (Velikhov 1959; Chandrasekhar 1960; see also Acheson 1973; Acheson & Hide 1973; Balbus & Hawley 1991; Knobloch 1992). For thin accretion discs and hence both the Rayleigh condition *d*(*r*^{2}*V* ^{2})/*dr*>0 for hydrodynamic stability and the condition *d*(*V* ^{2}/*r*^{2})/*dr*<0 for the MRI instability in the presence of a vertical magnetic field are satisfied.

Simulations of this instability in the nonlinear regime in the presence of dissipation have been carried out by several groups (Kersalé *et al.* 2006; Liu *et al.* 2006*b*). In its pure form, the standard MRI remains to be detected in a laboratory setting although an MRI-like instability of a turbulent flow in a spherical Couette–Taylor system has been observed (Sisan *et al.* 2004).

### (b) Helical magnetorotational instability

The case *V* _{z}=*const*.≠0, *V* _{ϕ}≠0 is also of interest since differential rotation always generates an azimuthal magnetic field from an imposed vertical field. In this case, the same procedure leads to the identity
2.2
Since *V* _{ϕ}*V* _{z}≠0, the form of the eigenvalue problem is dramatically altered. In particular, exponentially growing modes (*ω*=−i*λ*, *λ*>0) are only possible when
and hence when
Hence, unless either *V* _{ϕ} or *V* (or both) change sign somewhere in *a*<*r*<*b*, an exponentially growing instability is not possible (Knobloch 1992). Instead, the instability takes the form of growing oscillations (*ω* is complex). This is a consequence of the fact that the product *V* _{ϕ}*V* _{z}≠0 determines the sign of the helicity of the magnetic field. As a consequence, for fixed magnetic helicity the resulting perturbation equations are no longer invariant under the reflection (Knobloch 1996). As a result, the symmetry group of the problem with periodic boundary conditions in *z* is the group *SO*(2) of *proper* rotations of a circle. In contrast, when *V* _{ϕ}=0, the equations are invariant under the group *O*(2) of rotations *and* reflections of a circle. This distinction has profound consequences since the generic symmetry-breaking (i.e. *n*≠0) bifurcation in the former case is a Hopf bifurcation, while in the latter case both steady-state and Hopf bifurcations are of codimension one. Ecke *et al.* (1992) provide a simple illustration of these otherwise abstract concepts.

In the present case, the symmetry-breaking bifurcation introduces *z* dependence into the solution, and hence produces an axisymmetric *wave* that travels in either +*z* or −*z* directions, depending on the sign of the imposed magnetic helicity, i.e. the waves predicted by the theory have a preferred direction of propagation. Recent solution of the corresponding linear stability problem for non-ideal MHD in the Couette–Taylor geometry (Hollerbach & Rüdiger 2005) has confirmed these predictions and moreover shows that the critical magnetic Reynolds number, *Rm*^{helical}_{c}, required for the onset of the helical MRI is substantially lower than the corresponding critical magnetic Reynolds number, *Rm*^{non-helical}_{c}, for the non-helical case with imposed poloidal magnetic field only. It should be noted that, in periodic geometry in the *z* direction, the onset of helical MRI predicted by theory corresponds to the onset of *convective instability*, that is, instability that grows in the frame travelling with the group velocity of the waves. In the convectively unstable regime, a localized initial disturbance propagates in the direction of the group velocity faster than it spreads upstream. As a result at any given location the instability eventually dies out. In contrast, for *Rm*>*Rm*^{abs} the upstream growth of the disturbance exceeds the downstream advection at the group velocity, and the instability eventually grows at every fixed location (Huerre & Monkewitz 1990; Chomaz 1992; Tobias *et al.* 1998; Liu *et al.* 2006*a*). In particular, in bounded (i.e. non-periodic) domains, the presence of the boundaries ultimately leads to the decay of all disturbances in the convectively unstable regime, *Rm*^{helical}_{c}<*Rm*<*Rm*^{abs}, unless these are maintained by noise (Proctor *et al.* 2000). As a result in such domains, however large, the instability grows only when *Rm*>*Rm*^{abs}, up to corrections in inverse powers of the domain length, *L*_{z}. In accretion discs, this condition is likely to be satisfied, although this is not necessarily the case in the helical MRI experiments (Liu 2009; Priede & Gerbeth 2009).

This prediction has in turn led to a successful detection of the helical MRI mode in the PROMISE I experiment at Rossendorf (Stefani *et al.* 2006, 2007). This experiment uses the Couette–Taylor geometry, with an inner rotating cylinder to set up an angular velocity profile that decreases outwards. When the Ekman forcing at the lids is reduced with the help of split rings (as in the Rossendorf PROMISE II experiment or the Princeton experiment (Ji *et al.* 2006; Obabko *et al.* 2008)), the observed modes behave qualitatively much as predicted by the theory (Knobloch 1992, 1996; Tobias *et al.* 1998; Hollerbach & Rüdiger 2005).

It should be pointed out that when *V* _{z}=*V* _{z}(*r*), *V* _{ϕ}≠0 the structure of the problem is changed dramatically because of the presence of critical layers at locations where *ω*=±*nV* _{z}(*r*). Such critical layers can absorb and partially reflect incoming wavetrains and produce global instability modes. As a result, the properties of the eigenvalue problem are substantially altered. This case has not received attention in the astrophysics literature.

## 3. Evolution of the magnetorotational instability

For astrophysical applications we need to know what happens in a disc that is MRI-unstable. Specifically, we need to estimate the amplitude at which the MRI saturates in order to be able to estimate the efficiency of the resulting angular momentum transport, and for this we need to know whether the MRI can be treated as an ideal MHD problem throughout its evolution.

The evolution of the MRI has largely been studied using direct numerical simulations. These simulations have generally used the shearing box formulation of the problem, and until recently focused on ideal MHD. Figure 1*a* shows an example from the original paper by Hawley & Balbus (1991). The simulation is two-dimensional (i.e. axisymmetric), compressible and was started from an initial condition consisting of a localized slug of constant vertical field. Figure 1*b* illustrates the instability in the presence of an imposed radial field. The figure shows that the instability evolves into inward and outward radially moving sheets. These sheets are coupled by the magnetic field but numerical dissipation leads to reconnection. In the presence of any large-scale vertical field, the initially large vertical wavenumber of the instability gradually coarsens. The final state is unclear. Goodman & Xu (1994) suggested that the instability saturates as a result of a secondary shear instability between the inward and outward moving sheets. However, such instabilities are strongly affected by viscosity, and viscous processes are also typically omitted from the simulations.

Attempts to resolve the dissipative processes leading to reconnection were made by Sano *et al.* (1998) in the context of a compressible shearing sheet simulation. These authors found, based on integrations spanning 300 rotation periods, that the instability saturated when the Elsasser number, , but not otherwise. Subsequent non-axisymmetric computations suggest that saturation does occur even when *Λ*≥1 provided non-axisymmetric disturbances are permitted to evolve (Fleming *et al.* 2000).

These considerations have recently led to a re-examination of the role of dissipative processes in the saturation of the MRI. A careful study of the numerical dissipation in the ZEUS code for ideal MHD shows that the effective magnetic Prandtl number of the code, *Pm*≡*ν*_{num}/*η*_{num}>1 (Fromang & Papaloizou 2007); unfortunately in almost all circumstances prevailing in accretion discs, *Pm*<1 (Balbus & Henri 2008). Moreover, the efficiency of angular momentum transport appears to decrease with improved resolution (Fromang & Papaloizou 2007), suggesting that the transport rates estimated on the basis of ideal MHD simulations are affected by grid-scale dissipation and overestimate the efficiency of angular momentum extraction by the MRI. Related computations for an incompressible shearing box have been performed by Lesur & Longaretti (2007), who emphasize that the saturated turbulence level generated by the growth of the fastest MRI mode is not determined by the growth rate of this mode and in fact scales as *Pm*^{δ}, *δ*>0. The shearing box formulation is critically reviewed by Regev & Umurhan (2008).

The above results emphasize the need for an alternative approach (Knobloch & Julien 2005). This approach uses multi-scale asymptotics and takes explicit advantage of the separation between the short rotation period, the longer Alfvén crossing time and the even longer dissipative time scales. While this approach cannot treat accretion discs in detail, it is extremely helpful in identifying the various physical processes that are involved in the saturation of the (axisymmetric) MRI. The approach emphasizes the role played by these different time scales and indicates that complete saturation is only achieved on the slowest of these three time scales (Knobloch & Julien 2005) and thus beyond the reach of existing simulations.

The approach also emphasizes the important role played by the extraction of energy, by the growing MRI, from the background shear profile. In an astrophysical accretion disc this profile is determined, in steady state, by the gravitational field of the central object. However, depending on the growth rate of the instability, the instability may be fast enough to reduce the background shear, thereby quenching the instability. Since the gravitational field eventually restores the original shear profile, the growth rate of the instability will increase again, potentially leading to an oscillatory or otherwise time-dependent saturated state (Julien & Knobloch 2006). We believe that this outcome is very likely in global computations in which the *Ω*-effect quickly generates a strong azimuthal magnetic field, resulting in the helical MRI already mentioned. We believe that the existing shearing box formulations are deficient in that they do not permit any dynamics involving the background shear, which is taken as imposed and constant in both space and time—in addition to the concerns expressed by Regev & Umurhan (2008).

## 4. Theory

To allow for these possibilities, we describe here a formulation of the problem that does permit departures from Kepler flow associated with the instability, i.e. we recognize that the imposition of a difference in angular velocity across a (local part of the) disc represents a particular choice of the boundary conditions, and that a better choice might be to impose an *a priori* unknown angular velocity difference, and then determine its value self-consistently from the simulations using radial force balance. We find that such solutions do exist, and that the resulting disc may be supported radially by a combination of mechanical and magnetic pressure. The resulting one-parameter family of solutions allows us to examine the mechanisms responsible for the saturation of the instability. We find that saturation occurs at the expense of the angular velocity distribution and requires the presence of both resistivity *and* viscosity, potentially explaining the absence of saturation noted in two-dimensional (compressible) simulations when *Λ*>1 (Sano *et al.* 1998).

In order to extract the essence of the problem, we adopt the model considered originally by Balbus & Hawley (1991), i.e. we take a local model of a disc, with homogeneous structure in both radial and vertical directions, and constant imposed toroidal and poloidal magnetic fields. We assume, however, that the flow is incompressible. We do this for simplicity since compressibility effects are not fundamental to the instability (Balbus & Hawley 1991). The drawback of this assumption is that it makes a detailed comparison with the simulations impossible: the latter pose the problem in terms of the local *β* of the disc, i.e. the ratio of thermal to magnetic pressure, while in our formulation the absolute value of the (mechanical) pressure is arbitrary. We then identify a particular regime that permits us to study the evolution of the instability in the *strongly* nonlinear regime, i.e. far from threshold. The technique takes advantage of the small vertical scale of the instability, which is used as an expansion parameter and is thus ideally suited for fingering instabilities such as the MRI (Julien & Knobloch 2007). The approach generalizes to cylindrical geometry, but this aspect of the problem is omitted here.

### (a) The shearing sheet approximation

We consider a straight channel , filled with an electrically conducting incompressible fluid, and rotating about the *z* axis with constant angular velocity . We suppose that a linear shear flow is maintained in the channel, for example, by lateral boundaries that slide in the *y* direction with speeds , or through the shearing box approximation where a local expansion is performed about the global shear profile *Ω**(*r**) of a disc at with . In addition, we suppose that a constant magnetic field is present and consider *y*-independent perturbations of this state. In dimensionless form, these can be written as ** u**≡(

*u*,

*v*,

*w*)=(−

*ψ*

_{z},

*v*,

*ψ*

_{x}),

**≡(**

*b**a*,

*b*,

*c*)=(−

*ϕ*

_{z},

*b*,

*ϕ*

_{x}), and satisfy the equations (Julien & Knobloch 2006) 4.1 4.2 4.3and 4.4 where is proportional to the Alfvén speed associated with the imposed poloidal (vertical) magnetic field,

*J*(

*f*,

*g*)≡

*f*

_{x}

*g*

_{z}−

*f*

_{z}

*g*

_{x}, and , , and represent the dimensionless rotation rate, shear, kinematic viscosity and ohmic diffusivity, non-dimensionalized using a velocity scale

*U** and the channel width . The above equations are valid provided , where

*H** is the vertical pressure scale height and is the (horizontal) length scale associated with the rotational profile at . Note that drops out of these equations; as already mentioned, this is not so in an annulus, where hoop stresses are present, and the MRI becomes oscillatory (Knobloch 1992, 1996).

For normal modes of the form equations (4.1)–(4.4) yield the linear dispersion relation
4.5
where *p*^{2}≡*k*^{2}+*n*^{2}. This equation predicts a positive growth rate for sufficiently small wavenumbers *p*, whenever *σ*<0,*v*_{A}≠0, provided only that *ν* and *η* are sufficiently small. Moreover, for fixed *p*, the maximum growth rate is associated with the *x*-independent mode *k*=0.

Much of the work related to the MRI has focused on the astrophysical regime *Ω*,|*σ*|≫*v*_{A}≫*ν*,*η*. That is, the shear is the dominant source of energy for the instability, but the instability itself requires the presence of a (weaker) vertical magnetic field. Dissipative effects are weaker still and are often ignored entirely. This ideal limit (*ν*=*η*=0) is associated with the dispersion relation
4.6
In this case, the MRI occurs in the wavenumber range , where is the ideal MRI stability parameter and *n*>0. The maximum growth rate and associated wavenumber are given by
4.7
Equation (4.5) also yields an expression for the threshold for the instability, i.e. *λ*=0, which may be interpreted as defining a critical Elsasser number, *Λ*≡*v*^{2}_{A}/*Ωη*, namely,
4.8

Exact, fully nonlinear, but monochromatic solutions of equations (4.1)–(4.4) are known (Craik & Criminale 1986; Goodman & Xu 1994). For the gravest mode, *k*=0, these take the form
4.9
For this type of solution, all nonlinear terms vanish identically, and in the absence of three-dimensional parasitic instabilities, the solution grows exponentially without bound, i.e. (*Ψ*(*t*),*V* (*t*),*Φ*(*t*),*B*(*t*))∼e^{λt}. Thus, solutions of this type exhibit no saturation, a result supported by ideal numerical simulations of the MRI with constant imposed shear (Balbus & Hawley 1991, 1998). In these situations the volume-averaged turbulent shear stress,
also grows exponentially without bound. However, as shown below, if we relax the requirement that the background shear cannot change, the situation changes dramatically: the exponential growth is arrested by the change in the background shear and saturation occurs even in the absence of dissipation.

## 5. Asymptotic theory

Motivated by the astrophysical context, we suppose that the rotation time scale is relative to the Alfvén crossing time, while the dissipative processes act on the long time scale , where *ϵ*≪1, *δ*≪1. Specifically, we take *U**=*v**_{A} so that *v*_{A}=1. Then *η*=*S*^{−1}, where is the Lundquist number, and *ν*=*Pm* *S*^{−1}, where *Pm*≡*ν**/*η** is the magnetic Prandtl number. Likewise, |*σ*|=*Rm* *S*^{−1} and *Ω*=(*Ω**/|*σ**|)*Rm* *S*^{−1}, where *Rm*≡*σ***L*^{*2}/*η** is the magnetic Reynolds number. Thus, , and hence *ϵ*,*δ*≪1 implies that . Moreover, the stability parameter , implying a large spectrum of unstable modes, while . The results are illustrated in figure 2. For *ϵ*=*o*(*δ*), we find *Λ*≫1, and the role of dissipation in the MRI is subdominant. The dispersion relation therefore asymptotes to the ideal limit (4.6) (see solid curve in figure 2). In contrast, when , and the dispersion relation asymptotes to (4.5) (see dashed curve in figure 2).

We therefore write (Knobloch & Julien 2005; Jamroz *et al.* 2008*a*,*b*)
5.1
where *ϵ*≪1, *δ*≪1. We also pose a multiple scales expansion in the *x* direction, with ∂_{x} replaced by *δ*^{−1}∂_{x}+∂_{X}, where *X*=*δx*. In addition, vertical derivatives are large, , as are time derivatives . The inclusion of small vertical and horizontal scales allows us to fully retain all small-scale processes; the large scale *X* is required to allow averaged properties to vary across the channel. In parallel with the above assumptions, we need to make an assumption about the relative magnitude of the various fields. Here, the rapid shearing by the azimuthal flow suggests that we take and .

In what follows we separate the fields into mean and fluctuating components,
where the overbar denotes the average in both space 〈 〉_{V} and time 〈 〉_{t},
5.2
Averaging equations (4.2) and (4.4) yields the mean azimuthal equations
5.3and
5.4
where *J*_{X}(*f*,*g*)≡∂_{X}*f*∂_{z}*g*−∂_{z}*f*∂_{X}*g*. The fluctuating equations are
5.5
5.6
5.7and
5.8
where *J*_{x}(*f*,*g*)≡∂_{x}*f*∂_{z}*g*−∂_{z}*f*∂_{x}*g*, and .

We now expand the variables *ψ*,*v*,*ϕ*,*b* in terms of the small parameters *ϵ*,*δ* in the form , with similar expressions for the other three fields. From equations (5.5)–(5.8) we obtain *v*′_{0j}=*b*′_{0j}≡0 to all orders. Given this fact and that *v*′_{10}=0 from equation (5.5) at 𝒪(*δ*^{−1/2}) and *b*′_{10}=0 from equations (5.6) and (5.8) at , we find that the first non-trivial balances occur at in equations (5.5) and (5.7), and at 𝒪(*ϵ*^{1/2}*δ*^{1/2}) in equations (5.6) and (5.8).

The reduced fluctuating equations obtained are 5.9 5.10 5.11 and 5.12

Thus, nonlinear and dissipative terms are retained at leading order when *ϵ*=𝒪(*δ*) (case A) but these become subdominant when *ϵ*=*o*(*δ*) (case B), corresponding to *Λ*=𝒪(1) and *Λ*≫1, respectively. Both cases are of astrophysical interest.

The above asymptotic analysis captures the feedback mechanism arising through the appearance of local large-scale quantities, and . This type of feedback is also considered by Ebrahimi *et al.* (2009). The closure of the fluctuating equations (5.9)–(5.12) requires the computation of the Reynolds, Maxwell and mixed stresses. These enter the mean equations (5.3) and (5.4) at giving
5.13and
5.14
We are permitted to integrate once to obtain explicitly the quantities needed to close the system,
5.15and
5.16
where the constant *C* is determined by radial force balance across the channel in the saturated state (Julien & Knobloch 2006). Relations (5.15) and (5.16) indicate how the Reynolds, Maxwell and mixed stresses generated by the MRI feed back upon the imposed shear and toroidal magnetic field. Thus, they couple the dynamics of the small-scale instability to the large scale. This feedback is only present when dissipative effects are included, even at subdominant order; measures the decrease in the imposed shear arising from the MRI.

The fact that and are constant implies that and are linear functions in *X*. This property is identical to that possessed by the projection, via Taylor expansion, of the imposed background velocity shear profile *Ω*(*r*) onto the local geometry and emphasizes that such terms cannot be omitted. Indeed, on defining the imposed background shear stress, , the compensating viscous shear stress, , and the averaged turbulent shear generated by the MRI,
5.17
we see that . This states that the total shear stress within the shearing box is conserved and equal to *T*_{σ}. For outward transport of angular momentum, required of an accretion disc, , implying that and hence that the MRI *reduces* the imposed shear gradient. We can therefore evaluate the efficiency of angular momentum transport by examining the quantity . Similarly, for the mixed stress terms we define
so that
5.18

For numerical efficiency, in the computations that follow, only spatial averaging, denoted by is performed. For validity, this assumption demands that our integration of the reduced equations reaches a statistically steady state with . This is always the case.

## 6. Numerical simulations of the reduced equations: case A

### (a) Single-mode theory

In case A (*δ*=*ϵ*, *Λ*=𝒪(1)) equations (5.9)–(5.16) possess exact, saturated, single-mode solutions. These take the form
6.1
For these solutions, all Jacobians vanish identically and the only surviving nonlinearities arise through interactions between fluctuating and mean quantities. The resulting saturated state is given by the nonlinear dispersion relation (Knobloch & Julien 2005)
6.2
The resulting modifications to the shear and toroidal field are given by
6.3
where
6.4
Thus, when *C*=0. These results are reproduced in numerical simulations of equations (5.9)–(5.16) that start from initial conditions of the form (6.1) with periodic boundary conditions in the radial direction and chosen to maximize the linear theory growth rate (Jamroz *et al.* 2008*a*). The calculation shows that the analytically computed saturated state is stable (figure 3*a*).

### (b) Stress-free boundary conditions

In the presence of impenetrable but stress-free boundary conditions in the radial direction, boundary layers develop near the boundaries at *X*=0,*L*_{X}. In this case, a spatially modulated monochromatic initial condition of the form
where is the wavelength of the fastest growing mode and *L*_{X}≡*ϵ*^{−1}*L*=10*L*, evolves to the state shown in figure 4. Except in the boundary layers, the fields *ψ*′_{00} and *b*′_{11} have zero mean in the *z* direction, while the means of *v*′_{11} and *ϕ*′_{00} vary linearly with *X*. In this case the nonlinear dispersion relation (6.2) is no longer exact and the solution saturates at a slightly lower level (figure 3*b*). Figure 5 compares the single-mode analytical predictions with the results of numerical integration of the equations when *δ*=*ϵ*=0.1.

### (c) Random initial conditions

Figure 6 shows the evolution of *ϕ*′_{00} starting from a small-amplitude random perturbation *ψ*′_{00} and *ϕ*′_{00}=*v*′_{11}=*b*′_{11}=0. As before, the fastest growing linear wavenumber dominates the flow in the kinematic stage. As saturation is approached, the flow coarsens to larger and larger wavelengths in the vertical direction until the flow is dominated by the lowest vertical wavenumber supported on the computational domain. The single-mode theory (Knobloch & Julien 2005; Julien & Knobloch 2006) predicts this coarsening of the flow to the lowest wavenumber supported by the computational domain. Although an initially random field has no bias towards the single-mode theory and includes a mixture of many modes, the value of 〈∂_{X}*v*_{00}〉_{V} for the coarsened, saturated flow matches the single-mode theory for a mode with the largest wavelength supported by the computational domain. The saturated values of 〈∂_{X}*v*_{00}〉_{V} with random initial conditions with varying vertical box lengths are shown in figure 7 and follow the single-mode theory. Figure 8 shows that the Maxwell stresses dominate Reynolds stresses in the transport of angular momentum (cf. Sano *et al.* 1998). Moreover, while the values of the fluctuating quantities increase with , the back-reaction 〈∂_{X}*v*_{00}〉_{V} saturates (figure 8), as found in the single-mode approach.

## 7. Numerical simulations of the reduced equations: case B

We now present the corresponding results for case B (*ϵ*=*o*(*δ*), *Λ*≫1). All numerical simulations presented in this section were performed at parameter values together with the local Keplerian relation

### (a) Single-mode solutions

Figure 9*a* illustrates the evolution on a log–log scale of the single-mode solution (6.1) with real initial data, and in the reduced equations. For these solutions owing to their parity in the *z* direction; consequently (since *C*=0). Figure 9 reveals a subsequent transition to a regime of slower *algebraic* growth or decay of the mode amplitudes: *t*^{1/2} for the variables *Φ*_{0} and *V* _{0}, and *t*^{−1/2} for the variables *Ψ*_{0} and *B*_{0}. In addition to this algebraic behaviour, all variables have a superposed oscillation. This is evident in the algebraically decaying modes, where the oscillation swamps the decay. Figure 10*a* shows the corresponding evolution of 〈∂_{X}*v*_{00}〉_{V}. Remarkably, one can see that, although there is unbounded algebraic growth and decay in the solution for the MRI mode, the value of as . Thus 〈∂_{X}*v*_{00}〉_{V} saturates, as does the transport of angular momentum in the box. The saturated value is determined analytically in §8*b*.

### (b) Multi-mode solutions

Allowing for several simultaneous modes in the vertical yields more insight into the algebraically growing MRI mode. The computational domain selected has a length of 8*L* in the vertical, where is the wavelength of the fastest growing mode (Balbus & Hawley 1998), and *L*_{X}≡*ϵ*^{−1}*L*=10*L* in the horizontal. Similar results hold for smaller values of *ϵ*. We select a uniform distribution of wavelengths spanning a range of modes from to twice the cutoff wavelength . The mode amplitudes are sampled from a uniform distribution with upper bound 10^{−4}. Figure 11 shows the time evolution of the *z*-profile of *ϕ*_{00}, normalized with respect to its r.m.s. value. Because of the absence of quadratic Jacobian nonlinearities at *Λ*≫1, power is only retained in the range of unstable modes originally specified in the initial data. Moreover, although the power in these modes grows algebraically, the relative growth has an inverse relationship with the wavenumber Hence, at large times, the largest permissible wavelength in the box dominates. This can be seen in both figures 9(*b*) and 10(*b*): initially the flow evolution is dominated by modes with wavenumber near the maximum growth rate , but the solution coarsens with increasing time. Figure 10*b* shows that the exponential growth has again been curbed and that, in the algebraic regime (figure 9*b*), the perturbations produce a statistically steady 〈∂_{X}*v*_{00}〉_{V}. Also included in this figure are the theoretical predictions (8.5) of the saturated state based on and that based on the smallest wavenumber permitted by our initial data. The latter case is in good agreement with the simulations.

### (c) Random initial perturbations

Figure 12 illustrates the evolution of the flux function *ϕ*_{00} from an initial condition with random amplitudes of modes in *both* *x* and *z* in the computational domain used in §7*b*. The initial state is composed of modes within with a uniform sampling of amplitudes. As the flow evolves, all linearly stable spectral modes (satisfying ) are damped and the solution field is dominated by finger-like structures characteristic of the MRI. Figure 13*b* indicates that the perturbation fields once again ultimately evolve algebraically in time, and that the perturbations produce a statistically steady value of 〈∂_{X}*v*_{00}〉_{V}. Theoretical predictions (8.5) of the saturated state based on and the smallest permitted wavenumber, , are also indicated; the latter shows good agreement with the simulations. As in prior ideal (Balbus & Hawley 1998) and non-ideal (Julien & Knobloch 2006, 2007) simulations of the MRI where Jacobian nonlinear terms are retained, these simulations also exhibit coarsening of the solution with increasing time. However, since direct mode-to-mode coupling is absent, the flow only coarsens to the lowest wavenumber seeded by the initial condition, in the present case the box-filling mode and .

### (d) Dissipation and saturation

The above results indicate that, when , the kinematic regime characterized by exponential growth gives way with increasing time to a regime of persistent and weaker algebraic growth. Explicit inclusion of subdominant dissipation terms of type or in equations (5.9)–(5.12), where , cuts off this slow growth and results in the complete saturation of the instability. Computations indicate that saturation occurs primarily through the Reynolds and Maxwell stress terms, and not through viscous damping of the fluctuating fields.

## 8. Theory: case B

### (a) Exponential growth regime

The initial kinematic regime is characterized by exponential growth that is captured in detail by neglecting the nonlinear terms , in equations (5.9)–(5.12). In this regime, the instability evolves according to
8.1
(cf. Goodman & Xu 1994). This solution maximizes the velocity and magnetic correlations in the shear stress tensor , equation (5.17), corresponding to an outward transport of angular momentum. We have seen that, while this solution accurately matches the initial phase of the instability, it ignores the dominant saturation mechanism when *Λ*≫1: the coupling to the mean. The two-scale description given here shows that, as the perturbations grow, large-scale nonlinear terms become non-negligible because of the action of the MRI shear stress on the imposed background shear flow, and once this is the case the system can no longer be approximated by the linear system.

### (b) Algebraic growth and the saturation of angular momentum transport

The numerical simulations of the single-mode case using equations (5.9)–(5.12) show that the initial exponential growth gives way to algebraic growth/decay as time increases. Although the individual fields *ψ*_{00}, *v*_{11}, *ϕ*_{00}, *b*_{11} grow or decay algebraically, the quantity saturates at a constant value with small superposed oscillations (figure 9). This observation suggests the ansatz
8.2
where *α*>0 and *β*>0. Substitution into the reduced system (5.9)–(5.12) and (5.15)–(5.16) requires *α*=*β*=1/2 and yields, at leading order,
8.3and
8.4
together with the conditions
8.5and
8.6
The condition (8.5) is nothing but the requirement that the dispersion relation has a neutral mode; the frequency (8.6) is the frequency of the neutral mode when (8.5) holds and matches exactly the frequency visible in figure 9*a* (Jamroz *et al.* 2008*b*). Substituting the solution (8.2) into the turbulent stress balance (equation (5.15) with volume averaging) and using the form (8.3) and (8.4) of the eigensolutions yields
8.7
This expression does not contain secular terms proportional to or to *t*^{−1/2} sin *ωt*, let alone terms that grow in proportion to *t*, and hence saturates despite the algebraic growth of the contributing fields. This is a consequence of phase mixing: phase mixing introduces damping into the system that is analogous to Landau damping and this damping competes with energy injection by the MRI to produce a saturated state. To find the saturated amplitude we time-average the above result and compare it with equation (8.5):
8.8

### (c) Multi-mode saturation

Using the preceding single-mode analysis, we can now interpret the results obtained for more general initial conditions. For both multi-modes and the non-homogeneous multi-modes, we obtain a regime of exponential growth followed by a transition to algebraic growth. As in the single-mode case we see that the value of 〈∂_{X}*v*_{00}〉_{V} approaches a statistically steady value in the algebraic growth regime (figures 10*b* and 13*b*). We use equation (8.5) with and , where is the smallest wavenumber seeded by the initial condition, to define quantities we call and . Figures 10*b* and 13*b* show that, although the mode corresponding to has the largest growth rate and dominates the kinematic phase, it is the mode with wavenumber that dominates the long-time behaviour and is a good predictor of the saturated value of 〈∂_{X}*v*_{00}〉_{V}.

## 9. Conclusion

In this review we have described some recent developments connected with the MRI instability. We have focused on topics of theoretical interest at the expense of applications, be they astrophysical or laboratory. We have noted that understanding of the saturation of the instability is key to understanding the transport of angular momentum in an accretion disc, and pointed out the multiple time-scale nature of the saturation process. Thus, simulations must be run over much longer times, and all dissipative processes fully resolved.

We have pointed out that simulations that do not permit any modification, whether in space or time, of the background shear that feeds the instability miss a key saturation mechanism that operates even when dissipative processes are subdominant (case B), and have shown that the inclusion of this effect converts the exponential growth of single-mode solutions into much weaker algebraic growth. Despite this continued growth the correction to the shear profile does saturate, via phase mixing of the type familiar from Landau damping, and does so on a rotation time scale. We have seen that the saturated state can be predicted analytically, and that the result agrees well with the simulations starting from multi-mode initial conditions or from small-amplitude random initial conditions provided only that one uses the vertical wavenumber corresponding to the final coarsened state (Jamroz *et al.* 2008*b*).

In contrast, in case A, where *ϵ*=*δ*≪1 and , dissipative effects lead directly to saturation of the amplitude of the fingering modes, albeit on the dissipative time scale. Numerical simulations show that the single-mode theory of Knobloch & Julien (2005), summarized in §6*a*, captures the essence of the saturation process (Jamroz *et al.* 2008*a*).

Our formulation of the MRI employs a negative shear, *σ*<0. The sign implicitly defines the ‘outward’ direction, and our result implies that always. Thus, in our reduced model the MRI always acts to transfer angular momentum outwards. It is of interest to obtain a relation estimating the Shakura–Sunyaev *α* parameter (Shakura & Sunyaev 1973). Defining , where *c**_{s}=*Ω***H** is the local sound speed and is the *r*,*θ* component of the local stress tensor (), in our asymptotic regime
where *L**/*H**≪1 denotes the characteristic MRI length scale determined in terms of the vertical pressure scale height. This result is in qualitative agreement with recent three-dimensional shearing sheet simulations (Fromang & Papaloizou 2007), which suggest that angular momentum transport decreases with increasing resolution.

Of course, in three dimensions the MRI has additional saturation mechanisms available beyond modification of the background shear. Of these, the most promising are the parasitic (i.e. secondary) instabilities of the interpenetrating angular momentum sheets (Goodman & Xu 1994; Pessah & Goodman 2009). While there is vertical shear already associated with the radial motion, angular momentum conservation during initial stages of the instability also generates vertical shear owing to changes in the azimuthal velocity. These instabilities therefore require a full three-dimensional treatment. Simulations suggest that in shearing boxes of approximately unit aspect ratio these instabilities lead to strong intermittency in the associated Maxwell stress (Lesur & Longaretti 2007), while in boxes of large radial extent relative to height this intermittency is largely suppressed (Bodo *et al.* 2008). Our two-scale approach focuses on such horizontally extended boxes, but does not determine the stability properties of our solutions with respect to three-dimensional disturbances. We expect that these solutions are unstable to such perturbations but that their presence is reflected in the statistical properties of the resulting MRI turbulence, as occurs in other turbulent flows (Kawahara & Kida 2001; Viswanath 2007).

Existing simulations of the global problem in cylindrical geometry are not constrained by the constant imposed shear, and consequently are not encumbered by this highly restrictive assumption (Kersalé *et al.* 2006; Liu *et al.* 2006*b*; Ebrahimi *et al.* 2009), and in this sense are better able to capture the dominant saturation processes (Regev & Umurhan 2008). One possibility is that the energy extracted by the MRI goes into amplifying the magnetic field, although perhaps only on a resistive time scale (Ebrahimi *et al.* 2009). However, no simulations can reach the large Reynolds numbers for which the asymptotic approach described here is designed.

## Acknowledgements

We are grateful to B. Jamroz for assistance with numerical computations. This work was supported by NASA award NNG05GD37G (K.J.) and a University of Colorado SEED Grant (K.J.).

## Footnotes

↵1 In this paper, asterisks denote dimensional quantities.

One contribution of 13 to a Theme Issue ‘Turbulent mixing and beyond’.

- © 2010 The Royal Society