2.5.1. Sources of ionisation

As we focus on the outer part of PPDs ($R > 1\ \mathrm {AU}$), the gas is mostly cold with $T < 300\ \mathrm {K}$. This implies that thermal ionisation (owing to collision between molecules) is inefficient, and one has to rely on non-thermal ionisation processes. Here, we consider the following effects with their associated ionisation rate $\zeta$:

1. (i) X-ray ionisation owing to bremsstrahlung emission from an isothermal $T=5\ \mathrm {keV}$ corona localised around the central protostar (Igea & Glassgold Reference Igea and Glassgold1999; Bai & Goodman Reference Bai and Goodman2009, see their equation (21));

2. (ii) CR ionisation with $\zeta _\mathrm {CR}=\zeta _{\mathrm {CR},0} \exp (-\varSigma /96\ \mathrm {g}\,\mathrm {cm}^{-2})\,\mathrm {s}^{-1}$ (e.g. Umebayashi & Nakano Reference Umebayashi and Nakano1981) and $\zeta _{\mathrm {CR},0}=10^{-17}\ \mathrm {s}^{-1}$, corresponding to the interstellar value;

3. (iii) radioactive decay with $\zeta _\mathrm {rad}=10^{-19}\ \mathrm {s}^{-1}$ (Umebayashi & Nakano Reference Umebayashi and Nakano2009).

The amount of ionisation due to CRs is highly disputed. Some authors have proposed that owing to the wind coming from the young star, CRs are magnetically mirrored from the PPD, resulting in a significantly reduced ionisation rate due to CRs ($\zeta _{\mathrm {CR},0}\sim 10^{-20}\ \mathrm {s}^{-1}$, Cleeves, Adams & Bergin Reference Cleeves, Adams and Bergin2013). In contrast, it has been proposed that CRs could be accelerated in shocks produced in the protostellar jet by a Fermi process. This could result in ionisation rates as high as $\zeta _{\mathrm {CR},0}\sim 10^{-13}\ \mathrm {s}^{-1}$ (Padovani et al. Reference Padovani, Ivlev, Galli and Caselli2018). Observations of TW Hya tend to suggest a low ionisation rate owing to CRs ($\zeta _{\mathrm {CR},0}\lesssim 10^{-19}\ \mathrm {s}^{-1}$, Cleeves et al. Reference Cleeves, Bergin, Qi, Adams and Öberg2015), though this is still highly model dependent. Owing to these uncertainties, some authors (e.g.Ilgner & Nelson Reference Ilgner and Nelson2006) have simply omitted CR ionisation and consider only X-rays as the main source of ionisation. These difference and uncertainties in the treatment of the ionisation rate have to be kept in mind when comparing the results of different research groups.

We show in figure 6 the resulting ionisation rate following the disc structure presented in § 2.1. We find that CRs are shielded only in the innermost parts of the disc, where the column density goes above $100\ \mathrm {g}\,\mathrm {cm}^{-2}$. Most of the disc midplane up to $z\sim h$ has $\zeta \simeq \zeta _{\mathrm {CR},0}$, indicating that CRs are indeed the dominant source of ionisation in this region. Above $z\sim h$, X-rays start to penetrate the disc and the ionisation rate rises.

Figure 6. Ionisation rate $\log (\zeta )$ ($\mathrm {s}^{-1}$) as a function of radius and altitude (in disc scale height) resulting from X-rays, CRs and radioactive decay.

2.5.2. A simple chemical model

To illustrate the typical ionisation fractions expected in PPDs, we follow Oppenheimer & Dalgarno (Reference Oppenheimer and Dalgarno1974), Fromang, Terquem & Balbus (Reference Fromang, Terquem and Balbus2002) and Ilgner & Nelson (Reference Ilgner and Nelson2006) defining the following reaction network and rates with free electrons, neutral molecules $m$, molecular ions ${m^+}$ and metal atoms ${M}$:

(2.7)\begin{equation} \left.\begin{array}{c@{\quad}c@{}} {m} + \textrm{ionising radiation} \rightarrow {m}^{+} + {e}^{-} & \zeta,\\ {m}^{+} + e^{-} \rightarrow {m} & \delta,\\ {M}^{+} + e^{-} \rightarrow {M} & \delta_r, \\ {m}^{+} + {M} \rightarrow {m} + {M}^{+} & \delta_t, \end{array}\right\} \end{equation}

where $\zeta$ is the ionisation rate, $\delta$ is the dissociative recombination rate for molecular ions, $\delta _r$ the radiative recombination rate for metal atoms, and $\delta _t$ the rate of charge transfer from molecular ions to metal atoms. Following Fromang et al. (Reference Fromang, Terquem and Balbus2002), we take

(2.8)\begin{equation} \left.\begin{gathered} \delta_r=3\times 10^{-11} T^{-1/2}\ \mathrm{cm}^3\,\mathrm{s}^{-1},\\ \delta=3\times 10^{-6} T^{-1/2}\ \mathrm{cm}^3\,\mathrm{s}^{-1},\\ \delta_t=3\times 10^{-9}\ \mathrm{cm}^3\,\mathrm{s}^{-1}. \end{gathered}\right\} \end{equation}

In the absence of metals and grains, the rate equations admit a simple solution in steady state:

(2.9)\begin{equation} \xi=\sqrt{\frac{\zeta}{\delta n_n}}. \end{equation}

In the opposite metal-dominated limit, still without grains, one obtains

(2.10)\begin{equation} \xi=\sqrt{\frac{\zeta}{\delta_r n_n}}. \end{equation}

As $\delta _r\ll \delta$, one clearly sees that the absence of metals leads to a dramatic decrease in the ionisation fraction (Fromang et al. Reference Fromang, Terquem and Balbus2002).

2.5.3. Typical ionisation fraction profile

Grain-free case, Metal-free case: Combining (2.9) with the ionisation rate in § 2.5.1, one can obtain the ionisation fraction in the disc. However, this ionisation fraction depends not only on the disc chemistry one assumes but also on the disc structure. A lot of theoretical work has focused on the minimum mass solar nebula (MMSN) model, which assumes $\varSigma =1700\, R_{\mathrm {AU}}^{-3/2}\ \mathrm {g}\,\mathrm {cm}^{-2}$ (Wardle Reference Wardle2007; Bai & Stone Reference Bai and Stone2013b; Lesur, Kunz & Fromang Reference Lesur, Kunz and Fromang2014). This makes the disc much denser in the inner part, resulting in a stronger shielding of CRs and a lower ionisation fraction than less-dense discs. As an illustration, we show in figure 7 the resulting ionisation fraction with the disc structure presented in § 2.1 and with a MMSN disc model.

Figure 7. Ionisation fraction $\log (\xi )$ as a function of position in (a) our disc model (§ 2.1) and (b) in a MMSN. Note the difference in ionisation fraction close to the disc midplane for $R<10\ \mathrm {AU}$.

We observe that the lowest ionisation fraction reaches $10^{-14}$ in the MMSN case or $10^{-13}$ in our disc model. The lowest ionisation fractions are reached in the innermost parts of the disc, where the recombination is the fastest and $\textrm {CRs}+\text {X-rays}$ are efficiently shielded. The ionisation fraction progressively increases when X-rays start to penetrate, until one reach ionisation fractions as high as $10^{-6}$ at a few scale heights. Note that the differences between these models are only significant for $R < 10\ \mathrm {AU}$ because the column densities between the MMSN and our disc model are similar above this radius.

Inclusion of grains and metals: As demonstrated by Elmegreen (Reference Elmegreen1979) and Umebayashi & Nakano (Reference Umebayashi and Nakano1980) in the context of molecular clouds, and later applied to PPDs (Sano et al. Reference Sano, Miyama, Umebayashi and Nakano2000; Ilgner & Nelson Reference Ilgner and Nelson2006; Wardle Reference Wardle2007), grains tend to accelerate the recombination of electrons by removing them from the gas phase, resulting in a lower global ionisation fraction, which we have ignored here. To illustrate the effect of grains, let us add the following reactions to our simplified reaction network:

(2.11)\begin{equation} \left.\begin{gathered} \mathrm{grain} + {m}^{+} \rightarrow \mathrm{grain}^{+} + {m},\\ \mathrm{grain}^{-} + {m}^{+} \rightarrow \mathrm{grain} + {m},\\ \mathrm{grain} + {e}^{-} \rightarrow \mathrm{grain}^{-},\\ \mathrm{grain}^{+} + {e}^{-} \rightarrow \mathrm{grain},\\ \mathrm{grain} + {M}^{+} \rightarrow \mathrm{grain}^{+} + M,\\ \mathrm{grain}^{-} + {M}^{+} \rightarrow \mathrm{grain} + M,\\ \mathrm{grain}^{+} + \mathrm{grain}^{-} \rightarrow \mathrm{grain} + \mathrm{grain}. \end{gathered}\right\} \end{equation}

This reaction network only considers singly charged grains, whereas it is well known that grains can have many charges (Ilgner Reference Ilgner2012). We chose this approach to illustrate in the simplest model the effect of grains on the ionisation fraction, and later on the diffusivities because the abundance of multiply charged grains is usually lower than that of singly charged grains for $z < h$ (Wardle Reference Wardle2007).

The rates for these reactions are computed by assuming each species collides at its thermal velocity with a spherical grain of radius $a$ (see § 2.4 for more details). We assume a fixed sticking probability of electrons on grains, which corresponds to the probability of bouncing back from a grain.Footnote ³

The resulting ionisation fraction owing to electrons and charged grains are presented in figure 8. We observe essentially two trends. First, the smallest ionisation fraction is found when grains are present, whereas the highest ionisation fractions correspond to grain-free metal-rich cases, with variations owing to this composition effect of the order of three orders of magnitude. Second, the ionisation fraction increases with increasing radius. This effect is not only because the ionisation rate increases, but also because the recombination rate decreases owing to lower densities. Let us finally point out that when grains are present, they can become the main charge carrier, as is the case at $R=5\ \textrm {AU.}$

Figure 8. Ionisation fraction $\xi$ for three different compositions: row 1 (a,b), no grains, no metals; row 2 (c,d), no grains with $[{M}]=10^{-8}$; row 3 (e,f), with $a=0.1\ \mathrm {\mu }\mathrm {m}$ grains and metal atoms. The first column corresponds to $R=5\ \mathrm {AU}$ and the second column to $R=50\ \mathrm {AU}$.

In the following, we use the value $\xi =10^{-13}$ to evaluate several plasma parameters, keeping in mind this corresponds to a lower bound in our disc model.

3.3.1. Dynamical equation for the centre of mass

The set of equations (3.3) and (3.4) can, in principle, be solved simultaneously (O'Keeffe & Downes Reference O'Keeffe and Downes2014). However, it is numerically expensive because the numerical time steps are usually limited by $\tau _s$, which is much smaller than the timescales of interest (as described previously). Note, however, that there are situations where the multifluids approach cannot be avoided, such as when the timescale to reach the ionisation/recombination equilibrium becomes of the order of the timescales of interest (e.g. Ilgner & Nelson Reference Ilgner and Nelson2008), or when the neutral density is so low that the collision timescale $\tau _s$ becomes of the order of the timescales of interest, which can occur well above the disc in the early phases of star formation, when X-rays and UV are not yet produced by the central body.

However, if one focuses on disc dynamics and its immediate environment once the central star is formed, the single-fluid approximation is a perfectly reasonable approximation, as multi-fluid approaches tend to confirm (Rodgers-Lee, Ray & Downes Reference Rodgers-Lee, Ray and Downes2016). For this reason, I will focus here on the single-fluid approximation. To derive this single-fluid approximation, let us consider the dynamical equations for the centre of mass of the fluid, defining the total mass density $\rho =\sum _jn_j m_j$, the flow velocity $\boldsymbol {v}=\sum _j n_j m_j\boldsymbol {v}_j/\rho$ and the drift speed for each species $\boldsymbol {w}_j=\boldsymbol {v}_j-\boldsymbol {v}$ we sum equations (3.3) and (3.4) to obtain

(3.15)\begin{equation} \left.\begin{gathered} \dfrac{\partial \rho}{\partial t} +\boldsymbol{\nabla}\boldsymbol{\cdot } \rho \boldsymbol{v}=0,\\ \dfrac{\partial \rho \boldsymbol{v}}{\partial t}+\boldsymbol{\nabla}\boldsymbol{\cdot } (\rho\boldsymbol{v}\otimes\boldsymbol{v})=\boldsymbol{\nabla}\boldsymbol{\cdot} \left(\sum _j n_jm_j\boldsymbol{w}_j\otimes\boldsymbol{w}_j\right)- \boldsymbol{\nabla}P+\boldsymbol{f}+\dfrac{\boldsymbol{J}\boldsymbol{\times} \boldsymbol{B}}{c}+\sum _j{n_j q_j}\boldsymbol{E}, \end{gathered}\right\} \end{equation}

where we have introduced the total pressure and force $P$ and $\boldsymbol {f}$ as well as the total current $\boldsymbol {J}=\sum _j n_j q_j v_j$. These equations are exact. However, they do not correspond to the usual dynamical equations one is used to, and it is important to understand why each extra term can be neglected.

The first term on the right-hand side corresponds to the transport of momentum by the drift velocity. Physically, it can be interpreted as a diffusion of momentum owing to drifting particles. It can be neglected, provided that drift velocities are small, i.e. that $w_j < L \varOmega \sqrt {\rho /\rho _j}$ where $L$ is the typical length scale of interest and $\varOmega$ is the typical frequency.Footnote ⁶ The presence of the density ratio ensures that even for drift velocities comparable with $L\varOmega$, this term is negligible.

We also have a term involving the total charge of the flow $\sum _j n_j q_j$. As shown previously, this term is negligible provided that the timescale of interest is sufficiently long to recover charge neutrality, which is usually the case. We can therefore drop this term altogether to obtain the usual single-fluid equations

(3.16)\begin{equation} \left.\begin{gathered} \dfrac{\partial \rho}{\partial t} +\boldsymbol{\nabla}\boldsymbol{\cdot }\rho \boldsymbol{v}=0 ,\\ \dfrac{\partial \rho \boldsymbol{v}}{\partial t}+\boldsymbol{\nabla}\boldsymbol{\cdot }( \rho\boldsymbol{v}\otimes\boldsymbol{v})=-\boldsymbol{\nabla}P+\boldsymbol{f}+ \dfrac{\boldsymbol{J}\boldsymbol{\times}\boldsymbol{B}}{c}. \end{gathered}\right\} \end{equation}

3.3.2. Ohm's law

In the equation of motion for the centre of mass, we have left aside the fact that additional equations were required to obtain $\boldsymbol {B}$ and $\boldsymbol {J}$. Indeed, Maxwell's equations give us

(3.17)\begin{equation} \left.\begin{gathered} \dfrac{\partial \boldsymbol{B}}{\partial t}=-c\boldsymbol{\nabla}\boldsymbol{\times}\boldsymbol{E},\\ \boldsymbol{J}=\dfrac{c}{4{\rm \pi}}\boldsymbol{\nabla} \boldsymbol{\times} \boldsymbol{B}. \end{gathered}\right\} \end{equation}

