# Magnetic self-organization in Hall-dominated magnetorotational turbulence

Matthew W. Kunz, Geoffroy Lesur

► To cite this version:

Matthew W. Kunz, Geoffroy Lesur. Magnetic self-organization in Hall-dominated magnetorotational turbulence. Monthly Notices of the Royal Astronomical Society, 2013, 434, pp.2295-2312. 10.1093/mnras/stt1171. insu-03614296

HAL Id: insu-03614296

https://insu.hal.science/insu-03614296

Submitted on 20 Mar 2022

**HAL** is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire **HAL**, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Magnetic self-organization in Hall-dominated magnetorotational turbulence

Matthew W. Kunz<fnref n="1" />* and Geoffroy Lesur<fnref n="2" />†

Accepted 2013 June 24. Received 2013 June 21; in original form 2013 May 3

## ABSTRACT

The magnetorotational instability (MRI) is the most promising mechanism by which angular momentum is efficiently transported outwards in astrophysical discs. However, its application to protoplanetary discs remains problematic. These discs are so poorly ionized that they may not support magnetorotational turbulence in regions referred to as ‘dead zones’. It has recently been suggested that the Hall effect, a non-ideal magnetohydrodynamic (MHD) effect, could revive these dead zones by enhancing the magnetically active column density by an order of magnitude or more. We investigate this idea by performing local, three-dimensional, resistive Hall-MHD simulations of the MRI in situations where the Hall effect dominates over Ohmic dissipation. As expected from linear stability analysis, we find an exponentially growing instability in regimes otherwise linearly stable in resistive MHD. However, instead of vigorous and sustained magnetorotational turbulence, we find that the MRI saturates by producing large-scale, long-lived, axisymmetric structures in the magnetic and velocity fields. We refer to these structures as zonal fields and zonal flows, respectively. Their emergence causes a steep reduction in turbulent transport by at least two orders of magnitude from extrapolations based upon resistive MHD, a result that calls into question contemporary models of layered accretion. We construct a rigorous mean-field theory to explain this new behaviour and to predict when it should occur. Implications for protoplanetary disc structure and evolution, as well as for theories of planet formation, are briefly discussed.

**Key words:** accretion, accretion discs – instabilities – MHD – protoplanetary discs – stars: formation.

## 1 INTRODUCTION

Protoplanetary discs are poorly ionized. This fact casts doubt upon whether the most promising mechanism for enhanced angular-momentum transport in accretion discs, the magnetorotational instability (MRI; Balbus & Hawley 1991, 1998), is capable of driving the observationally inferred mass-accretion rates in these systems (e.g. Hartmann et al. 1998). Not only are the most potent sources of ionization (e.g. cosmic rays, stellar X-rays, UV radiation) shielded over significant portions of these discs, but also the presence of dust grains is anticipated to remove an appreciable fraction of the charge from the gas phase. All this makes it improbable that protoplanetary discs are magnetically coupled across their full radial and vertical extents (for a review, see Armitage 2011). This is particularly true in the dense mid-plane, where the principal source of ionization is likely to be weak radioactivity.

There are reasons to believe, however, that such pessimism is unwarranted. The MRI appears to be much more resilient in the face of diffusive losses than one may at first suspect. While a low degree of ionization is known to decouple the MRI-unstable charged species from the bulk neutral fluid in processes known as ambipolar diffusion (Blaes & Balbus 1994; Hawley & Stone 1998; Desch 2004; Kunz & Balbus 2004; Bai & Stone 2011; Simon et al. 2013) and Ohmic dissipation (Jin 1996; Sano & Miyama 1999; Fleming, Stone & Hawley 2000), it is somewhat surprising that the critical ionization fraction at ~1 au is just ~10⁻¹³ for typical protoplanetary discs (e.g. Balbus 2011). This is because the charged species have ~1 yr to communicate the magnetic field to the neutrals via collisions. This renders the innermost ($r \lesssim 0.1$ au), outermost ($r \gtrsim 30$ au) and surface-layer ($z \gtrsim 0.1$ au) regions magnetically active, either by thermal ionization of metals or by unshielded ionizing radiation.

Of course, these numbers come with large uncertainties, and much of the research concerning the MRI in protoplanetary discs has boiled down to determining the extent of magnetically active regions by coupling chemical networks of increasing complexity to

<footnote n="1">¹Department of Astrophysical Sciences, 4 Ivy Lane, Peyton Hall, Princeton University, Princeton, NJ 08544, USA</footnote>

<footnote n="2">²UJF-Grenoble I/CNRS-INSU, Institut de Planétologie et d'Astrophysique de Grenoble (IPAG) UMR 5274, F-38041 Grenoble, France</footnote>

<footnote>* NASA Einstein Postdoctoral Fellow.</footnote>

<footnote>† E-mail: geoffroy.lesur@obs.ujf-grenoble.fr</footnote>
either linear stability analyses (e.g. Gammie 1996; Igea & Glass-
gold 1999; Sano et al. 2000; Salmeron & Wardle 2003, 2005, 2008;
Wardle & Salmeron 2012) or some sort of non-linear criterion based
upon numerical experiments (e.g. Fromang, Terquem & Balbus
2002; Ilgner & Nelson 2006; Bai & Goodman 2009; Bai 2011).
However, while determining the chemical abundances and conse-
quent diffusivities in such discs is without a doubt essential to
improving our understanding of disc stability, structure and evolu-
tion, the results presented in this paper suggest that the philosophy
driving this approach can be misleading.

Here, we take an alternative route to understanding magnetoro-
tational turbulence in protoplanetary discs. We forego a detailed
study of disc chemistry and instead concentrate on the turbulent
disc dynamics themselves. Employing non-linear numerical simu-
lations and mean-field theory, we investigate the impact of Ohmic
dissipation and the Hall effect on magnetorotational turbulence.
We extend previous work by Sano & Stone (2002a,b, hereafter
SS02) into the Hall-dominated regime, and obtain qualitatively new
results. Instead of vigorous and sustained magnetorotational turbu-
lence, we find that the Hall-MRI saturates by producing large-scale,
long-lived, axisymmetric ('zonal') structures in the magnetic and
velocity fields. Their emergence – a result of the antidiffusive nature
of the Hall effect when the Maxwell stress increases with magnetic-
field strength – causes a reduction in turbulent transport by at least
two orders of magnitude from extrapolations based upon resistive
magnetohydrodynamics (MHD).

Our results suggest that existing estimates of the depth of magne-
torotationally active layers in protoplanetary discs based on damp-
ing by Ohmic dissipation and ambipolar diffusion are likely to be
in error. This conclusion has been reached before by other authors
(Wardle & Salmeron 2012), but for different reasons. Those au-
thors put emphasis on the fact that the Hall effect can render a
disc linearly unstable (Wardle 1999; Balbus & Terquem 2001) even
in the presence of strong Ohmic and ambipolar diffusion. In this
case, the critical magnetic Reynolds number Rm\_crit for magnetoro-
tational turbulence ought to be smaller than is often assumed. By
contrast, our results suggest that, even when a disc is deemed mag-
netically active from the perspective of linear analysis, the actual
turbulent transport that results may be much too small to be con-
sidered 'active'. Even for discs in which Rm \gtrsim 10³, an order of
magnitude larger than what is usually considered the critical value,
the Hall effect can cause a turbulent bifurcation to a low-transport
state (LTS). As a result, Hall-dominated regions of protoplanetary
discs (r ~ 5–10 au), while magnetically active, may nevertheless
exhibit prohibitively low accretion rates.

The paper is organized as follows. In Section 2.1, we present
the governing shearing-sheet equations of resistive Hall-MHD. In
Sections 2.2 and 2.3, we highlight two key physical concepts encap-
sulated by these equations – the conservation of canonical vorticity
and the close connection between the transport of magnetic flux
and the transport of angular momentum. Section 2.4 places these
considerations in the context of protoplanetary discs, from which
we obtain numerical estimates of the four dimensionless free pa-
rameters in our system (Section 2.5). We close Section 2 by proving
that MRI ‘channel’ modes remain exact non-linear solutions despite
the complicating features of the Hall effect (Section 2.6); we defer
to Appendix A an investigation of their stability to secondary ‘par-
asitic’ modes. Section 3 presents the numerical approach we have
adopted, the tests we have employed to verify its stability and accu-
racy (see also Appendix B), and the results of using this approach to
study the Hall-dominated MRI in the linear and non-linear regimes.
These results motivate the construction of a mean-field theory that

explains both the emergence of zonal structures and the transition to
an LTS observed in our simulations (Section 4). Finally, in Section 5
we summarize our results and briefly comment on their implications
for protoplanetary discs and planetesimal formation.

2 PRELIMINARIES

2.1 Shearing-sheet equations

We adopt the shearing-sheet approximation (Goldreich & Lynden-
Bell 1965), a useful framework for describing phenomena that vary
on lengthscales much less than the large-scale properties of the disc.
A small patch of the disc, co-orbiting with a fiducial point r₀ in the
mid-plane of the unperturbed disc at an angular velocity Ω = Ω₀êz,
is represented in Cartesian coordinates with the x and y directions
corresponding to the radial and azimuthal directions, respectively.
Differential rotation is accounted for by including the Coriolis force
and by imposing a background linear shear, v₀ = 2A₀xêy, where

$$
A_0 = \frac{r_0}{2} \left. \frac{d\Omega(r)}{dr} \right|_{r=r_0}
$$

is the Oort 'A' value; Keplerian rotation yields $A_0 = -(3/4)\Omega_0$. We
take the flow to be incompressible, a good assumption when the
magnetic pressure is much less than the gas pressure.

The equations of motion are then

$$
\begin{equation}
\begin{split}
\frac{\partial \mathbf{v}}{\partial t} = -\mathbf{v} \cdot \nabla \mathbf{v} - \frac{1}{\rho} \nabla P + \frac{\mathbf{J} \times \mathbf{B}}{c\rho} - 2\Omega_0 \times \mathbf{v} - 4A_0 \Omega_0 x \hat{\mathbf{e}}_x \\
\qquad + \nu \nabla^2 \mathbf{v},
\end{split}
\tag{1}
\end{equation}
$$

$$
\frac{\partial \mathbf{B}}{\partial t} = \nabla \times \left( \mathbf{v} \times \mathbf{B} - \frac{\mathbf{J} \times \mathbf{B}}{en_e} \right) + \eta \nabla^2 \mathbf{B}, \quad (2)
$$

subject to the constraints

$$
\nabla \cdot \mathbf{v} = 0, \tag{3}
$$

$$
\nabla \cdot \mathbf{B} = 0. \tag{4}
$$

Our notation is standard: $\rho$ is the (homogeneous) mass density, $v$ is the velocity, $P$ is the gas pressure, $\mathbf{B}$ is the magnetic field and

$$
\mathbf{J} = \frac{c}{4\pi} \nabla \times \mathbf{B}
$$

is the current density. The number density of electrons $n_e$ is taken
to be constant and uniform, as are the viscosity $\nu$ and resistivity
$\eta$. For future reference, we also introduce the total number density
$n = \rho/m$, where $m$ is the mean mass per particle, and the ion mass
density $\rho_i = m_i n_i$. Quasi-neutrality (i.e. $n_e = Z n_i$) is assumed.

Henceforth, the subscript ‘0’ on $A_0$ and $\Omega_0$ is dropped.

2.2 Lorentz force, Hall effect and canonical vorticity

In the incompressible approximation, the gas pressure *P* is determined not by an equation of state, but rather by satisfying the incompressibility condition (3). In fact, it is customary to eliminate the pressure by taking the curl of equation (1) to obtain an evolutionary equation for the flow vorticity *ω* ≡ ∇ × **v** + 2Ω:

$$
\frac{\partial \omega}{\partial t} = \nabla \times \left( \mathbf{v} \times \omega + \frac{\mathbf{J} \times \mathbf{B}}{c\rho} \right) + \nu \nabla^2 \omega. \quad (5)
$$
The form of equation (5) is very similar to that of equation (2). Just as the Lorentz force changes the number of vortex lines threading a fluid element, the Hall effect (represented by the penultimate term in equation 2) changes the number of magnetic-field lines threading a fluid element. Indeed, the origin of the Hall term is the differential motion between the electrons, to which the magnetic-field lines are tied (modulo Ohmic losses), and the drifting ions, which we take to be collisionally well coupled to the bulk neutral fluid.

Since the divergences of both the vorticity and the magnetic field are zero, any new vortex and magnetic-field lines that are made must be created as continuous curves that grow out of points or lines where the vorticity and magnetic field, respectively, vanish. Put simply, just as the effect of the Lorentz force on the vorticity is non-dissipative, so too is the Hall effect on the magnetic field; vorticity and magnetic flux can only be redistributed by these processes. We now prove that they must be redistributed in a specific way.

Consider the canonical momentum

$$ \rho_{\text{canonical}} \equiv m(\mathbf{v} + \mathbf{\Omega} \times \mathbf{r}) + \frac{eA}{c} \frac{n_e}{n}, \qquad (6) $$

and the associated canonical vorticity

$$ \omega_{\text{canonical}} \equiv \frac{1}{m} \nabla \times \rho_{\text{canonical}} = \omega + \frac{eB}{mc} \frac{n_e}{n}, \qquad (7) $$

where $A$ is the magnetic vector potential satisfying $B = \nabla \times A$. Students of plasma physics will recognize the final term in equation (7) as the vectorized Hall frequency $\omega_H \equiv (eB/mc)(n_e/n)$, at and above which small-wavelength circularly polarized waves with left-handed polarization cannot propagate. Combining equations (2) and (5), we find that the canonical vorticity satisfies

$$ \frac{\partial \omega_{\text{canonical}}}{\partial t} = \nabla \times (\mathbf{v} \times \boldsymbol{\omega}_{\text{canonical}}) + \nabla^2 (\mathbf{v} \boldsymbol{\omega} + \eta \boldsymbol{\omega}_H). \qquad (8) $$

This equation states that, in the absence of dissipative sinks, the canonical vorticity is frozen into the fluid. As a result, the combined number of vortex and magnetic-field lines threading a material surface is conserved; i.e. the canonical circulation

$$ \Gamma_{\text{canonical}} \equiv \oint_C \boldsymbol{\rho}_{\text{canonical}} \cdot d\boldsymbol{l} \quad \left( = \frac{1}{m} \int_S \boldsymbol{\omega}_{\text{canonical}} \cdot d\boldsymbol{S} \right) $$

around a simple closed contour $C$ bounding a material surface $S$ is a constant. This is simply Kelvin's (1869) circulation theorem generalized for Hall-MHD. An important consequence is that a local increase in magnetic flux must be accompanied by a local decrease in vorticity flux and vice versa.

Such behaviour is absent in ideal MHD, in which the magnetic flux is conserved for each fluid element independent of how the vorticity is advected. The difference is due to the fact that, in Hall-MHD, the ion-neutral fluid drifts relative to the field lines and, as such, has its momentum augmented by the magnetic field through which it travels. One may think of this as a consequence of Lenz's law. We refer the reader to the review by Polygiannakis & Moussas (2001) for further discussion of conserved quantities in Hall-MHD.

## 2.3 Hall electric field and Maxwell stress

Many of the results in this paper stem from the realization that the Hall electric field may be re-written in the following form:

$$ \frac{\mathbf{J} \times \mathbf{B}}{c e n_e} = \nabla \cdot \left( \frac{\mathbf{B} \mathbf{B}}{4 \pi e n_e} \right), \qquad (9) $$

dropping the extra $\nabla B^2$ term with impunity. This form is particularly useful, as it underscores the connection between the evolution of the magnetic flux and the Maxwell stress

$$ M_{ij} \equiv \frac{B_i B_j}{4\pi}, $$

whose $xy$-component plays the dominant role in transporting angular momentum in MRI-driven turbulence. In other words, the transport of magnetic flux in a partially ionized accretion disc is immediately tied to the efficiency and nature of the angular-momentum transport. This, along with the conservation of canonical vorticity, will turn out to be an extremely important property for understanding the subsequent analytical and numerical results.

## 2.4 Hall effect in protoplanetary discs

The Hall effect becomes important on lengthscales $\lesssim \ell_H$, where

$$ \ell_H \equiv \frac{v_A}{\omega_H} = \left( \frac{m_i c^2}{4\pi Z^2 e^2 n_i} \right)^{1/2} \left( \frac{\rho}{\rho_i} \right)^{1/2} \qquad (10) $$

and $v_A \equiv B/(4\pi\rho)^{1/2}$ is the Alfvén speed. The lengthscale $\ell_H$ is the ion skin depth divided by the square root of the mass-weighted ionization fraction, and is independent of magnetic-field strength. The smaller the degree of ionization, the broader is the range of scales that can be appreciably affected. This is why the cold, dense regions of protoplanetary discs are so easily susceptible to the Hall effect.

While the complexity involved in diagnosing the ionization rates, chemical abundances and consequent diffusivities in actual protoplanetary discs cannot be overstated (for a review, see Wardle 2007), it helps if we have at least some handle, if only rough, on the importance of the Hall effect in such discs. Assuming an equilibrium balance between cosmic ray ionization and dissociative recombination leads to the scaling $n_e \propto \sqrt{n}$, for which $\ell_H$ is a constant dependent only upon the ionization and recombination rates and the mean mass per particle. Taking their respective values to be $\zeta_{cr} = 10^{-17} \text{ cm}^3 \text{ s}^{-1}$, $\alpha_{dr} = 5 \times 10^{-7} \text{ cm}^3 \text{ s}^{-1}$ and $m = 2.33 m_p$ (e.g. Umebayashi & Nakano 1990), we find $\ell_H \simeq 0.4$ au. Assuming the standard model of the minimum-mass solar nebula (MMSN; Hayashi 1981), the disc scaleheight

$$ H \simeq 0.03 \left( \frac{r}{1 \text{ au}} \right)^{5/4} \text{ au} $$

is comparable to $\ell_H$ at a radius $r \sim 10$ au, inside of which the entire disc thickness becomes subject to the Hall effect. Further in around $\sim 1$ au, there is enough column density to effectively shield cosmic rays; there, the Hall effect is likely to be most important away from the mid-plane where Ohmic losses are less severe.

The presence of a small fraction of dust grains ($\sim 10^{-4}-10^{-2}$ by mass) can affect the extent of the Hall-dominated region by a substantial, though highly uncertain, amount by soaking up gas-phase charges and altering the equilibrium balance of chemical reactions. In this paper, we circumvent these complications by restricting ourselves to radii $\sim 5-10$ au where both Hall and Ohmic diffusion are considered to be dominant, and by allowing for a range of fields strengths and ionization fractions in a parameter study. Readers interested in the details of protoplanetary disc chemistry, the consequent values of the diffusivities and the implications for the (linear) MRI may consult, e.g., Salmeron & Wardle (2008) and Wardle & Salmeron (2012).
## 2.5 Dimensionless free parameters

Solutions to equations (1)–(4) are governed by four dimensionless free parameters: the plasma beta

$$ \beta \equiv \frac{(\Omega H)^2}{v_A^2}, $$

the viscous Elsasser number

$$ \Lambda_v \equiv \frac{v_A^2}{v\Omega}, $$

the (classical) Elsasser number

$$ \Lambda_{\eta} \equiv \frac{v_A^2}{\eta \Omega} $$

and the Hall Elsasser number

$$ \Lambda_H \equiv \frac{\omega_H}{\Omega}. $$

The Hall Elsasser number can be cast in the more traditional form

$$ \Lambda_H \equiv \frac{v_A^2}{\eta H \Omega}, $$

by introducing an effective Hall resistivity, $\Lambda_H \equiv v_A^2/\omega_H = v_A\ell_H$; however, this form is somewhat specious as the Hall effect is not dissipative. We also define the Reynolds number,

$$ \mathrm{Re} \equiv \frac{\Omega H^2}{\nu}, $$

and the magnetic Reynolds number,

$$ \mathrm{Rm} \equiv \frac{\Omega H^2}{\eta}, $$

which differ from their Elsasser counterparts by a factor of $\beta$.

The ratio of viscosity to resistivity is known as the magnetic Prandtl number, $Pm = \nu/\eta$. While $Pm$ does not appear explicitly in the equations, it is known to affect the saturated state of magnetototational turbulence (Fromang et al. 2007; Lesur & Longaretti 2007; Simon & Hawley 2009; Longaretti & Lesur 2010). In protoplanetary discs, our primary systems of interest here, $Pm \ll 1$ and thus we restrict our attention to low-$Pm$ flows. However, it should be noted that high-$Pm$ flows could exhibit significantly different behaviour than what is presented here (e.g. Balbus & Henri 2008).

Taking the density, temperature and rotation frequency at a radius of 10 au in the MMSN, and assuming a 10 mG magnetic field and µm-sized dust grains (see fig. 1 in Salmeron & Wardle 2008), we find typical values for these parameters of $\beta \approx 1000$, $\Lambda_{\eta} \approx 0.5$, $\Lambda_H \approx 0.01$, $\ell_H \approx 3H$ and $Rm \approx 500$; for all practical purposes, $\mathrm{Re}$ is infinite in protoplanetary discs. We adopt similar parameters in our numerical simulations (see Table 1) with the exception of $\mathrm{Re}$; numerical constraints demand that its value be $\lesssim 10^4$.

## 2.6 Hall-MRI channel modes

It is well known that, in the absence of the Hall effect, equations (1)–(4) admit exact non-linear solutions referred to as MRI ‘channel’ modes:

$$ \begin{aligned} \mathbf{B} &= \mathbf{B}_0 + \mathbf{B}_{\mathrm{ch}} \\ &= B_0 \hat{e}_z + b e^{\nu t} B_0 \cos K z (\hat{e}_x \sin \theta - \hat{e}_y \cos \theta), \end{aligned} \quad (11) $$

$$ \begin{aligned} \mathbf{v} &= \mathbf{v}_0 + \mathbf{v}_{\mathrm{ch}} \\ &= 2A x \hat{e}_y + b e^{\nu t} v_0 \sin K z (\hat{e}_x \cos \phi + \hat{e}_y \sin \phi), \end{aligned} \quad (12) $$

where $\gamma$ is the growth rate of the mode, $K$ is its vertical wavenumber, $b$ is a dimensionless measure of the channel amplitude, and $v_0$ and $B_0$ are constants (Goodman & Xu 1994). In the absence of dissipation, the two (constant) orientation angles $\phi$ and $\theta$ are equal. These solutions are exact because all non-linearities vanish (i.e. $\mathbf{v}_{\mathrm{ch}} \cdot \nabla \mathbf{v}_{\mathrm{ch}} = \mathbf{B}_{\mathrm{ch}} \cdot \nabla \mathbf{B}_{\mathrm{ch}} = \mathbf{v}_{\mathrm{ch}} \cdot \nabla \mathbf{B}_{\mathrm{ch}} = \mathbf{B}_{\mathrm{ch}} \cdot \nabla \mathbf{v}_{\mathrm{ch}} = 0$).

What seems to have gone unappreciated in the Hall-MRI literature is that these channel modes remain exact non-linear solutions in the presence of the Hall effect, since the channel current density

$$ \mathbf{J}_{\mathrm{ch}} = -\frac{K c}{4\pi} b e^{\nu t} B_0 \sin K z (\hat{e}_x \cos \theta + \hat{e}_y \sin \theta) \quad (13) $$

satisfies $\mathbf{J}_{\mathrm{ch}} \cdot \nabla \mathbf{B}_{\mathrm{ch}} = \mathbf{B}_{\mathrm{ch}} \cdot \nabla \mathbf{J}_{\mathrm{ch}} = 0$. This result also follows from equation (9), since there are no non-linear z-components of the Maxwell stress for a magnetic field described by equation (11).

Defining a dimensionless Hall parameter,<fnref n="1" />

$$ \mathrm{Ha} \equiv \frac{K^2 B_0 c}{8 \pi \epsilon n_e \Omega} = \frac{1}{2} \frac{\Omega}{\omega_{H,0}} \left( \frac{K v_{A,0}}{\Omega} \right)^2, $$

and an effective magnetic Reynolds number,

$$ \mathrm{Rm}_{\mathrm{eff}} \equiv \frac{\Omega}{\eta K^2} \left( \frac{1 - \mathrm{Ha}}{1 - \mathrm{Pm}} \right), $$

the properties of the channel solution may be found after some straightforward but tedious algebra:

$$ \gamma + \nu K^2 = -A \sin 2\theta - \frac{\Omega}{Rm_{\mathrm{eff}}} \left(1 + \frac{A}{\Omega} \sin^2 \theta\right), \quad (14) $$

$$ \begin{aligned} (K v_{A,0})^2 = & -4\Omega^2 \left( \mathrm{Ha} + \frac{A}{\Omega} \sin^2 \theta \right) \\ & \times \left[ 1 + \frac{1}{4} \mathrm{Rm}_{\mathrm{eff}}^{-2} + \frac{A}{\Omega} \cos^2 \theta \left( 1 + \frac{1}{2} \mathrm{Rm}_{\mathrm{eff}}^{-1} \tan \theta \right)^2 \right], \end{aligned} \quad (15) $$

$$ K v_0 = -2\Omega \left( \mathrm{Ha} + \frac{A}{\Omega} \sin^2 \theta \right) \left( 1 + \frac{1}{4} \mathrm{Rm}_{\mathrm{eff}}^{-2} \right)^{1/2}, \quad (16) $$

$$ \phi = \theta - \tan^{-1} \left( \frac{1}{2} \mathrm{Rm}_{\mathrm{eff}}^{-1} \right). \quad (17) $$

In the limit $Rm_{\mathrm{eff}} \to \infty$, equations (15) and (16) differ from equations 6 and 7 of Goodman & Xu (1994) only by the multiplicative factor

$$ \left[ 1 + \frac{2\Omega}{\omega_{H,0}} \left( 1 + \frac{A}{\Omega} \cos^2 \theta \right) \right]^{-1}; $$

the growth rate in this limit, $\gamma = -A \sin 2\theta$, is identical to that of the ideal MHD case (their equation 5). In other words, while the Hall effect cannot alter the growth rate of these modes, it can increase or decrease their wavenumbers and Mach numbers when the angular velocity and background magnetic field are antiparallel or parallel, respectively (Wardle 1999).<fnref n="2" /> Physically, this sensitivity to the polarity of the magnetic field appears because the magnetic

<footnote n="1">¹ This definition differs from that used in Balbus & Terquem (2001) by a factor of 4 and in Kunz (2008) by a factor of $|\Omega/A|$.</footnote>