The remaining unknown is, therefore, the electric field. Owing to our assumption of electro-neutrality, we cannot use Gauss's law to compute the electric field (since under our scheme of approximation, the total charge density is zero). However, we can use the dynamical equation for charged species to deduce the electric field that is consistent with quasi-neutrality.

Let us start with (3.3), and let us separate the velocity into a velocity for the centre of mass, and the drift velocity for species $j$:

(3.18)\begin{align} \rho_j \frac{\mathrm{d} \boldsymbol{w}_j}{\mathrm{d}t}&=- \boldsymbol{\nabla}P_j+\boldsymbol{f}_j+n_jq_j\left(\frac{\boldsymbol{w}_j\boldsymbol{\times}\boldsymbol{B}}{c}+\boldsymbol{E}_\mathrm{cm}\right)+\boldsymbol{R}_j\nonumber\\ & \quad -\rho_j\left[\boldsymbol{w}_j\boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{v}+\boldsymbol{v}\boldsymbol{\cdot} \boldsymbol{\nabla}\boldsymbol{w}_j+\frac{\boldsymbol{F}_{\mathrm{cm}}}{\rho}\right] \end{align}

where we have defined the electric field in the centre of mass frame $\boldsymbol {E}_\mathrm {cm}\equiv \boldsymbol {E}+\boldsymbol {v}\boldsymbol {\times } \boldsymbol {B}/c$ and the forces on the centre of mass $\boldsymbol {F}_\mathrm {cm}=-\boldsymbol {\nabla }P+\boldsymbol {f}+\boldsymbol {J}\boldsymbol {\times } \boldsymbol {B}/c$. Several terms can be neglected here assuming that the stopping time for the species is short compared to the other timescales of the problem.

1. (i) $\mathrm {d}_t \boldsymbol {w}_j$ can be neglected provided that $\varOmega \ll \tau _s^{-1}$ (i.e. this assumption is identical to the quasi-neutrality assumption discussed previously). In other words, the inertia of charged particles is negligible and they instantaneously reach their asymptotic velocity.

2. (ii) Similarly, the inertial term (second line) and external forces $\boldsymbol {f}_j$ can be neglected because they modify the impulsion on timescales long compared with $\tau _s$.

3. (iii) $\boldsymbol {\nabla } P_j\sim \rho _j c_{s,j}^2/\varLambda$ is negligible provided that $c_{s,j}\lesssim \varOmega \varLambda$.

The equations of motion for charged particles in the frame of the centre of mass therefore read

(3.19)\begin{equation} q_j\left(\frac{\boldsymbol{w}_j\boldsymbol{\times} \boldsymbol{B}}{c}+\boldsymbol{E}_{\mathrm{cm}} \right)-\gamma_{jn}m_j\rho \boldsymbol{w}_j=0, \end{equation}

where we have assumed that dominant collisions were due to neutrals. This is usually recast as

(3.20)\begin{equation} \boldsymbol{w}_j-\mu_j \boldsymbol{w}_j\boldsymbol{\times}\hat{\boldsymbol{b}}=\frac{c\mu_j}{B}\boldsymbol{E}_{\mathrm{cm}}, \end{equation}

where $\hat {\boldsymbol {b}}$ is a unit vector parallel to $\boldsymbol {B}$ and

(3.21)\begin{equation} \mu_j\equiv \frac{q_jB}{\gamma_{jn}\rho m_j c}, \end{equation}

is the Hall parameter (Wardle & Ng Reference Wardle and Ng1999). Equation (3.20) can be solved for $\boldsymbol {w}_j$, which gives the asymptotic velocity

(3.22)\begin{equation} \left.\begin{gathered} \boldsymbol{w}_{j,\parallel}=\dfrac{c\mu_j}{B}\boldsymbol{E}_{\mathrm{cm},\parallel},\\ \boldsymbol{w}_{j,\perp}=\dfrac{c\mu_j}{B(1+\mu_j^2)}\left[ \boldsymbol{E}_{\mathrm{cm},\perp}+\mu_j\boldsymbol{E}_{\mathrm{cm},\perp} \boldsymbol{\times}\hat{\boldsymbol{b}}\right]. \end{gathered}\right\} \end{equation}

We eventually obtain an expression closing our set of equations by relating the drift velocities to the current in the flow $\boldsymbol {J}=\sum _j n_j q_j \boldsymbol {w}_j$ and assuming quasi-neutrality $\sum _j n_j q_j=0$:

(3.23)\begin{equation} \left.\begin{gathered} \boldsymbol{J}_\parallel=\dfrac{c}{B} \left(\sum _j q_j n_j \mu_j\right) \boldsymbol{E}_\parallel, \\ \boldsymbol{J}_{\perp}=\dfrac{c}{B}\left(\sum _j \dfrac{q_jn_j\mu_j}{1+\mu_j^2}\right)\boldsymbol{E}_{\mathrm{cm},\perp}+ \dfrac{c}{B}\left(\sum _j \dfrac{q_j n_j}{1+\mu_j^2}\right) \hat{\boldsymbol{b}}\boldsymbol{\times}\boldsymbol{E}_{\mathrm{cm},\perp}. \end{gathered}\right\} \end{equation}

These expressions constitute the base of Ohm's law. We can identify three conductivity tensors, the Ohmic, Hall and Petersen conductivity tensors,

(3.24)\begin{equation} \left.\begin{gathered} \sigma_O= \dfrac{c}{B}\sum _j q_j n_j \mu_j,\\ \sigma_H=\dfrac{c}{B}\sum _j \dfrac{q_j n_j}{1+\mu_j^2},\\ \sigma_P=\dfrac{c}{B}\sum _j\dfrac{q_jn_j\mu_j}{1+\mu_j^2}, \end{gathered}\right\} \end{equation}

defined so that Ohm's law can be written in the more familiar form

(3.25)\begin{equation} \boldsymbol{J}=\sigma_\parallel \boldsymbol{E}_{\mathrm{cm},\parallel}+\sigma_H\hat{\boldsymbol{b}}\boldsymbol{\times} \boldsymbol{E}_{\mathrm{cm},\perp}+\sigma_P\boldsymbol{E}_{\mathrm{cm},\perp}. \end{equation}

This relation can be inverted one final time to obtain the electric field in the observer frame and write the induction equation as

(3.26)\begin{equation} \frac{\partial \boldsymbol{B}}{\partial t}=\boldsymbol{\nabla}\boldsymbol{\times} \left(\boldsymbol{v}\boldsymbol{\times}\boldsymbol{B}\right)-\boldsymbol{\nabla}\boldsymbol{\times} \left(\eta_O \boldsymbol{\nabla }\boldsymbol{\times}\boldsymbol{B}+ \eta_H(\boldsymbol{\nabla}\boldsymbol{\times}\boldsymbol{B}) \boldsymbol{\times}\hat{\boldsymbol{b}}+\eta_A(\boldsymbol{\nabla}\boldsymbol{\times} \boldsymbol{B})_\perp\right), \end{equation}

where the magnetic diffusivities are defined as

(3.27)\begin{gather} \eta_O=\frac{c^2}{4{\rm \pi}}\frac{1}{\sigma_O}, \end{gather}

(3.28)\begin{gather}\eta_H=\frac{c^2}{4{\rm \pi}}\frac{\sigma_H}{\sigma_H^2+\sigma_P^2}, \end{gather}

(3.29)\begin{gather}\eta_A= \frac{c^2}{4{\rm \pi}}\left(\frac{\sigma_P}{\sigma_H^2+\sigma_P^2}- \frac{1}{\sigma_O}\right), \end{gather}

where the subscripts $O$, $H$ and $A$ denotes Ohmic, Hall and ambipolar.

3.4.1. Simplified case of two charged species

In the simplest case of a plasma made of two singly charged species ($+$) and ($-$), we obtain the following simplified expressions from (3.27), (3.28) and (3.29):

(3.30)\begin{equation} \left.\begin{gathered} \eta_O=\dfrac{cB}{4{\rm \pi} e n_+}\left(\dfrac{1}{\mu_+-\mu_-}\right),\\ \eta_H=\dfrac{cB}{4{\rm \pi} e n_+}\left(\dfrac{\mu_++\mu_-}{\mu_--\mu_+}\right),\\ \eta_A=\dfrac{cB}{4{\rm \pi} e n_+}\left(\dfrac{\mu_+\mu_-}{\mu_--\mu_+}\right). \end{gathered}\right\} \end{equation}

First, all of these coefficients are proportional to $n_+^{-1}$, i.e. are inversely proportional to the ionisation fraction. Second, because $\mu _j\propto B$, we find that $\eta _O$ does not depend on $B$, whereas $\eta _H\propto B$ and $\eta _A\propto B^2$. Finally, we find that $\eta _H$ may have either sign. If $|\mu _-| > |\mu _+|$, we find $\eta _H > 0$, and $\eta _H < 0$ otherwise. As $\mu$ is essentially a measure of the collisionality and mass of the charge carrier, it indicates that the sign of the Hall effect depends on the nature of the charge carriers. In the case where the positive and negative charge carriers have identical masses and $\gamma$ so that $\mu _-=-\mu _+$, the Hall effect vanishes.

In the case of an electron–ion plasma, we have $|\mu _e|\simeq m_i/m_e|\mu _i|\gg |\mu _i|$. Hence, $\eta _O\propto |\mu _e|^{-1}$ and $\eta _A\propto |\mu _i|$, which justifies the usual statement that Ohmic diffusion is a result of electron–neutral collisions and ambipolar diffusion to ion–neutral collisions. We also have $\eta _H=|\mu _e|\eta _O$ and $\eta _A=|\mu _e\mu _i|\eta _O$. Hence, we can distinguish three regimes depending on the Hall parameter of the ion and electrons:

1. (i) $1 < \mu _i < |\mu _e|$ in which case $\eta _A > \eta _H > \eta _O$ and the regime is predominantly ambipolar;

2. (ii) $\mu _i < 1 < |\mu _e|$ in which case $\eta _H > (\eta _A,\eta _O)$ known as the Hall regime;

3. (iii) $\mu _i < |\mu _e| < 1$ where $\eta _O > \eta _H > \eta _A$ and which is dominated by Ohmic diffusion.

This allows us to delimit the Ohmic, Hall and ambipolar regime as a function of the neutral density and the field intensity (figure 9). As can be seen, the midplane of PPDs is expected to lie mostly in the Hall regime and possibly in the ambipolar regime in the outer-most parts of the disc.

Figure 9. Non-ideal regimes as a function of the neutral density and magnetic field intensity, computed for an electron–ion plasma and assuming $T<100\ \mathrm {K}$. The blue and green lines correspond to the typical values of a PPD midplane, for various plasma $\beta$ parameters.

A word of caution though: the physical nature of the Hall effect is different from the Ohmic and ambipolar counterparts (the Hall effect is dispersive, but not diffusive because $\boldsymbol {J}\boldsymbol {\times } \boldsymbol {B}\boldsymbol {\cdot }\boldsymbol {B}=\textbf {{0}}$). Being in the Ohmic- or ambipolar-dominated regime does not automatically imply that the Hall effect is dynamically unimportant.

Finally, we obtain the usual expressions for the diffusivities in the electron–ion case:

(3.31)\begin{align}\eta_O&=\frac{c^2\gamma_{en}m_n m_e}{4{\rm \pi} {e}^2}\frac{1}{\xi}\nonumber\\ &=2.3\times 10^{16} \left(\frac{\xi}{10^{-13}}\right)^{-1}\,\mathrm{cm}^2\,\mathrm{s}^{-1}\end{align}

(3.32)\begin{align}\eta_H&=\frac{cB}{4{\rm \pi} e n_e}\nonumber\\ &=5.0\times 10^{17}\left(\frac{\xi}{10^{-13}}\right)^{-1} \left(\frac{B}{1\ \mathrm{G}}\right)\left(\frac{n_n}{10^{14}\ \mathrm{cm}^{-3}}\right)^{-1}\,\mathrm{cm}^2\,\mathrm{s}^{-1} \end{align}

(3.33)\begin{align} \eta_A&=\frac{B^2}{4{\rm \pi}\gamma_{in}\rho \rho_i}\nonumber\\ &=1.6\times 10^{16} \left(\frac{\xi}{10^{-13}}\right)^{-1} \left(\frac{B}{1\ \mathrm{G}}\right)^2 \left(\frac{n_n}{10^{14}\ \mathrm{cm}^{-3}}\right)^{-2}\,\mathrm{cm}^2\,\mathrm{s}^{-1}. \end{align}

These values can be compared with diffusivities of everyday material such as iron ($\eta =8\times 10^2\ \mathrm {cm}^{2}\,\mathrm {s}^{-1}$), demineralised water ($\eta =1.4\times 10^{15}\ \mathrm {cm}^2\,\mathrm {s}^{-1}$) and dry air ($\eta =1.6\times 10^{24}\ \mathrm {cm}^2\,\mathrm {s}^{-1}$). Even though one might wrongfully conclude from this that MHD effects are irrelevant, the time scales (${\sim }1~\mathrm {year}$) and length scales (${\sim }1~\mathrm {AU}$) are also much larger than conventional everyday experiments. This illustrates the fact that dimensionless numbers should be compared and not dimensional quantities. As we show, one obtains magnetic Reynolds numbers of $O(1)$, which put these flows in a regime comparable with liquid sodium experiments on Earth.

3.4.2. Dimensionless numbers and application to disc models

It is customary to define dimensionless numbers in association with non-ideal effects in order to quantify their relative importance in the induction equation. First, one can define Elsasser numbers

(3.34)\begin{equation} \varLambda_{O,H,A}\equiv \frac{V_A^2}{\varOmega \eta_{O,H,A}}, \end{equation}

where $V_A$ is the Alfvén speed and $\varOmega$ is the rotation rate of the system. Note, however, that $\varLambda _O\propto B^2$ and $\varLambda _H\propto B$, which makes these numbers less useful when it comes to predicting the saturation of MHD instabilities because $B$ is a priori unknown. It is therefore useful to define two additional dimensionless numbers, the magnetic Reynolds number and the Hall Lundquist number

(3.35)\begin{equation} \left.\begin{gathered} {Rm}\equiv\dfrac{\varOmega H^2}{\eta_O},\\ \mathcal{L}_{H}\equiv \dfrac{V_A H}{\eta_H}. \end{gathered}\right\} \end{equation}

These two numbers do not depend on the field strength (at least in the two-species plasma case), and they turn out to be excellent saturation predictors in the non-linear regime of the MRI. We show in figure 10 the dimensionless numbers resulting from our grain-free metal-free ionisation model. As it can be seen, ${Rm} < 10^3$ only in the innermost regions of the disc. This is the region which was historically defined as the ‘dead zone’ (Gammie Reference Gammie1996). In addition, we find $10^{-1}<\mathcal {L}_{H}<10$ in most of the disc midplane with a sharp increase at the disc surface whereas $\varLambda _A\simeq 1$ in most of the disc.