<footnote n="2">² For $\Omega/\omega_{H,0} \to +\infty$, the Hall parameter $\mathrm{Ha} \to -(A/\Omega) \sin^2\theta$ and both the channel wavenumber (equation 15) and the Mach number (equation 16) tend to zero as the destabilizing shear is marginalized. Exact equality can only be achieved for $\mathrm{Ha} = \theta = 0$ or $\mathrm{Ha} = -(A/\Omega) \sin^2\theta = 1$, both of which are spurious solutions not satisfied by the original equations of motion.</footnote>
<table><caption>**Table 1.** List of the runs discussed in this paper, along with their defining dimensionless parameters. The viscous, resistive and Hall Elsasser numbers all refer to their values in the initial state, as does the plasma beta parameter; the Hall lengthscale $\ell_H$ is constant. The growth rate $\gamma$ refers to the most unstable pure-$k_z$ mode available in the simulation domain. The time-averaged turbulent transport $\bar{\alpha}$ is obtained using data between $t = 100$ and 630 (unless otherwise noted).</caption><thead><tr><th>Name</th><th>L<sub>x</sub> × L<sub>y</sub> × L<sub>z</sub></th><th>n<sub>x</sub> × n<sub>y</sub> × n<sub>z</sub></th><th>Λ<sub>ν</sub><sup>-1</sup></th><th>Λ<sub>η</sub><sup>-1</sup></th><th>Λ<sub>H</sub><sup>-1</sup></th><th>ℓ<sub>H</sub></th><th>β</th><th>γ</th><th>ᾱ</th></tr></thead><tbody><tr><td>ZB111</td><td>4 × 4 × 1</td><td>256 × 128 × 64</td><td>0.2</td><td>1</td><td>0</td><td>0</td><td>1000</td><td>0.40</td><td>3 × 10<sup>-2</sup></td></tr><tr><td>ZB1H1</td><td>4 × 4 × 1</td><td>256 × 128 × 64</td><td>0.2</td><td>1</td><td>17.4</td><td>0.55</td><td>1000</td><td>0.71</td><td>1.7 × 10<sup>-4</sup></td></tr><tr><td>ZB1H1L</td><td>8 × 8 × 1</td><td>512 × 256 × 64</td><td>0.2</td><td>1</td><td>17.4</td><td>0.55</td><td>1000</td><td>0.71</td><td>1.4 × 10<sup>-4</sup></td></tr><tr><td>ZTB1H1</td><td>4 × 4 × 1</td><td>256 × 128 × 64</td><td>0.2</td><td>1</td><td>17.4</td><td>0.55</td><td>1000</td><td>0.67</td><td>3.2 × 10<sup>-4</sup></td></tr><tr><td>ZB1H2</td><td>4 × 4 × 1</td><td>256 × 128 × 64</td><td>0.2</td><td>1</td><td>4.2</td><td>0.13</td><td>1000</td><td>0.60</td><td>1.4 × 10<sup>-1</sup></td></tr><tr><td>ZB1H3</td><td>4 × 4 × 1</td><td>256 × 128 × 64</td><td>0.2</td><td>1</td><td>8.6</td><td>0.27</td><td>1000</td><td>0.62</td><td>6.8 × 10<sup>-2</sup></td></tr><tr><td>ZB1H4</td><td>4 × 4 × 1</td><td>256 × 128 × 64</td><td>0.2</td><td>1</td><td>13</td><td>0.41</td><td>1000</td><td>0.68</td><td>4.7 × 10<sup>-4</sup></td></tr><tr><td>ZB1H5</td><td>4 × 4 × 1</td><td>256 × 128 × 64</td><td>0.2</td><td>1</td><td>21.8</td><td>0.69</td><td>1000</td><td>0.70</td><td>1.2 × 10<sup>-4</sup></td></tr><tr><td>ZB1H6</td><td>4 × 4 × 1</td><td>256 × 128 × 64</td><td>0.2</td><td>1</td><td>30.4</td><td>0.97</td><td>1000</td><td>0.53</td><td>1.4 × 10<sup>-5</sup></td></tr><tr><td>ZB311</td><td>4 × 4 × 1</td><td>256 × 128 × 64</td><td>0.32</td><td>1</td><td>0</td><td>0</td><td>3200</td><td>0.40</td><td>1.9 × 10<sup>-2</sup></td></tr><tr><td>ZB3H1</td><td>4 × 4 × 1</td><td>256 × 128 × 64</td><td>0.32</td><td>1</td><td>-1</td><td>0.018</td><td>3200</td><td>0.25</td><td>6.7 × 10<sup>-3</sup></td></tr><tr><td>ZB3H2</td><td>4 × 4 × 1</td><td>256 × 128 × 64</td><td>0.32</td><td>1</td><td>1</td><td>0.018</td><td>3200</td><td>0.50</td><td>4.2 × 10<sup>-2</sup></td></tr><tr><td>ZB3H3</td><td>4 × 4 × 1</td><td>256 × 128 × 64</td><td>0.32</td><td>1</td><td>2</td><td>0.035</td><td>3200</td><td>0.55</td><td>6.6 × 10<sup>-2</sup></td></tr><tr><td>ZB3H4</td><td>4 × 4 × 1</td><td>256 × 128 × 64</td><td>0.32</td><td>1</td><td>4</td><td>0.071</td><td>3200</td><td>0.62</td><td>1.0 × 10<sup>-1</sup></td></tr><tr><td>ZB3H5</td><td>4 × 4 × 1</td><td>256 × 128 × 64</td><td>0.32</td><td>1</td><td>8</td><td>0.14</td><td>3200</td><td>0.64</td><td>1.1 × 10<sup>-1</sup></td></tr><tr><td>ZB3H6</td><td>4 × 4 × 1</td><td>256 × 128 × 64</td><td>0.32</td><td>1</td><td>16</td><td>0.28</td><td>3200</td><td>0.70</td><td>1.1 × 10<sup>-2</sup></td></tr><tr><td>ZB3H7</td><td>4 × 4 × 1</td><td>256 × 128 × 64</td><td>0.32</td><td>1</td><td>32</td><td>0.57</td><td>3200</td><td>0.66</td><td>4.7 × 10<sup>-5</sup></td></tr><tr><td>ZB3H8</td><td>4 × 4 × 1</td><td>256 × 128 × 64</td><td>0.32</td><td>1</td><td>56.57</td><td>1.0</td><td>3200</td><td>0.74</td><td>1.8 × 10<sup>-5</sup></td></tr><tr><td>ZB3H9</td><td>4 × 4 × 1</td><td>256 × 128 × 64</td><td>0.32</td><td>1</td><td>100</td><td>1.8</td><td>3200</td><td>0.53</td><td>6.3 × 10<sup>-7</sup></td></tr><tr><td>ZB3I2<sup>a</sup></td><td>4 × 4 × 1</td><td>256 × 128 × 64</td><td>0.32</td><td>4</td><td>0</td><td>0</td><td>3200</td><td>0.17</td><td>5.1 × 10<sup>-3</sup></td></tr><tr><td>ZB3H10<sup>a</sup></td><td>4 × 4 × 1</td><td>256 × 128 × 64</td><td>0.32</td><td>4</td><td>-1</td><td>0.018</td><td>3200</td><td>0.08</td><td>2.7 × 10<sup>-4</sup></td></tr><tr><td>ZB3H11</td><td>4 × 4 × 1</td><td>256 × 128 × 64</td><td>0.32</td><td>4</td><td>1</td><td>0.018</td><td>3200</td><td>0.24</td><td>1.3 × 10<sup>-2</sup></td></tr><tr><td>ZB3H12</td><td>4 × 4 × 1</td><td>256 × 128 × 64</td><td>0.32</td><td>4</td><td>2</td><td>0.035</td><td>3200</td><td>0.30</td><td>2.4 × 10<sup>-2</sup></td></tr><tr><td>ZB3H13</td><td>4 × 4 × 1</td><td>256 × 128 × 64</td><td>0.32</td><td>4</td><td>4</td><td>0.071</td><td>3200</td><td>0.39</td><td>5.1 × 10<sup>-2</sup></td></tr><tr><td>ZB3H14</td><td>4 × 4 × 1</td><td>256 × 128 × 64</td><td>0.32</td><td>4</td><td>8</td><td>0.14</td><td>3200</td><td>0.50</td><td>9.5 × 10<sup>-2</sup></td></tr><tr><td>ZB3H15</td><td>4 × 4 × 1</td><td>256 × 128 × 64</td><td>0.32</td><td>4</td><td>16</td><td>0.28</td><td>3200</td><td>0.56</td><td>2.0 × 10<sup>-2</sup></td></tr><tr><td>ZB3H16</td><td>4 × 4 × 1</td><td>256 × 128 × 64</td><td>0.32</td><td>4</td><td>32</td><td>0.57</td><td>3200</td><td>0.62</td><td>4.9 × 10<sup>-4</sup></td></tr><tr><td>ZB10I1</td><td>4 × 4 × 1</td><td>256 × 128 × 64</td><td>0.2</td><td>50</td><td>0</td><td>0</td><td>10000</td><td>0.00</td><td>&lt;10<sup>-20</sup></td></tr><tr><td>ZB10I2</td><td>4 × 4 × 1</td><td>256 × 128 × 64</td><td>0.2</td><td>1</td><td>0</td><td>0</td><td>10000</td><td>0.41</td><td>8.7 × 10<sup>-3</sup></td></tr><tr><td>ZB10H1<sup>a, b</sup></td><td>4 × 4 × 1</td><td>256 × 128 × 64</td><td>0.2</td><td>50</td><td>25</td><td>0.25</td><td>10000</td><td>0.19</td><td>2.1 × 10<sup>-1</sup></td></tr><tr><td>ZB10H2</td><td>4 × 4 × 1</td><td>256 × 128 × 64</td><td>0.2</td><td>50</td><td>50</td><td>0.5</td><td>10000</td><td>0.32</td><td>1.8 × 10<sup>-3</sup></td></tr><tr><td>ZB10H3</td><td>4 × 4 × 1</td><td>256 × 128 × 64</td><td>0.2</td><td>50</td><td>100</td><td>1</td><td>10000</td><td>0.47</td><td>6.7 × 10<sup>-4</sup></td></tr><tr><td>ZB10H4</td><td>4 × 4 × 1</td><td>256 × 128 × 64</td><td>0.2</td><td>50</td><td>200</td><td>2</td><td>10000</td><td>0.55</td><td>4.8 × 10<sup>-5</sup></td></tr><tr><td>ZB10H5</td><td>4 × 4 × 1</td><td>256 × 128 × 64</td><td>0.2</td><td>50</td><td>300</td><td>3</td><td>10000</td><td>0.41</td><td>1.3 × 10<sup>-5</sup></td></tr></tbody></table>

<sup>a</sup>These runs exhibit a relatively small growth rate, and so the time-averaging procedure is performed between $t = 300$ and 630 in order to eliminate the influence of initial transients.

<sup>b</sup>This run shows strong bursts of turbulence associated with the breakup of channel modes. A longer time average should be used to obtain a properly converged $\bar{\alpha}$; consequently, the value of $\bar{\alpha}$ for this run is not used in the discussion.

‘epicycles’ induced by the Hall electric field introduce a handedness that enhances (if $\Omega/\omega_{H,0} < 0$) or reduces (if $\Omega/\omega_{H,0} > 0$) the effective magnitude of the destabilizing shear. For

$$ \frac{\omega_{H,0}}{\Omega} < -2 \left(1 + \frac{A}{\Omega}\right) \left( -\frac{1}{2} \text{ for Keplerian rotation} \right), \quad (18) $$

this multiplicative factor becomes negative at all orientation angles, and the channel modes become stable whistler waves. Note that the Hall-dominated regime is not accessible to the MRI when $B_z < 0$.

By differentiating equation (14) with respect to $\theta$ and setting the result to zero, one may show that the most vigorously growing channel takes an orientation of

$$ \theta_{\max} = \phi_{\max} = \frac{\pi}{4} + \frac{1}{2} \cot^{-1}(2Rm_{\text{eff}}) \quad (19) $$

and exhibits the growth rate

$$ \gamma_{\max} + \nu K^2 = -A \left(1 + \frac{1}{4} Rm_{\text{eff}}^{-2}\right)^{1/2} - \frac{\Omega}{Rm_{\text{eff}}} \left(1 + \frac{A}{2\Omega}\right). \quad (20) $$

In the limit of vanishing viscosity and resistivity, this reduces to the well-known result that the maximum growth rate of the MRI is given by the Oort ‘A’ value (Balbus & Hawley 1991).

Being non-linear solutions, Hall-MRI channel modes are analytically amenable to a stability analysis. Such an analysis indicates that these channels are subject to secondary ‘parasitic’ instabilities similar to those found by Goodman & Xu (1994) and further studied by Pessah & Goodman (2009), Latter, Lesaffre & Balbus (2009) and Pessah (2010). We defer this calculation to Appendix A, as our numerical simulations indicate that they do not play an important role in the saturation of the MRI in the Hall-dominated regime.

# 3 SHEARING-BOX SIMULATIONS

## 3.1 Numerical approach

### 3.1.1 The SNOOPY code

Equations (1)–(4) are solved using SNOOPY, a three-dimensional (3D) incompressible spectral code. SNOOPY is based on the FFTW 3
library to compute 3D Fourier transforms using domain decomposition. It uses a low-storage, third-order Runge-Kutta (RK3) scheme to compute all the terms except for the linear dissipation terms, which are integrated implicitly. All non-linearities are solved using a pseudo-spectral algorithm, avoiding aliasing errors with the 2/3 rule (Canuto et al. 1988). Shearing-sheet boundary conditions are implemented by solving the equations in a sheared frame comoving with the mean Keplerian flow. This algorithm is similar to the FARGO scheme (Masset 2000), although it can be extended to arbitrary order in space and time. This procedure implies a periodic remap in Fourier space, which is performed every $\Delta t_{\text{remap}} = |2A|L_y/L_x$, where $L_x$ and $L_y$ are, respectively, the radial and azimuthal extents of the shearing box (Umurhan & Regev 2004). SNOOPY is now a well-tested code, having been used for both hydrodynamical and MHD problems such as resistive MRI, dynamo and subcritical transitions to turbulence.

We have implemented a new module in SNOOPY to account for the Hall effect. This scheme integrates the Hall term using the same RK3 and pseudo-spectral algorithms employed elsewhere in the code. To guarantee stability, the integration time step must be smaller than all of the physical time-scales of the system. When the Hall effect dominates the dynamics, the shortest physical time-scale is given by the whistler-wave time-scale at the grid, $\tau_w \equiv \omega_H^{-1}(\Delta x/\pi \ell_H)^2$. This constraint is quite demanding, as it implies that very short time steps are needed to ensure stability ($\Delta t < \sqrt{3} \tau_w$; see Appendix B). The numerical results presented in Section 3.2 are therefore obtained at limited resolution (typically 64 points per H); even at this resolution, a simulation of the Hall-dominated MRI requires ~50 times more computational time than an ideal MHD simulation.