Figure 10. Magnetic Reynolds number ${Rm}$ (a), Hall Lundquist number $\mathcal {L}_H$ (b) and ambipolar Elsasser number $\varLambda _A$ (c) in our disc model (§ 2.1), using a simple ion–electron approximation with a metal-free chemistry. White values are ${ > }10^3$.

3.4.3. The role of grains

When it comes to the conductivity of PPDs, grains play essentially two roles.

1. (i) By capturing free electrons, they become predominantly negatively charged, and they increase the recombination rate with ions thanks to their large cross-section and reaction rates at the grain surface. The end product is generally a reduced ionisation fraction, possibly by several orders of magnitude (see § 2.5.3).

2. (ii) Owing to their high inertia, charged grains enter into the conductivity tensor as a very low Hall parameter species. In this case, the scaling laws obtained for the two species case do not hold anymore. The abundance of charged grains, therefore, changes the diffusion regime in which the system lies.

To illustrate how the diffusivity depends on the presence of grains, we present in figure 11 an example of diffusivity computation in plasmas of different compositions (these compositions are identical to those discussed in § 2.5.3). As can be seen, the addition of $0.1\ \mathrm {\mu }\mathrm {m}$ size grains to the system has a dramatic effect on the diffusivities: all of the dimensionless numbers decrease by several orders of magnitude close to the midplane, whereas ambipolar diffusion becomes stronger than Hall in the case with grains. Overall, ambipolar diffusion increases by $10^4$ whereas Hall and Ohmic increase by $10^2$ compared with the fiducial metal-free case. The Hall effect also changes sign at the disc surface. This arises because of the presence of negatively charged grains, which contribute to the disc conductivity tensor by reducing the ‘effective’ Hall parameter of negative charge carriers ($\text {electrons}+\text {grains}-$). In the end, the Hall conductivity becomes dominated by ions when they become sufficiently abundant: at the disc surface.

Figure 11. Dimensionless diffusivity for three different compositions: row 1 (a,b), no grains, no metals (identical to figure 10); row 2 (c,d), no grains with $[M]=10^{-8}$; row 3 (e,f), with $a=0.1\ \mathrm {\mu }\mathrm {m}$ grains and metal atoms. The first column corresponds to $R=5\ \mathrm {AU}$ and the second column $R=50\ \mathrm {AU.}$ Dashed lines correspond to negative diffusivities for the Hall effect.

In the grain-free case, the presence of metal atoms tends to decrease the diffusivities by typically one to two orders of magnitude. However, let us point out that the effect of metal atoms disappears once grains are sufficiently abundant (Ilgner & Nelson Reference Ilgner and Nelson2006). Our grain-free metal-rich model is, therefore, a best-case scenario for the ionisation fraction and the diffusivities.

The simplified grain model we have used is by no mean the final answer to this question. However, it demonstrates the strong effect of grains on the dynamics of the plasma. The intensity of this effect also depends on the grain size and grain abundance, a lower abundance or larger grain size leading to a smaller effect (Ilgner & Nelson Reference Ilgner and Nelson2006; Salmeron & Wardle Reference Salmeron and Wardle2008). Here, we have purposely chosen very small grains with an interstellar abundance to illustrate a worst-case scenario for the ionisation fraction. Finally, if one assumes polycyclic aromatic hydrocarbons (PAHs) are present in the gas phase, they then behave as very small grains, capturing all of the floating electrons and also affecting significantly the amplitude of non-ideal effects (Bai Reference Bai2011b).

5.2.1. Introduction

The shearing box model formally comes in two forms, the ‘large’ shearing box (LSB, also known as the stratified shearing box) model which is essentially a numerical equivalent of Hill's approximation in a finite size numerical domain, and the ‘small’ shearing box (SSB, also known as the unstratified shearing box), which is a local approximation in Hill's approximation. In essence, the SSB model assumes that one zooms on a region close to the disc midplane so that the box size $L\ll H$. In this case, vertical gravity can be neglected in (5.13) and because one expects $w\sim \varOmega L$, the flow is strongly subsonic so that an incompressible approximationFootnote ⁹ can be made in place of (5.12). The asymptotic of these two models are discussed in details by Umurhan & Regev (Reference Umurhan and Regev2004). Here, we mostly use the LSB, keeping in mind that an incompressible approximation is possible to study the basic effects of the shear.

5.2.2. Boundary conditions

The shearing box model makes use of Hill's approximation in the simplest possible numerical setup: a periodic box. Because Hill's approximation is local, the shearing box model is also local and should satisfy the asymptotic rules shown above. In particular, for a box of size $L$, we should have $L\sim H\ll R_0$. However, the presence of the comoving derivative makes things a bit more difficult. The fields $\rho$, $w$, etc. are advected everywhere at the azimuthal velocity $-q\varOmega x\boldsymbol {e}_y$. It is, therefore, physically inconsistent to assume $Q(-L_x/2,y,z)=Q(L_x/2,y,z)$ for any quantity $Q$ as one would do with $x$ periodic boundary conditions in a box of size $L_x$. In order to take into account the constant shear in the boundary conditions, one therefore enforces periodic boundary condition in a Lagrangian view known as ‘shearing sheet’: $Q(-L_x/2,y,z)=Q(L_x/2,y+q\varOmega L_x t,z)$, which can be represented graphically as in figure 14.

Figure 14. Radial boundary conditions in a shearing box. The background shear is represented in blue. At $t=0$ the boundary conditions are strictly periodic (red). At $t > 0$, the periodic boundary conditions are shifted in time (green), according to the advection by the mean flow.

In the azimuthal direction, periodic boundary conditions are the most natural choice. In the vertical direction, however, no obvious choice comes to mind. Depending on the problem at hand, one can use outflow, free-slip or even periodic boundary conditions.

5.2.3. Stress and accretion measurement

One of the goals of the shearing box approach is to measure directly the turbulent stress in the global dynamical equations (4.21). One therefore needs to quantify

(5.17)\begin{equation} W_{xy}=\overline{\rho w_x w_y}-\frac{\overline{B_x B_y}}{4{\rm \pi}} \end{equation}

which is Hill's equivalent of the radial stress term in (4.21). Despite a non-zero radial stress, the shearing box model does not exhibit any accretion owing to this component. This is because the shearing box model is a local asymptotic expansion. The terms leading to accretion are higher-order terms in $H/R_0$, which have been neglected, and which would break the in–out symmetry of Hill's potential.

5.2.4. Conserved quantities

The shearing box model is very useful in the sense that it allows one to have a tight control on the conserved quantities of the flow: vertical and radial magnetic flux, momentum and energy are all conserved thanks to the relatively simple boundary conditions used in the horizontal direction. We thus define

(5.18)\begin{equation} \langle \cdot \rangle\equiv \iiint\mathrm{d}x\,\mathrm{d}y\,\mathrm{d}z. \end{equation}

Mass, momentum and flux conservation equations eventually read

(5.19)\begin{gather} \partial_t\langle \rho\rangle +\left[\rho v_z\right]_{z=z_b}=0 \end{gather}

(5.20)\begin{align} \partial_t\langle \rho \boldsymbol{w}\rangle + \left[\rho \boldsymbol{w}w_z+\left(P+\frac{B^2}{8{\rm \pi}}\right) \boldsymbol{e}_z-\frac{1}{4{\rm \pi}}\boldsymbol{B}B_z\right]_{z=z_b}&=- 2\varOmega_0 \boldsymbol{e}_z\boldsymbol{\times}\langle \rho \boldsymbol{w}\rangle\nonumber\\ &\quad+q\varOmega_0\langle \rho w_x\rangle \boldsymbol{e}_y- \varOmega_0^2 \langle \rho z\rangle \boldsymbol{e}_z \end{align}

(5.21)\begin{gather} \partial_t\langle \boldsymbol{B}\rangle + \left[w_z\boldsymbol{B}-B_z\boldsymbol{w}\right]_{z=z_b} =-q\varOmega_0\langle B_x\rangle\boldsymbol{e}_y +\textrm{non-ideal terms}. \end{gather}

Interestingly, the momentum equation exhibits source terms connected to the effective potential and the Coriolis force. The vertical magnetic flux is found to be exactly conserved whereas the horizontal flux is not necessarily conserved. Horizontal flux can escape the box through the vertical boundary or, alternatively, can be amplified thanks to the $\varOmega$-effect that shows up as a source term of toroidal field. These simple conservation laws allows one to control very carefully the box physics. Moreover, because the vertical flux is conserved, one can classify shearing-box setup as a function of the average vertical flux. It is also possible to do this for the toroidal component in the $y$ direction, though conservation can be violated by flux escape at the boundary.

Energetics of the shearing box may be found by dotting the momentum equation with $\boldsymbol {w}$ and the induction equation with $\boldsymbol {B}$. One eventually obtains an equation for the mechanical energy (kinetic plus magnetic) in the box

(5.22)\begin{equation} \partial_t \mathcal{E}_\mathrm{Mech}+\boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{\mathcal{F}}_\mathrm{Mech}= P\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{w}+ \boldsymbol{E}_\mathrm{NI}\boldsymbol{\cdot} \boldsymbol{J}+q\varOmega_0 \left(\rho w_xw_y-\frac{B_xB_y}{4{\rm \pi}}\right)-\varOmega_0^2\rho w_z z, \end{equation}

where

(5.23)\begin{equation} \left.\begin{gathered} \mathcal{E}_\mathrm{Mech}=\dfrac{1}{2}\rho w^2+\dfrac{1}{8{\rm \pi}}B^2,\\ \boldsymbol{{\mathcal{F}}}_\mathrm{Mech}=\dfrac{1}{2}\rho w^2 \boldsymbol{w}+\dfrac{1}{4{\rm \pi}}\left(B^2\boldsymbol{w}-(\boldsymbol{w}\boldsymbol{\cdot} \boldsymbol{B})\boldsymbol{B}-c\boldsymbol{E}_\mathrm{NI}\boldsymbol{\times}\boldsymbol{B} \right)+P\boldsymbol{w}. \end{gathered}\right\} \end{equation}

This set of equations is the local equivalent of (4.17). Several comments can be made on the energetics. First, we find an energy flux made of three contributions, kinetic energy, Poynting flux (split into advective, magnetic and non-ideal contributions) and a pressure term. Second, we find several source/sink terms.

1. (i) $PdV$ work. This term is usually small when thermal effects are negligible, but can become important in thermally driven winds for example.

2. (ii) Non-ideal effects. For Ohmic and ambipolar diffusion, $\boldsymbol {E}_{\mathrm {NI}}\sim -\boldsymbol {J}$, hence for these two terms, $\boldsymbol {E}_{\mathrm {Ohm/ambipolar}}\boldsymbol {\cdot }\boldsymbol {J} < 0$. Meanwhile, the Hall effect has no contribution to this term because $\boldsymbol {J}\boldsymbol {\times }\boldsymbol {B}\boldsymbol {\cdot }\boldsymbol {J}=0$. Overall, the non-ideal term therefore appears as a sink of mechanical energy, as expected. Energy is dissipated into heat.

3. (iii) Radial stress. We recover the radial stress term of the global angular momentum equation (4.21), which appears as a source term here.

4. (iv) Potential work. The last term denotes the work done by the vertical gravity on the gas, which can become important in ejecting disc models. In this case, this term is negative.

5. (v) Viscous terms (not shown). When the flow is turbulent, viscous terms can become important at the small scale. In this case, they show up as an additional sink term in the mechanical energy equation.

Overall, unless strong thermal effects are present, the only source of mechanical energy in this system is the radial stress term, which is therefore positive definite, as already pointed out in the global version of this equation. The dynamics of the disc will therefore be dictated by how this source term is balanced by the various loss/flux terms in the energetics.

Let us finally point out that in a shearing box, the energy flux $\boldsymbol {{\mathcal {F}}}_\mathrm {Mech}$ is periodic/shear periodic. If we average the energy equation, we find that the only relevant flux component is that escaping through the vertical boundaries, as for the mass and momentum fluxes.

6.1.1. Linear hydrodynamic stability

Let us now consider a particle evolving in Hill's effective potential under the influence of the effective gravity and the Coriolis force. The particle is initially at rest at $(x,y)=0$. Magnetic fields are neglected in this first approach. The equation of motion for the fluid particle may be written

(6.1)\begin{gather} \frac{\mathrm{d}^2x}{\textrm{d}t^2}=2q\varOmega_0^2 x+2\varOmega_0\frac{\textrm{d}y}{\mathrm{d}t}, \end{gather}

(6.2)\begin{gather}\frac{\mathrm{d}^2y}{\mathrm{d}t^2}=-2\varOmega_0\frac{\mathrm{d}x}{\mathrm{d}t}, \end{gather}

(6.3)\begin{gather}\frac{\mathrm{d}^2z}{\mathrm{d}t^2}=-\varOmega_0^2 z. \end{gather}

We first note that the vertical and horizontal equations of motion are separable. In the vertical direction, it describes oscillations of the fluid particle around the midplane at frequency $\varOmega _0$.

In the horizontal direction, the equations describes epicycles. To show it, let us first integrate (6.2):

(6.4)\begin{equation} \mathcal{L}=\frac{\mathrm{d}y}{\mathrm{d}t}+2\varOmega_0 x, \end{equation}

where $\mathcal {L}$ is the local angular momentum of the particle, as defined in (5.16). Our particle being initially in equilibrium at $(x,y)=0$, it has $\mathcal {L}=0$ and we can write the radial equation of motion as

(6.5)\begin{equation} \frac{\textrm{d}^2x}{\textrm{d}t^2}=-2\varOmega_0^2(2-q) x \end{equation}

hence, our effective gravitational potential, which was initially unstable ($\partial _x\psi _{\mathrm {eff},\mathrm {Hill}} < 0$), is stabilised thanks to the conservation of angular momentum, provided that $q < 2$ ($q=3/2$ for astrophysical discs). The oscillations described by this particle have a frequency

(6.6)\begin{equation} \omega^2=2\varOmega_0^2(2-q)\equiv \kappa^2. \end{equation}

This characteristic frequency is named epicyclic frequency. In the particular case of a Keplerian disc ($q=3/2$), we find $\kappa ^2=\varOmega ^2$, i.e. the epicyclic frequency coincides with the orbital frequency. As a result, orbits are closed, a well-known property of the two-body problem (e.g. figure 15).

Figure 15. Epicyclic oscillations of a fluid particle orbiting a point mass resulting in an closed elliptic orbit.

This shows that at the linear level, pure Keplerian flows are stable. This does not mean that non-linear (or subcritical) instabilities cannot exist in these flows, given that the shear is a natural reservoir of free energy to trigger instabilities. Such non-linear instabilities are well known to develop in non-rotating sheared flows and in pipe flows. This question of subcritical instabilities is at the origin of many experiments, numerical simulations and theoretical developments. Today, there seems to be a consensus on the fact that pure Keplerian flows appear to be stable for Reynolds numbers up to a few million. However, thermal effects, such as heating and vertical stratification, are known to affect this picture, leading to thermally driven linear and non-linear instabilities. This is a whole field of research, which is not covered here. The interested reader may consult Fromang & Lesur (Reference Fromang and Lesur2019) for a more detailed overview of this topic.

6.1.2. Linear MHD stability

As we have shown, PPDs are hydrodynamically stable at the linear level. In MHD, however, things start to become a bit more interesting (see also Balbus & Hawley Reference Balbus and Hawley1998 for a similar treatment). Let us embed our disc in an external and constant magnetic field $\boldsymbol {B}_0$, which we assume is vertical. Assuming we still consider infinitesimal displacements around the equilibrium position of the fluid particles, the velocities are infinitely small, and the induction equation for magnetic fluctuations $\delta \boldsymbol {b}$ reads

(6.7)\begin{equation} \frac{\partial \delta \boldsymbol{b}}{\partial t}=B_0\frac{\partial \boldsymbol{v}}{\partial z}. \end{equation}

Clearly, the stability will now depend on how we move the particles with respect to each other. Let us consider a set of particles initially at $(x,y)=0$ and let us perturb these particles with a vertical harmonic perturbation

(6.8)\begin{equation} \boldsymbol{x}=\boldsymbol{x}_0\exp(\textrm{i}kz), \end{equation}

the resulting magnetic perturbation can be obtained by integrating the induction equation with respect to time:

(6.9)\begin{equation} \delta \boldsymbol{b}=\textrm{i}k B_0\boldsymbol{x}. \end{equation}

In order to model how the field affects the dynamics, we have to include the Lorentz force $\boldsymbol {F}_L$ in the equation of motion. In the horizontal direction, only the magnetic tension term $\boldsymbol {B}\boldsymbol {\cdot }\boldsymbol {\nabla } \boldsymbol {B}$ appears, so we have

(6.10)\begin{align} \frac{\boldsymbol{F}_L}{\rho}&= \frac{\boldsymbol{B}_0\boldsymbol{\cdot}\boldsymbol{\nabla} \delta \boldsymbol{b}}{4{\rm \pi}\rho}\nonumber\\ &=-\frac{k^2B_0^2}{4{\rm \pi}\rho}\boldsymbol{x}\nonumber\\ &=-V_A^2k^2\boldsymbol{x}, \end{align}

where $V_A$ is the Alfvén speed. The horizontal equations of motion are therefore reduced to

(6.11)\begin{equation} \left.\begin{gathered} \dfrac{\textrm{d}^2x}{\textrm{d}t^2}=2q\varOmega_0^2 x+2\varOmega_0\dfrac{\textrm{d}y}{\textrm{d}t}-V_A^2k^2x,\\ \dfrac{\textrm{d}^2y}{\textrm{d}t^2}=-2\varOmega_0\dfrac{\textrm{d}x}{\textrm{d}t}-V_A^2k^2y, \end{gathered}\right\} \end{equation}

where it is clear that the magnetic forces are acting as a restoring force (hence, the usual representation of a spring for the Lorentz force). Note also that angular momentum conservation is now broken by the azimuthal tension force. It is this effect that leads to an instability.

To show this, let us assume $\boldsymbol {x}=\boldsymbol {x}\exp (\sigma t)$. The equations of motion lead to the following eigenvalue problem

(6.12)\begin{equation} \left.\begin{gathered} (\sigma^2+V_A^2k^2)x=2q\varOmega_0^2x+2\varOmega_0\sigma y,\\ (\sigma^2+V_A^2k^2)y=-2\varOmega_0\sigma x, \end{gathered}\right\} \end{equation}

which allows us to obtain the dispersion relation

(6.13)\begin{equation} (\sigma^2+V_A^2k^2)^2-2q\varOmega_0^2(\sigma^2+V_A^2k^2)+4\varOmega_0^2\sigma^2=0, \end{equation}

where we recover epicyclic oscillations when $V_A=0$ with $\sigma ^2=-2\varOmega _0^2(2-q)=-\kappa ^2$ and pure Alfvénic oscillations when $\varOmega _0=0$ with $\sigma ^2=-V_A^2k^2$. Expanding this dispersion relation leads to

(6.14)\begin{equation} \sigma^4+\sigma^2\left(\kappa^2+2V_A^2k^2\right)+ V_A^2k^2\left(V_A^2k^2-2q\varOmega_0^2\right)=0. \end{equation}

This dispersion relation describes a linear instability when $\sigma ^2$ is positive, i.e. when

(6.15)\begin{equation} V_A^2k^2-2q\varOmega_0^2 < 0. \end{equation}

This instability is the MRI. It appears when the magnetic tension force is not too strong, as suggested by (6.15). It is possible to solve the full dispersion analytically to obtain the eigenvalues (see figure 16). When $V_Ak<\sqrt {2q}\varOmega _0$, positive eigenvalues are found,which are the signature of the MRI. The maximum growth rates are obtained for $V_Ak=\sqrt {2q}\varOmega _0/2$ with $\sigma _\mathrm {max}=q\varOmega /2$. Above the limit (6.15), the unstable branch becomes a stable Alfvén wave, which shows that the MRI is mostly an Alfvénic perturbation. In addition to this pair of branches, we find a pair of epicyclic modes which are stable for all $kV_A$ (figure 16), and have a non-zero frequency for $V_A=0$. If compressibility is added, the degeneracy between slow magnetosonic and Alfvén waves is lifted. In this case, it can be shown that the MRI arises from the slow magnetosonic mode (Balbus & Hawley Reference Balbus and Hawley1992).

Figure 16. Real part (black) and imaginary part (red dashed line) of the solutions of (6.14) with $q=3/2$. The MRI appears for weak enough fields $V_Ak < \sqrt {3}\varOmega _0$.

The physical interpretation of the MRI is straightforward: consider two fluid particles attached to a vertical field line and assume we slightly move these particles radially. First, they will start an epicyclic motion and drift azimuthally (figure 17). As they drift away, the azimuthal magnetic tension will act as a spring bringing back the particles together, slowing down the inner particle and accelerating the outer particle. This results in a loss of angular momentum for the inner particle, which falls further down, and conversely for the outer particle. This mechanism can only work if the radial magnetic tension is sufficiently weak; otherwise, the particles return to their initial point resulting in an Alfvénic oscillation. It is this radial component of the Lorentz force that is the stabilising agent of the MRI.

Figure 17. Physical representation of the MRI mechanism (see the text).

6.4.1. Historical background

The linear MRI in the non-ideal MHD regime has been explored by many researchers. After the discovery of the MRI in the disc context by Balbus & Hawley (Reference Balbus and Hawley1991), it was soon realised that PPDs, but also discs in cataclysmic variables could be in the non-ideal MHD regime, casting doubts on the applicability of this instability to these objects. Blaes & Balbus (Reference Blaes and Balbus1994) were the first to consider this problem, by working out the ambipolar-dominated MRI in the two-fluid approximation. Ohmic diffusion was first considered by Jin (Reference Jin1996) in unstratified models and Sano & Miyama (Reference Sano and Miyama1999) in stratified discs, which led to the dead zone model of PPDs (Gammie Reference Gammie1996). Later, Wardle (Reference Wardle1999) considered the three non-ideal MHD effects in simple axial geometry whereas Desch Reference Desch2004 considered also oblique modes and toroidal fields. Balbus & Terquem (Reference Balbus and Terquem2001) isolated the physics of the Hall-MRI and Kunz (Reference Kunz2008) demonstrated that one of the Hall-MRI branches was actually a new instability: the HSI. Finally, ambipolar diffusion was revisited in the single-fluid approximation by Kunz & Balbus (Reference Kunz and Balbus2004), demonstrating the existence of oblique ambipolar modes and their origin.

In the following, we revisit and discuss each of these effects with a unified notation and geometry. Note, however, that our approach and dispersion relation is formally identical to that of Desch (Reference Desch2004).

6.4.2. Linearised equations

As in the cases described previously, we start from (5.12)–(5.15) in which we assume the disc is threaded by a mean field having a vertical and azimuthal components: $\boldsymbol {B}=B_{0,y}\boldsymbol {e}_y+B_{0,z}\boldsymbol {e}_z$. In the following, we consider small axisymmetricFootnote ¹² perturbations of the equilibrium. We seek solutions of the form $\boldsymbol {w}=\boldsymbol {u}\exp [\boldsymbol {k}\boldsymbol {\cdot }\boldsymbol {x}+\sigma t]$ and $\boldsymbol {B}=\boldsymbol {B}_0+\boldsymbol {b}\exp [\boldsymbol {k}\boldsymbol {\cdot } \boldsymbol {x}+\sigma t]$, where $\boldsymbol {k}=k_x\boldsymbol {e}_x+k_z\boldsymbol {e}_z$ is the wavenumber and $\sigma$ is the linear growth rate of the instability. We moreover neglect vertical stratification and vertical gravity and assume the flow is incompressible, which implies that we consider the SSB approximation. As we show later, this is enough to capture most of the physics relevant to the problem. We explore in § 6.4.6 the effect of the vertical stratification on linear modes.

Under these assumptions, the linearised equations read

(6.33)\begin{equation} \sigma\boldsymbol{u}=-\textrm{i}\boldsymbol{k} \varPi+\textrm{i}\frac{\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{B}_0}{ 4{\rm \pi} \rho} \boldsymbol{b}+2\varOmega u_y \boldsymbol{e}_x-(2-q) \varOmega u_x\boldsymbol{e}_y, \end{equation}

(6.34)\begin{align} \sigma \boldsymbol{b}&=\textrm{i}(\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{B}_0) \boldsymbol{u} -q\varOmega b_x \boldsymbol{e}_{\boldsymbol{y}}-\eta_O k^2\boldsymbol{b}+\eta_H\left(\boldsymbol{k}\boldsymbol{\cdot} \hat{\boldsymbol{B}}_0\right) \boldsymbol{k}\boldsymbol{\times} \boldsymbol{b}, \nonumber\\ &\quad -\eta_A\left[\left(\boldsymbol{k}\boldsymbol{\cdot} \hat{\boldsymbol{B}}_0\right)^2\boldsymbol{b}-\left(\boldsymbol{b}\boldsymbol{\cdot} \hat{\boldsymbol{B}}_0\right)\left(\left[\boldsymbol{k}\boldsymbol{\cdot} \hat{\boldsymbol{B}}_0\right]\boldsymbol{k}-k^2\hat{\boldsymbol{B}}_0\right)\right], \end{align}

(6.35)\begin{gather} \boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{u}=0, \end{gather}

(6.36)\begin{gather} \boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{b}=0. \end{gather}

The expression of the non-ideal terms can be interpreted in the following way. First, Ohmic diffusion acts as a pure linear damping operator, as expected. The Hall term is proportional to $\boldsymbol {k}\boldsymbol {\times } \boldsymbol {b}$, which means it rotates the magnetic perturbation around the $\boldsymbol {k}$ direction, keeping its norm constant. Note that the direction of rotation is given by the sign of $\eta _H$, which shows that the handedness given by the Hall effect is connected directly to the microphysics of the plasma (see § 3.4.1). Finally, ambipolar diffusion involves an anisotropic diffusion term that we discuss in the following.

The solenoidal conditions can be used to eliminate $u_z$ and $b_z$ from the equations in favour of the horizontal components, leading to a fourth-order problem

(6.37)\begin{gather} \sigma u_x=\textrm{i}\frac{\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{B}_0}{4{\rm \pi} \rho} b_x+2\varOmega \frac{k_z^2}{k^2}u_y, \end{gather}

(6.38)\begin{gather} \sigma u_y=\textrm{i}\frac{\boldsymbol{k}\boldsymbol{\cdot} \boldsymbol{B}_0}{4{\rm \pi} \rho} b_y-(2-q)\varOmega u_x, \end{gather}

(6.39)\begin{gather} \sigma b_x=\textrm{i}(\boldsymbol{k}\boldsymbol{\cdot} \boldsymbol{B}_0) u_x-\eta_O k^2b_x-\eta_H\left(\boldsymbol{k}\boldsymbol{\cdot} \hat{\boldsymbol{B}}_0\right)k_z b_y-\frac{\eta_A}{B_0^2} \left(k^2B_{0,z}^2b_x-k_xB_{0,y}\left(\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{B}_0\right) b_y\right), \end{gather}

(6.40)\begin{align} \sigma b_y&=\textrm{i}(\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{B}_0) u_y-q\varOmega b_x-\eta_O k^2b_y+\eta_H\left(\boldsymbol{k}\boldsymbol{\cdot} \hat{\boldsymbol{B}}_0\right)\frac{k^2}{k_z}b_x\nonumber\\ &\quad-\frac{\eta_A}{B_0^2}\left(\left[(\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{B}_0) ^2+k^2B_{0,y}^2\right]b_y-\frac{k^2}{k_z^2}k_xB_{0,y}(\boldsymbol{k}\boldsymbol{\cdot} \boldsymbol{B}_0)b_x\right). \end{align}

The role played by ambipolar diffusion here is a bit more self-explanatory. We observe that the diagonal terms in the second pair of equations are always negative definite, hence ambipolar diffusion is really acting as a diffusion term on the diagonal components with an amplitude controlled by the magnitude, but also the orientation of $\boldsymbol {B}_0$. However, there are also off-diagonal terms proportional to $k_xB_{0,y}$. As we show, these terms can lead to oblique unstable modes (see also Kunz & Balbus Reference Kunz and Balbus2004).

We next solve the set of equations for $\sigma$ described previously and look for unstable eigenvalues. We follow Pandey & Wardle (Reference Pandey and Wardle2012) and first obtain an equation for the velocity fluctuations

(6.41)\begin{equation} \boldsymbol{u}=\frac{\textrm{i}\boldsymbol{k}\boldsymbol{\cdot} \boldsymbol{B}_0}{4{\rm \pi}\rho(\sigma^2+\kappa^2)} \left(\begin{array}{cc} \sigma & 2\varOmega \\ -(2-q)\varOmega & \sigma \end{array}\right) \boldsymbol{b}, \end{equation}

so that the induction equation becomes, in matrix form,

(6.42)\begin{align} &\left[\left(\begin{array}{cc} \sigma + \eta_O k^2+\tau_{{A}} k^2 V_{{A}z}^2 & \ell_{{H}} \omega_{{A}} k_z-\tau_{{A}} k_xV_{{A}y}\omega_{{A}}\\ q\varOmega - \dfrac{k^2}{k_z^2}\left(\ell_{{H}}\omega_{{A}} k_z+\tau_{{A}} k_xV_{{A}y}\omega_{{A}}\right) & \sigma+\eta_O k^2+\tau_{{A}}(\omega_{{A}}^2+k^2V_{{A}y}^2) \end{array}\right)\right.\nonumber\\ &\left.\quad +\frac{\omega_{{A}}^2}{\sigma^2+\dfrac{k_z^2}{k^2}\kappa^2} \left(\begin{array}{cc} \sigma & 2\varOmega\dfrac{k_z^2}{k^2} \\ -(2-q)\varOmega & \sigma \end{array}\right) \right]\boldsymbol{b}=0\end{align}

where we have introduced the Alvén speed $V_{A}=B_0/(4{\rm \pi} \rho )^{1/2}$, the Alfvén frequency $\omega _{{A}}\equiv \boldsymbol {k}\boldsymbol {\cdot } {\boldsymbol {V}}_{{A}}$, the Hall length $\ell _{{H}}\equiv \eta _H/V_{A}$, the ambipolar time $\tau _{{A}}\equiv \eta _A/V_{A}^2$ and the epicyclic frequency $\kappa ^2\equiv 2\varOmega ^2(2-q)$.