### 3.1.2 Testing the Hall-MHD module

We have assessed the stability and accuracy of the Hall-MHD module in SNOOPY using two tests. First, we verified that linear waves can propagate stably in all three spatial directions while satisfying the linear dispersion relation (equation 36 of Balbus & Terquem 2001)

$$ \frac{\omega}{\omega_H} = k \ell_H \left[ \sqrt{1 + \left( \frac{k \ell_H}{2} \right)^2} \pm \frac{k \ell_H}{2} \right] \quad (21) $$

across a range of wavenumbers $k$. At small wavenumbers (low frequencies), these waves are circularly polarized Alfvén waves; at large wavenumbers, right-handed waves (plus sign) go over to the high-frequency whistler-wave branch, whereas left-handed waves (minus sign) are cut off at $\omega_H$. In Fig. 1, the numerical eigenfrequencies (circles) are overlaid on the two solutions (solid lines) of equation (21). The former were obtained by exciting a small-amplitude velocity perturbation $\delta v_x = 2 \times 10^{-5} \cos(kz)$ along a mean magnetic field $B_0 \hat{e}_z$ and Fourier transforming $v_x(z=0, t)$ in time. This procedure gives two peaks in the spectra, which correspond to the eigenfrequencies of the right- and left-handed waves. The agreement between the analytical and numerical solutions is very good all the way down to the grid scale (the Nyquist frequency $k_N \ell_H = 3.2\pi$), a benefit of SNOOPY's spectral decomposition. Note that these tests were carried out with $\nu = \eta = 0$, demonstrating that the code can stably propagate whistler waves without the need for explicit dissipation. In Appendix B, we prove analytically why this is true for the RK3 integration scheme used in SNOOPY.

For our second test, we simulated the Hall-shear instability (HSI) described by Kunz (2008, hereafter K08). The basic ingredients of

<figure><img src="image_3.png" /><figcaption>**Figure 1.** Comparison between the analytical dispersion relation (equation 21; solid lines) and the numerical eigenfrequencies (circles) of linear waves in SNOOPY.</figcaption></figure>

this instability are a background shear, which generates a streamwise magnetic-field component from a transverse one, and the Hall effect, which induces a circular polarization that (for certain wavevectors) conservatively reorients stream-wise magnetic fields into the transverse direction. The induced transverse component is sheared further and an exponentially growing instability ensues. Since the interaction between shear and the Hall effect is strongest when the motions implied by the shear lie in the same plane as the magnetic ‘epicycles’ induced by the Hall effect, the instability is maximized when the vorticity, wavevector and magnetic field all share a mutual axis: $(k \cdot B)(k \cdot \omega)$ must be negative for instability. Note that this instability occurs even in non-rotating systems stable to the ‘classical’ MRI. In Fig. 2, we present the results of two numerical experiments with $\Omega = 0$, $A = -3/4$ and $(\Lambda_v^{-1}, \Lambda_\eta^{-1}, \beta) = (0.2, 1, 1000)$. We find a linear instability with a growth rate $\gamma = 0.69$ when $\Lambda_H^{-1} = 17.4$, matching the theoretical growth rate of the most unstable mode ($\gamma \simeq 0.70$) within 2 per cent. When $\Lambda_H^{-1} = 0$, no instability occurs; instead, the sheared fluctuations resistively decay. These results confirm the accuracy of the Hall-MHD module in SNOOPY.

<figure><img src="image_4.png" /><figcaption>**Figure 2.** Evolution of kinetic energy in a non-rotating shearing box demonstrating the HSI. Solid line: with Hall effect and Ohmic dissipation; dot-dashed line: with Ohmic dissipation only. The dashed line represents the theoretical growth rate.</figcaption></figure>
### 3.1.3 Units and runs

For our Hall-MRI simulations, the equations are put in a dimensionless form by choosing units natural to the system. The unit of time $\tau = \Omega^{-1}$ is the inverse of the rotation frequency. Lengthscales are measured in units of the vertical box size $L_z$, and the magnetic-field strength is measured in units of the initial Alfvén speed $v_{A,0}$. We also introduce a dimensionless measure of the turbulent transport,

$$ \alpha = \frac{\delta v_x \delta v_y - \delta B_x \delta B_y}{L_z^2 \tau^{-2}}. $$

This definition is made equivalent to the Shakura & Sunyaev (1973) definition of $\alpha$ by setting $L_z = H$, which we do for all of our simulations. Keplerian rotation (i.e. $A = -3/4$) is assumed.

In this work, we only consider horizontally extended ('slab') shearing boxes with $L_x, L_y > L_z$. As discussed by Bodo et al. (2008), a slab configuration exhibits better convergence properties than bar configurations (as used by SS02) due to the absence of recurrent channel modes. Horizontally extended domains have also been shown to impact quantitative measurements; Bodo et al. (2008) reported that the time-averaged turbulent transport in an $L_x = 8$ slab is a factor ~2 smaller than the transport measured in an $L_x = 1$ bar. One should bear this in mind when comparing our results with those of SS02. For most of our simulations, we consider boxes with aspect ratio $L_x \times L_y \times L_z = 4 \times 4 \times 1$ and resolution $n_x \times n_y \times n_z = 256 \times 128 \times 64$; we have also performed two versions of our fiducial simulation, one with an extended domain ($8 \times 8 \times 1$) and one with increased resolution ($384 \times 192 \times 96$).

All the runs discussed in this paper are summarized in Table 1. Unless otherwise stated, the initial conditions are white noise<fnref n="3" /> on all components of the velocity and magnetic fields with a typical rms amplitude equal to $7 \times 10^{-7}$. We add to the magnetic field a constant poloidal component whose magnitude is determined by a choice of initial $\beta$. These runs have been integrated up to $t = 630$, corresponding to 100 orbits in a Keplerian shearing box; our fiducial run was integrated to $t = 2000$. Finally, we introduce an averaging procedure $\overline{\cdot}$, which denotes a volume and time average. The time average is performed using data only from $t = 100$ (or, in some cases, $t = 300$) to the end of the run, in order to avoid transient effects from the (arbitrary) initial conditions and the initial growth of channel modes.

We close this section by providing a translation between our dimensionless parameters and those employed by SS02. Those authors used $X = 2\Lambda_H^{-1}$ and $Re_M = \Lambda_\eta$, and focused mainly on initial values of $X = 1, 2$ and $4$. For their fiducial initial $\beta = 3200$, this implies $\ell_H/H \sim$ a few $\times 10^{-2}$. Here we are concerned with somewhat larger values of $X$, which place $\ell_H \sim H$ and render the entire disc prone to the Hall effect.

## 3.2 Results

### 3.2.1 Averaged turbulent transport and the Hall effect

The relationship between the averaged turbulent transport and the intensity of the Hall effect was first explored numerically by SS02. In order to make a proper comparison with this earlier work, we have reproduced runs Z6–Z9 from SS02 (named ZB3H1–ZB3H4 here). We have also explored the regime of very strong Hall effect, up to $\Lambda_H^{-1} = 100$ (run ZB3H8), corresponding to a Hall term 16 times

<figure><img src="image_5.png" /><figcaption>Figure 3. Volume-averaged turbulent transport $\overline{\alpha}$ as a function of the inverse Hall Elsasser number $\Lambda_H^{-1}$ from runs ZB3(I1, H1–H9) (black circles) and runs Z6–Z9 from SS02 (red crosses).</figcaption></figure>

larger than the largest one considered by SS02. The results of this first exploration are presented in Fig. 3.

In the regime of weak Hall effect ($\Lambda_H^{-1} \lesssim 5$), we qualitatively recover the SS02 results: increasing $\Lambda_H^{-1}$ increases the amount of turbulent transport. This demonstrates the sensitivity of disc stability to the orientation of the background magnetic field (Wardle 1999). The actual values of $\alpha$ we find, however, are significantly smaller than those found by SS02. This is probably due to the different aspect ratios of the shearing boxes employed.

For larger $\Lambda_H^{-1}$, a totally new behaviour appears: the transport follows a steep decline down to $\overline{\alpha} = 4.7 \times 10^{-5}$ at $\Lambda_H^{-1} = 32$. At first sight, one may think that this behaviour is due to the stabilization of the MRI when the Hall effect is too dominant. To check this, we present in Fig. 4 the growth rate of the largest channel modes for the parameters used in runs ZB3H(1–4). In each of these runs, there exists at least one channel mode with $\gamma > 0.5$; i.e. a vigorous linear instability is initially present in all of these runs. The

<figure><img src="image_6.png" /><figcaption>Figure 4. Linear growth rate $\gamma$ as a function of the inverse Hall Elsasser number $\Lambda_H^{-1}$ for simulations ZB3H(1–9). Each line corresponds to a channel mode with $k_z = 2\pi i/L_z$.</figcaption></figure>

<footnote n="3">³ The non-linear outcome of a simulation does not depend significantly on the initial conditions.</footnote>
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Volume-averaged turbulent transport $\alpha$ as a function of time for different inverse Hall Elsasser numbers $\Lambda_H^{-1}$ (runs ZB3I1, ZB3H3, ZB3H6, ZB3H7, ZB3H9).</figcaption></figure>

behaviour for $\Lambda_H^{-1} > 5$ indicates that a new saturation mechanism is at work, one that is not related to the linear properties of the flow.

In order to isolate the cause of this new saturation mechanism, we focus on the evolution of the box-averaged turbulent transport at different $\Lambda_H^{-1}$ (Fig. 5). Simulations with $\Lambda_H^{-1} \lesssim 5$ show the same qualitative behaviour exemplified in ideal and resistive MHD simulations of the MRI – a time-steady $\alpha \sim$ few $\times 10^{-2}$ with moderate transport spikes. The intermediate case ($\Lambda_H^{-1} = 16$) shows oscillations between a low-transport state (LTS $\equiv \alpha < 10^{-4}$) and a high-transport state (HTS) qualitatively similar to that seen in the $\Lambda_H^{-1} \lesssim 5$ runs. For $\Lambda_H^{-1} = 32$ and 100, the turbulence stays in the LTS after an initial burst caused by the breakup of the channel mode. Note that the qualitative behaviour observed during the linear phase in all of these simulations (i.e. $t \lesssim 50$) is very similar.

Inspecting snapshots of the two extreme cases (Fig. 6) reveals the origin of the LTS. In the ideal case ($\Lambda_H^{-1} = 0$; top), the flow exhibits turbulent fluctuations of $B_z$. This is the ‘traditional’ saturated state of the MRI as described by Hawley, Gammie & Balbus (1995) and others. In the Hall-dominated case ($\Lambda_H^{-1} = 32$; bottom), we observe a coherent, axisymmetric, large-scale structure in $B_z$, which we refer to as a *zonal magnetic field*. In this zonal-field configuration, the vertical magnetic flux is accumulated in some radial ($x$) region, leaving most of the box with a very weak $B_z$ (typically $|B_z| < 10^{-3}$). In this configuration, almost no turbulent activity is observed. Note that the total vertical magnetic flux is conserved, indicating that this feature is due to a *redistribution* of magnetic flux.

### 3.2.2 Characterizing the LTS

In the previous section, we described a new saturated state of the MRI called the LTS. This state is characterized by very weak turbulent transport, despite the presence of a vigorous linear instability in the initial equilibrium. Here we conduct a dedicated study of the LTS by examining our fiducial run ZB1H1 in detail. This simulation exhibits the same kind of LTS as described above, with an averaged turbulent transport $\bar{\alpha} = 1.8 \times 10^{-4}$. Increasing the resolution to 384 $\times 192 \times 96$ does not change the outcome of the saturated

state: a similar LTS is observed, demonstrating that our simulations have converged.<fnref n="4" />

We quantify the presence of a zonal field by defining an averaging procedure,

$$ \langle \cdot \rangle \equiv \frac{1}{L_y L_z} \iint dy dz, $$

and computing the evolution of the vertically and azimuthally averaged vertical component of the magnetic field $\langle B_z \rangle$. In the top panel of Fig. 7, we present the resulting space-time diagram. This diagram clearly exhibits a strong zonal field with a typical radial thickness $\sim 1$. Outside of this zonal-field region, the averaged field is weak with $|\langle B_z \rangle| \lesssim 10^{-2}$. A closer inspection shows that the system initially exhibits two zonal-field regions, centred at $x \simeq 0.3$ and $x \simeq 1.8$. At $t \simeq 160$ a rapid reorganization occurs, and these two regions merge to produce one zonal field that survives for more than $1000 \Omega^{-1}$. While zonal fields are generally very long-lived structures in isolation, this demonstrates that they may become strongly unstable when another is nearby.

It was shown in Section 2.2 that, in the presence of the Hall effect, a new conserved quantity replaces the magnetic flux: the canonical vorticity. Since canonical vorticity is conserved without dissipative effects (e.g. viscosity, resistivity), we expect the formation of axisymmetric bands of vertical vorticity that are anticorrelated with the zonal fields. To check for this effect, we have computed the vertically and azimuthally averaged vertical component of the flow vorticity $\langle \omega_z \rangle = -\partial_x \langle v_y \rangle + 2A$ in run ZB1H1. The resulting space-time diagram is shown in Fig. 7(b), and clearly demonstrates the formation of a zonal-vorticity region akin to a zonal flow. In accordance with expectations from conservation of canonical vorticity, we also find that the flow vorticity is anticorrelated with the vertical magnetic field. However, the vorticity and magnetic field do not have exactly the same shape – the mean vorticity appears to be concentrated around the edges of the zonal-field region. This difference is due to the explicit dissipation, and in particular due to the fact that $Pm \ll 1$: magnetic-field lines diffuse quite rapidly through the bulk ion/neutral fluid, whereas vortex lines do not. Therefore, conservation of canonical vorticity is only approximately verified in our simulations, owing to the presence of non-negligible dissipative terms.

Since our box size is limited, one may suspect that the presence of only one zonal-field region in run ZB1H1 is an artefact of the boundary conditions. To check this, we have performed a simulation in a wider box ($8 \times 8 \times 1$) with the same physical parameters as run ZB1H1. The space-time diagram of this simulation (run ZB1H1L) is presented in Fig. 8. We observe the formation of two zonal-field regions of size $\approx 1.5$, which survive for the remainder of the simulation. This indicates that zonal-field regions have an intrinsic width independent of the radial and azimuthal box size (provided the latter is significantly larger than $H$).

In order to verify that the formation of a zonal field is not dependent upon the initial conditions of the simulation, we have also run a purely resistive ($\Lambda_H^{-1} = 0$) MRI simulation with the same parameters as run ZB1H1. This simulation was run up until $t = 630$, at which point the Hall effect was switched on with $\Lambda_H^{-1} = 17.4$. Within approximately four orbits, the fully developed 3D turbulence disappeared and was replaced by a large-scale zonal field with an averaged turbulent transport $\bar{\alpha} \sim 10^{-4}$. This demonstrates that the

<footnote n="4">⁴Due to the extreme cost of such a high-resolution simulation, this particular run was stopped at $t = 160$; it is therefore not listed in Table 1.</footnote>
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Snapshot of the vertical magnetic field $B_z$ in runs ZB3I1 (top) and ZB3H7 (bottom) at $t = 630$.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 7.** Space-time diagram of the vertically and azimuthally averaged vertical component of the magnetic field $\langle B_z \rangle$ and the vorticity $\langle \omega_z \rangle$ in run ZB1H1. The appearance of vorticity bands, anticorrelated with the zonal-field structures, is observed.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 8.** Space-time diagram of the vertically and azimuthally averaged vertical component of the magnetic field $\langle B_z \rangle$ in run ZB1H1L. Two stable zonal-field regions are produced.</figcaption></figure>

LTS and the zonal field associated with it are robust non-linear features of the Hall-dominated MRI.

To understand how this zonal-field structure is sustained, we return to the argument given in Section 2.3. In particular, it was shown that the Maxwell stress directly enters into the induction equation through the Hall effect. We therefore compute the mean stress and magnetic field in run ZB1H1, averaging these quantities in $y$, $z$ and time (from $t = 500$ to 600). The resulting profiles are presented in Fig. 9, and exhibit a very clear correlation between the averaged Maxwell stress $\langle M_{xy} \rangle$ and the magnetic-field profile $\langle B_z \rangle$. In particular, there are inflection points in the Maxwell stress at the boundaries of the zonal-field region. Between these points where the magnetic field is relatively strong, the MRI is magnetically quenched; elsewhere, the field is relatively weak and
<figure><img src="image_11.png" /><figcaption>**Figure 9.** $(y, z, t)$ averages of the turbulent stress (top) and magnetic field (bottom) for run ZB1H1. A clear correlation between the Maxwell stress and the mean vertical magnetic field is exhibited, a feature which is responsible for the formation and sustainment of a zonal magnetic field (see Section 3.2.2).</figcaption></figure>

turbulent fluctuations persist with a low level of transport. As we show explicitly in Section 4, such residual fluctuations are responsible for maintaining the integrity of the zonal-field structure in the face of resistive (both molecular and turbulent) diffusion.

### 3.2.3 A numerical criterion for the LTS

We have shown that an LTS exhibiting axisymmetric ('zonal') fields emerges in several simulations of the Hall-dominated MRI. To make any prediction about the saturation level, one must know when the system will choose the LTS instead of the 'classical' turbulent MRI state (HTS). To this end, we have systematically explored the parameter space $(\Lambda_{\eta}, \Lambda_{\nu}, \beta)$. In this parameter space, we include a regime which is stable without the Hall effect (runs ZB10XX). We present in Fig. 10 the growth rates of the most unstable Hall-MRI modes present in some representative runs from ZB3XX and ZB10XX. As expected, all the runs but ZB10I1 are linearly unstable with growth rates $\gamma > 0.1$. Note that the vertical wavelength of the most unstable mode increases with $\ell_H$ (Section 2.6).

In Fig. 11, we summarize all our results on a single plot exhibiting the mean turbulent stress $\bar{\alpha}$ as a function of $\ell_H$. Despite the

<figure><img src="image_12.png" /><figcaption>**Figure 10.** Location of some representative runs from ZB3XX (top) and ZB10XX (bottom) on a linear stability diagram. The colour contours denote the growth rates of the most unstable eigenmodes.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 11.** Mean turbulent stress versus Hall effect for $\beta = 1000, 3200$ and $10000$, and $\Lambda_{\eta}^{-1} = 1, 4$ and $50$. The turbulent stress decreases sharply when $\ell_H \gtrsim 0.2$, separating the LTS from the classical HTS.</figcaption></figure>

differing initial $\beta$, viscosities and resistivities, all of the values of $\bar{\alpha}$ tend to collapse on to a single curve dependent primarily upon $\ell_H$. We find that the system stays in the HTS up to $\ell_H \simeq 0.2$ for all our simulations, independent of the mean-field strength and of the resistivity; in this case, the typical turbulent stress $\bar{\alpha} \sim 10^{-2}–10^{-1}$. Beyond $\ell_H \sim 0.2$, the system transitions rapidly to an LTS
characterized by $\alpha \lesssim 10^{-3}$. A very important characteristic is that the presence of the LTS is not correlated to any linear property of the Hall-MRI. For instance, all of the situations where we observe an LTS are characterized by strongly unstable MRI modes with $\gamma \gtrsim 0.5$ (compare Figs 10 and 11). This illustrates the fact that a flow which is strongly unstable to the Hall-MRI does not necessarily evolve into fully developed turbulence.

### 3.2.4 Other magnetic-field configurations

One may suspect that the behaviour described above is very specific to situations with only a mean vertical magnetic field. We have therefore run a simulation with a magnetic field inclined by 45° in the azimuthal (y) direction of strength $\beta = 10^3$ (run ZTB1H1). We find results very similar to the pure-$B_z$ run ZB1H1. The transport is increased by a factor $\sim 2$, but the large-scale zonal-field structure characteristic of the LTS remains.

We have also performed a limited number of simulations either without a mean magnetic field or with only a mean azimuthal magnetic field, finding no sustained zonal-field configurations in either case. While the mean-field theory presented in the next section does suggest an explanation for this negative result, a more systematic exploration of these cases is needed to verify our conclusions. This is deferred to a future publication.

## 4 MEAN-FIELD THEORY OF ZONAL FIELDS AND FLOWS

In this section, we formulate a mean-field theory that explains the observed transport bifurcation from an HTS to an LTS and the appearance of zonal magnetic fields and flows. We start by separating the velocity and magnetic fields into fluctuating and non-fluctuating parts:

$$v = v_0 + \langle v \rangle + \delta v \quad \text{and} \quad B = \langle B \rangle + \delta B.$$

Upon averaging over azimuth and height, equations (3) and (4) become $\partial_x \langle v_x \rangle = 0$ and $\partial_x \langle B_z \rangle = 0$, respectively; it follows from equation (2) that $\langle B_x \rangle = 0$ if it is so initially. For clarity of presentation, we further assume that $\langle v_x \rangle = 0$. Not only is this assumption supported by our numerical results, but it also allows us to cancel global epicyclic oscillations non-essential for understanding the emergence of zonal fields and flows.

Introducing the Reynolds stress,

$$R_{ij} \equiv \rho \delta v_i \delta v_j,$$

and the Faraday tensor,

$$F_{ij} \equiv \delta v_i \delta B_j - \delta v_j \delta B_i,$$

the pertinent mean-field equations are the z-component of the averaged induction equation (2),

$$\frac{\partial \langle B_z \rangle}{\partial t} = - \frac{\partial \langle F_{xz} \rangle}{\partial x} - \frac{c}{en_e} \frac{\partial^2 \langle M_{xy} \rangle}{\partial x^2} + \eta \frac{\partial^2 \langle B_z \rangle}{\partial x^2}, \quad (22)$$

and the z-component of the averaged vorticity equation (5),

$$\frac{\partial \langle \omega_z \rangle}{\partial t} = -\frac{1}{\rho} \frac{\partial^2 \langle R_{xy} \rangle}{\partial x^2} + \frac{1}{\rho} \frac{\partial^2 \langle M_{xy} \rangle}{\partial x^2} + \nu \frac{\partial^2 \langle \omega_z \rangle}{\partial x^2}. \quad (23)$$

The mean vorticity,

$$\langle \omega_z \rangle = \frac{\partial \langle v_z \rangle}{\partial x} + 2A,$$

includes a contribution from the background shear (=2A, which is negative in Keplerian discs). Note that the averaged Maxwell stress $\langle M_{ij} \rangle$ comprises products of only the fluctuating magnetic fields.

In order to solve equations (22) and (23), we must construct models for $\langle F_{ij} \rangle$, $\langle R_{ij} \rangle$ and $\langle M_{ij} \rangle$. The Faraday tensor has been shown to be accurately modelled by a turbulent resistivity with coefficient $\eta_t$ (Lesur & Longaretti 2009), and we take this to be the case in what follows. We follow a similar approach for the Reynolds stress by modelling it as a turbulent viscosity with coefficient $v_t$. Adopting these simplifications, our mean-field equations become

$$\frac{\partial \langle B_z \rangle}{\partial t} \approx (\eta + \eta_t) \frac{\partial^2 \langle B_z \rangle}{\partial x^2} - \frac{c}{en_e} \frac{\partial^2 M}{\partial x^2}, \quad (24)$$

$$\frac{\partial \langle \omega_z \rangle}{\partial t} \approx (v + v_t) \frac{\partial^2 \langle \omega_z \rangle}{\partial x^2} + \frac{1}{\rho} \frac{\partial^2 M}{\partial x^2}. \quad (25)$$

We consider two models for the Maxwell stress, each of which will produce zonal behaviour very similar to that seen in our non-linear numerical simulations.

### 4.1 Case I: $\langle M_{xy} \rangle = M(\langle B_z \rangle)$

As a first approach, we take the xy-component of the Maxwell stress to be a function only of the local vertical magnetic flux,

$$\langle M_{xy} \rangle \equiv M(\langle B_z \rangle),$$

and we concentrate on the evolutionary evolution for the mean vertical magnetic field (equation 24). We suppose that there is some $\langle B_z^0 \rangle$ that satisfies this equation in steady state (e.g. $\langle B_z^0 \rangle$ constant) and we examine small deviations $\langle B_z^1 \rangle$ about that state. Linearizing equation (24), we find that such deviations satisfy

$$\frac{\partial \langle B_z^1 \rangle}{\partial t} \approx \left( \eta + \eta_t - \frac{c}{en_e} \frac{\mathrm{d}M}{\mathrm{d}\langle B_z \rangle} \right) \left. \frac{\partial^2 \langle B_z^1 \rangle}{\partial x^2} \right|_{\langle B_z^0 \rangle} \quad (26)$$

This equation has a simple interpretation. While resistivity acts diffusively on $\langle B_z^1 \rangle$, the Hall term may be diffusive or antidiffusive depending upon the local gradient of the Maxwell stress.

Fortunately, even without a specific model for $M$, progress can be made. For sufficiently large values of $\langle B_z \rangle$, we expect the MRI to be stable and $M \to 0$. We also expect $M \to 0$ for sufficiently small values of $\langle B_z \rangle$, since unstable modes exist only at small wavelengths where Ohmic dissipation becomes important and suppresses turbulent transport. In between these extremes, we know that the Maxwell stress is negative since the MRI transports angular momentum outwards. Therefore, $d^2M/d\langle B_z \rangle^2 < 0$, and so there must be a value of $\langle B_z^0 \rangle = B_{z,crit}$ above which $dM/d\langle B_z \rangle > 0$ and below which $dM/d\langle B_z \rangle < 0$.

First, let us consider $\langle B_z^0 \rangle < B_{z,crit}$. Then both the Ohmic and Hall contributions to equation (26) are positive and any deviations from steady state diffusively decay. Now let us consider the opposite case, $\langle B_z^0 \rangle > B_{z,crit}$. Then the Ohmic and Hall contributions to equation (26) have opposite signs and so the Hall effect acts antidiffusively. If the Hall effect can overcome diffusive processes, any increment in the local magnetic flux continues to grow and contract until $\langle B_z \rangle$ becomes large enough for $M \to 0$. By flux conservation, there must be accompanying patches of decreased magnetic flux, which we anticipate having low levels of turbulent transport as well. We associate this scenario with the transition to the LTS.

We now make these ideas concrete by specifying a simple model for the Maxwell stress,

$$M = -0.01(\gamma/\gamma_{\text{max}}), \quad (27)$$
<figure><img src="image_14.png" /><figcaption>Figure 12. Top: space-time diagram of the vertically and azimuthally averaged vertical component of the magnetic field ($\langle B_z \rangle$) in a toy model (Section 4.1) based upon our mean-field theory. The appearance of a zonal magnetic field is observed. Bottom: regions of increased field strength (red crosses) evolve antidiffusively towards $\langle M_{xy} \rangle = 0$, whereas regions of decreased field strength (blue circles) evolve diffusively towards an LTS.</figcaption></figure>