After a long but straightforward calculation, one eventually obtains the dispersion relation which can be written

(6.43)\begin{equation} \sigma^4+\mathcal{C}_3\sigma^3+\mathcal{C}_2\sigma^2+\mathcal{C}_1\sigma+\mathcal{C}_0=0 \end{equation}

with

(6.44)\begin{gather} \mathcal{C}_3=2\eta_Ok^2+\tau_{{A}}(k^2V_{{A}}^2+\omega_{{A}}^2), \end{gather}

(6.45)\begin{gather} \mathcal{C}_2=\frac{k_z^2}{k^2}\kappa^2+2\omega_{{A}}^2+\eta_O^2k^4+\tau_{{A}}^2\omega_{{A}}^2k^2V_{{A}}^2+q\varOmega\tau_{{A}}\omega_{{A}} k_xV_{{A}y}+ \ell_{{H}} k_z\omega_{{A}}\left(\frac{k^2}{k_z^2}\ell_{{H}} k_z\omega_{{A}} -q\varOmega\right), \end{gather}

(6.46)\begin{gather} \mathcal{C}_1=\mathcal{C}_1\left(\omega_{{A}}^2+\frac{k_z^2}{k^2}\kappa^2\right), \end{gather}

(6.47)\begin{align} \mathcal{C}_0\!&=\!\omega_{{A}}^2\left(\omega_{{A}}^2\underbrace{-2q\varOmega^2\frac{k_z^2}{k^2}}_\textrm{MRI} \right)\!+\!\kappa^2k_z^2k^2\eta_O^2\!+\!\ell_{{H}}\omega_{{A}} k_z\left(\underbrace{(4-q)\varOmega\omega_{{A}}^2}_\textrm{ion-cyclotron instability}\underbrace{-q\varOmega\frac{k_z^2}{k^2}\kappa^2}_\textrm{HSI}+\ell_{{H}} \omega_{{A}} k_z \kappa^2\right)\!, \nonumber\\ &\quad +\kappa^2\tau_{{A}}^2\omega_{{A}}^2k_z^2V_{{A}}^2+\underbrace{q\varOmega \tau_{{A}} k_x V_{{A}y}\omega_{{A}} \left(\frac{k_z^2}{k^2}\kappa^2+\omega_{{A}}^2\right)}_{\textrm{Oblique ambipolar modes}}. \end{align}

The stability of this linear system can be analysed in the vicinity of $\sigma =0$. In this case, a necessary and sufficient condition for instability is $\mathcal {C}_0<0$. This allows us to identify three sources of instability: the usual MRI, the ion-cyclotron instability, the Hall-shear instability (HSI) and the term at the origin of Oblique ambipolar modes, which is not a genuine instability branch.

Before exploring the non-ideal regime, let us point out that in the ideal MHD limit, this dispersion relation shows that modes with $k_x\ne 0$ have a lower growth rate than $k_x=0$ modes. As a result, $k_x=0$ modes are always the most unstable eigenmodes of the system. These modes are often called ‘channel modes’ as they do not have any horizontal spatial dependency. They are also exact non-linear solutions of the full MHD equations (Goodman & Xu Reference Goodman and Xu1994). For this reason, they are very robust and they often show up in the non-linear regime, as we show later.

6.4.3. Ohmic diffusion

The effect of Ohmic diffusion on the stability of the MRI is physically very intuitive: it stabilises MRI modes, starting from the largest $\omega _{{A}}$ of the system (see figure 19). The stability condition in the presence of Ohmic diffusion is deduced from the condition $\mathcal {C}_0=0$ and reads

(6.48)\begin{equation} \frac{2q\varOmega^2}{k^2V_{{A}z}^2}-1<\kappa^2\frac{k^2}{k_z^2}\frac{\eta_O^2}{V_{{A}z}^4} \rightarrow\textrm{Stability.} \end{equation}

Clearly, the modes with $k_x=0$ are the last to be stabilised when one increases Ohmic diffusion. Let us focus on this case and assume the flow is Keplerian so that $\kappa =\varOmega$ and $q=3/2$. The stability of the most unstable ideal mode is often considered as a proxy for the stability of the flow. This mode has $\omega _{{A}}=\sqrt {3}\varOmega /2$ so that the inequality reduces to

(6.49)\begin{equation} \varLambda_O^2<\frac{1}{3}\rightarrow \textrm{Stability of the most unstable ideal MRI mode.} \end{equation}

Note, however, that this does not imply that all of the modes available to the disc are stable, and that the disc is stable. A more constraining criterion results from this and requires that the mode with the largest length scale is MRI stable, i.e. that

(6.50)\begin{equation} \frac{3\varOmega^2}{k_{z,\mathrm{min}}^2V_{{A}z}^2}-1 < \varLambda_O^{-2} \rightarrow\textrm{General disc stability.} \end{equation}

Figure 19. MRI growth rate as a function of the Ohmic Elsasser number $\varLambda _O$ for a Keplerian disc ($q=3/2$) with $k_x=0$. Note the damping of the most unstable mode for $\varLambda _O=1/\sqrt {3}$ and the survival of low growth rate modes in the limit $\omega _{{A}}\rightarrow 0,~ \varLambda _O\rightarrow 0$.

6.4.4. Hall effect

Hall-driven linear waves: The Hall effect is known to be at the origin of new linear waves. These waves can be captured by letting $\varOmega \rightarrow 0$ and neglecting Ohmic and ambipolar diffusion. In this case the dispersion relation (6.43) gives

(6.51)\begin{align} 0&=\sigma^4+\sigma^2\left(2\omega_{{A}}^2+\ell_{{H}}^2\omega_{{A}}^2k^2\right)+\omega_{{A}}^4\nonumber\\ &=\left(\sigma^2+\textrm{i}\sigma\ell_{{H}} \omega_{{A}} k+\omega_{{A}}^2\right) \left(\sigma^2-\textrm{i}\sigma\ell_{{H}} \omega_{{A}} k+\omega_{{A}}^2\right). \end{align}

We recognise two waves with frequency $\omega \equiv \textrm {i}\sigma$ given by

(6.52)\begin{equation} \omega=\omega_{{A}}\left[\pm\frac{\ell_{{H}} k}{2}+ \sqrt{\frac{\ell_{{H}}^2k^2}{4}+1}\right], \end{equation}

where ‘$+$’ waves are known as whistlers or electron-cyclotron modes, whereas ‘$-$’ waves are ion-cyclotron modes. The whistler frequency $\omega _{{H}}=\omega _{{A}}\ell _{{H}} k$ increases as $k^2$ in the limit $k\rightarrow \infty$ whereas the ion-cyclotron frequency tends to a constant $\omega _{\mathrm {IC}}=\omega _{{A}}/(\ell _{{H}} k)$. Hence, these two waves behave very differently at the small scale. For positive $\ell _{{H}}$, whistlers are right-handed polarised wave whereas ion-cyclotron are left-handed. Physically, whistlers are essentially an oscillation of the electrons fluid (or of the lightest charged particle), leaving all of the other components of the plasma unaffected. Whistler and ion-cyclotron waves become standard right and left-handed circularly polarised Alfvén wave in the limit $k\rightarrow 0$.

HSI: The HSIFootnote ¹³ is a new branch of instability (Kunz Reference Kunz2008), which is often confused with the traditional MRI, despite its different physical origin. It is essentially an instability of whistler waves under the action of shear.

To capture the HSI, one can let $\omega _{{A}}\rightarrow 0$ while keeping $\ell _{{H}} \omega _{{A}}>0$. This ‘low magnetisation’ limit allows one to decouple the ions from the electrons, as is evident from (6.42). Neglecting Ohmic and ambipolar diffusion, one obtains the following dispersion relation

(6.53)\begin{equation} \left(\sigma^2+\frac{k_z^2}{k^2}\kappa^2\right)\left[\sigma^2+\ell_{{H}}\omega_{{A}} k_z\left(\frac{k^2}{k_z^2}\ell_{{H}} k_z\omega_{{A}}-q\varOmega\right)\right]=0, \end{equation}

which exhibits a linear instability when

(6.54)\begin{equation} \ell_{{H}} V_{{A}z} k_z^2\left(k^2\ell_{{H}} V_{{A}z} -q\varOmega\right) < 0\rightarrow\textrm{HSI unstable}. \end{equation}

Interestingly, the HSI shows up only when $q \ell _{{H}} V_{{A}z}>0$ or, in other words, when the vertical field points in the same direction as the rotation axis in Keplerian discs,Footnote ¹⁴ assuming $\ell _{{H}} > 0$. When the whistler frequency becomes too large ($k^2\ell _{{H}} V_{{A}z} > q\varOmega$), the instability disappears. For a given $k_z$, the most unstable mode has $k_x=0$, hence $k=k_z$. For this reason, the HSI often shows up as channel-like mode in simulations, in a similar way to the MRI. Last, the maximum growth rate is identical to the MRI $\sigma _\mathrm {max}=q\varOmega /2$ and is obtained for $k^2\ell _{{H}} V_{{A}z}=q\varOmega /2$.

Physically, this instability is a result of sheared whistler waves. If we look at the disc from the top with the vertical field pointing towards us, the magnetic perturbation of a whistler wave will tend to rotate counter-clockwise (figure 20). If the vertical field is positive (i.e. aligned with the rotation axis), the Keplerian shear is going to stretch the perturbation in the opposite direction, amplifying the toroidal field from the radial field and feeding the instability. In contrast, if the vertical field is anti-aligned with the rotation axis, the direction of rotation of whistler perturbations is that of the shear, resulting in a damped whistler wave.

Figure 20. Physical principle of the HSI. The magnetic perturbation (in green) is rotated (a) clockwise or (b) counter-clockwise by the Hall effect depending on the polarity of the mean field $B_0$. When $B_0 > 0$, the rotated perturbation is amplified by the shear (in blue) while it is damped by the shear when $B_0 < 0$.

Ion-cyclotron instability: The ion-cyclotron instabilityFootnote ¹⁵ is more difficult to isolate compared with the HSI because ion inertia cannot be neglected in this case. However, it is still possible to filter out whistler waves by letting $k_z\rightarrow \infty$ and keeping constant $k_x\lesssim k_z$. In this case, the whistler wave frequency $\omega _{{H}}$ becomes infinite whereas the cyclotron frequency remains finite. As we are looking for a finite growth rate, we assume $\sigma$ is finite when $k\rightarrow \infty$. Neglecting ambipolar and Ohmic diffusion and keeping only $O(k^4)$ terms, the dispersion relation (6.43) becomes (see also Simon et al. Reference Simon, Lesur, Kunz and Armitage2015)

(6.55)\begin{equation} \sigma^2 \ell_{{H}}^2k^2\omega_{{A}}^2+\left[\omega_{{A}}^2+(2-q)\varOmega \ell_{{H}} \omega_{{A}} k_z\right]\left[\omega_{{A}}^2+2\varOmega \ell_{{H}} \omega_{{A}} k_z\right]=0. \end{equation}

This relation clearly describes ion-cyclotron modes because in the limit $\varOmega \rightarrow 0$, we recover the ion-cyclotron frequency $\sigma =\pm \textrm {i}\omega _{\mathrm {CI}}$. It describes unstable modes, resulting from the interaction between ion-cyclotron waves and epicyclic motions, provided that

(6.56)\begin{equation} -\frac{V_{{A}z}^2}{(2-q)\varOmega}<\ell_{{H}}V_{{A}z}<-\frac{V_{{A}z}^2}{2\varOmega} \rightarrow\textrm{Ion-cyclotron unstable}. \end{equation}

This inequality brings up two remarks: (i) the ion-cyclotron instability appears for anti-aligned field configurations, i.e. field configuration opposite to the HSI; (ii) this instability does not vanish in the limit $k_z\rightarrow \infty$. Therefore, there is no small-scale quenching similar to the MRI and the HSI, but there is a field strength limit.

General case: In the general case, one cannot compute the growth rate of the Hall-MRI easily. To illustrate the general growth rate of the Hall-MRI, we present in figure 21 the growth resulting from (6.43) in the Hall-only case. To quantify the intensity of the Hall effect, we have defined a Lundquist number based on the vertical wavenumber

(6.57)\begin{equation} \mathcal{L}_\mathcal{H}^*=\frac{1}{\ell_{{H}} k_z}=\frac{V_A}{\eta_H k_z}. \end{equation}

We have chosen the most unstable modes $k_x=0$ and assumed $k_z>0$ so that $\omega _{{A}} < 0$ corresponds to $B_{0z}$ anti-aligned with $\varOmega$. As can be seen in figure 21, we recover the MRI in the limit $\mathcal {L}_{\mathcal {H}}\rightarrow \infty$, which gives identical growth rates under the symmetry $\omega _{{A}}\rightarrow -\omega _{{A}}$. As $\mathcal {L}_{\mathcal {H}}^*$ decreases and the Hall effect increases, this symmetry is broken.

Figure 21. Hall-MRI growth rate as a function of the modified Hall Lundquist number $\mathcal {L}_{\mathcal {H}}^*=(\ell _{{H}} k_z)^{-1}$ for a Keplerian disc ($q=3/2$) with $k_x=0$. We have assumed $k_z > 0$ so that $\omega _{{A}} < 0$ corresponds to the anti-aligned case $V_{{A}z}<0$. The white plain line corresponds to the HSI stability limit (6.54) and the white dashed lines to the ion-cyclotron stability limits (6.56).

For $\omega _{{A}} > 0$ (aligned field configuration), the MRI becomes the HSI, and the optimum growth rate starts to move to lower $\omega _{{A}}$. This is expected from our analysis because the maximum growth rate of the HSI is found for $\omega _{{A}}/\varOmega \simeq \mathcal {L}_{H}^* /2$. We show in white filled contours the stability limit of the HSI obtained from (6.54). This limit matches the full dispersion relation for $\mathcal {L}_{H}^*\lesssim 0.3$.

For $\omega _{{A}} < 0$ (anti-aligned field configuration), the MRI becomes the ion-cyclotron instability and the optimum growth rate moves to higher $|\omega _{{A}}|$. The stability contours of the ion-cyclotron instability (6.56) are shown as dashed white lines. As for the HSI, they match the full dispersion relation for $\mathcal {L}_{H}^*\lesssim 0.3$.

Effect of Ohmic diffusion on the Hall-MRI: As shown previously, Ohmic diffusion tends to damp MRI modes, starting from the largest $\omega _{{A}}$ and moving to lower $\omega _{{A}}$ as $\varLambda _O$ decreases. On the other hand, the Hall effect creates two distinct branches from the MRI, depending on the field alignment configuration. As Ohmic diffusion suppresses first high $|\omega _{{A}}|$ modes, the ion-cyclotron instability will be the first to be stabilised. On the other hand, because the HSI is living at lower $|\omega _{{A}}|$ compared with the MRI, it will be less affected by Ohmic diffusion than the MRI. This effect is illustrated in figure 22, demonstrating that HSI modes survive more easily to Ohmic diffusion, even when ideal MRI modes are suppressed.