where $\gamma$ is the growth rate obtained by solving the dispersion relation (equations 14–17) with $K = 2\pi H^{-1}$, $B_0 = \langle B_z \rangle$ and $A = -3/4$; its dependence on $\langle B_z \rangle$ is shown in the bottom panel of Fig. 12 as the solid line. While we do not advocate such a crude relationship between the non-linear Maxwell stress and the linear properties of the Hall-MRI, this model does satisfy all of the qualitative expectations for $M$ highlighted in the previous two paragraphs. (One could have equally well approximated $M$ by a non-positive function quadratic in $\langle B_z \rangle$ with upward concavity.)

Using equation (27), we solve equation (24) in a periodic domain of length $L_x = 4$, in which a uniform vertical magnetic field of strength $\langle B_z \rangle = 0.01$ is initially disturbed by low-amplitude white noise.<fnref n="5" /> The result is graphically presented in the top panel of Fig. 12 and bears a striking resemblance to the results of the fully non-linear numerical simulations (top panel of Fig. 7). A zonal magnetic field is produced, with accompanying regions of low magnetic-field strength.

To quantify the evolution of $\langle B_z \rangle$, data are collected at fixed intervals in time along lines of constant $x$ where the magnetic field achieves its global maximum and minimum; these points are marked in the top panel of Fig. 12 by the white crosses and white circles, respectively. These data are mapped on to the solid line (equation 27) in the bottom panel of Fig. 12. Regions of locally in-

creased magnetic-field strength (red crosses) find themselves associated with antidiffusive transport (i.e. $dM/d\langle B_z \rangle > 0$), further increasing their magnetic-field strength and moving to the right along the curve until $M = 0$. By contrast, regions of locally decreased magnetic-field strength (blue circles) find themselves associated with diffusive transport (i.e. $dM/d\langle B_z \rangle < 0$), further decreasing their magnetic-field strength and moving to the left along the curve as $M \to 0$. Note that $M \neq 0$ in the low-field region; otherwise, the necessary inflection point in $M(x)$ would vanish and the zonal field would diffuse away.

## 4.2 An analytical criterion for the LTS

From these results, we deduce that a transition from an HTS to an LTS can occur provided that

$$ \left. \frac{c}{en_e} \frac{dM}{d\langle B_z \rangle} \right|_{\langle B_z^0 \rangle} > \eta + \eta_t \quad (28) $$

for some $\langle B_z^0 \rangle$. To make this criterion more quantitative, we must estimate the size of each of these terms in the HTS. First, we neglect molecular resistivity, which we assume to be small compared to the turbulent resistivity. Secondly, we assume that the turbulent resistivity is related to the turbulent transport via $\eta_t \sim \bar{\alpha} \Omega H^2 / Pm_t$, where $Pm_t \sim 2$ is the turbulent Prandtl number estimated for this particular component of the turbulent resistivity by Lesur & Longaretti (2009).<fnref n="6" /> Finally, we approximate $dM/d\langle B_z \rangle$ by $\bar{\alpha} \rho (\Omega H)^2 / B_{z,\text{stab}}$, using the definition of $\bar{\alpha}$ and introducing $B_{z,\text{stab}}$ as the critical magnetic-field strength above which the longest wavelength MRI mode is stabilized by magnetic tension. Using these estimates, the bifurcation condition (28) reduces to

$$ \ell_H \gtrsim \frac{v_{A,\text{stab}}}{\Omega Pm_t} \quad (29) $$

where $v_{A,\text{stab}} \equiv B_{z,\text{stab}} / (4\pi\rho)^{1/2}$. Solving the linear dispersion relation for $v_{A,\text{stab}}$, our criterion becomes $\ell_H \gtrsim 0.2H$, tantalizingly close to the value deduced from Fig. 11.

## 4.3 Case II: $\langle M_{xy} \rangle = M(\langle B_z \rangle, \langle \omega_z \rangle)$

Equation (25) indicates that the gradient of the Maxwell stress also affects the mean vorticity. While the model for the Maxwell stress presented in Section 4.1 is successful at explaining the bifurcation from an HTS to an LTS, it does not take into account this effect, nor does it take into account the feedback of a vorticity-dependent Maxwell stress on the evolution of the magnetic field. Here, we generalize the form of the Maxwell stress to allow for this interplay:

$$ \langle M_{xy} \rangle = M(\langle B_z \rangle, \langle \omega_z \rangle). $$

To leading order in the perturbation amplitudes, equations (24) and (25) then become

$$ \frac{\partial \langle B_z^1 \rangle}{\partial t} \approx \left( \eta + \eta_t - \frac{c}{en_e} \frac{\partial M}{\partial \langle B_z \rangle} \right) \frac{\partial^2 \langle B_z^1 \rangle}{\partial x^2} - \frac{c}{en_e} \frac{\partial M}{\partial \langle \omega_z \rangle} \frac{\partial^2 \langle \omega_z^1 \rangle}{\partial x^2}, \quad (30) $$

$$ \frac{\partial \langle \omega_z^1 \rangle}{\partial t} \approx \left( \nu + v_t + \frac{1}{\rho} \frac{\partial M}{\partial \langle \omega_z \rangle} \right) \frac{\partial^2 \langle \omega_z^1 \rangle}{\partial x^2} + \frac{1}{\rho} \frac{\partial M}{\partial \langle B_z \rangle} \frac{\partial^2 \langle B_z^1 \rangle}{\partial x^2}, \quad (31) $$

<footnote n="5">⁵ In order to regularize the solutions of our mean-field equations, we have used a hyper-resistivity $\alpha \propto K^4$ to damp small-scale growing modes. This extra term mimics the role played by the $z$ direction, which would destabilize and damp structures with radial lengths $\lesssim H$.</footnote>

<footnote n="6">⁶ See also Guan & Gammie (2009) and Fromang & Stone (2009) for other magnetic-field configurations.</footnote>
<figure><img src="image_15.png" /><figcaption>Figure 13. Space-time diagrams of the vertically and azimuthally averaged vertical component of the magnetic field ($B_z$) (top) and the flow vorticity ($\omega_z$) (bottom) in a toy model (Section 4.3) based upon our mean-field theory. A zonal-vorticity band emerges anticorrelated with the zonal magnetic field, consistent with the behaviour found in our numerical simulations (see Fig. 7).</figcaption></figure>

where the partial derivatives of $M$ are evaluated at $\langle B_z^0 \rangle$ and $\langle \omega_z^0 \rangle$. In principle, the Lorentz and Hall terms may be diffusive or antidiffusive, depending upon the local gradient of the Maxwell stress. In practice, $\partial M/\partial\langle\omega_z\rangle$ is generally non-negative for MRI-driven turbulence (e.g. Pessah, Chan & Psaltis 2008), and so the Lorentz force acts diffusively on the mean vorticity. As in the previous section, $\partial M/\partial\langle B_z \rangle$ is positive above some critical field strength and negative below it. We capture this physics by taking the Maxwell stress to be proportional to the channel growth rate, allowing for its dependence on vorticity by promoting $2A \rightarrow \langle \omega_z \rangle$ in the dispersion relation. We then solve equations (24) and (25) as before, with $\langle \omega_z^0 \rangle = 2A$. The result is graphically presented in Fig. 13. A zonal magnetic field is produced, with the flow vorticity being anticorrelated with the vertical magnetic field. This model captures all of the salient features of the LTS.

## 5 DISCUSSION

In this paper, we have described the saturation of the MRI in a plasma dominated by the Hall effect. This situation is of particular importance to protoplanetary discs, which are poorly ionized and in which one expects a particularly strong Hall effect near the disc mid-plane (e.g. Wardle 2007). A strong Hall effect is considered to be particularly important to the disc dynamics, as it can destabilize plasmas which would otherwise be stabilized by Ohmic losses and could therefore eliminate (or at least significantly reduce) purported ‘dead zones’ in such discs (Wardle & Salmeron 2012).

Using the 3D spectral code SNOOPY, we recovered the linear properties of the Hall-dominated MRI. We also demonstrated analytically that MRI ‘channel’ modes remain non-linear exact solutions notwithstanding the complicating influence of the Hall effect, and we explored their stability with respect to secondary parasitic modes. We subsequently confirmed the numerical results of SS02, which correspond to the weak-Hall regime ($\ell_H \sim 10^{-2}H$). By extending those authors’ numerical experiments into the Hall-dominated regime ($\ell_H \gtrsim 0.2H$), we arrived at our most interesting and perhaps most important finding: the presence of a completely new saturation mechanism at work in MRI-driven turbulence. We have shown that this new saturation mechanism is associated with a reduction in turbulent transport by at least two orders of magnitude compared to ‘classical’ magnetorotational turbulence.

This LTS is produced by a redistribution of the vertical magnetic flux into long-lived, narrow ($\sim H$), axisymmetric regions, outside of which the mean field averages almost to zero. We refer to these regions as zonal magnetic fields. Coincident with these zonal-field regions are zonal flows, which we have shown to result from (approximate) conservation of canonical vorticity. The accumulation of vertical magnetic flux into zonal-field regions increases the local magnetic tension enough to stabilize their internal flow; outside of these regions, very weak magnetorotational turbulence persists. In fact, the surrounding weak turbulence actively reinforces the integrity of the zonal-field structure by contributing a confining azimuthal electric field. This reinforcement is necessary in the face of resistive (both molecular and turbulent) diffusion.

To explain this behaviour, we have developed a mean-field model based on $y-z$ averages of the full set of resistive Hall-MHD equations in the shearing sheet. The key ingredient of this model is the introduction of the $xy$-component of the Maxwell stress $M_{xy}$ into the induction equation by the Hall electric field. This term introduces a new type of radial magnetic diffusion into the averaged induction equation, which can be either diffusive or antidiffusive depending on the sign of $d\langle M_{xy} \rangle/d\langle B_z \rangle$. Since both very weak and very strong magnetic fields quench turbulent transport by the MRI (the former through Ohmic losses and the latter by strong tension), we conclude that there will always be regions of MRI-driven turbulence where $d\langle M_{xy} \rangle/d\langle B_z \rangle > 0$ and therefore the Hall effect will act antidiffusively. When $\ell_H \gtrsim 0.2H$, this antidiffusive behaviour appears to overtake resistive and turbulent diffusion. We have constructed a toy model based upon our mean-field theory that reproduces the zonal-flow structure revealed by our simulations and captures all the salient features of the LTS.

Our results indicate that, despite being strongly linearly unstable, Hall-dominated accretion discs could undergo a global bifurcation that substantially reduces the amount of MRI-driven turbulent transport. This contradicts the recent suggestion by Wardle & Salmeron (2012) that the Hall effect could revive turbulent ‘dead’ zones long suspected to exist in the poorly ionized interiors of protoplanetary discs (e.g. Gammie 1996). Our results also bring into question numerous attempts to estimate the size of the dead zone by coupling chemical networks with MRI stability criteria, whether they be obtained through linear analyses or by naively extrapolating results from resistive MHD simulations into the Hall-dominated regime. Since strong Hall diffusion can non-linearly render magnetically active regions magnetorotationally ‘dead’, existing estimates for the extent of the active and dead zones may require revision.

Despite the simplicity of this result and the very general physical principles upon which it is based (e.g. field-strength-dependent turbulent transport, conservation of canonical vorticity), it is prudent to remind the reader that these results have been obtained using
a crude representation (an incompressible and unstratified shearing box) of an actual protoplanetary disc. Since the magnetic quenching responsible for the stabilization of zonal fields is dependent upon the vertical scaleheight, which is artificially set in our simulations to be the vertical size of the shearing box, vertical stratification could affect the dynamics, extent and magnetic-field strength of the zonal structures. That being said, stratified discs do exhibit a natural lengthscale that stabilizes MRI modes (Gammie & Balbus 1994; Latter, Fromang & Gressel 2010; Lesur, Ferreira & Ogilvie 2013) so that, in principle, zonal fields and flows should also be produced and sustained in stratified shearing boxes. Compressibility, ambipolar diffusion and the presence of dust grains introduce additional complications by altering the ionization fraction and the turbulent response of the disc in height-dependent ways.<fnref n="7" /> Simulations of layered discs (without the Hall effect or ambipolar diffusion) have even shown that magnetorotational turbulence in active surface layers could drive a Reynolds stress (Fleming & Stone 2003; Oishi & Mac Low 2009) and/or a large-scale Maxwell stress (Turner & Sano 2008) in the dead zone. We intend to explore the impact of such complexities on the robustness of our results in subsequent work.

We close by remarking that the presence of zonal magnetic fields and zonal flows is of great interest in its own right. While zonal flows have been previously found to naturally emerge in both ideal and resistive simulations of MRI-driven turbulence, those structures are generally weak in amplitude and are believed to be generated by random contributions of the Maxwell stress (Johansen, Youdin & Klahr 2009). The zonal structures we find in Hall-dominated magnetorotational turbulence, on the other hand, are strong in amplitude and are driven by a coherent Maxwell stress acting in concert with conservation of canonical vorticity. A natural suspicion is that these dominant zonal structures may act as particle-trapping sites, enabling fast planetesimal formation through a gravitationally unstable dust layer. The enticing link between Hall-induced turbulent bifurcation to an LTS, the occurrence of zonal fields and flows, and the formation of planetesimals will be the subject of a forthcoming publication (Lesur & Kunz, in preparation).

## ACKNOWLEDGEMENTS