Figure 22. Same as figure 21 but including Ohmic diffusion with $\varLambda _O=0.1$. Note the complete stabilisation of the ion-cyclotron branch and the strong damping of ideal MRI modes. Only HSI modes subsist at low $\mathcal {L}_{H}^*$ with growth rates comparable with the ideal MRI case.

We can go further and estimate the stability limit of the HSI in the limit $\omega _{{A}}\rightarrow 0$ while keeping $\ell _{{H}} \omega _{{A}}>0$. In this limit, the dispersion relation of the HSI including Ohmic diffusion reads

(6.58)\begin{equation} \left(\sigma^2+\frac{k_z^2}{k^2}\kappa^2\right) \left[\left(\sigma+\eta_Ok^2\right)^2+\ell_{{H}}\omega_{{A}} k_z \left(\frac{k^2}{k_z^2}\ell_{{H}} k_z\omega_{{A}}-q\varOmega\right)\right]=0, \end{equation}

from which we can deduce a general instability criterion. For the most unstable HSI mode, this criterion reads

(6.59)\begin{equation} \varLambda_O^2 < \frac{1}{2}\varLambda_H^2 \rightarrow \textrm{Stability of the most unstable HSI mode}, \end{equation}

where $\varLambda _H$ is the Hall Elsasser number (see § 3.4.2). This expression can be compared directly with the Ohmic-only criterion (6.49), and shows that, as expected, the HSI can make the system unstable even when $\varLambda _O\ll 1$, provided that $\varLambda _H\ll 1$ as well.

For this reason, some researchers have proposed that the MRI could be ‘resuscitated’ in regions having a strong Ohmic diffusion thanks to the presence of the Hall effect, which could be dominant in some parts of the disc (Wardle & Salmeron Reference Wardle and Salmeron2012). The physical interpretation of this effect is relatively simple. In the case of the MRI, the maximum growth rate is found for $\omega _{{A}}\sim \varOmega$, in other words, the Alfvénic and rotation frequencies match. In the HSI case, it is the whistler and rotation frequencies $\omega _{{H}}\sim \varOmega$ that have to match. In the limit of strong Hall effect, $\omega _{{H}}= \ell _{{H}} k_z \omega _{{A}}$. Therefore, when $\ell _{{H}} kz\gg 1$, the scale at which the whistler frequency matches the rotation frequency is much larger than the scale at which the Alfvén frequency matches the rotation frequency. In other words, the optimum HSI mode has a much larger wavelength than its MRI counterpart. For this reason, the HSI is less sensitive to diffusion than the MRI, because diffusion first damps small-scale modes.

Although these statements are true in the linear regime, the non-linear saturation of the HSI is not guaranteed to be similar to that of the MRI. Effectively, the linear analysis is unable to predict the turbulent angular momentum transport one could obtain from the HSI. As we show in § 8.4, the HSI in the non-linear regime is indeed full of surprises.

6.4.5. Ambipolar diffusion

The linear ambipolar diffusion is a non-diagonal operator as shown in (6.42). The non-diagonal terms are non-zero only when $k_x V_{{A}y}\ne 0$, i.e. when both non-axial wavevectors and guide fields are considered. Otherwise, ambipolar diffusion acts as a usual diffusion operator by damping magnetic perturbation. Note, however, that even in this case, it does not act as a scalar diffusion, unless $k_x=0$ and $V_{{A}y}=0$.

Diagonal case ($k_x V_{{A}y}= 0$): In this case, the stability criterion is very similar to that of Ohmic diffusion. By requiring $\mathcal {C}_0 > 0$, one obtains

(6.60)\begin{equation} \frac{2q\varOmega^2}{k^2V_{{A}z}^2}-1<\kappa^2\tau_{{A}}^2\left(\frac{V_{{A}}}{V_{{A}z}}\right)^2\rightarrow\textrm{Stability.} \end{equation}

This criterion shows that, as for Ohmic diffusion, the most unstable modes have $k_x=0$. Moreover, it shows that the stability condition is affected by the toroidal field strength, a stronger $V_{{A}y}$ leading to a more stable system. As for Ohmic resistivity, we can derive a criterion for the stability of the most unstable ideal MRI mode in Keplerian flows, which reads

(6.61)\begin{equation} \varLambda_A^2<\frac{1}{3}\left(\frac{V_{{A}}}{V_{{A}z}}\right)^2 \rightarrow\textrm{Stability of the most unstable ideal MRI mode,} \end{equation}

and a criterion for the stability of all the modes available smaller than the minimum vertical wavenumber ${k_{z,\mathrm {min}}}$:

(6.62)\begin{equation} \frac{3\varOmega^2}{k_{z,\mathrm{min}}^2V_{{A}z}^2}-1 < \varLambda_A^{-2} \left(\frac{V_{{A}}}{V_{{A}z}}\right)^2\rightarrow\textrm{General disc stability.} \end{equation}

This shows that the MRI can be stabilised even for $\varLambda _A > 1$ provided that the toroidal field dominates over the poloidal field. The field topology is therefore a key point for the stability under the action of ambipolar diffusion. The physical reason for this is relatively simple. The MRI mechanism is only sensitive to the magnetic tension, which is due to $B_z$ in the absence of non-axisymmetric perturbations. In contrast, ambipolar diffusion increases with the total field strength. Hence, an increase in $V_{{A}y}$ results in an increase of the effective diffusion $\eta _A$, whereas the MRI feedback loop is left unperturbed.

Non-diagonal case ($k_xV_{{A}y}\ne 0$): In this case, the non-diagonal terms of ambipolar diffusion can act as a positive feedback loop in the induction equation. The criterion for this positive feedback is easily obtained from the dispersion relation (6.43): $q\varOmega k_x k_z V_{{A}y} V_{{A}z} < 0$. An illustration of the effect of this term on the general stability property is shown in figure 23 where we have assumed $V_{{A}y}=V_{{A}z}$, $q=3/2$ and $\varLambda _A=0.4$. We observe that oblique modes ($k_x\ne 0$) tend to be more unstable than axial modes when $\varLambda _A < 1$. It can be shown that oblique modes are always unstable, albeit with a vanishing growth rate, in the limit $\varLambda _A\rightarrow 0$ provided that $k_x/k_z$ is sufficiently large (Kunz & Balbus Reference Kunz and Balbus2004).

Figure 23. MRI growth rate with ambipolar diffusion and $\varLambda _A=0.4$. We have chosen $V_{{A}y}=V_{{A}z}$ and $k_x\ne 0$ to illustrate the effect of non-diagonal ambipolar terms. The most unstable mode in this case is found for $k_x < 0$ and a weakly unstable branch exists for $k_x/k_z\simeq 3.5$ in the limit $\omega _{{A}}\rightarrow \infty$. We name these modes ‘oblique ambipolar modes’. Note, however, their low relative growth rate.

6.4.6. Vertical stratification

All of the results described previously were computed ignoring vertical stratification. It should, however, be pointed out that because discs are vertically stratified, unstratified results should be taken with care. Let us emphasise here a few important results regarding the effect of stratification on MRI modes. In this section, we assume that the disc is vertically isothermal (see § 4.2) and is only threaded by a vertical field (no toroidal field). We quantify the intensity of the imposed field with the plasma beta parameter in the midplane:

(6.63)\begin{equation} \beta_{\mathrm{mid}} \equiv\frac{8{\rm \pi} P_{\mathrm{mid}}}{B^2}=\frac{2c_s^2}{V_{{A}z}^2 (z=0)}. \end{equation}

Ideal MRI: The ideal MRI case with an isothermal disc is presented in details by Latter, Fromang & Gressel (Reference Latter, Fromang and Gressel2010). We therefore just summarise here the main conclusions. First, stratified eigenvalues (growth rates) satisfy the same dispersion relation (6.14) as non-stratified modes with $V_{{A}}=V_{{A}}(z=0)$. Vertical wavenumbers $k$ are quantised, but eigenmodes are not simple harmonic functions in the $z$ direction. Latter et al. (Reference Latter, Fromang and Gressel2010) have shown that the first eigenmodes have $k_nH=1.1584,2.0796,2.9829,3.8798,\ldots$ for $n=1\ldots 4$. Taking $k_1$ as the lowest wavenumber accessible to the system, the MRI is stable for all possible eigenmodes in Keplerian discs ($q=3/2$) if $V_{{A}z}(z=0)>\sqrt {3}\varOmega /k_1$, i.e. if

(6.64)\begin{equation} \beta_{\mathrm{mid}}<\frac{2(k_1H)^2}{3} \simeq 0.89 \rightarrow\textrm{Stability.} \end{equation}

This value can vary slightly depending on the boundary conditions for the perturbation as $z\rightarrow \infty$. Overall, it is safe to assume that the MRI exists provided that $\beta _{\mathrm {mid}}\gtrsim 1$.

Second, eigenmodes come in two symmetries: odd symmetry modes (which corresponds to odd $n$) have odd $v_{x,y}(z)$, even $B_{x,y}(z)$ and exhibit a maximal magnetic perturbation at the midplane and no velocity perturbation at this location. Even symmetry modes (even $n$), on the other hand, have a velocity jet in the midplane and no magnetic perturbation at $z=0$ (see figure 24). The MRI does not choose specifically any symmetry in these local stratified models: both even and odd $n$ follow the same dispersion relation.

Figure 24. (a,b) MRI eigenmodes computed at $\beta _{\mathrm {mid}}=5$ for $n=1$ (left, odd mode, $\sigma =0.72\varOmega$) and $n=2$ (right, even mode $\sigma =0.67\varOmega$). (c,d) Schematic representation of the field perturbation owing to odd (left) and even (right) modes.

Third, the fastest growing modes tend to be more oscillatory in the disc midplane at higher $\beta _{\mathrm {mid}}$. This is because qualitatively, one expects $kV_{{A}z}\sim \varOmega$ for the fastest growing mode, hence in the midplane, one obtains $k_zH\sim \sqrt {\beta _{\mathrm {mid}}}$. An example of such a mode is shown in figure 25.

Figure 25. Fastest growing MRI eigenmode ($\sigma =0.75\varOmega$) in a stratified model for $\beta _{\mathrm {mid}}=10^3$. Note the strongly oscillating behaviour close to the midplane.

Ohmic and ambipolar diffusion: As shown in the non-stratified section, Ohmic and ambipolar diffusion tend to suppress the MRI when the Elsasser numbers are less than $1$. Typically, this situation occurs in the densest regions of the disc (see figure 10). We show in figure 26 the most unstable MRI eigenmode with $\beta =10^3$ resulting from our fiducial metal-free mode at $R=1\ \mathrm {AU}$. As expected, the MRI is strongly suppressed where $\varLambda _O < 1$. Ambipolar diffusion tends to reduce the overall growth rate but does not affect the shape of the eigenmode significantly. This is a perfect linear illustration of the historical ‘layered accretion’ paradigm where the MRI still survives at the surface of the disc, leaving the disc midplane up to one to two scale-heights essentially magnetically dead (Gammie Reference Gammie1996).

Figure 26. Fastest growing MRI eigenmode in a stratified model for $\beta _{\mathrm {mid}}=10^3$ using the diffusivity profiles from figure 10 at $R=1\ \mathrm {AU}$. (a) Including Ohmic diffusion only ($\sigma =0.72\varOmega$) and (b) including Ohmic and ambipolar diffusion ($\sigma =0.56\varOmega$). Note the lack of perturbations in the disc midplane due to Ohmic diffusion and partially to ambipolar diffusion. Equivalent eigenmodes with the same growth rates are found on the $z>0$ side of the disc.

Note that when the linear eigenmodes do not propagate in the disc, as is the case here, there are no well-defined ‘odd’ and ‘even’ modes. Instead, there are two families of modes localised either at the top or the bottom side of the disc, with identical growth rates. It is, in principle, possible to combine these top and bottom modes to reconstruct even and odd modes, but here we have chosen here to stress the lack of propagation through the disc of the perturbation.

Hall effect: As discussed previously, the Hall effect can potentially revive dead zones when the field is aligned with the rotation axis, thanks to the HSI branch of the Hall-MRI. This effect is illustrated in figure 27 where the MRI eigenmode now propagates in the midplane resulting in a fully active disc column at $R=1\ \mathrm {AU}$ when $V_{{A}z}\varOmega > 0$. In this case, however, the growing perturbation is mostly magnetic and velocity perturbations are mostly absent in the midplane of the eigenmode. This is because the HSI is an unstable whistler mode, which is essentially an electronic perturbation leaving the ions (and the neutrals) unperturbed. When the field is anti-aligned, we recover results similar to the purely diffusive case, albeit with a slightly reduced growth rate.

Figure 27. Fastest growing MRI eigenmode in a stratified model for $\beta _{\mathrm {mid}}=10^3$ using the diffusivity profiles from figure 10 at $R=1\ \mathrm {AU}$ including Ohmic, Hall and ambipolar diffusion: (a) assuming $V_{{A}z}\varOmega > 0$ ($\sigma =0.65\varOmega$) and (b) assuming $V_{{A}z}\varOmega < 0$ ($\sigma =0.51\varOmega$). The ‘dead zone’ is now subject to long-wavelength perturbations when $V_{{A}z}\varOmega > 0$, as a result of the HSI.

8.1.1. With a mean field

The first historical models of MRI turbulence (Hawley et al. Reference Hawley, Gammie and Balbus1995) tested both mean vertical and toroidal fields and were computed in the ideal MHD framework. MRI unstable modes grow, break up presumably thanks to secondary instabilities that are essentially Kelvin–Helmholtz-type modes (Goodman & Xu Reference Goodman and Xu1994; Latter, Lesaffre & Balbus Reference Latter, Lesaffre and Balbus2009) and end up in fully developed MHD turbulence.

In the mean vertical field case, Hawley et al. (Reference Hawley, Gammie and Balbus1995) found that and $\alpha \simeq 6.7 \beta _{\mathrm {mean}}^{-1/2}$. However, it was later pointed out by Bodo et al. (Reference Bodo, Mignone, Cattaneo, Rossi and Ferrari2008) that, in the presence of a mean vertical field, the box aspect ratio used by Hawley et al. (Reference Hawley, Gammie and Balbus1995) was prone to recurrent channel mode solutions, which was absent in wider boxes (more elongated in $x$). This implies that $\alpha$ and $\langle \beta \rangle$ are overestimated by a factor of two compared with wider boxes (Bodo et al. Reference Bodo, Mignone, Cattaneo, Rossi and Ferrari2008). We therefore correct (Hawley et al. Reference Hawley, Gammie and Balbus1995) for this and obtain the following scaling law in the range $400<\beta _{\mathrm {mean}} < 5\times 10^4$:

(8.5)\begin{equation} \left.\begin{aligned} & \alpha\simeq 0.61\langle\beta\rangle^{-1},\\ & \alpha\simeq 3.3 \beta_{\mathrm{mean}}^{-1/2}, \end{aligned}\right\} \end{equation}

a scaling which has also been verified in incompressible simulations (Longaretti & Lesur Reference Longaretti and Lesur2010). The extrapolation down to $\beta \rightarrow 1$ suggests that $\alpha \gtrsim 1$ could be reached. However, it is very difficult to explore this limit numerically because the MRI in this case becomes very strong and channel modes never break up into developed turbulence (Hawley et al. Reference Hawley, Gammie and Balbus1995). Lesur & Longaretti (Reference Lesur and Longaretti2007) explored this limit in the incompressible regime and found that turbulence was in this case very intermittent with long ‘quiet’ episodes of linear growth followed by strong and sudden bursts of turbulence. In this situation, a constant $\alpha$ is probably not a good model of the disc physics.

In the mean toroidal field case, the scaling deduced from Hawley et al. (Reference Hawley, Gammie and Balbus1995) in the range $2 < \beta _{\mathrm {mean}} < 1200$ gives lower $\alpha$ values:

(8.6)\begin{equation} \left.\begin{gathered} \alpha\simeq 0.51\langle\beta\rangle^{-1},\\ \alpha\simeq 0.35 \beta_{\mathrm{mean}}^{-1/2}. \end{gathered}\right\} \end{equation}

8.1.2. Zero mean field: the ‘MRI dynamo’

The existence of an MRI-driven dynamo was first demonstrated by Hawley, Gammie & Balbus (Reference Hawley, Gammie and Balbus1996). In this configuration, no external field is imposed and turbulence regenerates a field on which the MRI can grow, feeding back turbulent motions. Owing to this feedback loop, the system needs a finite-amplitude perturbation to sustain turbulence (Rincon et al. Reference Rincon, Ogilvie, Proctor and Cossu2008), implying that the instability is in this case subcritical. The dynamo feedback loop has been the subject of intense studies after the discovery of dynamo cycles (Lesur & Ogilvie Reference Lesur and Ogilvie2008b), both based on a quasi-linear theory of the toroidal MRI (Lesur & Ogilvie Reference Lesur and Ogilvie2008a) and on a dynamical system approach (Herault et al. Reference Herault, Rincon, Cossu, Lesur, Ogilvie and Longaretti2011; Riols et al. Reference Riols, Rincon, Cossu, Lesur, Longaretti, Ogilvie and Herault2013).

According to Hawley et al. (Reference Hawley, Gammie and Balbus1996), the MRI dynamo yields $\alpha \sim 0.01$. However, it was later realised that the value of $\alpha$ in ideal MHD simulations (where no physical dissipation is introduced) depends on the numerical resolution (Fromang & Papaloizou Reference Fromang and Papaloizou2007), with $\alpha \propto 1/N$ where $N$ is the number of grid points in one direction.Footnote ¹⁶ This problem of numerical convergence (because numerics do not seem to be converging as one increases the resolution) implies that the estimation for $\alpha$ from zero net flux simulation is intrinsically flawed. Obviously, this is most probably a result of numerical dissipation (the only dissipation channel of these simulations), which is not acting as a real physical dissipation operator. If one introduces explicit viscosity $\nu$ and resistivity $\eta _O$ in the system, $\alpha$ then converges to a finite value, which seems to depend only on ${Pm}\equiv \nu /\eta _O$ (Fromang Reference Fromang2010) and no longer depend on the resolution. This, however, leads to another complication known as the ‘Pm effect’.

8.1.3. The Pm effect

The Pm effect shows up in quasi-ideal simulations both with and without a mean vertical field. When one introduces only a small amount of viscosity $\nu$ and resistivity $\eta _O$, the large-scale linear MRI modes are largely unaffected. For this reason, this regime corresponds to a ‘quasi ideal-MHD’ regime. However, it turns out that the saturation level depends on the magnetic Prandtl number. A lot of literature has been devoted to this effect (Fromang et al. Reference Fromang, Papaloizou, Lesur and Heinemann2007; Lesur & Longaretti Reference Lesur and Longaretti2007). In simulations with a mean field (Lesur & Longaretti Reference Lesur and Longaretti2007; Simon & Hawley Reference Simon and Hawley2009), the Pm effect results in an increase of $\alpha$ when ${Pm}$ increases. In the zero mean field case, the effect is even stronger because the MRI dynamo, and therefore MHD turbulence, simply disappears when ${Pm}\lesssim 1$ (Fromang et al. Reference Fromang, Papaloizou, Lesur and Heinemann2007; Walker & Boldyrev Reference Walker and Boldyrev2017), giving $\alpha =0$ in that regime.

This effect shows up even in situations where the Reynolds and Elsasser numbers are much larger than one, indicating that linear or quasi-linear theory cannot explain it (Longaretti & Lesur Reference Longaretti and Lesur2010). Instead, it has been proposed that non-local energy transfers in spectral space (Lesur & Longaretti Reference Lesur and Longaretti2011) could be responsible for this effect. If true, because non-locality is necessarily bounded (Aluie & Eyink Reference Aluie and Eyink2010), the ${Pm}$ effect should disappear as $(\eta _O,\nu )\rightarrow 0$. An alternative viewpoint is to separate the small ‘turbulent’ scales from the large scales where the dynamo mechanism is presumably lying. By carefully computing the energy exchange between the scales, one finds that the small scales act as an effective viscosity when ${Pm} < 1$ and tend to damp the large-scale mechanism (Riols et al. Reference Riols, Rincon, Cossu, Lesur, Ogilvie and Longaretti2013). This, however, is not fully satisfactory as it does not prove that the MRI dynamo is non-existent in the limit ${Pm}\rightarrow 0$ and ${Rm}\rightarrow \infty$.

Despite numerous efforts, recent results indicate that the ${Pm}$ effect is still very much alive in simulations, up to ${Rm}=5\times 10^4$ and $\varLambda _O=O(100)$ (Potter & Balbus Reference Potter and Balbus2017). This effect is still a very open question regarding the MRI saturation level in the ‘quasi-ideal’ regime.

Note, finally, that the Pm effect is relevant only for the innermost parts of PPDs where ${Rm}\gg 1$ is expected. For the regions above $1\ \mathrm {AU}$, large-scale physics is dominated by magnetic diffusion and the $Pm$ effect is not relevant.

8.5.1. Ideal MHD

Hawley (Reference Hawley2001) and Steinacker & Papaloizou (Reference Steinacker and Papaloizou2002) were among the first to note the formation of ‘ring-like’ structures in MRI simulations with net vertical flux. The simulations are in these cases semi-global: vertical stratification is neglected whereas the radial curvature is retained, leading to a cylindrical setup. It is found that the net vertical flux is trapped in a low-density region, forming a gap. In these gaps, $\alpha$ can reach values as high as 1, consistently with the fact that these gaps correspond to low $\beta _{\mathrm {mean}}$ regions. Hawley (Reference Hawley2001) proposed that this could be the signature of a viscous instability: a local density minimum results in a local decrease of $\beta _{\mathrm {mean}}$ (assuming the mean field is kept at its initial value). As $\alpha \propto \beta _{\mathrm {mean}}^{-1/2}$ (see (8.5)), $\alpha$ increases in this region, which removes mass from the region because of angular momentum conservation. Unfortunately, this proposition has never been investigated further, leaving the origin of self-organisation in these simulations unexplained.

The same phenomenon was reported in local shearing box models with a mean vertical field by Bai & Stone (Reference Bai and Stone2014). They proposed that the non-diagonal components of the turbulent resistivity tensor (Lesur & Longaretti Reference Lesur and Longaretti2009) could be at the origin of the effect by acting as an ‘effective negative diffusivity’. This explanation is, however, dubious because several authors have measured the turbulent resistivity tensor of MRI turbulence (Fromang & Stone Reference Fromang and Stone2009; Guan & Gammie Reference Guan and Gammie2009; Lesur & Longaretti Reference Lesur and Longaretti2009) and found that the effective resistivity was always positive. Strikingly, Bai & Stone (Reference Bai and Stone2014) simulations clearly show that despite having a globally concentrated field $B_z(x)$, the turbulent electromotive force $\boldsymbol {{\mathcal {E}}}=\langle \boldsymbol {w}\boldsymbol {\times } \boldsymbol {B}\rangle$ does not depend on $x$ (Bai & Stone Reference Bai and Stone2014, figure 3). Hence, the turbulent resistivity prescription $\boldsymbol {{\mathcal {E}}}=\eta _\mathrm {turb}\boldsymbol {\nabla }\boldsymbol {\times }\langle \boldsymbol {B}\rangle$ likely breaks down altogether and should be replaced with a more elaborate closure scheme.

In simulations without net flux, Fromang & Nelson (Reference Fromang and Nelson2005) reported the spontaneous formation of giant anticyclones, though this could be a boundary condition artefact (Fromang, private communication). Johansen, Youdin & Klahr (Reference Johansen, Youdin and Klahr2009) also reported the formation of ‘pressure bumps’ that are caused by large-scale fluctuations of $\alpha$. In contrast to the simulations with a mean field vertical field, the features observed in the zero net field case are transient and only survive for a limited time, which depends on the box size. Johansen et al. (Reference Johansen, Youdin and Klahr2009) proposed a model based on a stochastic $\alpha$, which predicts long-lived axisymmetric structures in quasi-geostrophic equilibrium, as observed in their simulations.

8.5.2. Hall-MHD

The first mention of self-organisation in Hall-MHD appears in Kunz & Lesur (Reference Kunz and Lesur2013), where self-organisation has a dramatic effect on the saturation level of Hall-dominated MRI (see § 8.4). Self-organisation appears in an obvious manner by looking at the vertical field component of the flow (figure 32). Although in the ideal-MHD case, self-organisation appears as a ‘second-order’ effect on top of MRI turbulence, in the Hall-MHD case, self-organisation is the main saturation mechanism of the MRI. In other words, turbulence is mostly suppressed by self-organisation. Self-organisation shows up when $\mathcal {L}_H \lesssim 5$, and $\alpha$ essentially vanishes as a result of the lack of turbulence (figure 30). This surprising result also holds in the case with zero net flux or in the mixed case having both a mean azimuthal and vertical magnetic flux.

Figure 32. Vertical field component in snapshots of MRI turbulence. (a) Ideal-MHD simulation with $\beta _{\mathrm {mean}}=3200$. (b) Hall-MRI simulation with $\beta _{\mathrm {mean}}=3200$ and $\mathcal {L}_H=1.75$. Figure from Kunz & Lesur (Reference Kunz and Lesur2013).

The origin of Hall-driven self-organisation can be tracked down to the induction equation in the presence of a Hall effect. Indeed, When the Hall length is constant, the induction equation reads

(8.12)\begin{equation} \partial_t\boldsymbol{B}=\boldsymbol{\nabla}\boldsymbol{\times} (\boldsymbol{w}\boldsymbol{\times} \boldsymbol{B})+\ell_{{H}} \boldsymbol{\nabla}\boldsymbol{\times}\left(\boldsymbol{\nabla }\boldsymbol{\cdot} \frac{-{\boldsymbol{BB}}}{4{\rm \pi}}\right), \end{equation}

which highlights the role of the Maxwell stress in the induction equation. Guided by numerical simulations that show the appearance of a vertical magnetic field with variations in the $x$ direction, let us define an average

(8.13)\begin{equation} \bar{Q}=\iint \,\mathrm{d}y\,\mathrm{d}z Q \end{equation}

so that the induction equation for the ‘mean’ vertical field reads

(8.14)\begin{equation} \partial_t\overline{B_z}=\partial_x \left(\overline{w_z B_x-w_x B_z}\right)+\ell_{{H}} \partial_x^2 \frac{-\overline{B_xB_y}}{4{\rm \pi}}, \end{equation}

where we recognised the radial Maxwell stress term $\mathcal {M}_{xy}=-B_xB_y/4{\rm \pi}$, also present in the angular momentum conservation equation (4.21). This demonstrates that in Hall-MHD, the transport of magnetic flux is tightly linked to the transport of mass. Owing to energetic constraints, $\mathcal {M}_{xy} > 0$ (see § 5.2.4) in shear-driven instabilities/turbulence. Therefore, a concentration of magnetic field owing to the Hall effect is possible at local stress minimum. In a turbulent flow, short-lived stress minima occur randomly in the flow, and these minima tend to accumulate vertical magnetic flux according to the previous equation. When the Hall effect is strong enough, a local minimum can accumulate enough flux to become stable for the HSI. In this case, the flow becomes stable and the local turbulent stress vanishes, becoming a permanent minimum. This minimum continues to accumulate magnetic flux thanks to the remaining stress present on both sides until the flux outside of the minimum of stress becomes negative. At this point, the stress also vanishes in the regions $\overline {B_z} < 0$ and the systems settles down into a quasi-stationary state with very low-stress level (see also figure 33).

Figure 33. Hall self-organisation phenomenology. We start from a local fluctuation of the stress $\mathcal {M}_{xy}$ (a). This minimum creates a local maximum of $\overline {B_z}$, which overshoot the maximum $\overline {B_z}$ allowed by the HSI. Because the flow becomes locally stable, the stress vanishes (b). On the boundaries of this stable region, the Maxwell stress still transport magnetic flux towards the stable region, making it larger and emptying the rest of the domain (d). At some point, the total flux in the rest of the domain becomes negative, and it becomes HSI-stable. The stress therefore vanishes in this region as well, leaving only the interface with a minimal stress (c). Figure inspired by Kunz & Lesur (Reference Kunz and Lesur2013).

Although this process was initially identified in unstratified shearing box simulations, it was later unambiguously identified in cylindrical unstratified simulations (Béthune, Lesur & Ferreira Reference Béthune, Lesur and Ferreira2016). However, vertically stratified simulations do not seem to exhibit this process, for reasons not yet identified to date (Lesur et al. Reference Lesur, Kunz and Fromang2014).

8.5.3. Ambipolar diffusion

Self-organisation owing to ambipolar diffusion was first mentioned by Bai & Stone (Reference Bai and Stone2014) in unstratified shearing boxes. It was also observed in cylindrical simulations including ambipolar diffusion by Béthune et al. (Reference Béthune, Lesur and Ferreira2016), albeit at a very low level. However, whether ambipolar diffusion plays an active role in the self-organisation mechanism is an open question. Clearly, the ‘strength’ of self-organisation (quantified by the ratio $B_{z\mathrm {max}}/B_{z\mathrm {min}}$) in the ideal-MHD case is larger than the ambipolar diffusion case by a factor $10$ (figure 2 in Bai & Stone Reference Bai and Stone2014). In addition, the mechanism proposed for self-organisation only involves ideal-MHD terms (Bai & Stone Reference Bai and Stone2014). Finally, the existence of zonal flows in this configuration (non-stratified, ambipolar diffusion dominated) largely depends on the box aspect ratio. Bai & Stone (Reference Bai and Stone2014) largely explored the situation with $L_x=L_y=4 H$. However, a choice of box with $L_y>L_x$ tends to break zonal flows (figure 34). Overall, it is very possible that self-organisation in ambipolar-dominated unstratified shearing boxes is a mere numerical artefact.