Support for MWK was provided by NASA through Einstein Postdoctoral Fellowship Award Number PF1-120084, issued by the Chandra X-ray Observatory Center, which is operated by the Smithsonian Astrophysical Observatory for and on behalf of NASA under contract NAS8-03060. GL acknowledges support by the European Community via contract PCIG09-GA-2011-294110. This work was granted access to the HPC resources of IDRIS under allocation x2013042231 made by Grand Équipement National de Calcul Intensif (GENCI). Some of the computations presented in this paper were performed using the CIMENT infrastructure (https://ciment.ujf-grenoble.fr), which is supported by the Rhone-Alpes region (GRANT CPER07_13 CIRA: http://www.ci-ra.org). The Texas Advanced Computer Center at the University of Texas at Austin also provided HPC resources under grant number TG-AST090105. This work used the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by NSF grant OCI-1053575. The authors would like to thank Nuno Loureiro, Jake Simon and Jim Stone for useful conversations, as well as the expert referee for a prompt and constructive report.

## REFERENCES

* Armitage P. J., 2011, ARA&A, 49, 195

* Bai X.-N., 2011, ApJ, 739, 50

* Bai X.-N., Goodman J., 2009, ApJ, 701, 737

* Bai X.-N., Stone J. M., 2011, ApJ, 736, 144

* Balbus S. A., 2011, in Garcia P. J. V., ed., Physical Processes in Circumstellar Disks around Young Stars. University of Chicago Press, Chicago, p. 237

* Balbus S. A., Hawley J. F., 1991, ApJ, 376, 214

* Balbus S. A., Hawley J. F., 1998, Rev. Mod. Phys., 70, 1

* Balbus S. A., Henri P., 2008, ApJ, 674, 408

* Balbus S. A., Terquem C., 2001, ApJ, 552, 235

* Blaes O. M., Balbus S. A., 1994, ApJ, 421, 163

* Bodo G., Mignone A., Cattaneo F., Rossi P., Ferrari A., 2008, A&A, 487, 1

* Boyd J. P., 2001, Chebyshev and Fourier Spectral Methods. Dover Press, New York

* Canuto C., Hussaini M. Y., Quarteroni A., Zang T. A., 1988, Spectral Methods in Fluid Dynamics. Springer, Berlin

* Desch S. J., 2004, ApJ, 608, 509

* Falle S. A. E. G., 2003, MNRAS, 344, 1210

* Fleming T., Stone J. M., 2003, ApJ, 585, 908

* Fleming T. P., Stone J. M., Hawley J. F., 2000, ApJ, 530, 464

* Fromang S., Stone J. M., 2009, A&A, 507, 19

* Fromang S., Terquem C., Balbus S. A., 2002, MNRAS, 329, 18

* Fromang S., Papaloizou J., Lesur G., Heinemann T., 2007, A&A, 476, 1123

* Gammie C. F., 1996, ApJ, 457, 355

* Gammie C. F., Balbus S. A., 1994, MNRAS, 270, 138

* Goldreich P., Lynden-Bell D., 1965, MNRAS, 130, 125

* Goodman J., Xu G., 1994, ApJ, 432, 213

* Guan X., Gammie C. F., 2009, ApJ, 697, 1901

* Hartmann L., Calvet N., Gullbring E., D'Alessio P., 1998, ApJ, 495, 385

* Hawley J. F., Stone J. M., 1998, ApJ, 501, 758

* Hawley J. F., Gammie C. F., Balbus S. A., 1995, ApJ, 440, 742

* Hayashi C., 1981, Prog. Theor. Phys. Suppl., 70, 35

* Igea J., Glassgold A. E., 1999, ApJ, 518, 848

* Ilgner M., Nelson R. P., 2006, A&A, 445, 205

* Jin L., 1996, ApJ, 457, 798

* Johansen A., Youdin A., Klahr H., 2009, ApJ, 697, 1269

* Kunz M. W., 2008, MNRAS, 385, 1494 (K08)

* Kunz M. W., Balbus S. A., 2004, MNRAS, 348, 355

* Latter H. N., Lesaffre P., Balbus S. A., 2009, MNRAS, 394, 715

* Latter H. N., Fromang S., Gressel O., 2010, MNRAS, 406, 848

* Lesur G., Longaretti P.-Y., 2007, MNRAS, 378, 1471

* Lesur G., Longaretti P.-Y., 2009, A&A, 504, 309

* Lesur G., Ferreira J., Ogilvie G. I., 2013, A&A, 550, A61

* Longaretti P.-Y., Lesur G., 2010, A&A, 516, A51

* Masset F., 2000, A&AS, 141, 165

* Oishi J. S., Mac Low M.-M., 2009, ApJ, 704, 1239

* Pessah M. E., 2010, ApJ, 716, 1012

* Pessah M. E., Goodman J., 2009, ApJ, 698, L72

* Pessah M. E., Chan C.-K., Psaltis D., 2008, MNRAS, 383, 683

* Polygiannakis J. M., Moussas X., 2001, Plasma Phys. Control. Fusion, 43, 195

* Salmeron R., Wardle M., 2003, MNRAS, 345, 992

* Salmeron R., Wardle M., 2005, MNRAS, 361, 45

* Salmeron R., Wardle M., 2008, MNRAS, 388, 1223

* Sano T., Miyama S. M., 1999, ApJ, 515, 776

* Sano T., Stone J. M., 2002a, ApJ, 570, 314 (SS02)

* Sano T., Stone J. M., 2002b, ApJ, 577, 534 (SS02)

* Sano T., Miyama S. M., Umebayashi T., Nakano T., 2000, ApJ, 543, 486

* Shakura N. I., Sunyaev R. A., 1973, A&A, 24, 337

* Simon J. B., Hawley J. F., 2009, ApJ, 707, 833

* Simon J. B., Bai X.-N., Stone J. M., Armitage P. J., Beckwith K., 2013, ApJ, 764, 66

* Thomson W. (Lord Kelvin), 1869, Trans. R. Soc. Edinburgh, 25, 217

* Turner N. J., Sano T., 2008, ApJ, 679, L131

* Umebayashi T., Nakano T., 1990, MNRAS, 243, 103

* Umurhan O. M., Regev O., 2004, A&A, 427, 855

<footnote n="7">⁷ Preliminary simulations in unstratified boxes indicate that compressibility and strong ambipolar diffusion do not qualitatively change our results.</footnote>
Wardle M., 1999, MNRAS, 307, 849
Wardle M., 2007, Ap&SS, 311, 35
Wardle M., Salmeron R., 2012, MNRAS, 422, 2737

# APPENDIX A: CHANNEL STABILITY AND PARASITES

## A1 Formulation of the problem and method of solution

The fact that the Hall-MRI channel flows are exact non-linear solutions allows one to analytically examine their stability to parasitic modes. One simply considers the channel solution to be part of a time-dependent background state, upon which small-amplitude perturbations are applied. Unfortunately, this time dependence complicates matters somewhat, since one cannot Fourier decompose in time. We circumvent this difficulty by assuming that the growth rate $\sigma$ of the parasites is much greater than the growth rate $\gamma \sim \Omega$ of the 'background' channel solution (Goodman & Xu 1994). Since $\sigma \sim b\Omega$, this amounts to the assumption that the channel mode has grown to large amplitude ($b \gg 1$) and, from the perspective of the parasites, may be considered a stationary equilibrium. This ordering implies that the effects of rotation and shear on the parasites, as well as the uniform background magnetic field $B_0$, may be ignored. What results is an eighth-order boundary-value problem in $z$ with eigenvalue $\sigma$.

Following this prescription, we substitute

$$v = v_{ch} + \delta v, \quad B = B_{ch} + \delta B, \quad P = P_0 + \delta P$$

into equations (1)–(3), linearize in the perturbation amplitudes and search for Fourier modes with the space-time dependence $\delta \propto \exp(\sigma t + ik \cdot x)$. The perturbation wavevector

$$k = k (\hat{e}_x \sin \theta_k - \hat{e}_y \cos \theta_k)$$

is parallel to the channel magnetic (velocity) field when $\theta_k = \theta$ ($\theta_k = \phi \pm \pi/2$). The resulting set of equations is

$$\sigma \delta v = -i \mathbf{k} \cdot \mathbf{v}_{ch} \delta v - \delta v_z \frac{dv_{ch}}{dz} - \frac{1}{\rho} \left( i \mathbf{k} + \hat{e}_z \frac{d}{dz} \right) \delta \Pi + i \mathbf{k} \cdot \mathbf{B}_{ch} \frac{\delta \mathbf{B}}{4\pi\rho} + \frac{\delta B_z}{4\pi\rho} \frac{d\mathbf{B}_{ch}}{dz} + v \left( \frac{d^2}{dz^2} - k^2 \right) \delta v, \quad (A1)$$

$$\sigma \delta \mathbf{B} = -i \mathbf{k} \cdot \left( \mathbf{v}_{ch} - \frac{\mathbf{J}_{ch}}{en_e} \right) \delta \mathbf{B} - \left( \delta v_z - \frac{\delta \mathbf{J}_z}{en_e} \right) \frac{d\mathbf{B}_{ch}}{dz} \\ + i \mathbf{k} \cdot \mathbf{B}_{ch} \left( \delta \mathbf{v} - \frac{\delta \mathbf{J}}{en_e} \right) + \delta B_z \frac{d}{dz} \left( \mathbf{v}_{ch} - \frac{\mathbf{J}_{ch}}{en_e} \right) \\ + \eta \left( \frac{d^2}{dz^2} - k^2 \right) \delta \mathbf{B}, \quad (A2)$$

$$i \mathbf{k} \cdot \delta \mathbf{v} + \frac{d \delta v_z}{dz} = 0, \quad (A3)$$

where

$$\delta \Pi = \delta P + \frac{\mathbf{B}_{ch} \cdot \delta \mathbf{B}}{4\pi}$$

is the perturbation of the total pressure and $J_{ch}$ is the channel current density (equation 13).

The linear operator associated with equations (A1)–(A3) is periodic in $L \equiv 2\pi K^{-1}$, and we thus make the Floquet ansatz $\delta \propto \delta(z) \exp(ik_z z)$ where the function $\delta(z)$ is $L$-periodic. Following

Latter et al. (2009), we solve the resulting set of equations numerically via a pseudo-spectral technique: the $z$-domain is partitioned into $N (=256)$ grid points and the operator is discretized using Fourier cardinal functions (see Boyd 2001). This procedure leads to a $7N \times 7N$ generalized algebraic eigenvalue problem, which we solve using the QZ algorithm. The resulting (complex) eigenvalues are the growth rates $\sigma$ of the eigenmodes.

In the following sections, we give representative examples of three classes of parasitic eigenmodes. First, we take $k \perp B_{ch}$ and $\eta = 0$. The resulting 'kink' and 'kink-pinch' modes are essentially hydrodynamical disturbances of the channel jets, whose growth is afforded by the shear free energy stored in the channel mode. Secondly, we take $k \parallel B_{ch}$ and $\eta \neq 0$. The resulting 'pinch-tearing' mode extracts little to no free energy from the flow shear, and instead relies on Ohmic dissipation to facilitate the reconnection of neighbouring field lines of opposite polarity and thereby tap into the magnetic energy stored in the channel mode. Finally, we take $k \parallel B_{ch}$ and $\eta = 0$ and investigate whether the HSI can act as a parasitic mode by feeding off the shear free energy of the channel. In all cases, we set $\nu = 0$ and $\Omega/\omega_{H,0} = 17.4$, the latter of which corresponds to our fiducial simulation ZB1H1.

## A2 Kink and kink-pinch parasites

The kink (i.e. type I) and kink-pinch (i.e. type II) parasitic modes are described in detail in Goodman & Xu (1994, section 3.4) and Latter et al. (2009, section 2.3.3), to which we refer the reader. These parasites are most clearly identified when $\eta = 0$ and $k \perp B_{ch}$, conditions which simultaneously alleviate the stabilizing influence of magnetic tension and maximize the access to the shear free energy of the channel. In brief, the kink mode is a Kelvin-Helmholtz instability feeding upon the inflection points in the channel velocity profile. It is associated with $k_z = 0$ and, as such, is the dominant parasite for a two-stream channel. The kink-pinch mode is a hybrid mode, exhibiting both kink-like and pinch-like characteristics. Rather than exhibiting a phase velocity that matches the inflection points of the channel velocity profile, the kink-pinch mode appears to influence the flow most strongly near the points in $z$ where the channel magnetic field changes sign. This mode grows about an order of magnitude slower than the pure kink mode and, having $k_z \neq 0$, only afflicts channels with more than two streams.

In the top panel of Fig. A1, we plot the $\delta v_x$, $\delta v_z$ and $\delta \Pi$ components of the eigenfunction of a kink mode with $\theta = \phi = \pi/4$, $k/K = 0.5$, $\theta_k = -\pi/4$ and $k_z = 0$. The solid (dashed) lines denote the real (imaginary) parts. In the bottom panel, we display coloured isocontours of the real part of $v_x$ at $y = 0$ in the $(x, z)$ plane. The Hall effect does not significantly impact the structure of the mode since $k \cdot B_{ch} = 0$. Just as in ideal MHD, a pressure gradient develops across each jet and gives rise to a kink motion, deflecting the channels upwards or downwards. The Hall effect does, however, reduce the parasite growth rate from $\sigma/b\Omega = 0.1959$ to 0.008 616. This is because a strong Hall effect with $\Omega/\omega_{H,0} > 0$ substantially decreases the shear free energy $Kv_0$ stored in the channel mode (see equation 16).

We find a similar situation for the kink-pinch mode ($k_z = 0.5$), shown in Fig. A2. As in ideal MHD, the pressure perturbation $\delta \Pi$ changes sign at $Kz = \pi/2$ and kinks the jet centred there, while the vertical velocity perturbation $\delta v_z$ changes sign at $Kz = 3\pi/2$ and pinches the jet centred there. This gives rise to the alternate kinking and pinching of neighbouring jets, which can be observed
<figure><img src="image_16.png" /><figcaption>**Figure A1.** Top: the $\delta v_x$, $\delta v_z$ and $\delta \Pi$ components of the eigenfunction of a kink mode with $\theta = \phi = \pi/4$, $k/K = 0.5$, $\theta_k = -\pi/4$ and $k_z = 0$. The solid (dashed) lines denote the real (imaginary) parts. The total eigenfunction is normalized so that $\max|\delta v_x| = 1$. The growth rate $\sigma/b\Omega = 0.008616$. Bottom: coloured isocontours of the real part of $v_x$ at $y=0$ in the $(x, z)$ plane. The background is a four-stream Hall-MRI channel with jets centred at $Kz = n\pi/2$ with $n=1, 3, 5$ and $7$. The perturbation is normalized so that $\max|\delta v_x| = v_{ch}$.</figcaption></figure>

in the bottom panel of the figure. For the same reason as for the kink mode, the growth rate decreases from $\sigma/b\Omega = 0.0931 + 0.1319i$ to $0.004093 + 0.0058i$. The complex conjugate of this mode has the pressure perturbation kink at the $Kz = 3\pi/2$ jet and the vertical velocity perturbation pinch at the $Kz = \pi/2$ jet.

For both the kink and kink–pinch modes, the growth rates are much too small to be of relevance for the saturation of the Hall-dominated MRI. In order for $\sigma \sim \gamma$, the usual criterion for determining when secondary parasitic modes might overtake the primary channel mode and instigate a breakdown into 3D turbulence (Pessah & Goodman 2009), the channel must grow to amplitudes $b \sim 100$. By this stage, compressibility effects become important, a feature not captured in our parasite linear analysis nor in our non-linear simulations.

## A3 Pinch-tearing parasite

When resistivity is included, pinching motions are subject to tearing instabilities driven by magnetic reconnection. This not only modifies the kink–pinch mode by opening a pathway to extract the extra energy stored in the channel’s magnetic field, but also facilitates

<figure><img src="image_17.png" /><figcaption>**Figure A2.** Top: the $\delta v_x$, $\delta v_z$ and $\delta \Pi$ components of the eigenfunction of a kink–pinch mode with $\theta = \phi = \pi/4$, $k/K = 0.5$, $\theta_k = -\pi/4$ and $k_z = 0.5$. The solid (dashed) lines denote the real (imaginary) parts. The total eigenfunction is normalized so that $\max|\delta v_x| = 1$. The growth rate $\sigma/b\Omega = 0.004093 + 0.0058i$. Bottom: coloured isocontours of the real part of $v_x$ at $y=0$ in the $(x, z)$ plane. The background channel and normalization are as in Fig. A1. The entire pattern is moving to the left because $\sigma$ possesses a positive imaginary part.</figcaption></figure>

the growth of pinch-tearing modes (see section 2.3.4 of Latter et al. 2009). The latter manifest most readily when $\theta_k = \theta$, i.e. when the wavevector and the channel magnetic field are aligned. As $\theta_k$ deviates from $\theta$, the pinch-tearing mode bifurcates into two complex-conjugate modes, each one localized on one or the other jet. This is the situation we display in Fig. A3, which highlights the $\delta v_x$, $\delta B_x$ and $\delta B_z$ components of the eigenfunction of a pinch-tearing mode with $\theta_k = \pi/3$, $k/K = 0.4$, $k_z = 0$ and $b = 5$. The eigenmode localizes preferentially on the lower jet, with a growth rate $\sigma/b\Omega = 0.02868 - 0.03449i$; its complex conjugate localizes on the upper jet.

While the Hall effect with $\Omega/\omega_{H,0} > 0$ decreases the magnetic shear energy stored in the channel, it also generally increases the growth rate of the pinch-tearing mode. Indeed, this parasite grows faster than both the kink and kink–pinch modes described above by a factor of $\sim 10$. Nevertheless, its growth rate remains quite small; in simulations of the Hall-dominated MRI in which a two-stream $k_x = 0$ channel mode is deliberately excited, this parasite grows over such a long time-scale that the shearing of its wavevector (a feature not taken into consideration here) wraps the mode into an almost
<figure><img src="image_18.png" /><figcaption>**Figure A3.** Top: the $\delta v_x$, $\delta B_x$ and $\delta B_z$ components of the eigenfunction of a pinch-tearing mode with $\theta_k = \pi/3$, $k/K = 0.4$, $k_z = 0$ and $b = 5$. The channel orientation $\theta \simeq 0.255\pi$ is obtained by solving the dispersion relation for the fastest growing mode with $\beta$, $\Lambda_H$ and $\Lambda_\eta$ taken from our fiducial simulation ZB1H1. The solid (dashed) line denotes the real (imaginary) part. The total eigenfunction is normalized so that $\max|\delta v_x| = 1$. The growth rate $\sigma/b\Omega = 0.02868 - 0.03449i$. Bottom: coloured isocontours of the real parts of $v_x$ and $B_z$ at $y = 0$ in the $(x, z)$ plane. The background channel and normalization are as in Fig. A1. The entire pattern is moving to the right because $\sigma$ possesses a negative imaginary part.</figcaption></figure>

axisymmetric configuration with $k > 1$.<fnref n="8" /> In fact, the channel grows to such large amplitudes ($b > 300$) that the numerical time step drops precipitously and further time integration becomes untenable. In

simulations where a $k_x = 0$ channel mode is not exclusively excited by hand (e.g. all of those listed in Table 1), the linear phase is disrupted at $b \sim 10$, not by parasites, but rather by the interference of large-amplitude $k_x \neq 0$ channels.

### A4 HSI as a channel parasite

In Section 3.1.2 we employed the HSI as a test of our numerical algorithm. In the context of the channel mode, the background shear is provided by the sinusoidal variation of $v_{ch}$ and is strongest at the vertical locations satisfying $Kz = n\pi$ with $n = 0, \pm 1, \pm 2, \dots$. The magnitude of the channel magnetic field is also greatest at those heights. As a result, any local patch centred on $Kz = n\pi$ resembles a non-rotating shearing sheet with a vorticity-aligned magnetic field – precisely the setup considered by K08. While this analogy suggests that the HSI could act as a powerful channel parasite, this expectation is not borne out in our numerical simulations. Here we adapt the K08 calculation to determine why this is so.

The HSI is most easily examined in the limit of negligible resistivity, in which case the channel magnetic and velocity fields are mutually perpendicular ($\theta = \phi$; see equation 17). We can therefore erect an orthonormal coordinate system oriented with the channel mode: $\hat{e}_b = B_{ch}/B_{ch}$, $\hat{e}_v = v_{ch}/v_{ch}$ and $\hat{e}_z = \hat{e}_b \times \hat{e}_v$. In this geometry, wavevectors parallel to the channel magnetic field ($k = k\hat{e}_b$) have the greatest potential for growth. The $z$- and $v$-components of the linearized induction equation (A2) become<fnref n="9" />

$$ \sigma \delta B_z - b \cos Kz \frac{ck^2 B_0}{4\pi en_e} \delta B_v = ikB_0 b \cos Kz \delta v_z, \quad (A4) $$

$$ \begin{aligned} \sigma \delta B_v + b \cos Kz \left[ \frac{ck^2 B_0}{4\pi en_e} \left( 1 - \frac{1}{k^2} \frac{d^2}{dz^2} - \frac{K^2}{k^2} \right) - K v_0 \right] \delta B_z \\ = ikB_0 b \cos Kz \delta v_v. \end{aligned} \quad (A5) $$

It is clear from equation (A5) that the shear of the channel mode (represented by the final term in the brackets) uses $\delta B_z$ to generate $\delta B_v$. The Hall terms, on the other hand, generate $\delta B_z$ at the expense of $\delta B_v$. This effect is present even in the absence of shear and arises because the $v$-component of the perturbed electron velocity differs from the ion-neutral velocity by

$$ -\frac{\delta J_v}{en_e} = \frac{ick}{4\pi en_e} \left( 1 - \frac{1}{k^2} \frac{d^2}{dz^2} \right) \delta B_z. $$

The induced magnetic field is sheared further, and there is the potential for runaway.

It is a straightforward exercise to show from equations (A1), (A4) and (A5) that, whether $k \gg K$, $d/dz$ (the limit captured by the K08 analysis) or $d/dz \gg k$, $K$ (a WKBJ treatment), a necessary condition for instability is

$$ 1 < \frac{K v_0}{\omega_{H,0}}. \quad (A6) $$

Physically, this inequality states that the time required for an ion to execute one orbital gyration around a magnetic-field line must be longer (by a factor of $n_e/n$) than the time it takes for a magnetic perturbation to grow by shear. If this condition is not met, the ions are well coupled to the electrons (and thereby to the magnetic field), and we are left with simple linear-in-time growth due to shearing of

<footnote n="8">⁸ We refer the reader to appendix B of Latter et al. (2010) for a quantitative assessment of how shear retards parasite growth.</footnote>

<footnote n="9">⁹ These equations may be profitably compared with equations 46a and 46b of K08 after making the replacements $z \to x$ and $v \to y$.</footnote>
the magnetic-field perturbation by the channel flow. (This criterion is analogous to equation 43 of K08.)

What complicates matters beyond those investigated in K08 is that here the channel shear, which provides the free energy for growth, is itself a function of $\omega_{H,0}$. Plugging in our expression for the channel shear (equation 16 with $Rm_{eff} \to \infty$), our instability criterion (A6) for the HSI parasites becomes equation (18), precisely the stability criterion for the Hall-MRI channels themselves! In other words, if the Hall-MRI channels are active in the disc, then the HSI cannot act as a parasitic instability.

## APPENDIX B: NUMERICAL STABILITY IN HALL-MHD

Falle (2003) suggested that explicit schemes for numerically solving the equations of Hall-MHD are unconditionally unstable due to the existence of small-wavelength whistler waves. Although this conclusion is correct for the numerical schemes Falle (2003) considered, here we demonstrate that higher order time-explicit schemes, such as the one used in SNOOPY, are stable without the need for physical (e.g. Ohmic or ambipolar) or artificial (e.g. hyper-resistive) wave damping.

We start by considering the induction equation (equation 2) with the first (ideal) and third (Ohmic) terms on the right-hand side dropped. Decomposing the magnetic field into a fixed guide field $B_0$ and a small-amplitude fluctuation $\delta B(t) \exp(i\mathbf{k} \cdot \mathbf{x})$, we find that linear whistler waves are described by

$$ \frac{d\delta \mathbf{B}}{dt} = \frac{c\mathbf{k} \cdot \mathbf{B}_0}{4\pi e n_e} (\mathbf{k} \times \delta \mathbf{B}). \quad (B1) $$

In spectral codes such as SNOOPY the right-hand side of this equation is computed exactly using Fourier decomposition, and we adopt this scheme in what follows.

Without loss of generality we take the wavevector $\mathbf{k} = k\hat{\mathbf{e}}_z$ and magnetic-field perturbation $\delta \mathbf{B} = \delta B_x \hat{\mathbf{e}}_x + \delta B_y \hat{\mathbf{e}}_y$, ensuring $\mathbf{k} \cdot \delta \mathbf{B} = 0$. Equation (B1) can then be written as

$$ \frac{d\delta \mathbf{B}}{dt} = \mathbf{R} \delta \mathbf{B}, \quad \text{where} \quad \mathbf{R} = \frac{c k^2 B_{0,z}}{4\pi e n_e} \begin{pmatrix} 0 & -1 \\ 1 & 0 \end{pmatrix}. \quad (B2) $$

We integrate equation (B2) forward in time from $t^{(n)}$ to $t^{(n+1)}$ using an RK3 scheme similar to that used in SNOOPY. For a system of differential equations $\mathbf{y}' = f(\mathbf{y})$, this procedure reads

$$ \mathbf{q}_1 = f(\mathbf{y}^{(n)}) $$

$$ \mathbf{q}_2 = f\left(\mathbf{y}^{(n)} + \frac{h}{2}\mathbf{q}_1\right) $$

$$ \mathbf{q}_3 = f(\mathbf{y}^{(n)} - h\mathbf{q}_1 + 2h\mathbf{q}_2) $$

$$ \mathbf{y}^{(n+1)} = \mathbf{y}^{(n)} + \frac{h}{6}(\mathbf{q}_1 + 4\mathbf{q}_2 + \mathbf{q}_3), $$

where $h \equiv t^{(n+1)} - t^{(n)}$. Applying this algorithm to equation (B2), we find

$$ \delta \mathbf{B}^{(n+1)} = \mathbf{Q} \delta \mathbf{B}^{(n)} \quad (B3) $$

for

$$ \mathbf{Q} = \begin{pmatrix} 1 - \frac{1}{2}\epsilon^2 & -\epsilon + \frac{1}{6}\epsilon^3 \\ \epsilon - \frac{1}{6}\epsilon^3 & 1 - \frac{1}{2}\epsilon^2 \end{pmatrix} \quad \text{and} \quad \epsilon \equiv h \frac{c k^2 B_{0,z}}{4\pi e n_e}. $$

Note that the matrix $\mathbf{Q}$ is a third-order expansion of the formal solution $\delta \mathbf{B}^{(n+1)} = \exp(h\mathbf{R}) \delta \mathbf{B}^{(n)}$. Extensions to higher order are straightforward.

Stability is guaranteed if the eigenvalues of $\mathbf{Q}$,

$$ \lambda_{\pm} = 1 - \frac{\epsilon^2}{2} \mp i \left( \epsilon - \frac{\epsilon^3}{6} \right), $$

satisfy the inequality $|\lambda_{\pm}| < 1$. The numerical scheme is therefore stable provided $\epsilon < \sqrt{3}$; SNOOPY uses $\epsilon = 1.5$. It can easily be shown by this approach that similar schemes of first or second order in time, such as the ones considered by Falle (2003), are unconditionally unstable. The fourth-order Runge-Kutta scheme is stable for $\epsilon < 2\sqrt{2}$.

In conclusion, the third-order explicit time integrator employed in SNOOPY guarantees that linear whistler waves are stable, without the need for additional diffusion terms.

This paper has been typeset from a TeX/LA TeX file prepared by the author.