Figure 34. The $y$–$z$ average of $B_z$ as a function of time for non-stratified MRI simulations with $\varLambda _A=3$ and $\beta _{\mathrm {mean}}=1000$: (a) $L_x=L_y=4H$; (b) $L_x=4H$, $L_y=8H$; and (c) $L_x=4H$, $L_y=16H$. Note the disparity in zonal flows when $L_y>L_x$. All three simulations have been computed with a fixed resolution per scale height in the three spatial directions $n_{x,y,z}=64\ \mathrm {pts}/H$ using the Snoopy code (Lesur & Longaretti Reference Lesur and Longaretti2007).

9.1.1. Ideal MHD

The zero mean field stratified case was first explored by Brandenburg et al. (Reference Brandenburg, Nordlund, Stein and Torkelsson1995) and Stone et al. (Reference Stone, Hawley, Gammie and Balbus1996). They first noted the spontaneous appearance of a ‘butterfly’ diagram when looking at the space–time behaviour of the $x$–$y$ average toroidal magnetic field (figure 35; see also Davis et al. Reference Davis, Stone and Pessah2010 and Simon, Beckwith & Armitage Reference Simon, Beckwith and Armitage2012). This diagram shows quasi-periodic flip of the toroidal field, with a periodicity close to 10 local orbital periods. Whether this butterfly diagram leads to observational counterpart by modulating the turbulent transport is still debated (see, e.g., Hogg & Reynolds Reference Hogg and Reynolds2016 for an example). On the theoretical side, it is well reproduced by mean field models (Gressel Reference Gressel2010; Gressel & Pessah Reference Gressel and Pessah2015), but we are still lacking a first principle theory for this dynamo.

Figure 35. The $x$–$y$ average of $B_y$ ($\overline {B_y}$) as a function of time and $z$ for zero vertical net flux stratified MRI simulations. Note the presence of a butterfly pattern indicating a periodic reversal of the toroidal field followed by its increase away from the disc midplane.

In principle, once vertical stratification is included, one needs to specify both the vertical density and temperature profile as well as an equation of state. However, most of the literature has focused on isothermal simulations, in which the temperature is constant. In this case, it is found that

(9.3)\begin{equation} \alpha\simeq 0.02\pm 0.01\quad\textrm{(Stratified, isothermal, zero vertical flux)} \end{equation}

(e.g. Stone et al. Reference Stone, Hawley, Gammie and Balbus1996; Davis et al. Reference Davis, Stone and Pessah2010). In the outer region of PPDs, radiative transfer calculations tell us that the vertical temperature profile is approximately constant, so the isothermal approximation is probably a good model in this case.

Nevertheless, some researchers have explored non-isothermal models, such as Hirose and collaborators (Hirose, Krolik & Stone Reference Hirose, Krolik and Stone2006; Hirose, Blaes & Krolik Reference Hirose, Blaes and Krolik2009; Hirose et al. Reference Hirose, Blaes, Krolik, Coleman and Sano2014), Flaig and collaborators (Flaig, Kley & Kissmann Reference Flaig, Kley and Kissmann2010; Flaig et al. Reference Flaig, Ruoff, Kley and Kissmann2012) and Bodo and collaborators (Bodo et al. Reference Bodo, Cattaneo, Mignone and Rossi2013, Reference Bodo, Cattaneo, Mignone and Rossi2015). It is found that when the vertical profile becomes unstable for convection (i.e. it violates the Schwarzschild criterion), the turbulent transport of angular momentum can increase by up to an order of magnitude (Hirose et al. Reference Hirose, Blaes, Krolik, Coleman and Sano2014), thanks to a mechanism which is yet to be fully elucidated.

Following the numerical convergence issue pointed out by Fromang & Papaloizou (Reference Fromang and Papaloizou2007) in unstratified simulations (see § 8.1), numerical convergence in stratified setups has also been tested. Initial explorations were limited by numerical resources to 128 points per scale height (Davis et al. Reference Davis, Stone and Pessah2010) and showed the convergence of $\alpha$ as a function of the number of grid points, which led to the conclusion that ‘stratification saves the day’. However, this issue was tackled again with higher-resolution simulations: 200 points per scale height (Bodo et al. Reference Bodo, Cattaneo, Mignone and Rossi2014) and 256 points per scale height (Ryan et al. Reference Ryan, Gammie, Fromang and Kestener2017). These recent results show a weak dependence on the resolution with $\alpha \propto N^{-1/3}$ at the resolution explored, albeit with a different vertical boundary condition (Davis et al. (Reference Davis, Stone and Pessah2010) used periodic boundary conditions while Ryan et al. (Reference Ryan, Gammie, Fromang and Kestener2017) used outflow boundary conditions). Overall, these results indicate that numerical convergence is also an issue in stratified models.

9.1.2. Non-ideal MHD

Realistic Ohmic diffusion profiles were introduced in zero net flux stratified simulations for the first time by Fleming & Stone (Reference Fleming and Stone2003). These simulations exemplified (Gammie Reference Gammie1996) layered accretion model with a mid-plane dead zone and an active surface layer. The debate then crystallised on the thickness of the active layer, which, depending on the model, could lead to accretion rates compatible with observations.

This layered accretion paradigm, however, omits ambipolar diffusion and the Hall effect. Ambipolar diffusion is particularly important at low densities, right in the active layers of Gammie (Reference Gammie1996). It was quickly realised that ambipolar diffusion could dramatically reduce the strength of MRI turbulence in this layer, and even suppress it, because $\varLambda _A\lesssim 1$ in this region (Perez-Becker & Chiang Reference Perez-Becker and Chiang2011a,Reference Perez-Becker and Chiangb; Dzyurkevich et al. Reference Dzyurkevich, Turner, Henning and Kley2013). This is confirmed by direct numerical simulations including Ohmic and ambipolar diffusion but no Hall effect (Bai & Stone Reference Bai and Stone2013b; Simon et al. Reference Simon, Bai, Stone, Armitage and Beckwith2013b). These simulations suggest dramatically low values for $\alpha$, with $\alpha \simeq 3\times 10^{-6}$ at 1 AU (Bai & Stone Reference Bai and Stone2013b) to $\alpha \simeq 10^{-3}$ at 100 AU (Simon et al. Reference Simon, Bai, Stone, Armitage and Beckwith2013b).

These results imply accretion rates $\dot {M}\lesssim 10^{-10}\, M_\odot /\mathrm {yr}$ in the region 1–30 AU, which are too low by at least two orders of magnitude compared with observational constraints. However, the zero mean field configuration is rather artificial. Indeed, PPDs are expected to form from a collapsing cloud that has dragged some fraction of its initial magnetic flux. Therefore, a non-negligible poloidal flux is likely to be present in these objects, which motivated the need for simulations including a mean poloidal field.

9.2.1. Ideal MHD

The first viable shearing-box simulation of a vertically stratified disc with a mean vertical field was performed by Suzuki & Inutsuka (Reference Suzuki and Inutsuka2009) and Suzuki, Muto & Inutsuka (Reference Suzuki, Muto and Inutsuka2010), thanks to more robust numerical techniques compared with the initial attempts of Stone et al. (Reference Stone, Hawley, Gammie and Balbus1996) and carefully designed vertical boundary conditions. Their simulations have a relatively weak field $\beta _{\mathrm {mean}}\gtrsim 10^4$, but show significant deviations from the zero mean field scenario and, in particular, the presence of a strong outflow, quantified by the mass flux leaving the disc. In a box with $L_z=8H$, they find a vertical mass flux 100 times larger with $\beta _{\mathrm {mean}}=10^4$ compared with the zero net field case (Suzuki & Inutsuka Reference Suzuki and Inutsuka2009).

This problem was revisited by a large number of authors, both in the strong field limit $\beta _{\mathrm {mean}}\simeq 1$ (Moll Reference Moll2012; Ogilvie Reference Ogilvie2012; Lesur, Ferreira & Ogilvie Reference Lesur, Ferreira and Ogilvie2013, figure 36), and in the intermediate regime $\beta _{\mathrm {mean}}=10^3\text {--}10^4$ (Fromang et al. Reference Fromang, Latter, Lesur and Ogilvie2013; Bai & Stone Reference Bai and Stone2013a). It was quickly demonstrated that these outflows were a variation of Blandford & Payne (Reference Blandford and Payne1982) magneto-centrifugal paradigm (Fromang et al. Reference Fromang, Latter, Lesur and Ogilvie2013; Lesur et al. Reference Lesur, Ferreira and Ogilvie2013) with a strong time-dependency. Despite this connection to well-known launching mechanisms, some properties of shearing box outflows are intrinsically flawed.

Figure 36. Growth and saturation of the MRI in the presence of a strong vertical field ($\beta _{\mathrm {mean}}=10$) in ideal MHD. Tubes are magnetic field lines whereas gas density is represented with volume rendering. Note the initial growth of the linear mode in the disc midplane, which eventually saturates into a quasi-laminar outflow configuration similar to Blandford & Payne (Reference Blandford and Payne1982) paradigm. From Lesur et al. (Reference Lesur, Ferreira and Ogilvie2013).

First, the mass loss rate depends on the vertical box extension (Fromang et al. Reference Fromang, Latter, Lesur and Ogilvie2013), taller boxes leading to lower mass loss rates. This surprising result is actually expected from the shape of the gravitational potential in a shearing box (§ 5.1), which is unbounded when $z\rightarrow \infty$. Some authors have proposed to fix this issue by adding higher-order terms to the vertical gravity force, making it possible to escape at $z=\pm \infty$ with a finite amount of mechanical energy (Suzuki et al. Reference Suzuki, Muto and Inutsuka2010; McNally & Pessah Reference McNally and Pessah2015). This approach, however, violates the conservative nature of gravitation, because the resulting force does not derive from a potential anymore. A more rational approach would be to include all of the third-order terms in Hill's approximation. This, however, also introduces radial curvature terms, making the shear-periodic boundary conditions unadapted.

Second, the geometry of the outflow is problematic. In principle, shearing boxes allow both vertically even and odd magnetic configuration (figure 37). At the linear level (see § 6.4.6), the two symmetries show similar growth rates and properties. However, in the non-linear regime, the odd symmetry is problematic for the physical interpretation of the outflow. By simply looking at the poloidal field topology (figure 37), it is clear that the odd symmetry connects the bottom side of the disc to the central star, whereas the top side is connected to $R\rightarrow +\infty$. Physically, it means that no angular momentum is actually extracted from the disc: angular momentum comes from an unknown source at $z\rightarrow -\infty$ and flows through the disc up to $z\rightarrow +\infty$. The angular momentum conservation equation (4.21) clearly describes this situation: in the odd symmetry case, one has $B_y B_z(z)=B_yB_z(-z)$, so the surface stresses on both sides of the disc cancel out and the disc is not accreting any material.

Figure 37. Symmetries of the outflow configuration allowed in a shearing box. Poloidal field lines are represented in green whereas the toroidal field component $B_y < 0$ is shown in blue and $B_y > 0$ in red: (a) odd symmetry configuration; (b) even symmetry configuration. The usual (Blandford & Payne Reference Blandford and Payne1982) picture corresponds to the even symmetry in a shearing box model.

Naturally, this situation is rather unrealistic and the only physical configuration is the even symmetry one (or at least, a configuration where $B_y B_z$ have a different sign on both sides of the disc). However, it turns out that shearing boxes tend to settle into the odd configuration when $\beta _{\mathrm {mean}}\lesssim 10^3$ (Salvesen et al. Reference Salvesen, Simon, Armitage and Begelman2016). The reason for this unexpected result is still debated. It is certainly related to the fact that the shearing box has ‘too many symmetries’, and does not differentiate $x\rightarrow \pm \infty$. Some authors have enforced the even symmetry manually (Lesur et al. Reference Lesur, Ferreira and Ogilvie2013), which leads to physical outflow configurations. However, this is not satisfactory as the flow symmetry should be enforced by the global geometry of the disc and its surrounding, which is not captured in shearing box models.

Despite these difficulties with outflows, it is possible to measure the turbulent stress in these simulation. A systematic exploration of shearing box models with $\beta _{\mathrm {mean}}\in [10,\infty ]$ (Salvesen et al. Reference Salvesen, Simon, Armitage and Begelman2016) shows that

(9.4)\begin{equation} \alpha=10.1 \beta_{\mathrm{mean}}^{-0.53} \quad\textrm{for}\ \beta_{\mathrm{mean}} < 10^5 \end{equation}

with $\alpha$ recovering the zero net flux value when $\beta _{\mathrm {mean}} > 10^5$. This scaling leads to $\alpha$ values about 3 times larger than the unstratified estimate (§ 8.1). As in the unstratified case, the level of turbulent stress depends on the magnetic Prandtl number, $\alpha$ increasing with increasing Pm (Fromang et al. Reference Fromang, Latter, Lesur and Ogilvie2013), so caution is still needed when using these scalings in phenomenological models.

The characterisation of winds in a shearing box is always a bit difficult because of the symmetry issues noted previously. In addition, the usual MHD outflow invariants (§ 11) are only defined in global geometry, making them inaccessible in local models (but see Lesur et al. (Reference Lesur, Ferreira and Ogilvie2013) for local equivalents of the global invariants). Despite these difficulties, one usually defines an outflow rate $\zeta$ and a torque parameter $\upsilon$, which measure the mass and angular momentum evacuated by the magnetised wind, respectively.

It is customary to define the outflow rate as

(9.5)\begin{equation} \zeta\equiv \frac{\overline{\rho w_z}|_\mathrm{top}-\overline{\rho w_z}|_\mathrm{bottom}}{\rho_\mathrm{mid}c_s}, \end{equation}

where ‘top’ and ‘bottom’ subscripts denote the top and bottom of the disc (which, in principle, can be chosen freely) and $\overline {\quad }$ denotes a horizontal averaging procedure. In a real (global) disc, the local domain is emptied by the wind, but it is also replenished by the divergence of the accretion flow, so that a steady state is locally achieve. As there is technically no accretion flow in shearing boxes, this replenishment is absent, leading to boxes that slowly lose mass. The typical timescale over which a box loses a significant fraction of it mass is given by $\tau _\mathrm {loss}\sim (\zeta \varOmega )^{-1}$. This implies that the shearing box model is really valid in the limit $\zeta \ll 1$.

The outflow rate is known to depend on several key parameters. First, it should be noted that the outflow rate depends on the box extension, both horizontal and vertical. In the horizontal plane, it seems that convergence is reached for $L_x,L_y\gtrsim 4H$ (Fromang et al. Reference Fromang, Latter, Lesur and Ogilvie2013). In the vertical direction, however, as the box gets higher, the outflow rate decreases. Fromang et al. (Reference Fromang, Latter, Lesur and Ogilvie2013) and Lesur et al. (Reference Lesur, Ferreira and Ogilvie2013) showed that doubling the vertical extension of the box divides $\zeta$ by a factor of three for $\beta _{\mathrm {mean}}=10^4$ and by a factor of 1.6 for $\beta _{\mathrm {mean}}=10$. In any case, this indicates that the outflow rate does not converge in shearing boxes. This is because it is not possible to escape to infinity in the gravitational potential of Hill's approximation (see the discussion in § 5.1). Hence, the ‘border’ of the gravitational potential is set artificially by the location of the physical boundary in the z direction of the shearing box.