Magnetohydrodynamics of protoplanetary discs

Abstract

Protoplanetary discs are made of gas and dust orbiting a young star. They are also the birth place of planetary systems, which motivates a large amount of observational and theoretical research. In these lecture notes, I present a review of the magnetic mechanisms applied to the outer regions ($R\gtrsim 1\ \mathrm {AU}$) of these discs, which are the planet-formation regions. In contrast to usual astrophysical plasmas, the gas in these regions is noticeably cold ($T < 300\ \mathrm {K}$) and dense, which implies a very low ionisation fraction close to the disc midplane. In these notes, I deliberately ignore the innermost $(R\sim 0.1\ \mathrm {AU})$ region, which is influenced by the star–disc interaction and various radiative effects. I start by presenting a short overview of the observational evidence for the dynamics of these objects. I then introduce the methods and approximations used to model these plasmas, including non-ideal magnetohydrodynamics, and the uncertainties associated with this approach. In this framework, I explain how the global dynamics of these discs is modelled, and I present a stability analysis of this plasma in the local approximation, introducing the non-ideal magneto-rotational instability. Following this mostly analytical part, I discuss numerical models that have been used to describe the saturation mechanisms of this instability, and the formation of large-scale structures by various saturation mechanisms. Finally, I show that local numerical models are insufficient because magnetised winds are also emitted from the surface of these objects. After a short introduction on wind physics, I present global models of protoplanetary discs, including both a large-scale wind and the non-ideal dynamics of the disc.

PART ONE: Observations and physical description

1. Observational context

Recent years have seen a dramatic change in our understanding of protoplanetary discs (PPDs), both from an observational and a theoretical point of view. Observations are now able to resolve the outer regions (radii larger than 1 astronomical unit [AU]) and show the existence of many unexpected features: spiral arms, rings and crescent-like structures. Although these observations mostly probe the distribution of dust grains, they indicate that the gaseous structure of PPDs is much more complex and rich than initially anticipated. In this part, we review the most recent evidence for PPD structure and evolution, which can be used to constrain the most recent theoretical models.

1.1. Observational diagnostics

Today observations probe different regions of the disc. In order to interpret these observations and constrain theoretical models, it is essential to clearly understand the quantities and limitations of each kind of observation. A typical PPD can be separated into two parts: an inner dust-free disc (from a few stellar radii to the dust sublimation radius) made of hot gas (typically more than 1000 K) and an outer disc of gas and dust (figure 1). The disc outer edge can range from 100 AU to more than 1000 AU depending on the object under consideration.

Figure 1. PPD diagram showing the various observational diagnostics. Disc winds have been omitted for clarity.

Observations typically probe the following regions:

1. (i) The UV excess is a signature of the accretion shock at the foot of accretion columns. It is very often the only way to deduce the accretion rate in a specific disc.

2. (ii) The near and mid-infrared continuum (also known as infrared excess) is a result of stellar photons scattered by small dust grains (typically less than $1~\mathrm {\mu }\mathrm {m}$ in size). Scattered light probes the very surface of the dust layer as the dust disc is very optically thick at these wavelengths. For this reason, the intensity of scattered light is not related to the column density but to the amount of stellar light received by the layer. It therefore characterizes the disc geometry.

3. (iii) The (sub-)millimetre continuum probes the thermal emission of bigger dust grains (typically with a size of the order of 1 mm). If the dust layer is optically thin at these wavelengths (as usually assumed), the emissivity is related to the column density of dust, but also to its temperature.

4. (iv) Spectral lines, both in the infrared and at radio wavelengths, probe specific gas tracers such as gas or molecular transitions. These lines are usually optically thick, which implies that they only probe the surface of the gas layer. For this reason, direct estimates of the gas mass in the disc is very difficult, and one has to rely on proxies.

These observational properties are then used to derive several useful dynamical quantities.

1.2. Accretion

Because the thermal equilibrium of PPDs for $R\gtrsim 1\ \mathrm {AU}$ is dominated by the illumination of the central star (D'Alessio et al. 1998), a direct measurement of the accretion rate through viscous heating is not possible. For this reason, observational evidence of accretion in these regions are scarce and plagued by uncertainties. There are mainly two classes of accretion signature, which are all indirect.

The first is the observational signature of accretion columns at the stellar surface. These accretion columns are formed when the disc material is lifted and accreted by the stellar magnetic field. The gas then ends up in a nearly free-fall speed and hits the stellar surface, forming an accretion shock. The luminosity of this accretion shock observed in UV bands is directly related to the accretion rate in the accretion columns and, therefore, in the innermost disc. It should be kept in mind that accretion rates deduced by this method are not necessarily accretion rates in the entire disc, which can in principle vary with radius if the disc is not in steady state, or if the disc is losing mass from a wind. Typical results show accretion rates of the order of $10^{-8}\,M_\odot /\mathrm {year}$ with uncertainties of the order of an order of magnitude depending on the object under consideration¹ (e.g.figure 2a). These accretion rates tend to decrease over timescales of a few million years.

Figure 2. (a) Measurement of the accretion rate as a function of stellar age in NGC 2264 using the excess UV due to accretion columns. From Venuti et al. (2014). (b) Fraction of disc signature (accretion) and dust signature (infrared excess) as a function of the cluster age. Both show that discs have an average lifetime of a few million years. From Fedele et al. (2010).

The second observational evidence lies in the proportion of stars showing disc features (accretion on the stellar surface, or infrared excess signifying the presence of dust around the star) as a function of the stellar age. The disappearance of these signatures in older stars allows one to evaluate the typical gas and dusty disc lifetimes. These two time scales do not necessarily match as the gas disc could, for instance, disappear before the dusty disc. However, they both show the same trend: disc tends to disappear on a timescale of a few million years (figure 2b).

By combining this information, and assuming that accretion is approximately constant during the lifetime of these objects, one deduces that typical PPD masses range from $10^{-3}\,M_\odot$ to $10^{-1}\,M_\odot$, which is consistent with mass inferred from the total dust content of the disc (Andrews et al. 2013).

1.3. Ejection: winds and jets

PPDs are often observed in association with large-scale winds and jets. Jets are often seen in forbidden emission lines and correspond to fast collimated flow ($v > 100\ \mathrm {km}\,\mathrm {s}^{-1}$). Their high velocity suggests they are launched from the inner few AU of the disc (Frank et al. 2014). The typical outflow rate is estimated to be of the order of 10 % of the accretion rate in classical T-tauri stars (Frank et al. 2014).

In addition to these jets, a slower component is also observed in molecular lines. This ‘molecular outflow’ is denser and reach velocities $v\sim 1\text {--}10 \ \mathrm {km}\,\mathrm {s}^{-1}$ (figure 3). They could be a result of the interaction of the jet with its environment, or they could be a genuine outflow component, emitted from the disc at $R\gtrsim 1\ \mathrm {AU}$.

Figure 3. Observation of a disc and an atomic jet seen by the Hubble Space Telescope (Burrows et al. 1996) and a molecular wind observed in CO(2-1) by ALMA (Louvet et al. 2018) in HH30, a PPD seen edge-on. Courtesy of F.Louvet.

1.4. Structures

The progress in observational techniques (adaptative optics, interferometry) now allows astronomers to resolve the disc and look for signatures of planet formation, accretion or other unexpected processes. The first class of observations relies on polarimetric differential imaging (PDI) of scattered light emission in the near infrared. This technique allows one to obtain only the light scattered by dust grains (which is naturally polarised) and not the light of the central object. They have been used to probe the disc surface of various disc (mainly transitional discs). Stunning structures such as spiral and rings were found² in several objects (figure 4).

Figure 4. Scattered light images in the near infrared using PDI: (a) spiral structures observed in MWC758, from Benisty et al. (2015); (b) multiple ring structures observed in HD97048, from Ginski et al. (2016).

The second class of observation is based on interferometry at millimetric and sub-millimetric wavelengths. The ALMA observatory has been very successful at probing the very structure of PPDs with incredible resolution and unexpected results (figure 5, Andrews et al. 2018).

Figure 5. (a) Ring-like structures observed in TW Hydra. From Andrews et al. (2016). (b) Multiple ring structure in a deprojected image of HL-Tau from Partnership et al. (2015). (c) Horsehoe-like structure observed in Oph IRS 48 at sub-millimetre wavelengths (green, tracing millimetre-sized dust) and corresponding scattered light infrared emission (yellow, tracing $\mathrm {\mu }\mathrm {m}$ size dust) from van der Marel et al. (2013). (d) Spiral structures seen at sub-millimetre wavelengths in the young and massive disc of Elias 2-27, from Pérez et al. (2016).

Although these observations probe the dust distribution in the disc, they also tell us about the gas distribution and dynamics, because the grains that are observed are tightly coupled to the gas through a drag force. Such a direct connection has been recently confirmed observationally by simultaneously looking at the continuum (dust) and line emissions (probing the gas kinematics) (Teague, Bae & Bergin 2019).

All these observations indicate that discs are not smooth and symmetrical. They are instead structured on length scales comparable with our solar system. Structures are categorised in spirals, rings and horseshoes, which can be associated with specific physical processes in the disc. It should be noted some of these structures are found in transitional discs, i.e. truncated discs that are presumably in the final evolution stage of PPDs. All of these structures could be the signature of embedded planets perturbing the disc structure by gravitational interaction. However, other processes have been proposed that do not assume planets. One of the key questions is, therefore, whether or not these structures are necessarily a signature of embedded planets.

1.5. Turbulence

Turbulence is likely one of the key elements of any dynamical theory for the evolution of discs. Theoretical arguments (see § 4.4) show that turbulence should be subsonic in these systems, i.e., that chaotic motions of the gas are slower than the sound speed. This implies that turbulence is difficult to detect because the turbulent broadening of spectral lines is comparable with the thermal spreading of the molecules constituting the gas. For this reason, heavy molecules such as CN and CO tend to be preferred to detect turbulence, because their thermal velocity is lower compared with lighter molecules at a given equilibrium temperature. High-resolution spectra obtained from ALMA for CO lines indicates that turbulence is very weak, or non-existent (Flaherty et al. 2015, 2017). Spectral broadening smaller than 3 % of the local sound speed are found as best fits to observational data at large distances (typically more than 30 AU). This turbulent broadening is way smaller than the typical values expected from ideal magnetohydrodynamics (MHD) turbulence which typically predicts $\delta v\gtrsim 0.1 c_s$.

Another signature of turbulence (or, more precisely, the lack of turbulence) lies in the dust vertical distribution. Indeed, dust grains naturally tend to settle towards the midplane, unless turbulence stirs them up into the disc atmosphere. Direct measurements of the thickness of the dust layer allow one to deduce the level of hydrodynamical turbulence in the disc. Such a measurement has been done in the case of HL-tau, where the thickness of the rings is used as a tracer for the disc thickness (Pinte et al. 2016). The result is that $100\,\mathrm {\mu } \mathrm {m}$ grains have settled towards the midplane, with a vertical dust scale height about 10 times smaller than the gas scale height. This implies a very low level of turbulence in the disc, with typically $\delta v\sim 10^{-2}\,c_s$ ($\alpha \sim 10^{-4}$, see § 4.4).

1.6. Magnetic fields

Evidence for magnetic fields in PPDs is scarce. Typical values are expected to be of the order of a Gauss at 1 AU down to a few milli-Gauss at a few tens of astronomical units (Wardle 2007), although these theoretical values could vary by several orders of magnitude. For this reason, measurement through Zeeman effect is unfeasible except in the very inner disc. In this region, toroidal magnetic fields of a few kilo-Gauss have been measured, although it is not clear whether this field belongs to the host star or to the disc itself (Donati et al. 2005). At larger distances (tens of astronomical units), attempts at measuring the field strength through Zeeman splitting in molecular lines have only led to upper limits, with $B_z < 0.8\ \mathrm {mG}$ and $B < 30\ \mathrm {mG}$ (Vlemmings et al. 2019).

Topological information on the field is also accessible through polarisation in the continuum (i.e. dust thermal emission). It is assumed that dust grains tend to align perpendicularly to magnetic field lines, thereby emitting thermal radiation with a preferred polarisation, perpendicular to the local field orientation (Cho & Lazarian 2007; Stephens et al. 2014). However, polarisation in sub-millimetric radiation can also be due to self-scattering by dust grains (Kataoka et al. 2015; Yang et al. 2016). Campaigns using multiple wavelengths observations have attempted to disentangle these two effects (Stephens et al. 2017), but the interpretation of the results in terms of magnetic topology remains very uncertain.

Finally, magnetic field intensities can be deduced from meteoritic and cometary evidence in our own solar system, assuming that the field gets frozen in the body during its formation in the parent disc. Field strength of the order of $0.1\ \mathrm {G}$ around 1 AU are inferred from remnant magnetisation in meteorites following this idea (Fu et al. 2014), whereas upper limits with $B < 30\ \mathrm {mG}$ in the region around 15–45 AU is deduced from the magnetisation of Comet 67P/Churyumov-Gerasimenko (Biersteker et al. 2019).

2. Disc prototype

2.1. Fluid properties

PPDs are rather cold objects, with temperatures ranging from 1000 K in the inner ($0.1$ AU) disc down to 10 K in the outer (100 AU) disc. In order to characterise these discs, It is important to quantify the typical length scales and time scales relevant to the problem. Let us start with a typical disc model which matches disc observations (Andrews et al. 2009):

(2.1)\begin{equation} \left.\begin{gathered} \varSigma=300\,R_{\mathrm{AU}}^{-1}\ \mathrm{g}\,\mathrm{cm}^{-2},\\ T=280\,R_{\mathrm{AU}}^{-1/2}\ \mathrm{K},\\ \varOmega=2\times 10^{-7}\,R_{\mathrm{AU}}^{-3/2}\,\mathrm{s}^{-1}. \end{gathered}\right\} \end{equation}

Here, we have defined the main physical properties of a disc: its surface density $\varSigma$, which correspond to the usual mass density integrated in the direction perpendicular to the disc plane, its temperature $T$, and its angular velocity $\varOmega$ around the central object. We also define for convenience a dimensionless distance from the central object, in astronomical units: $R_{\mathrm {AU}}\equiv R/1\ \mathrm {AU}$.

This simple model leads to a $0.04\,M_\odot$ mass disc, extending from 0.07 to 200 AU, rotating around a solar mass star, typical of discs which have been observed. We can deduce some useful dynamical parameters associated from this simplified models. Defining the isothermal sound speed as $c_s\equiv \sqrt {P/\rho }$ and using the vertical hydrostatic equilibrium to define the disc vertical scale height (§ 4.2) $H=c_s/\varOmega$, one obtains

(2.2)\begin{equation} \left.\begin{gathered} c_s=10^5\,R_{\mathrm{AU}}^{-1/4}\,\mathrm{cm}\,\mathrm{s}^{-1},\\ H=5\times 10^{11}\,R_{\mathrm{AU}}^{5/4}\ \mathrm{cm},\\ \dfrac{H}{R}=0.03\,R_{\mathrm{AU}}^{1/4},\\ \rho_\mathrm{mid}=6\times 10^{-10}\,R_{\mathrm{AU}}^{-9/4}\ \mathrm{g}\,\mathrm{cm}^{-3},\\ n_\mathrm{mid}=1.5\times 10^{14}\,R_{\mathrm{AU}}^{-9/4}\ \mathrm{cm}^{-3},\\ P_\mathrm{mid}=6\,R_{\mathrm{AU}}^{-11/4}\ \mathrm{dyn}\,\mathrm{cm}^{-2}. \end{gathered}\right\} \end{equation}

2.2. Magnetic fields

Magnetic fields in PPDs are poorly constrained (§ 1.6). It is widely believed that fields are largely sub-thermal: the thermal pressure of the fluid dominates over the magnetic pressure (this requirement follows from the fact that the discs are approximately in Keplerian rotation). This translates into a plasma $\beta$ parameter

(2.3)\begin{align} \beta&\equiv\frac{P_{\mathrm{th}}}{P_{\mathrm{mag}}}\nonumber\\ &=\frac{8{\rm \pi} P}{B^2}\gg 1. \end{align}

In practice, $\beta \simeq 1$ constitutes a lower limit for the MRI to operate in geometrically thin discs (see § 6.4.6). Note also that if dynamo action is generating a field (both ordered or disordered), then $\beta$ cannot reach a value lower than $\beta \sim 1$, hence this value is actually a lower limit for the typical plasma $\beta$ expected in these discs. It is possible to connect the field strength to the plasma $\beta$ using the properties defined previously and obtain

(2.4)\begin{equation} B=12\,R_{\mathrm{AU}}^{-11/8}\beta^{-1/2} \ \mathrm{G}. \end{equation}

The upper bound $B\lesssim 10\ \mathrm {mG}$ for $R\sim 10\ \mathrm {AU}$ mentioned in § 1.6 tend to suggest $\beta \gtrsim 10^4$ in these regions, which confirms that the field strength is expected to be strongly sub-thermal.

2.3. Fluid approximation

PPDs are mostly constituted of neutral gas. In order to describe this gas, it is tempting to use the fluid approximation. For this approximation to be valid, the gas under consideration needs to be collisional, i.e. gas particles need to be subject to many collisions during one dynamical timescale. This ensures that at the microphysical level, the velocity distribution of the gas phase can be approximated by a Maxwellian distribution, allowing us to use a scalar pressure field.

Assuming the gas is mainly made of $H_2$ molecules of radius $10^{-8}\ \mathrm {cm}$, we can estimate the cross section of neutral molecules as $\sigma _{nn}=3\times 10^{-16}\ \mathrm {cm}^2$. This gives us an approximate mean free path $\ell _\mathrm {mfp}$ and collision frequency $\omega _\mathrm {coll}$

(2.5)\begin{equation} \left.\begin{gathered} \ell_\mathrm{mfp}=\dfrac{1}{n\sigma_{nn}}=22\,R_{\mathrm{AU}}^{9/4}\ \mathrm{cm},\\ \omega_\mathrm{coll}=\dfrac{v_\mathrm{th}}{\ell_\mathrm{mfp}}=5\times 10^3\,R_{\mathrm{AU}}^{-5/2}\ \mathrm{s}^{-1}. \end{gathered}\right\} \end{equation}

We therefore have $\ell _\mathrm {mfp}\ll R$ and $\omega _\mathrm {coll}\gg \varOmega$, which validate the fluid approximation to describe the dynamics of PPDs to a very good approximation. It should be noted that these quantities are evaluated at the disc midplane. If one looks at regions well above the disc, as in the case of outflows, $\ell _\mathrm {mfp}$ increases significantly. One finds that $\ell _\mathrm {mfp}\gtrsim H$ when $n\lesssim 10^4\ \mathrm {cm}^{-3}$, i.e. when the atmosphere is $10^{10}$ times less dense than the midplane at 1 AU. Such a strong density contrast is almost never reached in outflow models, where one finds density contrasts between $10^4$ and $10^7$ (e.g. figure 46). Nevertheless, it should be kept in mind that very weak outflows in the outermost parts of the disc can be close to the collisionless regime.

2.4. Grain population

The question of grains is of importance in PPDs. As is usually assumed, we consider a constant dust to gas mass fraction, equal to that of the interstellar medium (1/100). We further assume that grains are spherical with a radius $a$ and made of olivine with a density $\rho _o=3\ \mathrm {g}\,\mathrm {cm}^{-3}$. The density of grains is therefore

(2.6)\begin{equation} \left.\begin{gathered} \rho_\mathrm{grain}=6\times 10^{-12}\,R_{\mathrm{AU}}^{-9/4}\ \mathrm{g}\,\mathrm{cm}^{-3},\\ n_\mathrm{grain}=1.4\,R_{\mathrm{AU}}^{-9/4}a_{\mathrm{\mu}\mathrm{m}}^{-3}\,\mathrm{cm}^{-3}. \end{gathered}\right\} \end{equation}

In this last estimate, we have assumed that all the grains had the same size. This is an over-estimation because the sizes are actually distributed over a wide range of scales. In addition, the grain size distribution is expected to evolve with time as grains are known to be growing in PPDs. However, this order of magnitude estimate points to an important fact: the abundance of grains $n_{\mathrm {grain}}/n\sim 10^{-14}a_{\mathrm {\mu }\mathrm {m}}^{-3}$. Hence, if grains are smaller than $1\ \mathrm {\mu }\mathrm {m}$, the typical ionisation fraction of PPDs ($10^{-14}$) suggest that grains are more abundant than free charge carriers. As we show in § 3.4.3, this has a huge effect on the plasma conductivity tensor as grains can become the main charge carriers.

2.5. Ionisation fraction

The ionisation fraction $\xi \equiv n_{-}/n_n$, where $n_{-}$ is the number of free negative charge carriers, is a highly uncertain quantity, with very little constraints coming from observations. The ionisation fraction typically range from $10^{-16}\text {--}10^{-13}$ at 1 AU to $10^{-13}\text {--}10^{-10}$ at 100 AU. However, the resulting plasma is not necessarily a plasma made of electrons and molecular ions. Indeed, if dust grains are present and sufficiently abundant, they tend to suck electrons and ions in the gas phase, leading to a plasma made of positively and negatively charged grains (Sano et al. 2000).

Here, we illustrate how each physical process affects the ionisation fraction by considering a simple chemical network which includes singly charged grains. We combine this network with ionisation rate prescriptions for the various ionisation sources (X-rays, UV, cosmic rays (CRs) and radioactive decay).

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 1999; Bai & Goodman 2009, 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 1981) 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 2009).

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 2013). 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. 2018). 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. 2015), though this is still highly model dependent. Owing to these uncertainties, some authors (e.g.Ilgner & Nelson 2006) 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 (1974), Fromang, Terquem & Balbus (2002) and Ilgner & Nelson (2006) 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. (2002), 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. 2002).

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 2007; Bai & Stone 2013b; Lesur, Kunz & Fromang 2014). 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 (1979) and Umebayashi & Nakano (1980) in the context of molecular clouds, and later applied to PPDs (Sano et al. 2000; Ilgner & Nelson 2006; Wardle 2007), 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 2012). 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 2007).

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.³

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. Plasma description in PPDs

In this section, we explore the properties of the plasma constituting PPDs and ask whether they can be described using non-ideal MHD. For this limit to be valid, we have to satisfy the following three criteria.

1. (i) Binary Coulomb interactions should be negligible. This implies that the plasma parameter (defined in the following) is much larger than one.

2. (ii) Electro-neutrality is satisfied on timescales of interest, i.e. any charge separation is quickly eliminated by electrostatic interactions.

3. (iii) The behaviour of each fluid component (electrons, ions, neutrals, charged grains) can be described using a single-fluid approximation.

3.1. Plasma parameter

Several quantities allow one to characterise a plasma, the first being the nature of the electromagnetic interaction. The most fundamental quantity characterising a plasma is the Debye length that may be written in an electron–ion plasma

(3.1)\begin{align} \lambda_\mathrm{D}&\equiv\sqrt{\frac{k_BT_e}{4{\rm \pi} (1+Z) n_e e^2}}\nonumber\\ &=30\left(\frac{\xi}{10^{-13}}\right)^{-1/2} R_{\mathrm{AU}}^{7/8}(1+Z)^{-1/2}\ \mathrm{cm}, \end{align}

where $Z$ is the averaged number of charges on the ions and we have assumed electro-neutrality so that $n_i=n_e/Z$. The Debye length is clearly below the scales of interest in PPDs. Even if one considers charged grains, the same Debye length can be derived because it does not depend on the particle mass. In addition to this characteristic length, a ‘good’ plasma should have many particles in a Debye sphere, ensuring the screening of short-range Coulomb interaction. This is quantified by the plasma parameter $\varUpsilon$, equal to the number of charge carriers in a Debye sphere

(3.2)\begin{align} \varUpsilon &= 4{\rm \pi} n_e\lambda_\mathrm{D}^3 \nonumber\\ &=4\times 10^5 (1+Z)^{-3/2}R_{\mathrm{AU}}^{3/8}\left(\frac{\xi}{10^{-13}}\right)^{-1/2}. \end{align}

Hence, despite the low ionisation fraction and low temperature of these objects, they are still very much in the plasma regime where short-range Coulomb interactions can be neglected. Note, however, that reducing the ionisation fraction and, at the same time, increasing the number of charges could change this picture, breaking the plasma approximation altogether. However, this would require $Z\gtrsim 10^3$ in PPDs, a value that is never encountered, even in chemical models including grains (e.g. Wardle 2007).

3.2. Electro-neutrality and drag

PPDs are weakly ionised objects. This implies that the dynamical equations describing the flow and the approximations underlying their derivation should be clearly stated. In this section, we derive these equations, starting from the multi-fluid plasma description. We assume the gas is made of neutral and charged ‘particles’ (particle could mean electron, ion, or charged grain, indifferently). The multi-fluid approximation is valid because the collision timescales are short, as demonstrated previously. We therefore start from the following dynamical equations:

(3.3)\begin{gather} \frac{\partial n_j}{\partial t}+\boldsymbol{\nabla}\boldsymbol{\cdot}n_j\boldsymbol{v}_j=0, \end{gather}

(3.4)\begin{gather}\frac{\partial n_jm_j \boldsymbol{v}_j}{\partial t}+\boldsymbol{\nabla}\boldsymbol{\cdot } (n_jm_j\boldsymbol{v}_j\otimes\boldsymbol{v}_j)=-\boldsymbol{\nabla}P_j+\boldsymbol{f}_j+ n_jq_j\left(\frac{\boldsymbol{v}_j\boldsymbol{\times}\boldsymbol{B}}{c}+\boldsymbol{E}\right)+\boldsymbol{R}_j, \end{gather}

where $n_j$, $m_j$, $\boldsymbol {v}_j$, $P_j$, $q_j$, and $\boldsymbol {f}_j$ are the number density, mass, velocity, pressure, charge, and additional forces (gravity, etc.) on species $j$. We have also included a drag force $\boldsymbol {R}_j$ between this species and all of the other species. This force is a result of to inter-species collisions and can be written as

(3.5)\begin{equation} \boldsymbol{R}_j=\sum_k \gamma_{jk}\rho_j\rho_k(\boldsymbol{v}_k-\boldsymbol{v}_j), \end{equation}

because each fluid component is collisional and, therefore, has a Maxwellian velocity distribution. Here, $\gamma =\langle \sigma v\rangle _{jk}/(m_j+m_k)$ and $\langle \sigma v\rangle _{jk}$ is the momentum exchange rate between species $j$ and $k$. As expected from momentum conservation, we have $\sum _j \boldsymbol {R}_j=0$.

It is usually assumed that electro-neutrality follows from the fact that the plasma frequency $\omega _p$ is much larger than any frequency of interest. Although this is indeed a good criterion for a fully ionised plasma, it is not necessarily true for a weakly ionised plasma. Let us therefore revisit this criterion, starting from the linearised multi-fluid equations. We perturb only one species along the $x$-axis, leaving the others unperturbed. We moreover assume that the fluid pressure and other external forces are negligible compared with electromagnetic forces. The linearised equation of motion reads

(3.6)\begin{equation} \left.\begin{gathered} \dfrac{\partial \delta n}{\partial t}+n_0\partial_x v_x =0,\\ n_0m\dfrac{\partial v_x}{\partial t}=n_0qE_x-\gamma m n_0\rho v_x. \end{gathered}\right\} \end{equation}

Solving these equations requires an equation for the electric field, which is obtained from one of Maxwell's equation

(3.7)\begin{equation} \partial_x E=4{\rm \pi} q \delta n. \end{equation}

We can combine these equations to obtain a second-order relation on the density fluctuation

(3.8)\begin{equation} -\frac{\partial^2\delta n}{\partial t^2}=\omega_p^2\delta n+\frac{1}{\tau_s}\frac{\partial \delta n}{\partial t}, \end{equation}

where we have introduced the plasma frequency $\omega _p$ and the stopping time $\tau _s$ as

(3.9)\begin{equation} \left.\begin{gathered} \omega_p\equiv \left(\dfrac{4{\rm \pi} n q^2}{m}\right)^{1/2},\\ \tau_s\equiv \dfrac{1}{\gamma \rho}. \end{gathered}\right\} \end{equation}

Dynamically, this equation describes damped plasma oscillations with frequencies

(3.10)\begin{equation} \omega_{\pm}=\frac{\textrm{i}\tau_s^{-1}\pm\sqrt{4\omega_p^2-\tau_s^{-2}}}{2}, \end{equation}

for which we can distinguish two physical limits.

1. (i) $\omega _p\gg \tau _s^{-1}$ in which case the plasma is subject to plasma oscillations at frequency $\omega _p$ with a damping timescale equal to $\tau _s$. If we consider phenomena on frequencies much lower than $\omega _p$, we can average out the highest-order time derivative and obtain a simple closure relation between $v_x$ and $E_x$: $v_x=qE_x/\gamma m \rho$, which constitutes the base of Ohm's law. Once these oscillations are time-averaged, the plasma can be assumed to be electrically neutral.

2. (ii) $\omega _p\ll \tau _s^{-1}$ in which case the plasma is subject to over-damped oscillations with two imaginary frequencies $\omega _+=\textrm {i}\tau _s^{-1}$ and $\omega _-=\textrm {i}\omega _p^2\tau _s\ll \omega _+$ associated with two damping timescales $\tau _\pm =(\omega _\pm )^{-1}$. To interpret physically these timescales, let us consider a plasma at rest in which we introduce a localised charge deficit. First, the plasma is going to start moving to ‘fill’ the charge deficit. Owing to the drag, however, it very rapidly reaches an asymptotic velocity, given by $v_x=qE_x/\gamma m \rho$. Here $\tau _+$ is the time needed by the system to be put in motion and reach this quasi-stationary velocity. This velocity fluctuation, however, is smaller than that which would be obtained in a pure plasma oscillation, because the drag prevents the plasma from reaching high velocities. Hence, it takes a time $\tau _-$ to actually fill the charge deficit. This implies that Ohm's law, given by the asymptotic velocity, is valid on timescales longer than $\tau _+$, and that charge inertia can be neglected in that limit. However, charge neutrality is restored on the much longer timescale $\tau _-$.

To summarise, it is possible to neglect inertia for the charged species in the momentum equation provided that the timescales under consideration are larger than $\mathrm {max}(\tau _s,\omega _p^{-1})$, and recover Ohm's law without time derivative. Note that this condition is different from electroneutrality, which requires timescales longer than $\mathrm {max}(\omega _p^{-1},(\omega _p^2\tau _s)^{-1})$, which are significantly longer than $\tau _s$ when $\omega _p\tau _s\ll 1$. It should be pointed out that this analysis was done for a single species, whereas plasmas in PPDs can be made of many different species. Hence, the condition for electroneutrality needs to be satisfied only by the most mobile species of the plasma, which can then compensate for charge fluctuations, and not necessarily by all of the species present.

In PPDs, we obtain the following values for the plasma frequency, depending on the type of charge carrier

(3.11)\begin{equation} \left.\begin{gathered} \omega_{p,e}=2.2\times 10^{5} \left(\dfrac{\xi}{10^{-13}}\right)^{1/2}R_{\mathrm{AU}}^{-9/8} \ \mathrm{s}^{-1},\\ \omega_{p,i}=9.3\times 10^{2} \left(\dfrac{\xi}{10^{-13}}\right)^{1/2}R_{\mathrm{AU}}^{-9/8} \ \mathrm{s}^{-1},\\ \omega_{p,g}=1.8\times 10^{-3} \left(\dfrac{\xi}{10^{-13}}\right)^{1/2}R_{\mathrm{AU}}^{-9/8}a_{\mathrm{\mu}\mathrm{m}}^{-3/2}\ \mathrm{s}^{-1}, \end{gathered}\right\} \end{equation}

where $e$, $i$ and $g$ denote electrons, ions and grains, respectively. As can be seen, this frequency is always short compared with the timescales of interest, but grains tend to have significantly lower frequencies owing to their higher inertia.

The stopping times can be estimated starting from the momentum exchange rates $\langle \sigma v\rangle _{ij}$. As we are interested only in weakly ionised plasmas, collisions between charged species will be extremely rare. We therefore only consider neutral-charge collisions.

The ‘collision’ between electron/ions and neutrals are mainly a result of the electrostatic interaction between the approaching charge and the dipole induced on the neutral by the charge. This is estimated by

(3.12)\begin{equation} \left.\begin{gathered} \langle \sigma v\rangle_e=8.3\times 10^{-9}\times \mathrm{max} \left[1,\left(\dfrac{T}{100\ \textrm{K}}\right)^{1/2}\right]\mathrm{cm}^3\,\textrm{s}^{-1},\\ \langle \sigma v\rangle_i=2.4\times 10^{-9}\left(m_H/m_n\right)^{1/2}\ \mathrm{cm}^3\, \mathrm{s}^{-1}, \end{gathered}\right\} \end{equation}

where $\langle \sigma v\rangle _e$ is deduced from Draine, Roberge & Dalgarno (1983) and $\langle \sigma v\rangle _i$ is obtained from Draine (2011), following Bai (2011a).⁴ For grains above a size of a few $10^{-2}\ \mathrm {\mu }\mathrm {m}$, collisions mainly behave as billiard balls. In other words, $\sigma v$ is roughly equal to the velocity of the incident neutral times the cross-section of the grain. For spherical grains, this leads to⁵

(3.13)\begin{align} \langle \sigma v\rangle_g&={\rm \pi} a^2 \sqrt{\frac{2k_B T}{m_n}}\nonumber\\ &=2.6\times 10^{-3} a_{\mathrm{\mu}\mathrm{m}} \left(\frac{T}{100\ \textrm{K}}\right)^{1/2}\ \mathrm{cm}^3\, \mathrm{s}^{-1}. \end{align}

These rates allow us to compute stopping times for each species following the previous definition:

(3.14)\begin{equation} \left.\begin{gathered} \tau_{s,e}=6.7\times 10^{-7}\,R_{\mathrm{AU}}^{9/4}\ \mathrm{s},\\ \tau_{s,i}=4.9\times 10^{-5}\,R_{\mathrm{AU}}^{9/4}\ \mathrm{s},\\ \tau_{s,g}=8.1\times 10^4\,R_{\mathrm{AU}}^{9/4}a_{\mathrm{\mu}\mathrm{m}}^2\left(\dfrac{100\ \textrm{K}}{T}\right)\mathrm{s} , \end{gathered}\right\} \end{equation}

which shows that because of the low ionisation fraction and the neutral drag, $\omega _{p,j}\tau _{s,j}<1$ for ions and electrons, whereas it is greater than one for grains. This means that plasma oscillations are over-damped for ions and electrons (case (ii)) and are not directly relevant for quasi-neutrality. Nevertheless, $\omega _p\tau _s>10^{-2}$, so even in this case, electroneutrality is recovered on timescales shorter than a second. Grains, on the other hand, are usually in regime (i), with a relatively low plasma frequency (period of a few days for $1\ \mathrm {\mu }\textrm {m}$ size grains), decreasing rapidly with increasing grain size. Grains are usually not the only charge carrier in discs, so electro-neutrality is guaranteed by ions and electrons, but it should be kept in mind that, in a hypothetical situation where grains would be the only charge carrier, electro-neutrality could be violated, leading to phenomena similar to lightning. This, however, is not explored here, and we only consider situations where ions and electrons are still present in the system.

3.3. Single-fluid approximation

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 2014). 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 2008), 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 2016). 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.⁶ 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 1999). 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. Non-ideal diffusivities

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 1996). 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 2006). 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 2006; Salmeron & Wardle 2008). 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 2011b).

3.4.4. Conclusion on non-ideal MHD effects

Overall, there is no general consensus on the quantitative strength of non-ideal MHD effects in the outer part ($R > 1\ \mathrm {AU}$) of PPDs. It is clear that these effects are qualitatively very important though, and that these objects are far from the ideal MHD regime. Let us summarise here the source of uncertainty and their implication for the strength of non-ideal effects.

Ionisation rate: Because CRs can be both shielded by the stellar wind or ‘locally’ produced in shocks surrounding the forming star, there is tremendous uncertainty of six orders of magnitude on the ionisation rate due to CRs (see the discussion in § 2.5.1). As this mechanism is the main source of ionisation in the disc below two scale heights, this implies a three orders of magnitude uncertainty in the ionisation fraction $\xi$ and similarly on the diffusivity coefficients. Ionisation as a result of X-rays is also subject to caution because the X-ray flux coming from the star is largely variable, leading to order of magnitude fluctuations of the ionisation fraction close to the disc surface.

Disc structure: The disc structure is a fundamental parameter that determines the penetration depth of ionising radiations, but also the recombination rate. Denser disc models, such as the MMSN, tend to have lower ionisation fractions and larger diffusion. We have shown that by comparing a theoretical MMSN disc model with a model favoured by observations (§ 2.1), one can change the ionisation fraction by two orders of magnitude (see § 2.5.3). As the gas column density profile is largely unknown for $R\sim 10\ \mathrm {AU}$, one is forced to use the gas column density as a free parameter.

Grains: As shown previously, grains affect both the ionisation fraction and the dependence of the diffusivities on $\xi$. Overall, grains tend to reduce the ionisation fraction by several orders of magnitude (typically two to three). Owing to the change in composition (grains become the dominant charge carrier close to the midplane), diffusivities increase by two to four orders of magnitude in our worst-case scenario, compared with the metal-free case. In addition, the Hall diffusivity $\eta _H$ can be reversed. The effect of grains naturally depends on the assumed grain size and abundance. It is usually found that grain size $a > 1\ \mathrm {\mu }\mathrm {m}$ does not affect the conductivity tensor too much (Salmeron & Wardle 2008) and that a significant depletion of small grains also reduces their effect (Ilgner & Nelson 2006). All of these calculations assume all of the grains have the same size, which most presumably overestimates the abundance of grains and their effect on the conductivity tensor. More realistic grain size distribution, including more complex chemical reaction networks (e.g. Thi et al. 2019) tend to obtain diffusivities that are within an order of magnitude of the diffusivities discussed in our grain-free metal-free scenarios.

Finally, the physics at the grain surface is poorly understood, which gives a lot of freedom to chemical models. For instance, some authors assume a fixed sticking coefficients of electrons and ions on dust grains as we do (Sano et al. 2000; Wardle 2007), whereas others include dependences as a function of the grain size, charge and temperature (Ilgner & Nelson 2006; Bai 2011a); note, however, that the dependency of the sticking coefficient on the grain charge differs significantly between these authors. This can have an additional order of magnitude effect on the resulting diffusivities.

Overall, one is forced to conclude that the conductivity tensor of PPDs is plagued by uncertainties and that no chemical/ionisation/grain model is better than the other. Given the previous discussion, the uncertainty on the diffusion coefficient is at least ${\pm }3$ orders of magnitude, which has dramatic effects on the dynamical behaviour of these objects. Until more constraints are obtained for these coefficients, theoreticians are forced to explore in a more or less systematic manner the parameter space of the conductivity tensor.

Owing to these uncertainties, we focus in the following on the ‘intermediate’ case of a metal-free grain-free case that we discussed in § 3.4.2 and for which diffusivities are given in figure 10.

PART TWO: Disc dynamics: global and local views

4. Introduction

4.1. Motivations

Explaining accretion in discs is a long-standing problem of modern astrophysics. Even though angular momentum transport equations have been known for a long time, the road to quantifying the level of stress in various discs has been paved with unforeseen difficulties. The main idea that has been followed since the pioneering work of Shakura & Sunyaev (1973) is that accretion discs are somehow turbulent, and this turbulence generates a radial stress. The key is then to relate this radial stress to the other large-scale quantities such as the disc surface density $\varSigma$ and thickness $H$, the rotation rate $\varOmega$, the diffusivities and the magnetic field strength. This approach is, in essence, very similar to the mixing length theory of convection, except that in the disc case, one does not transport heat but angular momentum. In discs, it is called the $\alpha$-disc theory.

Here, we present the basic concepts behind accretion and the $\alpha$-disc theory. Then, we introduce the magneto-rotational instability (MRI), which is probably the most promising instability to explain the origin of accretion in astrophysical discs. Finally, we apply the MRI in the context of PPDs, taking into account non-ideal MHD effects.

4.2. Disc equilibrium

A PPD is typically made of gas (and possibly dust) orbiting a young stellar object of mass $M$. Here, we assume that the gravity of the orbiting gas onto itself (self-gravity) is negligible. This is not necessarily true in very massive discs or in the outer parts of young class 0 objects. Under these assumptions, the gravitational potential is simply that of the central object and the equilibrium may simply be written as

(4.1)\begin{equation} \left.\begin{gathered} 0=-\dfrac{1}{\rho}\dfrac{\partial P}{\partial R}-\partial_R\psi+\varOmega^2R,\\ 0= -\dfrac{1}{\rho}\dfrac{\partial P}{\partial z}-\partial_z\psi, \end{gathered}\right\} \end{equation}

where $(R,z)$ are cylindrical coordinates and $\varOmega$ is the angular velocity of the flow, which we assume only depends on $R$ and $\psi =-GM/(R^2+z^2)^{1/2}$ is the cylindrical potential. A useful quantity will be the Keplerian frequency, which corresponds to the orbital frequency of a test particle on a circular orbit at radius $R$:

(4.2)\begin{equation} \varOmega_K(R)=\sqrt{\frac{GM}{R^3}}. \end{equation}

In order to simplify the computation, let us assume that the disc is locally isothermal:⁷ $T(R)$. Under these assumptions, the sound speed may be written

(4.3)\begin{equation} c_s\equiv\sqrt{\frac{P}{\rho}}=\sqrt{\frac{kT}{\mu}}, \end{equation}

where $k$ is Boltzmann's constant and $\mu$ is the mean molecular mass. As the disc is locally isothermal, $c_s$ only depends on $R$, as the temperature does.

We start with the vertical equilibrium, which we consider close to the disc midplane ($z\ll R$) because we assume the disc is thin:

(4.4)\begin{align} c_s^2\partial_z\log \rho&=-\frac{GMz}{(R^2+z^2)^{3/2}}\nonumber\\ &\simeq z\varOmega_K ^2+O(z^3), \end{align}

where we have assumed $z\ll R$. We deduce from this the vertical density profile

(4.5)\begin{equation} \rho=\rho_0(r)\exp\left(-\frac{z^2}{2H^2}\right), \end{equation}

where we have defined the disc scale height

(4.6)\begin{equation} H\equiv c_s/\varOmega_K. \end{equation}

The thin disc approximation $H\ll R$ implies that the disc is cold, or in other words that $c_s\ll R\varOmega _K$.

In the radial direction, we first have to compare the radial pressure gradient with the gravitational potential

(4.7)\begin{equation} 0=\underbrace{-\frac{1}{\rho}\frac{\partial P}{\partial R}}_{\sim c_s^2/R}-\underbrace{\partial_R\psi}_{\sim \varOmega_K^2R}+\varOmega^2R. \end{equation}

The pressure gradient is $(H/R)^2$ smaller than the gravitational potential and can be neglected in the thin disc approximation. This means that the disc is, to a very good approximation, a Keplerian disc $\varOmega =\varOmega _K$. Note, however, that local (i.e. on radial length scales of the order of $H$) pressure variations may exist leading to measurable deviation from the Keplerian rotation. These variations are typically responsible for zonal flows and local pressure maxima.

4.3. Accretion theory

The energetics of MHD-driven discs has been discussed extensively by Balbus, Gammie & Hawley (1994), Balbus & Hawley (1998) and Balbus & Papaloizou (1999). Here, we revisit this question, and include the possibility of wind-driven accretion in the system. The accretion of mass in astrophysical discs is described by the equation of mass, angular momentum and mechanical energy conservation equations:

(4.8)\begin{gather} \frac{\partial \rho}{\partial t}+\boldsymbol{\nabla}\boldsymbol{\cdot}\rho\boldsymbol{u}=0 \end{gather}

(4.9)\begin{gather} \frac{\partial R\rho u_\phi}{\partial t}+ \boldsymbol{\nabla}\boldsymbol{\cdot}\left[R\rho u_\phi\boldsymbol{u}-R\frac{B_\phi \boldsymbol{B}}{4{\rm \pi}} \right]=0 \end{gather}

(4.10)\begin{align} &\frac{\partial \left(\dfrac{1}{2}\rho u^2+\rho \psi +\dfrac{B^2}{8{\rm \pi}}\right)}{\partial t} +\boldsymbol{\nabla}\boldsymbol{\cdot}\left[\left(\frac{1}{2}\rho u^2+\rho \psi+P +\frac{B^2}{4{\rm \pi}}\right)\boldsymbol{u}-\frac{\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{B}}{ 4{\rm \pi}}\boldsymbol{B}-\frac{\boldsymbol{\mathcal{E}}_{\mathrm{NI}}\times \boldsymbol{B}}{4{\rm \pi}}\right]\nonumber\\ &\quad = P\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}+ \frac{\boldsymbol{\mathcal{E}}_{\mathrm{NI}}\boldsymbol{\cdot} \boldsymbol{J}}{c}, \end{align}

where $\boldsymbol {\mathcal {E}}_{\mathrm {NI}}$ are electromotive forces owing to non-ideal effects. Note that molecular viscosity is usually negligible in these equations as it is several orders of magnitude smaller than non-ideal MHD effects. One notable exception is naturally when non-ideal MHD effects are absent, such as in ideal-MHD flows or in purely hydrodynamic flows subject to turbulence and/or spiral density waves. In these cases, viscosity becomes non-negligible in the energy equation because of the formation of small-scale structures, either through a direct turbulent cascade, or thanks to shocks. In any case, this viscosity then leads to an additional definite negative source term in the energy equation, which transforms mechanical energy into heat. The energy flux and angular momentum flux terms associated with viscosity are always negligible for practical applications.

In order to capture the dynamics of the disc, we separate the gravitational potential $\psi$ as a midplane potential $\varPsi$ and a deviation as one moves away from the disc midplane $\varPhi$:

(4.11)\begin{equation} \psi=\varPsi(R)+\varPhi(R,z). \end{equation}

We also separate the mean rotational motion of the disc from its deviations (not necessarily small):

(4.12a––c)\begin{equation} u_r=v_r;\quad u_\phi=\varOmega R+v_\phi;\quad u_z=v_z, \end{equation}

where we only assume that $\varOmega$ satisfies the radial equilibrium in the disc midplane

(4.13)\begin{equation} \varOmega^2R=\partial_R \varPsi. \end{equation}

Under these assumptions, it is possible to rewrite the angular momentum conservation as

(4.14)\begin{equation} \frac{\partial R\rho v_\phi}{\partial t}+\rho \boldsymbol{u} \boldsymbol{\cdot}\boldsymbol{\nabla}(\varOmega R^2)+\boldsymbol{\nabla} \boldsymbol{\cdot}\left[R\rho v_\phi\boldsymbol{u}-R \frac{B_\phi\boldsymbol{B}}{4{\rm \pi}} \right]=0, \end{equation}

where we have used the continuity equation to eliminate the terms proportional to $\varOmega R^2$. A similar procedure can be followed for the energy equation, which can be written as

(4.15)\begin{align} &\left(\frac{1}{2}\varOmega^2R^2+\varPsi\right)\left[ \frac{\partial \rho }{\partial t}+\boldsymbol{\nabla}\boldsymbol{\cdot}\rho \boldsymbol{u}\right]+\rho \boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla} \left(\frac{1}{2}\varOmega^2R^2+\varPsi\right)\nonumber\\ &\quad +\varOmega\left(\frac{\partial R\rho v_\phi}{\partial t}+\boldsymbol{\nabla}\boldsymbol{\cdot} \left[R\rho v_\phi\boldsymbol{u}-R\frac{B_\phi\boldsymbol{B}}{4{\rm \pi}} \right]\right)+ \left[R\rho v_\phi\boldsymbol{u}-R\frac{B_\phi\boldsymbol{B}}{4{\rm \pi}} \right] \boldsymbol{\cdot}\boldsymbol{\nabla} \varOmega\nonumber\\ &\quad+\frac{\partial \left(\dfrac{1}{2}\rho v^2+\rho \varPhi + \dfrac{B^2}{8{\rm \pi}}\right)}{\partial t}+\boldsymbol{\nabla}\boldsymbol{\cdot} \left[\left(\frac{1}{2}\rho v^2+\rho \varPhi+P +\frac{B^2}{4{\rm \pi}}\right)\boldsymbol{v}\right.\nonumber\\ &\quad\left. -\frac{\boldsymbol{v}\boldsymbol{\cdot}\boldsymbol{B}}{4{\rm \pi}} \boldsymbol{B}-\frac{\boldsymbol{\mathcal{E}}_{\mathrm{NI}}\times \boldsymbol{B}}{4{\rm \pi}}\right]=P \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}+\frac{\boldsymbol{\mathcal{E}}_{\mathrm{NI}}\boldsymbol{\cdot} \boldsymbol{J}}{c}. \end{align}

We recognise the mass conservation equation in the first line, and the angular momentum conservation equation in the second line. Substituting (4.8) and (4.14) into the previous equation allows us to recast energy conservation as

(4.16)\begin{align} &\rho \boldsymbol{u }\boldsymbol{\cdot}[\nabla \varPsi- \varOmega^2R\boldsymbol{\nabla}R]+\left[R\rho v_\phi\boldsymbol{u}-R \frac{B_\phi\boldsymbol{B}}{4{\rm \pi}} \right]\boldsymbol{\cdot}\boldsymbol{\nabla} \varOmega\nonumber\\ &\quad +\frac{\partial \dfrac{1}{2}\rho v^2+\rho \varPhi +\dfrac{B^2}{8{\rm \pi}}}{\partial t}+ \boldsymbol{\nabla}\boldsymbol{\cdot}\left[\left(\frac{1}{2}\rho v^2+\rho \varPhi+P +\frac{B^2}{4{\rm \pi}}\right)\boldsymbol{v}\right.\nonumber\\ &\left.\quad -\frac{\boldsymbol{v}\boldsymbol{\cdot}\boldsymbol{B}}{4{\rm \pi}}\boldsymbol{B}-\frac{\boldsymbol{\mathcal{E}}_{\mathrm{NI}}\times \boldsymbol{B}}{4{\rm \pi}}\right]=P \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}+\frac{\boldsymbol{\mathcal{E}}_{\mathrm{NI}}\boldsymbol{\cdot} \boldsymbol{J}}{c}, \end{align}

where we recognise the radial equilibrium in the first term, which can be cancelled out. Hence, we obtain an energy equation for the velocity fluctuations, which reads

(4.17)\begin{align} &\frac{\partial \left(\dfrac{1}{2}\rho v^2+\rho \varPhi +\dfrac{B^2}{8{\rm \pi}}\right)}{\partial t}+\boldsymbol{\nabla}\boldsymbol{\cdot}\left[\left(\frac{1}{2}\rho v^2+\rho \varPhi+P +\frac{B^2}{4{\rm \pi}}\right)\boldsymbol{v}-\frac{\boldsymbol{v}\boldsymbol{\cdot} \boldsymbol{B}}{4{\rm \pi}}\boldsymbol{B}-\frac{\boldsymbol{\mathcal{E}}_{\mathrm{NI}}\times \boldsymbol{B}}{4{\rm \pi}}\right] \nonumber\\ &\quad =P\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}-\left[R\rho v_\phi\boldsymbol{u}-R\frac{B_\phi\boldsymbol{B}}{4{\rm \pi}} \right]\boldsymbol{\cdot}\boldsymbol{\nabla} \varOmega+ \frac{\boldsymbol{\mathcal{E}}_{\mathrm{NI}}\boldsymbol{\cdot} \boldsymbol{J}}{c}. \end{align}

4.3.1. Averaged equations

In order to compute averaged conservation equations, we define an azimuthal average as

(4.18)\begin{equation} \langle Q\rangle =\frac{1}{2{\rm \pi}}\int \textrm{d}\phi Q \end{equation}

and a vertical integration of the azimuthal average

(4.19)\begin{equation} \bar{Q}=\int_{z=-h}^{z=+h}\textrm{d}z \langle Q\rangle, \end{equation}

so that the continuity equation (4.8) reads

(4.20)\begin{equation} \frac{\partial \varSigma}{\partial t}+\frac{1}{R}\frac{ \partial}{\partial R}R \overline{\rho u_r}+\left[\langle \rho v_z\rangle\right]_{z=-h}^{+h}=0, \end{equation}

where $\varSigma \equiv \bar {\rho }$ is the gas surface density.

The equation of angular momentum conservation (4.14) can be recast using the same averaging procedure (4.19) defined above to obtain an equation relating the mass accretion rate $\overline {\rho v_r}$ as to the radial and surface stresses

(4.21)\begin{equation} \overline{\rho v_r}\frac{\partial}{\partial R}\varOmega R^2+\frac{1}{R}\frac{\partial}{\partial R} R^2\left[\underbrace{\overline{\rho v_\phi v_r}-\frac{\overline{B_\phi B_r}}{4{\rm \pi}}}_{\textrm{Radial stress}}\right]+\underbrace{\left[R\langle \rho v_\phi v_z\rangle-R\frac{\langle B_\phi B_z\rangle}{4{\rm \pi}}\right]_{z=-h}^{+h}}_{\textrm{Surface stress}}=0, \end{equation}

where we have assumed that $v\ll \varOmega R$, which allows us to neglect the remaining time derivative.

This demonstrates the close relationship between the accretion rate and the transport of angular momentum by the stresses. Angular momentum can be transported outward in the disc by the radial stress, or evacuated from the disc by a torque applied at the disc surface, as for example when a magnetised wind is present.

This link between accretion and stress can also be seen by averaging of the mechanical energy equation (4.17):

(4.22)\begin{equation} \partial_t\overline{\mathcal{E}_m}+\frac{1}{R}\frac{\partial}{\partial R} R \overline{\mathcal{F}_{m,R}} +\left[\langle \mathcal{F}_{m,z}\rangle \right]_{z=-h}^{+h}= \overline{P\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{v}}- \underbrace{\left[\overline{\rho v_\phi v_R}-\frac{\overline{B_\phi B_R}}{4{\rm \pi}} \right]\frac{\textrm{d}\varOmega}{\textrm{d}\log R}}_{\textrm{Radial stress source term}}+ \overline{\frac{\boldsymbol{\mathcal{E}}_{\mathrm{NI}}\boldsymbol{\cdot}\boldsymbol{J}}{c}}, \end{equation}

where we have the mechanical energy of the fluctuations

(4.23)\begin{equation} \mathcal{E}_m=\frac{1}{2}\rho v^2+\rho \varPhi +\frac{B^2}{8{\rm \pi}} \end{equation}

and its associated energy flux

(4.24)\begin{equation} \left(\mathcal{E}_m+P+\frac{B^2}{8{\rm \pi}}\right) \boldsymbol{v}-\frac{\boldsymbol{v}\boldsymbol{\cdot}\boldsymbol{B}}{4{\rm \pi}} \boldsymbol{B}-\frac{\boldsymbol{\mathcal{E}}_{\mathrm{NI}}\times \boldsymbol{B}}{4{\rm \pi}}. \end{equation}

This energy equation demonstrates a very important fact: unless one assumes that the energy flux locally deposits energy (which implies that a source of energy is externally provided to the disc), then the only term that can balance diffusive (and viscous, when applicable) losses is the radial stress source term, which appears as a source term in the conservation of mechanical energy. Diffusive (and viscous) source terms being necessarily negative definite, we have

(4.25)\begin{equation} \left[\overline{\rho v_\phi v_R}-\frac{\overline{B_\phi B_R}}{4{\rm \pi}} \right]\frac{\textrm{d}\varOmega}{\textrm{d}\log R}<0. \end{equation}

As this term balances losses (which convert mechanical energy into heat), it is also equal to the local heating rate of the disc is we assume the fluctuations are statistically steady (as in a saturated turbulent state) and no energy escapes via the vertical energy flux.⁸ Note that the surface stress does not appear as a source term, as it does not lead to any local heating, despite driving accretion. That's one of the key difference between radially-driven and vertically-driven accretion.

4.4. $\alpha$ disc theory

This theory assumes no wind is present at the disc surface. In order to solve the long-term evolution of the disc, one needs to express the radial stress

(4.26)\begin{equation} \overline{W_{r\phi}}=\overline{\rho v_\phi v_r}-\frac{\overline{B_\phi B_r}}{4{\rm \pi}}, \end{equation}

as a function of vertically averaged quantities such as $\varSigma$ or $\bar {P}$. Historically, and based on a purely dimensional argument (Shakura & Sunyaev 1973), it is usually assumed that

(4.27)\begin{equation} \overline{W_{r\phi}}=\alpha \bar{P}, \end{equation}

where $\alpha$ is a dimensionless constant. Physically, it can, however, be justified as a mixing length theory: let us consider turbulent velocity fluctuations $v$ in a thin disc. The fluctuations are confined in the disc thickness $H$ with a forcing frequency $\varOmega _K$ (these two quantities are the only length and frequency accessible to an ideal system). Hence, we expect $v=\theta H\varOmega _K$ where $\theta$ is a dimensionless constant, of order unity. Therefore, $\overline {W_{r\phi }}=\theta ^2 \overline {\rho v^2}=\theta ^2 \overline {\rho H^2\varOmega _K^2}$. Using (4.6), one obtains $W_{r\phi }=\theta ^2 \overline {\rho c_s^2}=\theta ^2 \bar {P}$. Hence, thanks to the vertical equilibrium of a thin disc, the prescription of Shakura & Sunyaev (1973) shows up as a mixing length theory with a length $H$, a frequency $\varOmega _K$ and $\alpha =\theta ^2$.

Interestingly, because $H=c_s/\varOmega _K$, $\theta$ is actually a measure of the Mach number of the flow $\theta =v/c_s$. If the turbulence was strongly supersonic, then strong shocks would appear, dissipating rapidly turbulent fluctuations until they become subsonic. For this reason, and in the absence of any supersonic excitation, turbulence is expected to be essentially subsonic with $\theta \lesssim 1$ and, therefore, $\alpha < 1$.

This prescription may be seen as a viscous theory. Indeed, the $\alpha$-disc prescription leads to $\overline {W_{r\phi }}=\alpha \overline {P}=\alpha \varSigma c_s H \varOmega _K$. As $R\,\mathrm {d}\varOmega _K/\mathrm {d}R=-3/2\varOmega _K$, the stress can be recast as

(4.28)\begin{equation} \overline{W_{r\phi}}=-\frac{2}{3}\nu_t\varSigma \frac{\mathrm{d}\varOmega}{\mathrm{d}\log R}, \end{equation}

where we have defined an effective viscosity $\nu _t=\alpha c_s H$. Here, we clearly recognise the usual $R-\phi$ component of the viscous stress in the Navier–Stokes equations.

Plugging the $\alpha$ prescription into (4.21) and neglecting surface (wind) contribution leads to

(4.29)\begin{equation} \overline{\rho v_r}=-\frac{1}{R\partial_R(\varOmega_KR^2)}\frac{\partial}{\partial R} R^2\alpha c_s^2\varSigma. \end{equation}

This allows us to express the mass accretion rate $\dot {M}\equiv -2{\rm \pi} R \overline {\rho v_r}$ as

(4.30)\begin{equation} \dot{M}=\frac{4{\rm \pi}}{R\varOmega_K } \frac{\partial}{\partial R} R^2\alpha c_s^2\varSigma. \end{equation}

We can then use mass conservation (4.20) to obtain an equation for $\varSigma$

(4.31)\begin{equation} \frac{\partial \varSigma}{\partial t}=\frac{1}{R}\frac{\partial}{\partial R}\left[ \frac{1}{\partial_R(\varOmega_KR^2)}\frac{\partial}{\partial R} R^2\alpha c_s^2\varSigma \right], \end{equation}

which essentially constitutes a diffusion equation for the surface density. The diffusion timescale associated to accretion can be estimated using $c_s=\varOmega _K H$. One finds

(4.32)\begin{equation} \tau_{\mathrm{visc}}^{-1}\sim \alpha\varOmega_K\left(\frac{H}{R}\right)^2\ll\varOmega_K. \end{equation}

Accretion therefore occurs on timescales much longer than the orbital timescale in thin discs. This is usually a problem for simulations trying to capture the phenomenon of accretion. However, it allows us to separate accretion from dynamics occurring at the local orbital frequency, by stating that accretion is essentially non-existent on this timescale.

4.5. $\alpha$–$\upsilon$ disc theory

This theory is identical to the alpha disc theory for the radial stress part, but it also includes a contribution from the surface term, owing to a hypothetical wind. To do so, let us define

(4.33)\begin{equation} W_{z\phi}= \rho v_\phi v_z-\frac{B_\phi B_z}{4{\rm \pi}}, \end{equation}

and in a way similar to the $\alpha$ prescription, we assume

(4.34)\begin{equation} \left[\langle W_{z\phi} \rangle\right]_{z=-h}^{+h}=\upsilon P_{\mathrm{mid}}, \end{equation}

where $P_\mathrm {mid}$ is the midplane pressure of the disc. Using the same procedure as for the $\alpha$ disc, we can express the mass accretion rate as a function of $\alpha$ and $\upsilon$

(4.35)\begin{equation} \dot{M}=\frac{4{\rm \pi}}{R\varOmega_K } \left[\underbrace{\frac{\partial}{\partial R} R^2\alpha \overline{P}}_\mathrm{radial}+\underbrace{R^2\upsilon P_\mathrm{mid}}_\mathrm{vertical}\right]. \end{equation}

The comparison between the $\alpha$ term and the $\upsilon$ term is revealing as it compares the role played by the radial and vertical stresses. One can assume that in first approximation $\bar {P}\simeq P_{\mathrm {mid}}H$ so that the vertical contribution is $R/H (\upsilon /\alpha )$ times larger than the radial one. This implies, in particular in thin discs where $R/H\gg 1$, that magnetised winds can easily be the dominant source of accretion.

In addition, using (4.35) in the continuity equation, the vertical stress term shows up as a first-order radial derivative of $\varSigma$ (${=}$advection) whereas the radial term appears as a second-order derivative as in the usual alpha disc theory. For this reason, wind-driven discs cannot be treated as viscous discs, because the wind component appears as an advective term in the surface density evolution.

4.6. Beyond the $\alpha$ prescription

The $\alpha$ disc model is useful as a starting point to characterise the evolution of discs. However, it is not based on first principles, and it would be desirable to compute the turbulent stress $W_{r\phi }$ directly from the equations of motion for the gas.

This is, however, a rather complicated task that often implies using numerical tools, as the equations of motion cannot, in general, be solved analytically. As the disc is thin, and turbulence, in the $\alpha$ disc theory, is supposed to be confined by the scale height $H\ll R$, one can start by using this scale separation to look only at what is happening at the scale $H$, leaving the global scale ($R$) apart. This is the idea of local models, often called ‘shearing box’ models, following Hawley, Gammie & Balbus (1995).

5. Local models

5.1. The Hill's approximation

The Hill's approximation is a local view of the dynamics of an orbiting system, which was initially used by Hill (1878) to model the libration motions of the moon along its orbit. It has been used more recently as an efficient tool to model the dynamics of gas or stars in gravitating systems (e.g. Goldreich & Lynden-Bell 1965) and it was later implemented numerically in the so-called ‘shearing-box’ by Hawley et al. (1995). In this model, one considers the dynamics of the flow around an equilibrium point $R_0$, which is rotating with the disc at the angular velocity $\varOmega _O\equiv \varOmega _K(R_0)$. We define a Cartesian frame $(x,y,z)$, attached to this point so that $x$ is aligned with the radius, $y$ with the azimuth and $z$ is aligned with the vertical direction (figure 12).

Figure 12. Rotating frame on a circular orbit at $R_0$.

In this frame, the system follows the usual single-fluid equations of motion (MHD). As it is rotating, we have, in addition, a Coriolis force and a centrifugal force, so that the equations of motion read

(5.1)\begin{gather} \partial_t\rho+\boldsymbol{\nabla}\boldsymbol{\cdot}\rho \boldsymbol{v}=0, \end{gather}

(5.2)\begin{gather}\partial_t \boldsymbol{v}+\boldsymbol{v}\boldsymbol{\cdot}\boldsymbol{\nabla} \boldsymbol{v}=-\frac{1}{\rho}\boldsymbol{\nabla}P+ \frac{\boldsymbol{J}\boldsymbol{\times} \boldsymbol{B}}{\rho c}-2\varOmega_0 \boldsymbol{e}_z\boldsymbol{\times} \boldsymbol{v}+\varOmega_0 R^2\boldsymbol{e}_R-\boldsymbol{\nabla}\psi, \end{gather}

(5.3)\begin{gather}\partial_t P+\boldsymbol{v}\boldsymbol{\cdot}\boldsymbol{\nabla}P=- \gamma \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{v}, \end{gather}

(5.4)\begin{gather}\partial_t\boldsymbol{B}=\boldsymbol{\nabla}\boldsymbol{\times} \left(\boldsymbol{v}\boldsymbol{\times}\boldsymbol{B}+c\boldsymbol{E}_{\mathrm{NI}}\right), \end{gather}

where $\psi$ is the gravitational potential, $\boldsymbol {E}_\mathrm {NI}$ is the non-ideal electromotive force and $R\equiv \sqrt {(R_0+x)^2+y^2}$ is the cylindrical radius. We also assume the gas follows an ideal equation of state with first adiabatic exponent $\gamma$. As is well known, the centrifugal force derives from a potential of the form $\psi _c=-\varOmega _0^2 R^2 /2$. The effective potential (gravitational plus centrifugal) in the corotating frame therefore reads

(5.5)\begin{equation} \psi_{\mathrm{eff}}=-\frac{GM}{\left((R_0+x)^2+y^2+z^2\right)^{1/2}}- \frac{1}{2}\varOmega_0^2 \left((R_0+x)^2+y^2\right). \end{equation}

Hill's model focuses on a ‘small’ (i.e. of the order of the disc scale height in the case of a gaseous disc) region around the fiducial point $R_0$. We therefore expand the effective potential around this point, assuming $x\sim y\sim z\lesssim H$ to obtain Hill's potential

(5.6)\begin{equation} \psi_{\mathrm{eff,Hill}}=\varOmega_0^2\left[- \frac{3R_0^2}{2}-\frac{3}{2}x^2+\frac{1}{2}z^2+{O}\left(\frac{H^3}{R_0}\right)\right]. \end{equation}

This effective potential has been truncated at the first non-trivial order. It is, however, interesting to note that it does not depend on $y$ and that is does not contain any cross term such as $xy$ or $xz$. It is also independent from $R_0$ (apart from the constant term). This simplicity in the effective potential is what makes this model so useful for analytical and numerical computation. Any higher-order expansion will include curvature terms such as $x/R_0$ dependences and cross-dependences, making calculations much more tedious.

Let us emphasise already at this stage that Hill's potential is not adapted to global phenomenon. This can be seen by comparing the iso-potentials of $\psi _{\mathrm {eff}}$ and $\psi _{\mathrm {eff},\mathrm {Hill}}$ (figure 13). Hill's approximation is found to be symmetrical in $x\rightarrow -x$, implying that one does not know where the centre of attraction is located (both $x\rightarrow -\infty$ and $x\rightarrow +\infty$ are technically valid). Moreover, the neutral iso-potential $\psi (x,z)=-3\varOmega _0^2R_0^2/2$ has an asymptote for $z\rightarrow +\infty$ at $x=R_0(\sqrt {3}-1)$, which is absent in Hill's approximation. This asymptote is key for outflows to be ejected to $z\rightarrow \infty$, and is the main reason why local models always produce outflows which depend on the location of the $z$ boundary conditions (see § 9.2).

Figure 13. Effective potential in (a) the rotating frame $\psi _{\mathrm {eff}}$ and (b) its Hill's approximation $\psi _{\mathrm {eff,Hill}}$. Note the $x\rightarrow -x$ symmetry of Hill's potential, as well as the different asymptotic behaviour.

Finally, it can be shown that in the more general case of a central gravity of the form $\psi _{\mathrm {grav}}=r^{-2(q-1)}$, the equilibrium rotation profile is $\varOmega \propto R^{-q}$ and Hill's potential reads $\psi _{\mathrm {eff,Hill}}=\varOmega _0^2R_0^2(qx^2+z^2)/2+\mathrm {constant}$. In the case of a gravitational potential from a central point mass, we simply have $q=3/2$.

The equations of motion in Hill's approximation finally reads

(5.7)\begin{gather} \partial_t\rho+\boldsymbol{\nabla}\boldsymbol{\cdot}\rho \boldsymbol{v}=0, \end{gather}

(5.8)\begin{gather}\partial_t \boldsymbol{v}+\boldsymbol{v}\boldsymbol{\cdot} \boldsymbol{\nabla}\boldsymbol{v}=-\frac{1}{\rho}\boldsymbol{\nabla}P+ \frac{\boldsymbol{J}\boldsymbol{\times}\boldsymbol{B}}{\rho c}-2 \varOmega_0 \boldsymbol{e}_z\boldsymbol{\times}\boldsymbol{v}+ \varOmega_0^2(2q x\boldsymbol{e}_x-z\boldsymbol{e}_z) \end{gather}

(5.9)\begin{gather}\partial_t P+\boldsymbol{v}\boldsymbol{\cdot}\boldsymbol{\nabla}P=- \gamma \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{v} \end{gather}

(5.10)\begin{gather}\partial_t\boldsymbol{B}=\boldsymbol{\nabla}\boldsymbol{\times} \left(\boldsymbol{v}\boldsymbol{\times}\boldsymbol{B}+c \boldsymbol{E}_{\mathrm{NI}}\right), \end{gather}

where we have differentiated Hill's potential to obtain the tidal acceleration $\varOmega _0^2(q x\boldsymbol {e}_x-z\boldsymbol {e}_z)$. This set of equation admits a simple steady solution $\boldsymbol {V}_0$ for the velocity field by simply balancing the Coriolis and tidal forces in the $x$ direction:

(5.11)\begin{equation} \boldsymbol{V}_0= -q\varOmega_0 x \boldsymbol{e}_y. \end{equation}

This velocity field represents a constant radial shear. It corresponds to the local representation of the Keplerian flow, which is not in solid body rotation. As Hill's approximation only retains the first terms of the effective potential, the shear in Hill's model is linear with $x$ and does not depend on $z$. It is sometimes useful to work with the deviations from this velocity field. Let us define $\boldsymbol {w}\equiv \boldsymbol {v}-\boldsymbol {V}_0$ (where $w$ is not necessarily small compared with $v$), for which the equations of motion read

(5.12)\begin{gather} D_t\rho+\boldsymbol{\nabla}\boldsymbol{\cdot }\rho \boldsymbol{w}=0, \end{gather}

(5.13)\begin{gather}D_t \boldsymbol{w}+\boldsymbol{w}\boldsymbol{\cdot}\boldsymbol{\nabla} \boldsymbol{w}=-\frac{1}{\rho}\boldsymbol{\nabla}P+\frac{\boldsymbol{J} \boldsymbol{\times}\boldsymbol{B}}{\rho c}-2\varOmega_0 \boldsymbol{e}_z\boldsymbol{\times}\boldsymbol{w}+q\varOmega_0w_x \boldsymbol{e}_y- \varOmega_0^2 z\boldsymbol{e}_z \end{gather}

(5.14)\begin{gather}D_t P+\boldsymbol{w}\boldsymbol{\cdot}\boldsymbol{\nabla}P=- \gamma \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{w} \end{gather}

(5.15)\begin{gather}D_t\boldsymbol{B}=-q\varOmega_0B_x\boldsymbol{e}_y+ \boldsymbol{\nabla}\boldsymbol{\times}\left(\boldsymbol{w}\boldsymbol{\times} \boldsymbol{B}+c\boldsymbol{E}_\mathrm{NI}\right), \end{gather}

where we have defined the comoving derivative $D_t\equiv \partial _t-q\varOmega x\partial _y$. It is worth noting that this last set of equations does not present any $x$ dependency, except in the comoving derivative. This property is the key allowing the definition of a numerical ‘shearing box’, which we present later. In exchange, new source terms have appeared when going from $v$ to $w$: the ‘lift up’ effect $q\varOmega _0w_x\boldsymbol {e}_y$, which is essentially the advection of the mean Keplerian flow by radial motions, and the ‘$\varOmega$ effect’ of dynamo theory (actually due to shear, not rotation) $-q\varOmega _0B_x\boldsymbol {e}_y$. As we show later, these terms are key elements to the local physics of accretion discs. A local equivalent of the angular momentum equation can also be derived from (5.13). By defining $\mathcal {L}=w_y+(2-q)\varOmega _0 x$, we obtain the conservation equation

(5.16)\begin{equation} D_t \mathcal{L}+ \boldsymbol{w}\boldsymbol{\cdot} \boldsymbol{\nabla}\mathcal{L}=-\frac{1}{\rho}\partial_y P+\frac{[\boldsymbol{J} \boldsymbol{\times}\boldsymbol{B}]_y}{\rho c}, \end{equation}

where the Coriolis and tidal terms have vanished in the definition of the local angular momentum.

5.2. The shearing box model

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 approximation⁹ can be made in place of (5.12). The asymptotic of these two models are discussed in details by Umurhan & Regev (2004). 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. The linear MRI in local models

6.1. Lagrangian analysis

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 (2019) 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 1998 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 1992).

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.2. Application of the MRI to local models

In real disc models, the disc is characterised by a mean vertical and potentially a mean azimuthal magnetic field, as in the linear analysis. However, the vertical wavelength of the perturbation cannot be larger than the disc scale height.¹⁰ In other words, $k_{z,\mathrm {min}}\simeq 1/H$, which implies that any disc model threaded by a vertical field has a minimal Alfvén frequency $\omega _{{A}}\equiv k_z V_{{A}z}$. In addition, the vertical wavenumbers accessible to a disc are quantised because of the limited vertical extension (see § 6.4.6), hence $k_z=nk_{z,\mathrm {min}}$ and $\omega _{{A}}$ is also quantised. As an example, we show in figure 18 the modes accessible to a disc threaded by a vertical field with $V_{{A}z}=0.2\ \varOmega H$. In this particular example, only the four smallest $n$ are MRI-unstable, the most unstable mode being $n=2$.

Figure 18. MRI growth rate as a function of the Alfvén frequency $\omega _{{A}}$ (black). Quantised modes accessible to a disc model with $V_{{A}z}=0.2\ \varOmega H$ are shown in red with increasing $n$ from left to right. Only four modes are MRI-unstable in this example (see the text).

This quantisation of MRI modes also indicates that the MRI can be stabilised for sufficiently strong fields. Indeed, as the field strength increases, quantised unstable modes drift to the right of figure 18 because $k_{z,\mathrm {min}}$ increases. The MRI is entirely stabilised when the lowest $n$ enters the stable regime, i.e. when $V_{{A}z} > \sqrt {2q}\varOmega / k_{z,\mathrm {min}}$. With $k_{z,\mathrm {min}}\simeq 1/H$ this implies in Keplerian discs

(6.16)\begin{equation} V_{{A}z}\gtrsim \sqrt{3}\varOmega H\rightarrow\textrm{Stability}. \end{equation}

For this reason, the MRI is often seen as a ‘weak field’ instability, even though technically, this is more a result of the geometrically thin disc approximation. Note also that this criterion is not an effect of compressibility, as is sometimes thought. Indeed, because $\varOmega H=c_s$, the previous criterion implies that the disc stabilises for vertical field strength above equipartition. This is, however, just a coincidence resulting from the vertical equilibrium in a geometrically thin disc. The MRI also exists well above equipartition in a slightly modified form, provided that large enough wavelengths are allowed (Kim & Ostriker 2000).

6.3. Non-axisymmetric MRI

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 this subsection, we relax the axisymmetry hypothesis used previously, but we still assume the flow is incompressible, and neglect vertical gravity and stratification. We then consider the following set of equations:

(6.17)\begin{gather} \partial_t\boldsymbol{w}+\boldsymbol{w}\boldsymbol{\cdot} \boldsymbol{\nabla}\boldsymbol{w}-q\varOmega x \partial_y \boldsymbol{w}=-\boldsymbol{\nabla}\varPi+\frac{\boldsymbol{B}\boldsymbol{\cdot} \boldsymbol{\nabla}\boldsymbol{B}}{4{\rm \pi}\rho}+2\varOmega w_y \boldsymbol{e}_{\boldsymbol{x}}- (2-q)\varOmega w_x\boldsymbol{e}_{\boldsymbol{y}}, \end{gather}

(6.18)\begin{gather}\partial_t \boldsymbol{B}+\boldsymbol{w}\boldsymbol{\cdot}\boldsymbol{\nabla} \boldsymbol{B}-q\varOmega x \partial_y \boldsymbol{B}= \boldsymbol{B}\boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{w}- q\varOmega B_x \boldsymbol{e}_{\boldsymbol{y}}, \end{gather}

(6.19)\begin{gather}\boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{w}=0. \end{gather}

We are going to linearize this system, assuming the field can be decomposed as $\boldsymbol {B}=\boldsymbol {B}_0+\boldsymbol {b}$. The presence of the advection term $q\varOmega x\partial _y$, however, leads to some technical difficulties as it involves an explicit spatial dependency whenever modes are non-axisymmetric. We therefore follow Kelvin (1880) and Craik & Criminale (1986) using decomposition into time-dependent ‘waves’,¹¹ which, for any quantity $\tilde {X}$, assumes as spatial decomposition

(6.20)\begin{equation} \tilde{\boldsymbol{X}}=\boldsymbol{X}(t)\exp(\textrm{i}\boldsymbol{k}(\boldsymbol{t})\boldsymbol{\cdot} \boldsymbol{x}). \end{equation}

Using this decomposition, it is easy to show that

(6.21)\begin{equation} \partial_t \tilde{{\boldsymbol{X}}}-q\varOmega x\partial_y \tilde{{\boldsymbol{X}}}=\dot{{\boldsymbol{X}}}+\textrm{i}{\boldsymbol{X}} \left(\dot{{\boldsymbol{k}}}\boldsymbol{\cdot} \boldsymbol{x}-q \varOmega x k_y\right). \end{equation}

As it is the only term that shows this explicit $x$ dependency in the equations of motion, and because these equations are assumed to be valid for all $x$, we are forced to conclude that the term in parentheses cancels out:

(6.22a––c)\begin{equation} \dot{k}_x-q\varOmega k_y=0;\quad \dot{k}_y=0;\quad \dot{k}_z=0. \end{equation}

Without loss of generality, we can therefore assume a decomposition into ‘shearing waves’, defined as

(6.23)\begin{equation} \boldsymbol{k}(\boldsymbol{t})=\boldsymbol{k}_\textbf{{0}}+q \varOmega t k_y\boldsymbol{e}_{\boldsymbol{x}} \end{equation}

which is solution to the set of equations described previously. Using this shearing wave decomposition, we obtain

(6.24)\begin{gather} \dot{\boldsymbol{w}}=-\textrm{i}\boldsymbol{k}\varPi+\textrm{i} \frac{\boldsymbol{k} \boldsymbol{\cdot} \boldsymbol{B}_0}{4\pi\rho} \boldsymbol{b}+2\varOmega w_y \boldsymbol{e}_{\boldsymbol{x}}-(2-q) \varOmega w_x\boldsymbol{e}_{\boldsymbol{y}}, \end{gather}

(6.25)\begin{gather}\dot{\boldsymbol{b}}=\textrm{i}(\boldsymbol{k} \boldsymbol{\cdot} \boldsymbol{B}_0) \boldsymbol{w}-q\varOmega b_x \boldsymbol{e}_{\boldsymbol{y}}, \end{gather}

(6.26)\begin{gather}\boldsymbol{k} \boldsymbol{\cdot} \boldsymbol{w}=0, \end{gather}

where we have dropped the explicit time dependency of $k$ and the $\tilde {\,}$ symbols for simplicity. We next take the time derivative of (6.26):

(6.27)\begin{equation} q\varOmega k_y w_x+\boldsymbol{k}\boldsymbol{\cdot} \dot{\boldsymbol{w}}=0. \end{equation}

This allows us to express the generalised pressure

(6.28)\begin{equation} \varPi=-\frac{\textrm{i}}{k^2}\left(2\varOmega w_y k_x-2\varOmega(1-q)w_xk_y\right). \end{equation}

Thus, finally, the equations of motion read

(6.29)\begin{align} \dot{\boldsymbol{w}}&=\textrm{i}\frac{\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{B}_0}{4 {\rm \pi}\rho}\boldsymbol{b}+2\varOmega w_y (1-g_{xx})\boldsymbol{e}_{\boldsymbol{x}} +2(1-q)\varOmega w_x g_{xy}\boldsymbol{e}_{\boldsymbol{x}}\nonumber\\ &\quad -q\varOmega w_x g_{yy}\boldsymbol{e}_{\boldsymbol{y}}-(2-q)\varOmega w_x(1-g_{yy})\boldsymbol{e}_{\boldsymbol{y}}-2\varOmega w_y g_{xy}\boldsymbol{e}_{\boldsymbol{y}} \nonumber\\ &\quad +2(1-q)\varOmega w_x g_{yz}\boldsymbol{e}_{\boldsymbol{z}} -2\varOmega w_y g_{xz}\boldsymbol{e}_{\boldsymbol{z}}, \end{align}

where we have introduced $g_{ij}=k_ik_j/k^2$.

It is not possible to go any further without making any approximation. Indeed, although by construction $\boldsymbol {k}\boldsymbol {\cdot } \boldsymbol {B}_0$ does not have any time dependency, the pressure factors $g_{ij}$ do have one, so that a standard normal mode decomposition is prone to failure. It is possible to numerically integrate these equations as a function of time (see, e.g., Balbus & Hawley 1992). This always leads to transiently growing solutions, i.e. perturbations that only grow for a finite time. To understand why this is always the case, let us continue our analysis using a first-order Wentzel–Kramers–Brillouin (WKB) approximation.

To obtain a dispersion relation, one needs to assume that $k$ is ‘almost’ steady, i.e. $k_y\ll k_x$ so that $\mathrm {d} \log k/\mathrm {d}t\ll \varOmega$. This limit is often described as a strongly leading or strongly trailing wave. This limit implies that we can neglect all the $g_{yj}$ terms in the previous expansion. Assuming $X=\tilde {X}\exp [\sigma t+\textrm {i} \boldsymbol {k}(\boldsymbol {t})\boldsymbol {\cdot }\boldsymbol {x})]$, we then obtain

(6.30)\begin{equation} \sigma \boldsymbol{w}=\textrm{i}\frac{\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{B}_0}{4 {\rm \pi}\rho}\boldsymbol{b}+2\varOmega v_y (1-g_{xx})\boldsymbol{e}_{\boldsymbol{x}}-(2-q) \varOmega v_x\boldsymbol{e}_{\boldsymbol{y}}-2\varOmega v_y g_{xz} \boldsymbol{e}_{\boldsymbol{z}}, \end{equation}

(6.31)\begin{equation}\sigma \boldsymbol{b}=\textrm{i}(\boldsymbol{k} \boldsymbol{\cdot} \boldsymbol{B}_0) \boldsymbol{v}-q\varOmega b_x \boldsymbol{e}_{\boldsymbol{y}}. \end{equation}

These equations clearly exhibit an Alfvén mode in the $z$ direction, with $\sigma =\pm \textrm {i} (\boldsymbol {k}\boldsymbol {\cdot } \boldsymbol {V}_{{A},0})$ and ${\boldsymbol {V}}_{{A},0}\equiv \boldsymbol {B}_0/\sqrt {4{\rm \pi} \rho }$ as one would expect. We can then solve independently the horizontal problem to obtain the dispersion relation:

(6.32)\begin{equation} \sigma^4+\sigma^2(2\omega_A^2+2(2-q)\varOmega^2g_{zz})+\omega_A^2(\omega_A^2-2q\varOmega^2g_{zz}) \end{equation}

where we have defined the Alfvén frequency $\omega _A=\boldsymbol {k}\boldsymbol {\cdot } \boldsymbol {V}_{{A},0}$. This equation describes both the traditional MRI mode and the non-axisymmetric MRI, as its close resemblance with (6.14) suggests. It should, however, be noted that in the $k_z/k\rightarrow 0$ limit, the instability is lost because $g_{zz}\rightarrow 0$ and the last term of the dispersion relation, responsible for the MRI, vanishes. Physically, this happens because the pressure gradient balances the Coriolis force in the $x$ component of the equation of motion.

Now, it should be kept in mind that this dispersion relation is derived for a shearing wave whose $\boldsymbol {k}$ is slowly evolving with time. As $t\rightarrow \infty$, one then expects $|k_x|\rightarrow \infty$ and therefore $g_{zz}\rightarrow 0$. Hence, by nature, shearing waves automatically quench the growth of the MRI as they evolve. For this reason, non-axisymmetric structures, even without any dissipation, are necessarily transiently growing solution, and therefore never lead to a genuine linear instability. That being said, non-axisymmetry is required when only a large-scale toroidal field is available in the system, because one needs $\boldsymbol {k}\boldsymbol {\cdot } \boldsymbol {B}_0\ne 0$. Therefore, there is no ‘toroidal field MRI’ as is sometimes presented in the literature. There is just a transient growth in the linear phase, which can be described as a temporary MRI in the WKB approximation, but only a non-linear feedback can re-excite new shearing waves to keep increasing the energy of the fluctuations.

It should be stressed that the absence of any linear non-axisymmetric instability is observed only in the local limit (i.e. shearing box). If one considers a global disc, including curvature and radial boundaries, then a genuine linear instability can be recovered (e.g. Curry & Pudritz 1996), with properties similar to that found in the WKB approximation presented previously.

6.4. MRI in non-ideal MHD

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 (1991), 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 (1994) 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 (1996) in unstratified models and Sano & Miyama (1999) in stratified discs, which led to the dead zone model of PPDs (Gammie 1996). Later, Wardle (1999) considered the three non-ideal MHD effects in simple axial geometry whereas Desch 2004 considered also oblique modes and toroidal fields. Balbus & Terquem (2001) isolated the physics of the Hall-MRI and Kunz (2008) 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 (2004), 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 (2004).

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 axisymmetric¹² 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 2004).

We next solve the set of equations for $\sigma$ described previously and look for unstable eigenvalues. We follow Pandey & Wardle (2012) 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 1994). 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 HSI¹³ is a new branch of instability (Kunz 2008), 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,¹⁴ 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 instability¹⁵ 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. 2015)

(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 2012). 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 2004).

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 (2010). 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. (2010) 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 1996).

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.

7. The helicoidal MRI

The helicoidal MRI (HMRI) was first identified by Hollerbach & Rüdiger (2005) using a spectral analysis of a rotating Taylor–Couette flow. This instability is known to work for arbitrarily low magnetic Reynolds numbers in rotating sheared flows, implying that it could work both in liquid sodium experiments and in weakly ionised astrophysical discs such as PPDs. For the sake of completeness, let us show here the origin of this instability and discuss its application to Keplerian discs.

In contrast to the non-axisymmetric MRI (which is fully local), the HMRI is tightly linked to the presence of curvature in the physics of the system, hence it cannot be captured in Hill's approximation (even though a WKB analysis can capture it in the global geometry; see Kirillov, Stefani & Fukumoto (2014) and the following discussion). Therefore, let us return to the evolution equations in global geometry.

7.1. Full set of equations in cylindrical geometry

We consider the motion of a conductive fluid in cylindrical coordinates. We denote the velocity of the fluid $\boldsymbol {u}$ and the magnetic field $\boldsymbol {B}$. We consider the inviscid, incompressible and ideal equations of MHD, which reads, component by component,

(7.1)\begin{align} &\partial_t u_R+u_R\partial_R u_R+\frac{u_\phi}{R}\partial_\phi u_R-\frac{u_\phi^2}{R}+u_z\partial_z u_R\nonumber\\ &\quad =-\frac{1}{\rho_0}\partial_R\left(P+\frac{B^2}{8{\rm \pi}}\right)+\frac{1}{4{\rm \pi}\rho_0}B_R\partial_R B_R+\frac{1}{4{\rm \pi}\rho_0}\frac{B_\phi}{R}\partial_\phi B_R-\frac{1}{4{\rm \pi}\rho_0}\frac{B_\phi^2}{R}+\frac{1}{4{\rm \pi}\rho_0}B_z\partial_z B_R \end{align}

(7.2)\begin{align} &\partial_t u_\phi+u_R\partial_R u_\phi+\frac{u_\phi}{R}\partial_\phi u_\phi+\frac{u_Ru_\phi}{R}+u_z\partial_z u_\phi \nonumber\\ &\quad =-\frac{1}{\rho_0}\frac{1}{R}\partial_\phi\left(P+\frac{B^2}{8{\rm \pi}}\right)+\frac{1}{4{\rm \pi}\rho_0}B_R\partial_R B_\phi+\frac{1}{4{\rm \pi}\rho_0}\frac{B_\phi}{R}\partial_\phi B_\phi+\frac{1}{4{\rm \pi}\rho_0}\frac{B_RB_\phi}{R}+B_z\partial_z B_\phi \end{align}

(7.3)\begin{align} &\partial_t u_z+u_R\partial_R u_z+\frac{u_\phi}{R}\partial_\phi u_z+u_z\partial_z u_z \nonumber\\ &\quad =-\frac{1}{\rho_0}\partial_z\left(P+\frac{B^2}{8{\rm \pi}}\right)+\frac{1}{4 {\rm \pi}\rho_0}B_R\partial_R B_z+\frac{1}{4{\rm \pi}\rho_0}\frac{B_\phi}{R}\partial_\phi B_z+ \frac{1}{4{\rm \pi}\rho_0}B_z\partial_z B_z \end{align}

for the equation of motion, which we combine into the induction equation,

(7.4)\begin{equation} \partial_t B_R+u_R\partial_R B_R+\frac{u_\phi}{R}\partial_\phi B_R+u_z\partial_z B_R=B_R\partial_R u_R+\frac{B_\phi}{R}\partial_\phi u_R+B_z\partial_z u_R \end{equation}

(7.5)\begin{equation}\partial_t B_\phi +u_R\partial_R B_\phi+\frac{u_\phi}{R}\partial_\phi B_\phi+u_z\partial_z B_\phi+\frac{u_\phi B_R}{R}=B_R\partial_R u_\phi+\frac{B_\phi}{R}\partial_\phi u_\phi+B_z\partial_z u_\phi+\frac{u_R B_\phi}{R} \end{equation}

(7.6)\begin{equation}\partial_t B_z +u_R\partial_R B_z+\frac{u_\phi}{R}\partial_\phi B_z+u_z\partial_z B_z=B_R\partial_R u_z+\frac{B_\phi}{R}\partial_\phi u_z+B_z\partial_z u_z, \end{equation}

and the continuity equation,

(7.7)\begin{equation} \frac{1}{R}\partial_R R u_R+\frac{1}{R}\partial_\phi u_\phi+\partial_z u_z =0. \end{equation}

7.2. Linearisation

We linearise the equations with respect to a background flow

(7.8)\begin{equation} \boldsymbol{u}_0\equiv R\varOmega(R)\boldsymbol{e}_\phi, \end{equation}

where $\varOmega (R)$ is the angular velocity profile of the mean flow (arbitrary, for the moment). In addition, we consider the case with a mean magnetic field defined as

(7.9)\begin{equation} \boldsymbol{B}_0=B_{\phi,0}(R)\boldsymbol{e}_\phi+B_{z,0}\boldsymbol{e}_z \end{equation}

so the only spatial dependency is in the $R$ direction for the toroidal component of the field. In addition, we are going to assume that the flow is axisymmetric, so that we can cancel $\partial _\phi$ derivatives.

The velocity and magnetic fields are then expanded as

(7.10)\begin{equation} \left.\begin{gathered} \boldsymbol{u}=\boldsymbol{u}_0+\boldsymbol{v},\\ \boldsymbol{B}=\boldsymbol{B}_0+\boldsymbol{b}, \end{gathered}\right\} \end{equation}

where deviations are assumed to be infinitely small compared with the means. The linearised equations of motion eventually read

(7.11)\begin{equation} \left.\begin{gathered} \partial_t v_R=-\partial_R \varPi -\dfrac{2B_{\phi,0}b_\phi}{4{\rm \pi}\rho_0 R}+\dfrac{1}{4{\rm \pi}\rho_0}B_{z,0}\partial_z b_R+2\varOmega v_\phi\\ \partial_t v_\phi=\dfrac{1}{4{\rm \pi}\rho_0}b_R\dfrac{1}{R}\partial_RRB_{\phi,0}+\dfrac{1}{4{\rm \pi}\rho_0}B_{z,0}\partial_z b_\phi-(2\varOmega+R\partial_R\varOmega) v_R\\ \partial_t v_z=-\partial_z \varPi+\dfrac{1}{4{\rm \pi}\rho_0}B_{z,0}\partial_z b_z \end{gathered}\right\} \end{equation}

whereas the induction equation reads

(7.12)\begin{equation} \left.\begin{gathered} \partial_t b_R=B_{z,0}\partial_z v_R,\\ \partial_t b_\phi=b_RR\partial_R\varOmega+B_{z,0}\partial_z v_\phi-v_RR\partial_R\dfrac{B_{\phi,0}}{R},\\ \partial_t b_z=B_{0,z}\partial_z v_z. \end{gathered}\right\} \end{equation}

In order to make the notations more concise and consistent with the usual MRI derivation, we define the following coefficients

(7.13)\begin{equation} \left.\begin{gathered} q\equiv-\dfrac{\textrm{d}\log\varOmega}{\textrm{d}\log R},\\ p\equiv -\dfrac{\textrm{d}\log B_{\phi,0}}{\textrm{d}\log R}. \end{gathered}\right\} \end{equation}

It should be noted that the particular case $p=1$ corresponds to a case where no axial current is present in the system. This is the reference case considered by Hollerbach & Rüdiger (2005)

Using this notation, the equations eventually read

(7.14)\begin{equation} \left.\begin{gathered} \partial_t v_R=-\partial_R \varPi -\dfrac{2B_{\phi,0}b_\phi}{4{\rm \pi}\rho_0 R}+\dfrac{1}{4{\rm \pi}\rho_0}B_{z,0}\partial_z b_R+2\varOmega v_\phi\\ \partial_t v_\phi=-(p-1)\dfrac{B_{\phi,0}}{4{\rm \pi}\rho_0 R}b_R+\dfrac{1}{4{\rm \pi}\rho_0}B_{z,0}\partial_z b_\phi-\varOmega(2-q) v_R\\ \partial_t v_z=-\partial_z \varPi+\dfrac{1}{4{\rm \pi}\rho_0}B_{z,0}\partial_z b_z \end{gathered}\right\} \end{equation}

whereas the induction equation reads

(7.15)\begin{equation} \left.\begin{gathered} \partial_t b_R=B_{z,0}\partial_z v_R,\\ \partial_t b_{\phi}=-q\varOmega b_R+B_{z,0}\partial_z v_\phi+(p+1)\dfrac{B_{\phi,0}}{R}v_R,\\ \partial_t b_z=B_{0,z}\partial_z v_z. \end{gathered}\right\} \end{equation}

7.3. WKB approximation and pressure

We are going to look for a local solution, i.e. a solution with fast variation compared with $R$. We therefore focus on a tiny region around a fiducial radius $R_0$ defining $x\equiv R-R_0$, and expand the solution as

(7.16)\begin{equation} \boldsymbol{v},\boldsymbol{b}\propto\exp\left[\sigma t+\textrm{i}(k_x x+k_z z)\right], \end{equation}

thus we obtain

(7.17)\begin{equation} \left.\begin{gathered} \sigma v_R=-\textrm{i}k_R \varPi -2\omega_{A\phi}b_\phi+\textrm{i}\omega_{Az} b_R+2\varOmega v_\phi,\\ \sigma v_\phi=-(p-1)\omega_{A\phi}b_R+\textrm{i}\omega_{Az} b_\phi-\varOmega(2-q) v_R,\\ \sigma v_z=-\textrm{i}k_z \varPi+\textrm{i}\omega_{Az} b_z, \end{gathered}\right\} \end{equation}

whereas the induction equation reads

(7.18)\begin{equation} \left.\begin{gathered} \sigma b_R=\textrm{i}\omega_{Az} v_R,\\ \sigma b_\phi=-q\varOmega b_R+\textrm{i}\omega_{A\phi} v_\phi+(p+1)\omega_{A\phi}v_R,\\ \sigma b_z=\textrm{i}\omega_{Az} v_z, \end{gathered}\right\} \end{equation}

where we have defined the Alfvén frequency $\omega _{Az}=k_z B_{0,z}/\sqrt {4{\rm \pi} \rho _0}$ and the toroidal Alfvén frequency $\omega _{A\phi }=B_{\phi ,0}/(R_0\sqrt {4{\rm \pi} \rho _0})$. It is obvious from here that the vertical equation of motion and induction are just fed by the horizontal problem, but they do not have any feedback in the horizontal plane. We therefore consider only the horizontal equations without loss of generality. In order to solve for the pressure, we dot the equation of motion by $k$ to obtain an equation for $\varPi$. This gives

(7.19)\begin{equation} \varPi = \frac{2\textrm{i}k_R\omega_{A\phi}b_\phi}{k^2}-2\textrm{i}\varOmega v_\phi\frac{k_R}{k^2}. \end{equation}

This eventually leads to the following linear problem

(7.20)\begin{equation} \begin{pmatrix} -\sigma & 2\varOmega g_{zz} & \textrm{i}\omega_{Az} & -2\omega_{A\phi}g_{zz} \\ -(2-q)\varOmega & -\sigma & -(p-1)\omega_{A\phi} & \textrm{i}\omega_{Az}\\ \textrm{i}\omega_{Az} & 0 & -\sigma & 0\\ (p+1)\omega_{A\phi} & \textrm{i}\omega_{Az} & -q\varOmega & -\sigma \end{pmatrix} \begin{pmatrix} v_r\\v_\phi \\ b_r\\ b_\phi \end{pmatrix}=0, \end{equation}

where $g_{zz}=k_z^2/k^2$. This has non-trivial roots provided that the matrix determinant cancels out. This condition leads to a dispersion relation on $\sigma$:

(7.21)\begin{align} &\sigma^4+\sigma^2\left[2\omega_{Az}^2+\kappa^2 g_{zz}+2(p+1)\omega_{A\phi}^2g_{zz} \right]-8\textrm{i}\sigma\omega_{Az}\omega_{A\phi}\varOmega g_{zz} \nonumber\\ &\quad +\omega_{Az}^2\left[\omega_{Az}^2-2q\varOmega^2g_{zz}+2(p-1)\omega_{A\phi}^ 2g_{zz}\right]=0. \end{align}

We recognise the usual MRI-related form of the dispersion relation in the last term. The case $p=1$, which corresponds to the HMRI initially derived by Hollerbach & Rüdiger (2005), therefore does not affect the original MRI criterion. As we show later, the HMRI driving term is actually the linear term in $\sigma$.

7.4. HMRI with resistivity

Adding resistivity $\eta$ means that the linearised induction equations are modified as follows:

(7.22)\begin{equation} \left.\begin{gathered} \sigma b_R=\textrm{i}\omega_{Az} v_R-\eta k^2 b_R,\\ \sigma b_\phi=-q\varOmega b_R+\textrm{i}\omega_{A\phi} v_\phi+(p+1)\omega_{A\phi}v_R-\eta k^2 b_\phi,\\ \sigma b_z=\textrm{i}\omega_{Az} v_z-\eta k^2b_z. \end{gathered}\right\} \end{equation}

The resulting linear system is then transformed into

(7.23)\begin{equation} \begin{pmatrix} -\sigma & 2\varOmega g_{zz} & \textrm{i}\omega_{Az} & -2\omega_{A\phi}g_{zz} \\ -(2-q)\varOmega & -\sigma & -(p-1)\omega_{A\phi} & \textrm{i}\omega_{Az}\\ \textrm{i}\omega_{Az} & 0 & -\sigma-\eta k^2 & 0\\ (p+1)\omega_{A\phi} & \textrm{i}\omega_{Az} & -q\varOmega & -\sigma-\eta k^2 \end{pmatrix} \begin{pmatrix} v_r\\v_\phi \\ b_r\\ b_\phi \end{pmatrix}=0, \end{equation}

which results in the following dispersion relation

(7.24)\begin{align} &\sigma^2(\sigma+\eta k^2)^2+\sigma(\sigma+\eta k^2)\left[2\omega_{Az}^2+2(p+1)\omega_{A\phi}^2g_{zz}\right]+(\sigma+\eta k^2)^2\kappa ^2g_{zz} \nonumber\\ &\quad -\textrm{i}\omega_{Az}\omega_{A\phi}\varOmega g_{zz}\left[(8-2q)(\sigma+\eta k^2)+2q\sigma\right]\nonumber\\ &\quad +\omega_{Az}^2\left[\omega_{Az}^2-2q \varOmega^2g_{zz}+2(p-1)\omega_{A\phi}^2g_{zz}\right]=0. \end{align}

This can be combined into a standard fourth-order polynomial in $\sigma$:

(7.25)\begin{align} &\sigma ^4+2\sigma^3\eta k^2+\sigma^2\left[2\omega_{Az}^2+\kappa^2g_{zz}+2(p+1)\omega_{A\phi}^2g_{zz}+\eta^2k^4 \right]\nonumber\\ &\quad -\sigma\left[8\textrm{i}\omega_{Az}\omega_{A\phi} \varOmega g_{zz}-2\eta k^2\left(\omega_{Az}^2+\kappa^2g_{zz}\right)\right]\nonumber\\ &\quad +\omega_{Az}^2\left[\omega_{Az}^2-2q\varOmega^2g_{zz}+2(p-1)\omega_{A\phi}^2 g_{zz}\right]\nonumber\\ &\quad -(8-2q)\textrm{i}\omega_{Az}\omega_{A\phi}\varOmega g_{zz} \eta k^2+\eta^2k^4\kappa^2g_{zz}=0, \end{align}

which is formally equivalent to Kirillov et al. (2014, appendix B). Although this dispersion relation contains both the MRI and the HMRI in the resistive regime, it is difficult to disentangle the two instabilities. We therefore follow Liu et al. (2006) and consider the limit where the resistivity is much larger than the other terms.

7.5. Inductionless limit

In the inductionless limit, we neglect magnetic induction by considering roots $\sigma \ll \eta k^2$. We therefore introduce a small parameter $\varepsilon =1/\eta k^2$, and expand the dispersion relation at first order in $\varepsilon$. This suppresses the strongly damped magnetic modes, having $\sigma \sim -\eta k^2$, which are primarily damped Alfvén waves and on which the usual MRI lives. For this reason, the inductionless limit allows us to suppress (stabilised) MRI modes. Following this limit, we get a second-order dispersion relation:

(7.26)\begin{equation} \sigma^2+\sigma\varepsilon \left[2\omega_{Az}^2+2(p+1)\omega_{A\phi}^2g_{zz}\right]+\kappa ^2g_{zz}-2\textrm{i}\varepsilon \omega_{Az}\omega_{A\phi}\varOmega g_{zz}(4-q)=0. \end{equation}

The roots are simple to obtain and give at first order in $\varepsilon$:

(7.27)\begin{equation} \sigma=\pm \textrm{i}\kappa g_{zz}^{1/2}+\varepsilon\left[\pm \frac{1}{\kappa}\omega_{Az}\omega_{A\phi}\varOmega g_{zz}^{1/2}(4-q)-\omega_{Az}^2-(p+1)\omega_{A\phi}^2g_{zz}\right]. \end{equation}

As can be seen, this growth rate describes a small deviation from pure epicyclic oscillations. Therefore, in contrast to the MRI, which is an instability of the (slow) Alvén branch, the HMRI is an overstability of the epicyclic branch. This explains, in part, its survival in the limit of small magnetic Reynolds numbers.

The instability clearly arises for the $+$ sign of the roots, when

(7.28)\begin{equation} \frac{1}{\kappa}\omega_{Az}\omega_{A\phi}\varOmega g_{zz}^{1/2}(4-q)-\omega_{Az}^2-(p+1)\omega_{A\phi}^2g_{zz} > 0\rightarrow\textrm{Instability}. \end{equation}

As this condition is on a second-order polynomial on $\omega _{A\phi }$, it is strictly equivalent to requiring that the discriminant of the polynomial is positive, i.e. that

(7.29)\begin{equation} \frac{(4-q)^2}{2(2-q)}-4(p+1) > 0. \end{equation}

Assuming $\kappa ^2>0$, this is again equivalent to asking that

(7.30)\begin{equation} q^2+8qp-16p>0 \end{equation}

which is (again) a second-order polynomial in $q$, which we ought to resolve. The stability condition is then simply

(7.31)\begin{equation} q > 4(\sqrt{p(1+p)}-p)\quad\textrm{or}\quad q < -4(\sqrt{p(1+p)}+p)\rightarrow\textrm{Instability}. \end{equation}

In the current-free configuration ($p=1$), we recover the so-called Liu limit (Liu et al. 2006) $q > 4(\sqrt {2}-1)\simeq 1.657$, i.e. that a current-free toroidal field is stable in the Keplerian regime. It is the main reason why this instability has been mostly neglected in the astrophysical context.

However, if one allows for an axial current in the system ($p\ne 1$), it is possible to recover the instability for Keplerian rotation profiles. This can be deduced from (7.30):

(7.32)\begin{equation} p < \frac{q^2}{8(2-q)}\rightarrow\textrm{Instability}. \end{equation}

In the Keplerian regime, we therefore require $p < 9/16=0.5625$.

The maximum growth rates are obtained from the mean of the two roots of (7.28), which is simply

(7.33)\begin{equation} \omega_{A\phi,\mathrm{max}}= \frac{1}{2\kappa}\omega_{Az}\varOmega g_{zz}^{1/2}(4-q). \end{equation}

The maximum growth rate is then given by

(7.34)\begin{equation} \sigma_{\mathrm{max}}=\textrm{i}\kappa g_{zz}^{1/2}+\frac{V_{Az}^2}{\eta}\left[ g_{zz}\frac{(4-q)^2}{2(2-q)(p+1)}-1\right], \end{equation}

which turns out to be spatially scale-free. From these results, we can deduce a few important results for Keplerian discs with $q=3/2$. First, one finds that the growth rate in units of the local orbital frequency scales like the Ohmic Elsasser number, i.e. $\textrm {Im}(\sigma _\mathrm {max})/\varOmega \simeq \varLambda _O$. This implies relatively low growth rates in PPDs, unless they are hosting a dynamically strong vertical field (with $V_{Az}\sim \varOmega H$ where $H$ is the disc thickness). In addition, these ‘optimum’ growth rates can only be reached for $\omega _{A\phi }\sim \omega _{Az}$ (from (7.33)). Combining these constraints, we find that $V_{A\phi }\sim R\varOmega Hk$. Given that, by construction, $Hk > 1$ (because the vertical wavelength needs to fit in the disc), this shows that the azimuthal Alfvén velocity needs to be of the order of the Keplerian velocity (or larger). This is a very strong azimuthal field, well above equipartition. Among other things, such a strong field implies that the disc is no longer Keplerian as the radial tension and magnetic pressure forces are of the order of the central gravity.

For these reasons, the role played by the HMRI in the context of PPDs has been mostly ignored, as it lives in a parameter regime probably distant from that of real systems, which are known to be in Keplerian rotation, and hence for which $\omega _{A\phi }\ll \varOmega$.

PART THREE: Non-linear saturation of the MRI

The saturation of the MRI has been mostly studied in local shearing box simulations. As the seminal paper by Hawley et al. (1995), several studies have been dedicated to non-ideal effects and the role they play in the MRI saturation.

Let us first distinguish unstratified and vertically stratified shearing box models. In the unstratified model, the box is periodic in the vertical direction, making the system much simpler to analyse. It corresponds to the SSB model of Umurhan & Regev (2004) (see also § 5.2). When vertical stratification is included, the flow is usually allowed to escape through the vertical boundary conditions, potentially leading to outflows. Here, we first focus on unstratified models before moving to stratified models.

It is also important to separate simulations with and without a mean field. Although the MRI (the linear instability) does require a mean field to exist, it has been shown that in the non-linear regime, MHD turbulence exists without any mean field. We call this case the ‘MRI dynamo’. Although this was initially a mere curiosity, it turned out to be the fiducial configuration for many stratified and even global simulations owing to the technical difficulties associated with simulations with a mean field.

In all of these models, one uses box averages, defined for a quantity $Q$ by

(7.35)\begin{equation} \langle Q\rangle=\iiint\,\mathrm{d}^3 \boldsymbol{x} Q. \end{equation}

One of the key elements is then to quantify the $\alpha$ parameter from the turbulent stress, which is defined as

(7.36)\begin{equation} \alpha\equiv \frac{1}{\langle P\rangle}\left\langle \rho w_xw_y-\frac{B_xB_y}{4{\rm \pi}}\right\rangle, \end{equation}

which is also sometimes averaged in time.

8. Unstratified models

Except in rare circumstances, unstratified models are periodic in the vertical direction. In this case, the energy conservation equation (5.22) reads

(8.1)\begin{equation} \partial_t \langle \mathcal{E}_\mathrm{Mech}\rangle =\langle P\boldsymbol{\nabla} \boldsymbol{\cdot}\boldsymbol{w}\rangle +\langle\boldsymbol{E}_{\mathrm{NI}} \boldsymbol{\cdot}\boldsymbol{J}\rangle +q\varOmega_0 \left\langle \rho w_xw_y-\frac{B_xB_y}{4{\rm \pi}}\right\rangle. \end{equation}

In numerical experiments, because MRI turbulence is subsonic, the contribution of the $P\,\mathrm {d}V$ term is small whereas the $x$–$y$ stress term is definite positive (it is the term driving accretion). Therefore, a quasi-steady turbulent state for which $\mathcal {E}_\mathrm {mech}\sim \mathrm {constant}$ necessarily implies that non-ideal effects (or viscous effects, which have been ignored here) are non-negligible. This statement can be understood in terms of Kolmogorov's turbulent cascade argument. In the energy equation, the stress is a source of mechanical energy. In the cascade argument, it corresponds to the energy injection rate. For the cascade to be in a steady state, some dissipative effect must enter the picture at some scale to dissipate what has been injected at the beginning of the cascade. One concludes from this argument that there is no such thing as ideal MRI turbulence.

Here, we keep the terminology ‘ideal MRI’ in cases where dissipation is sufficiently small to be negligible on large scales $\varLambda _{O,H,A}\gg 1$, or when no physical dissipation is considered explicitly in the numerical method (an approach referred to as ‘implicit large Eddy simulations’ (ILES)). Nevertheless, numerical dissipation is always actively playing a role in these models, which is not without consequences, as we shall show.

Configurations with a mean field are often characterised by the plasma $\beta$ parameter of the mean field, as defined by

(8.2)\begin{equation} \beta_{\mathrm{mean}}\equiv\frac{8{\rm \pi} \langle P\rangle }{\langle B\rangle ^2}. \end{equation}

This should not be confused with the turbulent plasma $\beta$ parameter

(8.3)\begin{equation} \langle \beta\rangle \equiv\left\langle \frac{8{\rm \pi} P }{ B ^2}\right\rangle , \end{equation}

which is not a control parameter of the physical system because usually $B\gg \langle B\rangle$.

In principle, the vertical box extension $L_z$ does not necessarily match the pressure scale height $H=c_s/\varOmega$ because vertical stratification is ignored in these models. However, in most of the simulations published today, it is the case, so that the two scales can be identified. This allows us to obtain an alternative expression for $\beta _{\mathrm {mean}}$, which can be useful to interpret numerical simulations

(8.4)\begin{equation} \beta_{\mathrm{mean}}\equiv \frac {2\varOmega^2L_z^2}{V_{{A}z}^2}. \end{equation}

It is this expression for $\beta$ that is used in incompressible simulations.

8.1. Ideal MHD

8.1.1. With a mean field

The first historical models of MRI turbulence (Hawley et al. 1995) 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 1994; Latter, Lesaffre & Balbus 2009) and end up in fully developed MHD turbulence.

In the mean vertical field case, Hawley et al. (1995) found that and $\alpha \simeq 6.7 \beta _{\mathrm {mean}}^{-1/2}$. However, it was later pointed out by Bodo et al. (2008) that, in the presence of a mean vertical field, the box aspect ratio used by Hawley et al. (1995) 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. 2008). We therefore correct (Hawley et al. 1995) 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 2010). 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. 1995). Lesur & Longaretti (2007) 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. (1995) 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 (1996). 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. 2008), 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 2008b), both based on a quasi-linear theory of the toroidal MRI (Lesur & Ogilvie 2008a) and on a dynamical system approach (Herault et al. 2011; Riols et al. 2013).

According to Hawley et al. (1996), 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 2007), with $\alpha \propto 1/N$ where $N$ is the number of grid points in one direction.¹⁶ 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 2010) 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. 2007; Lesur & Longaretti 2007). In simulations with a mean field (Lesur & Longaretti 2007; Simon & Hawley 2009), 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. 2007; Walker & Boldyrev 2017), 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 2010). Instead, it has been proposed that non-local energy transfers in spectral space (Lesur & Longaretti 2011) could be responsible for this effect. If true, because non-locality is necessarily bounded (Aluie & Eyink 2010), 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. 2013). 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 2017). 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.2. Ohmic diffusion

In non-stratified shearing box simulations with a net vertical flux, the presence of at least one MRI unstable mode is given by (6.50). Using the expression for $\beta _{\mathrm {mean}}$, and assuming that $k_{z,\mathrm {min}}=2{\rm \pi} /L_z$, the stability condition of a non-stratified shearing box reads

(8.7)\begin{equation} {Rm}>\frac{\beta_{\mathrm{mean}}}{\sqrt{\dfrac{3}{2{\rm \pi}^2}\beta_{\mathrm{mean}}-4}}\rightarrow\textrm{Instability}. \end{equation}

The effect of a strong Ohmic diffusion on the non-linear saturation of the MRI was first explored by Sano, Inutsuka & Miyama (1998) in two dimensions and Fleming, Stone & Hawley (2000) in three dimensions. In the case with a mean field, it is found that when $\varLambda _O\lesssim 1$ (or, equivalently, $Rm\lesssim \beta _{\mathrm {mean}}$) and (8.7) is verified so that the system is MRI-unstable, MRI turbulence is affected: $\alpha$ becomes lower than in the ideal case, and turbulence becomes intermittent, with periods of linear growth followed by rapid decay owing to reconnection events. As in the low-$\beta$ ideal case, it is unclear whether this regime can be modelled with a constant $\alpha$ coefficient.

In the zero mean field case, Fleming et al. (2000) found that the MRI was disappearing below a critical Reynolds number $10^3<{Rm}_{c} < 10^4$. Whether ${Rm}_{c}$ depends on resolution, as the saturation level of the MRI dynamo does, is an open question.

One can combine these results into a global map showing how turbulence saturates in simulations with and without a net flux (figure 28). We have assumed that in the limit $\beta _{\mathrm {mean}} \rightarrow \infty$, the zero net flux subcritical threshold was to be considered instead of the linear stability limit. As a result, the region $\beta _{\mathrm {mean}}>10^5$ follows the zero mean field criterion, whereas $\beta _{\mathrm {mean}} < 10^5$ shows a transition region between fully developed turbulence with $\varLambda _O > 1$ and a linearly stable flow, which we have named intermittent turbulence, in reference to Fleming et al. (2000). As the data is sparse, we do not have any clear estimate of the values $\alpha$ in the intermittent turbulence region.

Figure 28. MRI turbulence regions as a function of Rm and $\beta _{\mathrm {mean}}$ in the Ohmic diffusion case. The red line corresponds to the linear stability criterion (8.7), the blue line to the limit $\varLambda _O=1$ and the green line to the zero net field limit ${Rm}_{c}=10^3$.

8.3. Ambipolar diffusion

The role played by ambipolar diffusion on the saturation level of the MRI was first studied by Hawley & Stone (1998) in the two-fluid limit. They found that MRI turbulence was unaffected by ambipolar diffusion for $\varLambda _A\gtrsim 100$. Because PPDs are strongly collisional, the two-fluid approach is not very efficient because numerical time steps are limited by the collision timescale. For this reason, this problem was revisited using the single-fluid approach by Bai & Stone (2011). They found that MRI is progressively suppressed by ambipolar diffusion as $\varLambda _A$ decreases.

Starting from Bai & Stone (2011) results, it is possible to create a simple phenomenology for MRI saturation under the effect of ambipolar diffusion. First, let us recall from (6.62) that the linear stability criterion with $k_{z,\mathrm {min}}=2{\rm \pi} /L_z$ is given by

(8.8)\begin{equation} \left.\begin{aligned} & \varLambda_A>\varLambda_{A,\mathrm{crit}}\frac{V_{{A}}}{V_{{A}z}}\rightarrow\mathrm{Instability},\\ & \mathrm{with}\quad \varLambda_{A,\mathrm{crit}}\equiv \frac{1}{\sqrt{\dfrac{3}{8{\rm \pi}^2}\beta_{\mathrm{mean}}-1}}. \end{aligned}\right\} \end{equation}

The linear stability of a box with a pure vertical field is therefore very similar to the Ohmic case. However, as the instability grows, $\langle V_{{A}}\rangle$ increases, which eventually leads to the violation of the stability criterion. This self-suppression effect of MRI-turbulence allows us to deduce the saturation level by assuming that at saturation

(8.9)\begin{equation} \varLambda_A= \varLambda_{A,\mathrm{crit}}\left(\frac{V_{{A}}}{V_{{A}z}}\right)_{\mathrm{sat}}. \end{equation}

The ratio of Alfvén speeds can be evaluated by $\langle \beta \rangle$:

(8.10)\begin{equation} \left(\frac{V_{{A}}}{V_{{A}z}}\right)_{\mathrm{sat}}\simeq \left(\frac{\beta_{\mathrm{mean}}}{\langle \beta\rangle}+1\right)^{\delta}. \end{equation}

In principle, one would naïvely expect $\delta =0.5$. However, we find that this estimate does not fit the numerical results of Bai & Stone (2011). This is because the large-scale field $V_{{A}}$ is not necessarily proportional to the instantaneous $\langle \beta \rangle$ computed from the fluctuations at all scales. Instead, we therefore choose $\delta =1$, which presents a better correlation to the available data. Noting that, similarly to the ideal case of Hawley et al. (1995), Bai & Stone (2011) found that $\langle \beta \rangle \simeq (2\alpha )^{-1}$, we obtain a simple estimate for $\alpha$ by combining the previous expression:

(8.11)\begin{equation} \alpha\simeq\frac{1}{2\beta_{\mathrm{mean}}}\left(\frac{\varLambda_A}{ \varLambda_{A,\mathrm{crit}}}-1\right). \end{equation}

Of course, this estimates diverges as $\varLambda _A\rightarrow \infty$ because the saturation mechanism owing to ambipolar diffusion becomes non-existent. In this case, the ideal-MHD estimate (8.5) should be used instead. A sample of predicted $\alpha$ from this saturation estimation is given in figure 29. These estimates are within 50 % of the calculated value with a pure vertical field of Bai & Stone (2011) and can be used as a proxy to estimate the transport in ambipolar-dominated discs. This estimate also successfully recovers the ideal regime when $\varLambda _A\gtrsim 50$, as reported by numerical simulations.

Figure 29. MRI turbulent transport deduced from (8.11) and (8.5) in the ambipolar-dominated regime with a mean vertical field. These estimates match the numerical values of Bai & Stone (2011) at ${\pm }50\,\%$. The white dashed line corresponds to the marginal stability limit $\varLambda _A=\varLambda _{A,\mathrm {crit}}$. The region below this line has $\alpha =0$ in the net vertical field case, but can reach $\alpha \sim 10^{-4}$ when a mean toroidal field component is introduced, thanks to the presence of unstable oblique modes (see the text).

In the case of a pure azimuthal field, Bai & Stone (2011) have shown that no turbulence is sustained below $\varLambda _A\lesssim 3$, and $\alpha$ progressively drops to 0 from the ideal MHD value at $\varLambda _A=100$.

Finally, in the case of a mixed mean vertical and azimuthal field, the subsistence of oblique modes (see § 6.4.5) even at $\varLambda _A\lesssim 0.1$ creates a weak transport with $\alpha \sim 3\times 10^{-4}$ for $V_{{A}y}\sim V_{{A}z}$ (Bai & Stone 2011). Above $\varLambda _A\sim 1$, the scaling (8.11) is approximately recovered as oblique modes become unimportant.

8.4. Hall effect

The effect of the Hall effect on the saturation level of MRI turbulence was first explored by Sano & Stone (2002) with a relatively weak Hall effect ($\mathcal {L}_H\gtrsim 20$). They found that with a mean vertical field, $\alpha$ increases with decreasing $\mathcal {L}_H$ in the aligned case, whereas $\alpha$ decreases in the anti-aligned case. However, PPDs are likely to have lower $\mathcal {L}_H$ than those studied by Sano & Stone (2002) (typically $\mathcal {L}_H\sim 1$ as in figure 10). For this reason, this problem was revisited by Kunz & Lesur (2013) with simulations in the $\mathcal {L}_H\sim 1$ regime. In the case with a mean vertical field aligned with the rotation axis, it is found that despite being violently unstable from the linear point of view owing to the HSI, the flow settles down into a quasi-laminar state for $\mathcal {L}_H\lesssim 5$, with negligible turbulent transport (figure 30). This ‘low transport state’ is characterised by a self-organised flow where the vertical field is concentrated in a narrow region in the $x$ direction (which can be identified as a ring in global geometry). The mechanism behind self-organisation in Hall-MHD is detailed in § 8.5.

Figure 30. Evolution of the turbulent transport as a function of the intensity of the Hall effect in the mean vertical field case. Data from Sano & Stone (2002) (SS02) and Kunz & Lesur (2013) (KL13). Here $\mathcal {L}_H < 0$ corresponds to anti-aligned field configuration. Note that the KL13 $\beta _{\mathrm {mean}}=10^4$ case is linearly stable for $\mathcal {L}_H^{-1}=0$ because of Ohmic diffusion, and exemplify the reactivation of the linear MRI under the action of Hall (see § 6.4.4). Note that the $\alpha$ values from Sano & Stone (2002) have been renormalised to match the definition of $\alpha$ in Kunz & Lesur (2013).

In the zero net-flux case, Sano & Stone (2002) found evidence that ${Rm}_c$ decreases from around $10^4$ in the ideal MHD regime possibly down to a few times $10^3$ for $\mathcal {L}_H\simeq 20$. At the same time, $\alpha$ increases by a factor of a few compared with the ideal case. For stronger Hall effects ($\mathcal {L}_H\lesssim 5$), a low transport state similar to the case with net flux is observed. However, in this case, the system switches back to a turbulent state periodically, resulting in short bursts of turbulence in the system (figure 31).

Figure 31. Turbulent transport $\alpha$ as a function of time in zero net flux simulations including the Hall effect. For a weak Hall effect ($\mathcal {L}_H\gtrsim 10$), $\alpha$ is larger than in the case without Hall effect by a factor ${\sim }3$. When the Hall effect increases, the system enters the low transport state ($\mathcal {L}_H\lesssim 5$). In contrast to the net field case, it periodically switches back to a high transport state, resulting in bursts of $\alpha$.

Overall, it is found that despite a powerful large-scale instability, the MRI in the Hall-dominated regime does not result into an efficient turbulent transport of angular momentum in unstratified boxes. This surprising result, however, does not hold in stratified boxes, where the Hall effect effectively leads to an enhanced radial stress (see § 9.2.2). Hence, the low transport state of Kunz & Lesur (2013) is really a peculiarity of the unstratified setup.

8.5. Self-organisation

Self-organisation is a process by which a disorganised (i.e. turbulent) flow creates large-scale and long-lived structures. There are several examples of self-organisation in nature, the best-known being probably the great red spot of Jupiter, resulting from small-scale turbulent motions which cascade to large scales forming a giant anticyclone. Self-organisation is a spontaneous symmetry-breaking process: the system starts from a statistically homogenous state and ends up in a heterogeneous state with well-identified structures. As such, self-organisation should be distinguished from local instabilities such as Kelvin–Helmholtz or the Rossby wave instability (RWI; Lovelace et al. 1999) that result from a special location in the flow (e.g. vortensity extremum).

Although self-organisation was clearly pointed out as a key phenomenon in Hall-MHD by Kunz & Lesur (2013), this phenomenon (or a weaker version of it) has been observed in ideal MHD simulations of MRI turbulence by several authors since early 2000. It is also a very promising mechanism to explain some of the structures observed in the sub-millimetric range (§ 1.4). Let us overview the different mechanisms that have been proposed to explain this phenomenon.

8.5.1. Ideal MHD

Hawley (2001) and Steinacker & Papaloizou (2002) 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 (2001) 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 (2014). They proposed that the non-diagonal components of the turbulent resistivity tensor (Lesur & Longaretti 2009) 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 2009; Guan & Gammie 2009; Lesur & Longaretti 2009) and found that the effective resistivity was always positive. Strikingly, Bai & Stone (2014) 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 2014, 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 (2005) reported the spontaneous formation of giant anticyclones, though this could be a boundary condition artefact (Fromang, private communication). Johansen, Youdin & Klahr (2009) 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. (2009) 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 (2013), 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 (2013).

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 (2013).

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

8.5.3. Ambipolar diffusion

Self-organisation owing to ambipolar diffusion was first mentioned by Bai & Stone (2014) in unstratified shearing boxes. It was also observed in cylindrical simulations including ambipolar diffusion by Béthune et al. (2016), 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 2014). In addition, the mechanism proposed for self-organisation only involves ideal-MHD terms (Bai & Stone 2014). Finally, the existence of zonal flows in this configuration (non-stratified, ambipolar diffusion dominated) largely depends on the box aspect ratio. Bai & Stone (2014) 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 2007).

9. Stratified shearing boxes

Stratified models have been mostly explored in the zero mean field configuration. When a mean vertical field is included, simulations lead to ‘high magnetic pressures that disrupt the vertical structure of the disk before the flow makes the transition to MHD turbulence’ (quoting Stone et al. 1996 p. 659, §  3). Although the situation is less dramatic when more adapted boundary conditions are used (see below), this statement explains the lack of simulations with a net vertical field until recently. The situation with a mean toroidal field is less interesting since the toroidal component can usually escape through the vertical boundaries, eventually leading to a situation identical to that without a mean field.

As in the non-stratified case, we define the plasma $\beta$ parameter of the mean vertical field threading the disc

(9.1)\begin{equation} \beta_{\mathrm{mean}}\equiv\frac{8{\rm \pi} \langle P\rangle_{z=0} }{\langle B_z\rangle_{z=0} ^2}, \end{equation}

where averages are taken in the disc midplane. In the following, we use a definition of $\alpha$ using the box averaged pressure, i.e.

(9.2)\begin{equation} \alpha\equiv \frac{1}{\langle P\rangle}\left\langle \rho w_xw_y-\frac{B_xB_y}{4{\rm \pi}}\right\rangle, \end{equation}

which does not depend on the vertical extension of the simulation domain (provided that all of the significant stress is contained in the box). However, some authors such as Stone and coworkers usually defines $\alpha$ from the vertically averaged stress and the midplane pressure $P_0\equiv P(z=0)$, which leads to predicted $\alpha$ values two or three times smaller than those obtained with the previous definition (Davis, Stone & Pessah 2010). Note also that Stone's definition yields $\alpha$, which decreases as the box size increases. Indeed, the stress being concentrated in the region $z\lesssim 2H$, the vertically averaged stress decreases as the box size increases above $2H$. These differences in the definition of $\alpha$ should be kept in mind when comparing results from different groups.

9.1. Zero mean field

9.1.1. Ideal MHD

The zero mean field stratified case was first explored by Brandenburg et al. (1995) and Stone et al. (1996). 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. 2010 and Simon, Beckwith & Armitage 2012). 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 2016 for an example). On the theoretical side, it is well reproduced by mean field models (Gressel 2010; Gressel & Pessah 2015), 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. 1996; Davis et al. 2010). 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 2006; Hirose, Blaes & Krolik 2009; Hirose et al. 2014), Flaig and collaborators (Flaig, Kley & Kissmann 2010; Flaig et al. 2012) and Bodo and collaborators (Bodo et al. 2013, 2015). 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. 2014), thanks to a mechanism which is yet to be fully elucidated.

Following the numerical convergence issue pointed out by Fromang & Papaloizou (2007) 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. 2010) 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. 2014) and 256 points per scale height (Ryan et al. 2017). 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. (2010) used periodic boundary conditions while Ryan et al. (2017) 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 (2003). These simulations exemplified (Gammie 1996) 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 (1996). 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 2011a,b; Dzyurkevich et al. 2013). This is confirmed by direct numerical simulations including Ohmic and ambipolar diffusion but no Hall effect (Bai & Stone 2013b; Simon et al. 2013b). These simulations suggest dramatically low values for $\alpha$, with $\alpha \simeq 3\times 10^{-6}$ at 1 AU (Bai & Stone 2013b) to $\alpha \simeq 10^{-3}$ at 100 AU (Simon et al. 2013b).

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. Mean field and outflows

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 (2009) and Suzuki, Muto & Inutsuka (2010), thanks to more robust numerical techniques compared with the initial attempts of Stone et al. (1996) 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 2009).

This problem was revisited by a large number of authors, both in the strong field limit $\beta _{\mathrm {mean}}\simeq 1$ (Moll 2012; Ogilvie 2012; Lesur, Ferreira & Ogilvie 2013, figure 36), and in the intermediate regime $\beta _{\mathrm {mean}}=10^3\text {--}10^4$ (Fromang et al. 2013; Bai & Stone 2013a). It was quickly demonstrated that these outflows were a variation of Blandford & Payne (1982) magneto-centrifugal paradigm (Fromang et al. 2013; Lesur et al. 2013) 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 (1982) paradigm. From Lesur et al. (2013).

First, the mass loss rate depends on the vertical box extension (Fromang et al. 2013), 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. 2010; McNally & Pessah 2015). 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 1982) 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. 2016). 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. 2013), 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. 2016) 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. 2013), 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. (2013) 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. 2013). In the vertical direction, however, as the box gets higher, the outflow rate decreases. Fromang et al. (2013) and Lesur et al. (2013) 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.

Assuming $L_z\simeq 10H$, the combination of the data published by Suzuki & Inutsuka (2009), Bai & Stone (2013a) and Fromang et al. (2013) leads to

(9.6)\begin{equation} \zeta\simeq 5\times 10^{-5}+\frac{10}{\beta_{\mathrm{mean}}}\quad\mathrm{for}\ \beta\gtrsim 10^2. \end{equation}

This relation gets shallower for $\beta _{\mathrm {mean}} < 10^2$ (Bai & Stone 2013a; Lesur et al. 2013), and eventually leads to a sharp decrease of $\zeta$ when $\beta \lesssim 1$ (see figure 9 in Lesur et al. 2013). It is also expected that the prefactor is a decreasing function of $L_z$ and could trace the geometrical aspect ratio $H/R$ of the disc, which is not captured in shearing boxes (Bai & Stone 2013a). However, this dependence is complex because it likely depends on the magnetisation. It should also be pointed out that this relation is close to that found in RMHD simulation of the MRI in the context of dwarf novae (Scepi et al. 2018). Hence, the disc thermodynamics do not seem to greatly affect the outflow rate. Nevertheless, this scaling should be taken with caution, and is probably only an upper bound to the real outflow rates one would obtain solving for the global problem.

The second parameter characterising the outflow is the angular momentum extracted by the wind, $\upsilon$, defined as

(9.7)\begin{equation} \upsilon \equiv \frac{\overline{T_{yz}}|_\mathrm{top} -\overline{T_{yz}}|_\mathrm{bottom} }{P_\mathrm{mid}}, \end{equation}

where we have defined the stress tensor $T_{yz}=\rho w_y w_z-B_yB_z/4{\rm \pi}$. This quantity directly enters the angular momentum conservation law (see § 4.5) so that, once it is known, one can automatically compute the accretion rate associated with the wind. Evaluating $\upsilon$ is, however, notoriously difficult in shearing boxes because its sign is not well-defined (see the symmetry discussion given previously). A naïve temporal averaging therefore leads to a negligible value of the angular momentum extracted by the wind. In order to circumvent this difficulty, several strategies have been used: Fromang et al. (2013) computed $\upsilon$ on a short time period, during which the polarity remains fixed, whereas Bai & Stone (2013a) and Scepi et al. (2018) computed the time-averaged absolute value of $\upsilon$. Another difficulty lies in the altitude at which this quantity is evaluated. Bai & Stone (2013a) evaluated it at the box boundary, whereas Scepi et al. (2018) computed its maximum as a function of $z$. In all these cases, the results behave like the scaling inspired from Scepi et al. (2018):

(9.8)\begin{equation} \upsilon \simeq \frac{(4\pm 3)\times 10}{\beta_{\mathrm{mean}}}\left[\left (\frac{\beta_{\mathrm{mean}}}{4.7\times 10^4}\right)^2+1\right]^{0.3}, \end{equation}

where the uncertainty is a result of the different methods of measurement found in the literature. As for $\zeta$, this scaling probably gets shallower for $\beta _{\mathrm {mean}}<10^2$. However, this has never been properly evaluated in the literature.

It should be noted that for $\beta _{\mathrm {mean}} > 1$, $\upsilon <\alpha$. Hence, the vertical stress is always weaker than the radial stress, by a factor $O(\beta _{\mathrm {mean}}^{1/2})$. However, the respective contribution of these two terms to the mass accretion rate in the disc is proportional to $(R/H) (\upsilon /\alpha )$ (see § 4.5), implying that the vertical stress is the dominant mass accretion mechanism whenever $\beta _{\mathrm {mean}}\lesssim (R/H)^2$ (Fromang et al. 2013; Bai & Stone 2013a).

9.2.2. Non-ideal MHD

Simulations with a mean vertical field first focused on the effect of Ohmic and ambipolar diffusion (Bai & Stone 2013b; Simon et al. 2013a), neglecting the Hall effect for technical reasons. It was found that the presence of a mean field allows the formation of a magnetised outflow at the disc surface ionised by far UV radiation, leading to accretion because of the angular momentum extracted by the outflow. Complete models including also the Hall effect (Bai 2014; Lesur et al. 2014; Bai 2015; Simon et al. 2015) showed similar outflows, with, in addition, the presence of a midplane laminar stress owing to the Hall effect (Lesur et al. 2014).

The difference between Ohmic only, $\text {Ohmic}+\text {ambipolar}$, and $\text {Ohmic}+\text {ambipolar}+\text {Hall}$ effect is demonstrated in figure 38 for an MMSN disc at 1 AU. In the Ohmic-only case, one recovers the turbulent surface layer and a dead midplane. Adding ambipolar diffusion leads to a different picture where the disc surface becomes mostly laminar, with a weakly magnetised outflow. This might be surprising at first sight because the linear analysis shows little difference between the Ohmic and Ohmic+ambipolar cases regarding the localisation and growth rates of the eigenmodes (see figure 26). The difference between Ohmic and ambipolar cases is a result of the fact that $\eta _A\propto B^2$ whereas $\eta _O$ does not depend on $B$. As eigenmodes grow, the horizontal field grows as well and $\eta _A$ increases. This rapidly leads to the saturation of the eigenmode by diffusion, in a way similar to § 8.3. In the end, the perturbation never grows to the large-enough amplitudes required to break up in developed turbulence, as in the Ohmic case.

Figure 38. Maxwell stress $M_{xy}$ averaged horizontally as a function of $t$ and $z$ in simulations with Ohmic diffusion only (a), Ohmic and ambipolar diffusion (b) and Ohmic, ambipolar and Hall effect (c). Simulations computed at 1 AU assuming a MMSN, with a mean vertical field $\beta _{\mathrm {mean}}=10^5$. From Lesur et al. (2014).

When the Hall effect is eventually added, the disc midplane is subject to a laminar stress in the region 1–10 AU, in addition to the weak surface outflow that subsists. The midplane stress is tightly linked to the field polarity: a negative mean field ($\boldsymbol {\varOmega }\boldsymbol {\cdot }\boldsymbol {B} < 0$) leads to its disappearance. This laminar stress is a result of the HSI (Lesur et al. 2014), which only shows up for positive polarities (see the discussion in §§ 6.4.4 and 6.4.6). Hence, the Hall effect is indeed able to revive the MRI in dead zones provided that the mean field has a positive polarity, as already deduced from the linear theory. However, it saturates as a laminar stress. This effect is recovered up to $R\sim 30\ \mathrm {AU}$ but disappears at 100 AU in MMSN models (Simon et al. 2015) because the Hall effect is much weaker in the outermost regions of the disc. The presence of this laminar stress suggests that the correlation lengths of the midplane structures are larger than the horizontal box size, which contradicts the spirit of the shearing box approximation and calls for global simulations to properly characterise this transport process.

The top/down symmetry in non-ideal MHD models is also problematic in simulations including non-ideal MHD effects. When $\beta _{\mathrm {mean}} < 10^5$, these simulations exhibit most of the time an odd symmetry with respect to the disc midplane, which implies that the outflow is not extracting any net angular momentum from the disc. There is today no clear explanation for this trend. It could be that the midplane current layer required by the odd symmetry configuration is expelled by the strong Ohmic or ambipolar diffusion (Bai & Stone 2013b), or that the HSI spontaneously saturates in an even configuration (Lesur et al. 2014), or a combination of these effects. Using similar techniques to that used in ideal-MHD shearing box, several groups have managed to measure the transport coefficient in non-ideal MHD simulations, computing the vertical stress on one side, or computing its absolute value. A representative summary of the results is given in figure 39.

Figure 39. Measurements of transport coefficients in the literature, assuming ionisation structures at 1 AU and at 30 AU. The ideal MHD relations (9.4), (9.6) and (9.8) are shown in black dashed lines. We have used data from Lesur et al. (2014)$=$L14, Bai (2014)$=$B14, Bai (2015)$=$B15 and Simon et al. (2015)$=$S15. Simulation with the vertical field aligned with the rotation axis are in red, simulations with vertical field anti-aligned are in blue and simulations without Hall effect are in green. Note that Bai (2014) has two chemical models at 1 AU, with and without grains, hence the presence of two sets of points. Note that points on the same $\beta _{\mathrm {mean}}$ have been slightly shifted horizontally to improve readability.

Several trends can be seen in this figure. First and foremost, all of the transport coefficients are reduced by the non-ideal MHD effects (with the notable exception of one run from Lesur et al. (2014) with a very strong Hall effect at 1 AU). The most reduced coefficient is the radial angular momentum transport $\alpha$, whereas the vertical (wind) transport $\upsilon$ appears to be the less affected. This difference in behaviour is the main reason why winds are today favoured in PPDs: the efficiency of wind-driven transport is less affected by non-ideal MHD effects than radial angular momentum transport. In addition to these remarks, one should note that the sensitivity to the field polarity is more pronounced at 1 AU than at 30 AU, as expected from the profile of dimensionless Hall number (figure 10). As already guessed in the previous discussion, aligned fields tend to have a larger $\alpha$, $\zeta$ and $\upsilon$. It is interesting to also note that, for simulations at 30 AU, the scaling with $\beta _{\mathrm {mean}}$ is similar to that found in ideal MHD, apart from a constant offset.

9.2.3. Outflow-induced self-organisation

There is much evidence of self-organisation in stratified shearing boxes. The first piece of evidence can be seen in the space–time diagrams of Simon & Armitage (2014). Bai (2015) then explicitly reported the spontaneous formation of self-organised flows, which is particularly strong for $\beta _{\mathrm {mean}}\lesssim 10^4$. The same process was later identified by Simon et al. (2018) and Riols & Lesur (2019). Self-organisation is clearly dissociated from the Hall effect in stratified flows, as it appears also in simulations where the Hall effect is absent (Simon & Armitage 2014; Simon et al. 2018; Riols & Lesur 2019).

This self-organisation shows again a very tight intricacy between the field and the gas: magnetic field lines tends to concentrate into ‘gaps’, leaving regions with a lot of mass but no poloidal field (figure 40). Several scenarios have been proposed to explain this effect. Bai & Stone (2014) proposed that this is a result of the anisotropy of the turbulent diffusivity tensor, following the unstratified box argument (§ 8.5.1).

Figure 40. Shearing box simulation exhibiting self-organisation with ambipolar diffusion only. (a) Magnetic configuration, with poloidal field lines and colours showing the toroidal field component amplitude. (b) Density map (colour) and poloidal velocity streamlines. Note the poloidal field concentration in the region $(x,z)\simeq 0$, associated with a minimum of density (${=}$‘gap’). Figure from Riols & Lesur (2019).

An alternative viewpoint was proposed by Riols & Lesur (2019) who noted that, in the presence of vertical stratification, the average radial flow was converging towards the region of flux accumulation, and hence in the gap. This counter-intuitive finding implies that gaps are necessarily emptied by the outflow, and the radial flow is simply trying to replenish the gaps. This radial flow drags the poloidal field towards the gap (cf. figure 41), and is therefore a pure advection process, not a turbulent anti-diffusion effect as proposed by Bai & Stone (2014). The stronger field in the gap leads to a more efficient ejection, as observed in the ejection ‘plume’ (cf.figure 40), leading to a quasi-steady state from the disc viewpoint . This feedback loop turns out to be a linear instability of the wind-emitting disc, with predictable growth rate and optimum disturbance. Let us note that the linear instability criterion is simply

(9.9)\begin{equation} -\frac{\mathrm{d}\log \zeta}{\mathrm{d}\log \beta_{\mathrm{mean}}} > -\frac{\mathrm{d}\log \alpha}{\mathrm{d}\log \beta_{\mathrm{mean}}}\rightarrow \textrm{Instability}, \end{equation}

which essentially compares the time needed to empty the disc by the outflow to the time needed to refill the disc by ‘viscous’ diffusion (Riols & Lesur 2019). As shown previously, this criterion is generally satisfied in stratified shearing box models, so this ‘wind instability’ mechanism could potentially explain the self-organisation observed in stratified models.

Figure 41. Self-organisation feedback loop proposed by Riols & Lesur (2019). Consider a small density deficit (a), this deficit induces a radially converging flow (b), which drags the poloidal field line inwards (c). The stronger field, leads to a more efficient local ejection index that empties the region even more (d).

10. Conclusion

In ideal MHD, the MRI is found to be a robust angular momentum transport mechanism, but in the absence of an external poloidal field, it fails to account for observed accretion rates if one includes the relevant non-ideal MHD effects.

When a mean field is added, magnetised outflows are found as a result of the saturation of the MRI. In ideal MHD, outflows are expected to be the dominant angular momentum transport mechanism when $\beta _{\mathrm {mean}}\lesssim (R/H)^2$. When non-ideal effects are introduced, however, the radial transport of angular momentum is reduced significantly, whereas outflows are much less affected. As a result, when non-ideal MHD effects are strong, magnetised outflows are expected to become the dominant mechanism of angular momentum transport, even for relatively weak fields.

We note, however, that the shearing box approach cannot model these outflows properly because of several inconsistencies (top/down symmetry of the outflow, mass loss rate depends on boundary conditions). In addition, the presence of a laminar midplane stress owing to the HSI and the existence of a self-organisation mechanism associated with the outflow make the shearing box ill-suited to study the dynamics of these objects. The only way to connect our local understanding of the physics to the global dynamics of these objects is, therefore, to perform global simulations, in which outflows are modelled properly.

PART FOUR: Global models of PPDs

In this section, we revisit the physical concepts underlying ejection in accretion discs. We then use these concepts to interpret the most recent numerical models of PPDs that exhibit accretion and outflows.

11. Some outflow definitions and properties

Before exploring the connection between a weakly ionised disc and an outflow, let us underline a few important definitions and properties on outflows. To this end, let us work first in the ideal MHD approximation and assume for the moment that the outflow is stationary and axisymmetric. In a cylindrical frame $(R,\phi ,z)$, the conservation of mass and momentum read

(11.1)\begin{gather} \boldsymbol{\nabla} \boldsymbol{\cdot}\rho\boldsymbol{u}_p=0 \end{gather}

(11.2)\begin{gather}\rho \boldsymbol{u}_p\boldsymbol{\cdot}\boldsymbol{\nabla} u_R= \rho\varOmega^2 r-\partial_R P+\frac{J_\phi B_z}{c}-\frac{J_z B_\phi}{c}-\rho \partial_R\psi \end{gather}

(11.3)\begin{gather}\frac{1}{R}\rho \boldsymbol{u}_p\boldsymbol{\cdot }\boldsymbol{\nabla} \varOmega R^2=\frac{1}{R}\boldsymbol{\nabla}\boldsymbol{\cdot } \left(R\frac{\boldsymbol{B}_pB_\phi}{4{\rm \pi}}\right) \end{gather}

(11.4)\begin{gather}\rho \boldsymbol{u}_p\boldsymbol{\cdot}\boldsymbol{\nabla} u_z= -\partial_z P-\partial_z\left(\frac{B_\phi^2+B_R^2}{8{\rm \pi}}\right)+ \frac{B_R\partial_R B_z}{4{\rm \pi}}-\rho\partial_z\psi \end{gather}

where the index $p$ denotes the poloidal components $(R,z)$ of the vector fields, $\varOmega \equiv u_\phi /R$ ($\varOmega$ is not necessarily Keplerian) and $\psi =GM/(R^2+z^2)^{1/2}$ is the gravitational potential, assumed to be solely a result of the central star. Equation (11.3) can be recast in conservative form to obtain an angular momentum conservation equation

(11.5)\begin{equation} \boldsymbol{\nabla}\boldsymbol{\cdot }\left(\rho \boldsymbol{u}_p\varOmega R^2-R\frac{\boldsymbol{B}_p B_\phi}{4{\rm \pi}}\right)=0, \end{equation}

which shows that the magnetic field can carry angular momentum in the Maxwell stress $\boldsymbol {B}_pB_\phi$.

In addition to the equation of motion, one needs to solve the induction equation in the ideal regime in the poloidal and azimuthal directions

(11.6)\begin{equation} \boldsymbol{\nabla}\boldsymbol{\times}\left(\boldsymbol{u}_p\boldsymbol{\times} \boldsymbol{B}_p\right)=0, \end{equation}

(11.7)\begin{equation}\boldsymbol{\nabla}\boldsymbol{\cdot}\frac{1}{R}\left(\varOmega R \boldsymbol{B}_p-B_\phi\boldsymbol{u}_p\right)=0. \end{equation}

Finally, the magnetic field has to satisfy the solenoidal condition $\boldsymbol {\nabla }\boldsymbol {\cdot }\boldsymbol {B}=0$. In an axisymmetric configuration this can be used to express the poloidal field components as a function of a scalar streamfunction $a$:

(11.8)\begin{equation} \boldsymbol{B}_p=\frac{1}{R}\boldsymbol{\nabla}a\boldsymbol{\times}\boldsymbol{e}_\phi. \end{equation}

By construction, the value of $a$ is constant along a poloidal field line. Hence, we can label each poloidal field line with the value of its streamfunction. From these equations, it is possible to derive several physical constants of motion that are useful to interpret outflow solutions.

11.1. Frozen-in condition

The frozen-in condition is derived from mass conservation (11.1) and the poloidal induction equation (11.6). Let us first start with the vertical component of the induction equation, which reads

(11.9)\begin{equation} \frac{1}{R}\frac{\partial}{\partial R} R\left(u_zB_R-u_RB_z\right)=0, \end{equation}

from which one deduces that

(11.10)\begin{equation} u_zB_R-u_RB_z=\frac{\beta(z)}{R}, \end{equation}

where $\beta$ is an unknown scalar function (physically, $\beta /R$ is the $\phi$ component of the electromotive force $E_\phi$). Using the radial induction equation, we may see that $\partial _z\beta =0$ so that $\beta$ is a constant. To avoid any singularity for $E_\phi$ at $R=0$, we are then forced to have $\beta =0$ and hence $\boldsymbol {B}_p$ and $\boldsymbol {u}_p$ are parallel, i.e.

(11.11)\begin{equation} \boldsymbol{u}_p= \mu(R,z)\boldsymbol{B}_p. \end{equation}

We can inject this relation in the mass conservation equation to obtain

(11.12)\begin{equation} \boldsymbol{B}_p\boldsymbol{\cdot}\boldsymbol{\nabla} (\rho\mu) =0, \end{equation}

which indicates that $\eta \equiv \rho \mu$ is a constant along magnetic field lines, i.e. $\eta =\eta (a)$. We can then recast (11.11) into

(11.13)\begin{equation} \boldsymbol{u}_p=\eta(a)\frac{\boldsymbol{B}_p}{4{\rm \pi}\rho}, \end{equation}

which constitutes the frozen-in condition where $\eta$ describes the amount of mass loaded along a poloidal field line. Here $\eta$ is a direct measure of the gas density $\rho _{A}$ at the Alfvén surface, which is the region defined by $\boldsymbol {u}_p=\boldsymbol {V}_{{A},p}=\boldsymbol {B}_p/\sqrt {4{\rm \pi} \rho }$:

(11.14)\begin{equation} \eta(a)=\sqrt{4{\rm \pi} \rho_{A}(a)}. \end{equation}

Note regarding the applicability of the frozen-in condition in shearing boxes: when deriving the frozen-in condition, we have used a regularity condition at $R=0$ for the electromotive force. Such a condition does not exist in the shearing box approximation so that $E_\phi$ can be non-zero in principle. This implies that in a stationary shearing box solution, the poloidal field and velocity are not necessarily parallel (Lesur et al. 2013). Physically, it means that field lines can be indefinitely dragged radially inward or outward, without affecting the stationarity condition, thanks to the assumed radial periodicity. This is physically impossible from a global point of view, but it illustrates once again the numerous drawbacks of the shearing box model when dealing with outflows.

11.2. Magnetic surface rotation

The rotation of magnetic surfaces is obtained by substituting the frozen-in condition (11.13) into the azimuthal component of the induction equation (11.7):

(11.15)\begin{equation} \boldsymbol{B}_p\boldsymbol{\cdot}\boldsymbol{\nabla}\left( \varOmega-\frac{B_\phi\eta(a)}{4{\rm \pi} \rho R}\right)=0, \end{equation}

which allows us to define a new invariant along field lines

(11.16)\begin{equation} \varOmega^*(a)\equiv \varOmega-\frac{\eta(a)}{4{\rm \pi} \rho R}B_\phi=\varOmega_{A}-\frac{V_{{A}\phi}(R_{A})}{R_{A}}, \end{equation}

where the $A$ indices denote quantities evaluated at the Alfvén surface and $R_{A}$ is the Alfvén radius, where the outflow becomes super-Alfvénic. It can easily be checked that in a frame rotating at $\varOmega ^*(a)$, the total field is parallel to the total velocity. In other words, there is no induced electromotive force in the frame rotating at $\varOmega ^*(a)$. Therefore, $\varOmega ^*$ can be interpreted as the rotation speed of magnetic surfaces. By combining (11.13) and (11.16), we can express the total velocity as a function of the magnetic field,

(11.17)\begin{equation} \boldsymbol{u}=\frac{\eta(a)}{4{\rm \pi} \rho}\boldsymbol{B}+R\varOmega^*(a)\boldsymbol{e}_\phi . \end{equation}

This illustrates the fact that only the poloidal components of the field and the velocity are parallel in general. In the particular case where $\eta =0$ (no motion along the poloidal field lines), we recover Ferraro's iso-rotation law $\varOmega =\varOmega ^*(a)$.

11.3. Angular momentum

The angular momentum invariant is easily deduced from the angular momentum conservation equation (11.5):

(11.18)\begin{equation} \boldsymbol{B}_p\boldsymbol{\cdot}\boldsymbol{\nabla} \left(\varOmega R^2-\frac{RB_\phi}{\eta(a)}\right)=0 \end{equation}

from which we deduce the angular momentum invariant $\mathcal {L}(a)$:

(11.19)\begin{align} \mathcal{L}(a)&\equiv \varOmega R^2-\frac{RB_\phi}{\eta(a)}, \end{align}

(11.20)\begin{align} &=\varOmega^*(a)R_A^2, \end{align}

where the last equality has been obtained at the Alfvén surface using (11.14). Hence, the amount of angular momentum transported by the outflow is an invariant made of two parts: the classical kinetic contribution and a magnetic part stored in $B_\phi$.

11.4. Bernoulli invariant

The Bernoulli invariant is obtained from the scalar product of the poloidal equations of motion (11.2)–(11.4) with $\boldsymbol {u}$. One obtains

(11.21)\begin{equation} \boldsymbol{u}_p\boldsymbol{\cdot}\boldsymbol{\nabla} \left[{\frac{u^2}{2}+\psi_G}\right]=-\boldsymbol{u}_p\boldsymbol{\cdot} \frac{\boldsymbol{\nabla}P}{\rho}+\boldsymbol{u}\boldsymbol{\cdot} \frac{\boldsymbol{J}\boldsymbol{\times}\boldsymbol{B}}{\rho c}, \end{equation}

where $\psi _G=GM/(R^2+z^2)^{1/2}$ is the gravitational potential. We then use (11.17) to obtain the work of the Lorentz force:

(11.22)\begin{align} \boldsymbol{u}\boldsymbol{\cdot}\frac{\boldsymbol{J}\boldsymbol{\times}\boldsymbol{B}}{c} &=R\varOmega^*(a)\frac{\boldsymbol{J}_p\boldsymbol{\times}\boldsymbol{B}_p}{\rho c}\nonumber\\ &=\varOmega^*(a)\frac{\boldsymbol{B}_p\boldsymbol{\cdot}\boldsymbol{\nabla} \left(RB_\phi\right)}{4{\rm \pi}\rho }\nonumber\\ &=\boldsymbol{u}_p\boldsymbol{\cdot}\boldsymbol{\nabla} \left(\frac{R\varOmega^*(a)B_\phi}{\eta(a)}\right). \end{align}

In addition, we may express the work of pressure forces using the enthalpy per unit mass $\mathcal {H}$

(11.23)\begin{equation} \boldsymbol{u}_p\boldsymbol{\cdot}\frac{\boldsymbol{\nabla}P}{\rho}= \boldsymbol{u}_p\boldsymbol{\cdot}\boldsymbol{\nabla}\mathcal{H}+\delta \mathcal{Q} \end{equation}

in which we have also considered the effect of an additional heating term owing to radiative heating/cooling $\delta \mathcal {Q}$. Putting all the terms together and integrating along one particular streamline $a$, we find that the quantity

(11.24)\begin{equation} \mathcal{B}\equiv \frac{u^2}{2}+\psi_G+\mathcal{H}-\frac{R\varOmega^*(a) B_\phi}{\eta(a)}+\int_{\mathcal{S}(a)} \delta\mathcal{Q}\,\mathrm{d}s \end{equation}

is conserved along poloidal field lines and streamlines. Note that because the heating term is not a proper differential, it depends on the integral of the heating term along the chosen streamline $\mathcal {S}(a)$. Of course, for an outflow to escape up to $z\rightarrow \infty$, one needs $\mathcal {B} > 0$ on the streamlines (assuming $\mathcal {H} > 0$). In the midplane of a Keplerian disc, one has $u^2/2+\psi _G=-v_K^2/2$ where $v_K$ is the Keplerian velocity, so that additional ingredients are required to produce an outflow. Two extreme situations can be identified.

Thermal winds: Here magnetic effects are neglected and ejection is possible because the disc is hot (large initial enthalpy) or because a lot of heating is applied along the streamlines (named photo-evaporation in the PPDs community).

Cold MHD winds: Here thermal effects are neglected and the toroidal field acts as an energy reservoir to launch the outflow.

11.5. Dimensionless numbers characterising an outflow

Based on the MHD invariants, it is possible to define a series of dimensionless numbers that characterise an outflow streamline unambiguously. These numbers use the physical properties at the base of the outflow. Let us therefore write $\varOmega _0$ as the rotation rate at the base of the outflow and $R_0$ as its cylindrical radius (one typically has $\varOmega _0=\varOmega _K(R_0)$), and $B_0$ as the poloidal field strength threading the disc at the location where the outflow is launched. We then define

(11.25)\begin{equation} \left.\begin{gathered} \kappa\equiv \eta \dfrac{\varOmega_0 R_0}{B_0},\\ \omega^*\equiv \dfrac{\varOmega^*}{\varOmega_0},\\ \lambda\equiv \dfrac{\mathcal{L}}{\varOmega_0 R_0^2}=\omega \left(\frac{R_A}{R_0}\right)^2,\\ e\equiv\dfrac{\mathcal{B}}{\varOmega_0^2R_0^2/2}. \end{gathered}\right\} \end{equation}

The Bernoulli invariant can be easily expressed in terms of the other invariants, using the dimensionless rotation rate $\omega \equiv \varOmega /\varOmega _0$:

(11.26)\begin{equation} e=\frac{u_p^2}{\varOmega_0^2R_0^2}+\omega^2\frac{R^2}{R_0^2}-\frac{2R_0}{\sqrt{R^2+z^2}}+2\omega^*\left(\lambda-\omega\frac{R^2}{R_0^2}\right)+\theta, \end{equation}

where $\theta =2\left [\mathcal {H}+\int _{\mathrm {a}=\mathrm {cste}} \delta \mathcal {Q}\,\mathrm {d}s\right ]/\varOmega _0^2R_0^2$ is the dimensionless thermal energy content of the flow. This expression clearly demonstrates the contribution of $\lambda$ to the energy content of the flow. This parameter, often called the magnetic level arm, is of key importance, as can be seen if one computes its value at the outflow base, assuming $u_p\ll \varOmega _0R_0$, $\omega =1$ and $R=R_0$:

(11.27)\begin{equation} e\simeq 2\omega^*(\lambda-1)+\theta -1. \end{equation}

Obviously, for the outflow to be able to propagate up to infinity, one needs $e > 0$. We recover here the two extreme limits discussed previously: purely thermal winds, which have $\lambda =1$ and require $\theta > 1$; or cold MHD winds with $\theta =0$, which need

(11.28)\begin{equation} 2\omega^*(\lambda-1) > 1. \end{equation}

For all practical applications of MHD outflows, one has $\omega ^*\simeq 1$ to a very good approximation. This implies that an outflow can exist only if

(11.29)\begin{equation} \lambda>\tfrac{3}{2}. \end{equation}

11.6. Ejection efficiency

The existence of an outflow is tightly linked to the process of accretion happening inside the disc because the energy of the wind is obtained from the accretion power of the disc. To understand this connection, let us first define the accretion rate of the disc

(11.30)\begin{equation} \dot{M}_\mathrm{acc}(R)\equiv-\int_{z^-}^{z^+}\mathrm{d}z \int_0^{2{\rm \pi}} R\, \mathrm{d}\phi \rho u_r, \end{equation}

where the integration is performed on the box vertical extension, defined by $z^-$ and $z^+$. It is also useful to define the outflow rate of the wind between the inner radius of the disc $R_\mathrm {in}$ and the radius under consideration

(11.31)\begin{equation} \dot{M}_\mathrm{wind}(R)\equiv \int_{R_\mathrm{in}}^R \,\mathrm{d}R \int_0^{2{\rm \pi}} R\,\mathrm{d}\phi [\rho u_z]_{z^-}^{z^+}. \end{equation}

These two quantities are connected by the continuity equation (4.20):

(11.32)\begin{equation} 2{\rm \pi}\frac{\partial \varSigma}{\partial t}+\frac{1}{R}\frac{\partial}{\partial R}\left(\dot{M}_\mathrm{wind}- \dot{M}_\mathrm{acc}\right)=0. \end{equation}

At this stage, it is useful to introduce the ejection efficiency index

(11.33)\begin{align} \xi&\equiv \frac{\mathrm{d}\log \dot{M}_\mathrm{acc}}{\mathrm{d}\log R} \end{align}

(11.34)\begin{align} &=\frac{1}{\dot{M}_\mathrm{acc}}\frac{\mathrm{d} \dot{M}_\mathrm{wind}}{\mathrm{d}\log R}, \end{align}

which quantifies what fraction of the mass is being lost in the wind, the second line being obtained from the continuity equation, assuming stationarity. As expected, $\xi =0$ corresponds to a situation without any wind.

One can also relate the accretion rate to the ejection rate using the angular momentum conservation equation (4.21) as

(11.35)\begin{equation} \dot{M}_\mathrm{acc} =\frac{4{\rm \pi}}{\varOmega_K}\left[\underbrace{\frac{1}{R}\frac{\partial}{\partial R} R^2\overline{W_{R\phi}}}_{\tau_{R}} + \underbrace{R\langle W_{z\phi}\rangle_{z^-}^{z^+}}_{\tau_{z}}\right], \end{equation}

where we have defined the stress tensor $W_{i\phi }\equiv \rho v_iv_\phi -B_iB_\phi /4{\rm \pi}$. This allows us to define the radial and vertical contribution to the accretion of the disc $\tau _{R\phi }$ and $\tau _{z\phi }$. It is then useful to introduce the ratio of these two quantities

(11.36)\begin{equation} \varLambda\equiv \frac{\tau_{z}}{\tau_{r}} \end{equation}

so that the accretion rate is simply

(11.37)\begin{equation} \dot{M}_\mathrm{acc} =\frac{4{\rm \pi}}{\varOmega_K}\tau_{z}\frac{1+\varLambda}{\varLambda}. \end{equation}

It is then possible to relate the vertical torque $\tau _{z\phi }$ to the mass ejection rate by noting that the torque is evaluated high above the disc midplane, so that the kinetic contribution to the stress is negligible, which implies $\tau _{z\phi }\simeq - [RB_\phi B_z/4{\rm \pi} ]_{z^-}^{z^+}\simeq -RB_0 B_\phi (z^+)/2{\rm \pi}$ where the second equality assumes the outflow is top/down symmetric and that the poloidal field strength does not vary much between the midplane and the disc surface. It is then simple to show that the vertical torque is directly connected to the MHD invariants

(11.38)\begin{equation} \tau_{z}\simeq 2 (\varOmega^*R_A^2-\varOmega_0 R_0^2)[\rho u_z](z_+). \end{equation}

The vertical mass flux being directly related to the radial derivative of the outflow rate, we can express the accretion rate as

(11.39)\begin{equation} \dot{M}_\mathrm{acc} \simeq\frac{2}{R_0^2\varOmega_0}\frac{1+\varLambda}{\varLambda}(\varOmega^*R_A^2-\varOmega_0 R_0^2)\frac{\mathrm{d}\dot{M}_\mathrm{wind}}{\mathrm{d}\log R}. \end{equation}

From which we obtain an expression for the mass ejection index

(11.40)\begin{equation} \xi\simeq \frac{\varLambda}{2(\varLambda+1)}\frac{1}{\lambda-1}. \end{equation}

This relation reveals several key features of MHD outflows. First, it shows that large level arms $\lambda$ are associated with small ejection indices. Interestingly, this result does not depend on the radial contribution to the angular momentum budget $\varLambda$. Second, the energy constraint (11.29) imposes an upper bound on $\xi$ in cold MHD winds: $\xi \lesssim \varLambda /(\varLambda +1)\lesssim 1$. Note that this relation allows for outflows approximately as massive as the mass accretion rate, but it does not allow for outflows vastly more massive than this. Outflows with $\xi \gg 1$ therefore necessarily require some thermal energy driving to escape the potential well, as one would expect.

11.7. Connection to shearing box simulations

In the previous part, we had to introduce several local quantities to characterise outflows in local shearing box models. As pointed out, however, shearing boxes lack global constraints, which implies that some of the solutions are likely unphysical. As a first step, it is therefore useful to relate these local quantities to global MHD invariants to test the domain of validity of shearing box solutions. In this subsection, we make the assumption that the global outflow is top/down symmetric, so that the invariants are the same on both sides of the disc.

We first start with the outflow rate $\zeta$ (defined in (9.5)), which can be easily connected to the mass loading parameter $\kappa$

(11.41)\begin{equation} \kappa=\frac{1}{4}\frac{R}{H}\beta_{\mathrm{mean}} \zeta . \end{equation}

The magnetic level arm can also be related to $\upsilon$ (defined in (9.7)), provided that we neglect the kinetic contribution to the vertical stress, which is valid if we choose the disc upper boundary to lie high enough above the midplane

(11.42)\begin{equation} \lambda=1+\frac{H}{R}\frac{\upsilon}{\zeta}. \end{equation}

The energetic constraint (11.29) therefore leads to a new constraint in shearing boxes: $2\zeta /\upsilon < H/R$. This constraint cannot be satisfied by the shearing box scalings (9.6) and (9.8) when $H/R\lesssim 0.2$. Hence, shearing box wind solutions always eject too much mass with too little energy to escape the global potential well of discs with realistic aspect ratios. In other words, if one puts a shearing box wind solution in a global disc configuration, the ejected material should fall back onto the disc.

Finally, we can relate the stress rate $\varLambda$ to $\alpha$ and $\upsilon$. For this, let us assume that the disc surface density follows a power law: $\varSigma =\varSigma _0(R/R_0)^{-p}$, that $\alpha$ and $H/R$ are constant with radius, and that the disc is vertically isothermal. Under these assumptions, the contributions to the mass accretion are

(11.43)\begin{equation} \left.\begin{gathered} \tau_R=\alpha \varOmega_0^2R_0^2\left(\dfrac{H}{R}\right)^2(1-p)\varSigma\\ \tau_z=\frac{1}{\sqrt{2{\rm \pi}}}\upsilon\varOmega_0^2R_0^2\dfrac{H}{R}\varSigma \end{gathered}\right\} \end{equation}

so that the ratio simply reads

(11.44)\begin{equation} \varLambda=\frac{1}{\sqrt{2{\rm \pi}}(1-p)}\frac{\upsilon}{\alpha}\frac{R}{H}. \end{equation}

The scalings (9.6) and (9.8) then suggest that $\varLambda \rightarrow 0$ when $\beta _{\mathrm {mean}}\rightarrow \infty$. Combining these results with (11.40), this implies that the ejection index $\xi$ tends to decrease as $\beta _{\mathrm {mean}}\rightarrow \infty$, as one would expect.

12. Outflow phenomenology

The launching of an outflow is tightly linked to physics of the accreted material because the outflow energy eventually comes from the accretion power of the disc. Let us divide the overall structure into a ‘disc region’ and an ‘outflow’ region (figure 42) and describe the physical processes in each region.

Figure 42. Global disc–wind interaction scheme. We distinguish a disc region in dashed lines from the outflow region. The poloidal field line is represented in red and the poloidal streamlines in green. In addition, the toroidal field component is shown in blue in a frame corotating at $\varOmega _K(R_0)$.

12.1. Disc region

In the disc region, the flow is mostly accreted, whereas the field lines are stationary. The fact that poloidal field and streamlines are not parallel implies that non-ideal MHD effects are necessarily present in the disc. Historically, these non-ideal effects have been treated in two ways: (i) assuming that the disc ionisation fraction is extremely low as in PPDs, so that ambipolar and Ohmic diffusion are very large (Konigl 1989; Wardle & Konigl 1993); or (ii) assuming that the disc was turbulent, the turbulence leading to a ‘turbulent resistivity’ modelled as a non-isotropic diffusivity tensor (Ferreira & Pelletier 1993).

The toroidal component of the magnetic field $B_\phi$ is the key ingredient of the interaction between the disc and the outflow. In the disc, the toroidal field is produced by the shearing of the radial field by the differential rotation of the disc $\partial _tB_\phi \simeq B_RR\partial _R \varOmega$. This effect is actually the main energy source for the outflow, which converts shear energy into magnetic energy stored in $B_\phi$ at the disc surface. It is often assumed that outflows are top-down symmetric to satisfy the global symmetries of the system, so that $B_R(z=0)=0$ and $B_\phi (z=0)=0$ follows. The bending of poloidal field lines in the disc implies $\partial _z B_R > 0$ for $z > 0$ and hence $\partial _zB_\phi <0$ up to the disc surface. The actual value of $B_\phi$ depends on the competition between the shearing of $B_R$ and the magnetic diffusion which damps the shear amplification. Overall, one can estimate $B_\phi ^+$ at the disc surface by balancing shear and diffusion

(12.1)\begin{equation} B_\phi^+\sim B_z\tan(\theta)\frac{\varOmega h^2}{\eta}, \end{equation}

where $h$ is the disc scale height, $\eta$ is the magnetic diffusivity and $\theta$ is the inclination of the poloidal field at the disc surface (see figure 42). Accretion of the disc material naturally follows from the profile of $B_\phi (z)$. The toroidal field is responsible for an azimuthal magnetic tension force $B_z\partial _z B_\phi /4{\rm \pi}$, which slows down the rotating material and leads to accretion. We can deduce the accretion rate from (4.21):

(12.2)\begin{equation} \rho u_R \sim \frac{B_\phi^+ B_z}{4{\rm \pi} \varOmega h}\sim \frac{B_z^2\tan(\theta) h}{4{\rm \pi} \eta}, \end{equation}

where $\rho$ is measured in the middle of the accretion flow (which, in general, is the disc midplane, but can also be off-midplane for dissymmetric outflows). This accretion is physically a result of a transfer of angular momentum from the accreted material to the toroidal field. In the end, the angular momentum is stored in $B_\phi ^+$ and is eventually used to launch the outflow.

The outflow base has to be replenished by the disc. Hence, a net positive vertical acceleration is required in the disc region to push material upward, even if at modest velocities. A careful examination of (11.4) allows us to isolate the role played by each term in the vertical acceleration: the magnetic pressure is necessarily directed downward because $B_R^2$ and $B_\phi ^2$ are both increasing functions of $z$. Gravity is also directed downward. The magnetic tension term is usually small (it involves a radial derivative whereas the other terms involve vertical derivatives, which are larger by a factor of $R_0/h$), and typically directed downward if we assume the most natural situation with $B_z$ decreasing with radius. Hence, the only term leading to an upward acceleration is the vertical thermal pressure gradient. The role played by the thermal pressure at the base of the outflow has often been missed. However, it has some important consequences. For instance, $B_\phi ^+$ cannot reach arbitrarily large values that would otherwise prevent thermal pressure from pushing materials to the wind base and vertically squeeze the disc. Ejection therefore requires $(B_\phi ^+)^2/8{\rm \pi} \lesssim P^+$ (Ferreira 1997).

A summary of the physics of the disc region is as follows.

1. (i) The region is non-ideal as the gas has to be allowed to stream through poloidal field lines which are stationary.

2. (ii) The toroidal field at the disc surface stores the energy and angular momentum taken from the accreted material.

3. (iii) The upward motion needed to replenish the outflow base is a result of the vertical thermal pressure gradient. This sets a limit to the amount of toroidal field that can be stored at the disc surface because magnetic pressure prevents this upward motion.

12.2. Outflow region

In order to understand the dynamics of the outflow from its base, it is convenient to introduce the Alfvénic Mach number $\xi$ defined as

(12.3)\begin{equation} \xi\equiv \frac{u_p}{V_{{A},p}}=\eta \sqrt{\frac{4{\rm \pi} }{\rho}}, \end{equation}

where the second equality is deduced from (11.13). In the outflow, we can expect a continuous acceleration of the flow, so that $\xi$ is an increasing function as one moves along a poloidal streamline, with $\xi \sim 0$ at the outflow base. We can eliminate $B_\phi$ in favour of the angular velocity by combining the magnetic surface rotation invariant (11.16) and the angular momentum invariant (11.19) to obtain

(12.4)\begin{equation} \varOmega=\varOmega^*\left(\frac{1-\xi^2(R_A/R)^2}{1-\xi^2}\right). \end{equation}

This expression looks singular for $\xi =1$, which corresponds to the Alfvén point. However, at this particular point, $R=R_A$ so that $\varOmega$ is actually smooth across this point. Second, in the limit of low Mach numbers $\xi \ll 1$, which corresponds to the base of the outflow, we can expand this expression to obtain

(12.5)\begin{equation} \varOmega\simeq \varOmega^*\left[1-\xi^2\left(\frac{R_A^2}{R^2}-1\right)+{O}(\xi^4)\right]. \end{equation}

Hence, the outflow is rotating at a constant angular velocity up to the point where $\xi ^2\simeq 1/(R_A^2/R^2-1)$. Physically, the angular momentum stored in $B_\phi ^+$ by the disc is progressively used to accelerate the flow via the azimuthal magnetic tension force, leading to the apparent solid rotation profile. This works until most of the toroidal field has been used and the angular momentum is all in kinetic form. On the opposite limit $\xi \rightarrow \infty$, we find $\varOmega \simeq \varOmega ^*R_A/R^2$, i.e. a constant angular momentum rotation profile, as expected.

The vertical acceleration of the outflow results from angular momentum conservation. As demonstrated previously, angular momentum stored in $B_\phi$ is converted into kinetic angular momentum. Hence, $B_\phi$ decreases along the streamline. This leads to a magnetic pressure force $\partial _z B_\phi ^2$ directed upward in (11.4) and, hence, a vertical acceleration of the outflow. As the outflow bends toward the vertical axis, $B_R$ decreases as well, leading to an additional magnetic pressure term $\partial _z B_R^2$ accelerating the outflow vertically. Overall, the vertical acceleration is a magnetic pressure effect owing to the decrease of $B_\phi$ and $B_R$ along the streamlines.

One can use the Bernoulli invariant to characterise the topology of the outflow close to the wind base. For simplicity and following Blandford & Payne (1982), let us assume that the wind base is located at $(R_0,z=0)$. We consider a fluid particle, initially following a Keplerian rotation orbit at $(R_0,z=0)$ and we assume this particle follows the streamline of the outflow so that at a later time, the particle is located at $(R_0+\delta R,\delta z)$. As we focus on the launching region of the outflow, we assume that $\varOmega \simeq \varOmega ^*$. During this displacement, the Bernoulli invariant should be conserved, so we should have

(12.6)\begin{align} -\frac{1}{2}(\varOmega^*R_0)^2+\psi_G(R_0,0)& > -\frac{1}{2}[\varOmega^*(R_0+\delta R)]^2+\psi_G(R_0+\delta R,\delta z)\nonumber\\ &\rightarrow \frac{1}{2}(\varOmega^*R_0)^2 \left(-3\left(\frac{\delta R}{R_0}\right)^2+\left(\frac{\delta z}{R_0}\right)^2\right) < 0, \end{align}

where the inequality comes from the assumption that $u_p(R_0+\delta R, \delta z)>u_p(R_0,0)$ (i.e. the flow accelerates). It can be converted into a criterion on $\theta$: $\tan \theta > 1/\sqrt {3}$ or $\theta >30^\circ$. This well-known criterion, initially derived by Blandford & Payne (1982), is valid provided that thermal effects (enthalpy, heating) are negligible, i.e. in cold winds. This critical angle is usually interpreted in the framework of the ‘magneto-centrifugal acceleration’ of Blandford & Payne (1982): magnetic field lines are assumed to behave as rigid poloidal wires on which fluid particles are drifting (‘bead on a wire’ analogy). If the field lines are inclined sufficiently, the particles can overcome the vertical gravity and are accelerated by the centrifugal force.

This simple physical interpretation based on the Bernoulli invariant close to the wind base is very useful as a first approach to outflow physics, but it should not be taken too strictly as it leads to several misconceptions about the very nature of outflows. Let us underline the main physical differences between the magneto-centrifugal acceleration picture and the processes at work in an outflow.

1. (i) Magnetic field lines do not behave as rigid poloidal wires. As shown previously, the magnetic configuration has a strong toroidal field at the wind base $B_\phi ^+$, which indicates that the field is wound up in this region.

2. (ii) The constant angular velocity approximation is only valid close to the wind base (and far from the Alfvén radius). It is a result of the toroidal field tension, which accelerates the flow azimuthally, hence the necessity of a wound field configuration. Winds having $R_A/R_0\sim 1$ will therefore experience a very small region of solid rotation (if any).

3. (iii) As discussed previously, vertical acceleration is always a magnetic pressure effect that does not rely on gravity, centrifugal acceleration or magnetic tension.

As proposed by Contopoulos & Lovelace (1994) and Ferreira (1997), it is therefore preferable to qualify outflows as ‘magnetically driven’ instead of ‘centrifugally driven’.

Let us summarise here the physics of a cold outflow close to the launching region up to the Alfvén point.

1. (i) The flow is accelerated azimuthally by the magnetic tension owing to $B_\phi$. This results in a transfer of angular momentum from $B_\phi$ to the gas as it accelerates.

2. (ii) As $B_\phi ^2$ (and, possibly, $B_R^2$) decreases with $z$, the flow is accelerated vertically by magnetic pressure.

3. (iii) Close to the launching point and far from the Alfvén radius, the outflow is approximately in solid rotation.

4. (iv) In a cold wind, energy conservation at the wind base implies $\theta > 30^\circ$.

13. Global numerical models

Many numerical models of PPDs have been published in the literature. Here, we focus on models tackling the effect of non-ideal MHD effects and winds on the dynamics of a disc, because we focus on the outer parts of PPDs. We therefore exclude ‘ideal MHD’ models and simulations without a large-scale magnetic field.

The first model to investigate this regime was published by Gressel et al. (2015). However, the limited vertical extension of this work (typically four scale heights) makes the interpretation of outflow properties difficult. Therefore, we focus on models with a larger vertical extension such as those published by Béthune, Lesur & Ferreira (2017) and Bai (2017).

13.1. Global topology

One of the main problems with shearing-box models is the presence of an odd symmetry for the solution, leading to difficulties to interpret the role played by the outflow. The use of global models avoids this problem because, in this case, the gravitational potential is not symmetrical with respect to $r$. Many models seem to converge towards dissymmetric outflows, i.e. outflows that are neither even or odd, but are essentially odd for $-3H < z < 3H$ and exhibit a strong current sheet on one side of the disc. In this configuration, the outflow is dissymmetric, and more mass ejected from one side of the disc than the other (with ratios of $\dot {M}_w$ reaching a factor of a few). This kind of solution was found both with only Ohmic and ambipolar diffusion (Gressel et al. 2015) and with all three non-ideal effects (Bai 2017; Béthune et al. 2017).

The presence of an odd symmetry in the midplane region is reminiscent of the shearing box solutions described in § 9.2. However, this symmetry is not verified far away from the midplane as global effects enter the scene. Gressel et al. (2015) have proposed that the global field structure was playing a significant role in shaping the outflow. Indeed, by choosing initially a field configuration with $\partial _RB_z > 0$, it is possible to produce outflows directed inwards (towards the star). It has therefore been proposed that the global radial magnetic pressure gradient was responsible for the symmetry breaking observed in these models. In addition, the vertical rotation profile (and therefore the vertical temperature structure) might be playing a role by shearing the poloidal field (Gressel et al. 2015). The question of the origin of the global outflow configuration therefore remains open (figure 43).

Figure 43. Global simulation for a wind in a PPD which exhibits a dissymmetric outflow. Black lines are poloidal magnetic field lines, green arrows represent the poloidal velocity and the background colour traces the azimuthal field $B_\phi$. Close to the midplane, the configuration has an odd symmetry, as found in shearing box models and a current sheet is found at $z\sim 3H$. Figure from Béthune et al. (2017).

Outflows are not always dissymmetric. Indeed, both Béthune et al. (2017) and Bai (2017) have reported symmetric (even) outflow configurations. These symmetric configurations seem to be found mostly when $B_z$ is weak enough ($\beta _\mathrm {mid}\gtrsim 10^4$) (Bai 2017; Béthune et al. 2017) and for anti-aligned cases ($B_z\varOmega < 0$). The sensitivity of the outflow configuration on the field polarity suggests that the HSI is partly responsible for the outflow configuration. However, there is no one-to-one correspondence between the field alignment/strength and the outflow configuration, and the choice of initial conditions also seems to be playing a non-negligible role. If true, it might be desirable to start from a magnetic configuration as close as possible to the configuration expected from core collapse calculations (Bai 2017).

13.2. Accretion

The question of the engine driving accretion can be directly addressed in global simulations. Indeed, one can measure individually each term in the angular momentum conservation equation

(13.1)\begin{equation} \dot{M}_\mathrm{acc}=\frac{4{\rm \pi}}{R\varOmega_K}\left\{\underbrace{\frac{\partial}{\partial R} R^2\left[\overline{\rho v_\phi v_r}-\frac{\overline{B_\phi B_r}}{4{\rm \pi}}\right]}_{\tau_r}+\underbrace{\left[R^2\rho v_\phi v_z-R^2\frac{B_\phi B_z}{4{\rm \pi}}\right]_{z=-h}^{+h}}_{\tau_z}\right\}, \end{equation}

where we have defined the mass accretion rate $\dot {M}_\mathrm {acc}=-2{\rm \pi} R\overline {\rho v_r}$ and the radial and vertical torques $\tau _{r,z}$. An example of such a measure is given in figure 44. We find that accretion is mostly a result of the wind, the surface torque being the main contribution to angular momentum extraction in the disc.

Figure 44. Measured accretion rate as a function of time in a non-ideal model. The accretion rate and torque contributions have been average radially. Most of the accretion is a result of the wind ($\tau _z$ term) whereas the radial stress does not seem to contribute to the angular momentum budget. Figure from Béthune et al. 2017.

Note, however, that the fact that the radial torque is negligible does not imply that the radial stress is also negligible. As can be seen from (13.1), one can cancel the radial torque if the radial stress is proportional to $1/R^2$. In Béthune et al. (2017), a strong laminar radial stress is indeed present in the disc, with effective $\alpha$ values reaching a few times $10^{-2}$. However, because of the surface density profile, the net torque exerted on the disc mostly cancels out. Similar behaviours have been obtained by Bai (2017), with disc regions exhibiting positive and negative torques.

Overall, it is clear that in these models, the wind torque is playing a significant if not dominant role in the mass accretion rate. The radial torque, on the other hand, is less straightforward as it depends on the initial conditions chosen for the model. Accretion, decretion or both can be obtained from the radial torque, despite the presence of a relatively strong positive laminar radial stress when the field is aligned with the vertical rotation axis, as in shearing box models (see § 9.2.2).

Despite the uncertainties, these models predict mass accretion rates $\dot {M}_\mathrm {acc}\sim 10^{-8}-10^{-7}\,M_\odot /\mathrm {yr}$. However, even these values should be interpreted with care, as they are usually measured in the middle of the simulation domain, typically 5–20 AU from the central star. If an outflow is indeed present and carrying mass away, the mass accretion rate onto the star can be significantly smaller than that derived in the bulk of the disc. From the definition of the ejection efficiency (11.33), we have

(13.2)\begin{equation} \dot{M}_\mathrm{acc}(R)=\dot{M}_\mathrm{acc}(R_0)\left(\frac{R}{R_0}\right)^{\xi} . \end{equation}

Large ejection efficiencies ($\xi =O(1)$) such as that found in recent global models therefore lead to dramatically reduced mass accretion rates at the inner radius of the disc.

Unless one assumes an ionisation rate much higher than that expected in these objects (see § 3.4), the flow is mostly laminar, with a very low time-dependency. Hence, the radial stress measured in these models is not the usual turbulent stress found in ideal simulations, but really a purely magnetic term with no velocity counterpart. This implies that dust grains present in the disc will be less subject to turbulent fluctuations and will therefore settle towards the disc midplane more rapidly.

13.3. Ejection and mass loss rate

The outflow is not only responsible for carrying angular momentum away from the disc, but it also contributes significantly to the mass loss of the disc. All of the simulations published up to now find that the mass loss rate in the outflow is, broadly speaking, comparable with or even larger than the mass accretion rate in the disc ($\dot {M}_\mathrm {w}\gtrsim 10^{-8}\text {--}10^{-7}\,M_\odot /\mathrm {yr}$), which implies $\xi \sim 1$.

The mass outflow rate is tightly connected to the amount of flux threading the disc, with $\dot {M}_\mathrm {w}\propto \beta _\mathrm {mid}^{-1/2}$ (Béthune et al. 2017), indicating that the mass flux is proportional to the magnetic flux threading the disc. Interestingly, similar scalings are obtained in non-ideal shearing box models (Bai & Stone 2013b; Lesur et al. 2014) whereas steeper dependences are found in ideal shearing box models (Suzuki & Inutsuka 2009; Bai & Stone 2013a).

The engine driving ejection can be isolated first by looking at the magnetic level arm of these outflows. Most of the models published up to now find level arms $\lambda < 2$. This is coherent with the very high mass loss rates found in these simulations (high $\xi$, see (11.40)). Although it is, in principle, possible to obtain cold outflows with low $\lambda$ and high $\xi$ in discs threaded by a weak field ($\beta _{\mathrm {mean}}\gg 1$, see for instance Jacquemin-Ide, Ferreira & Lesur 2019), some of the outflows published to date in PPDs have $\lambda < 3/2$ (e.g. Béthune et al. 2017), which violate the cold MHD wind constraint (11.29). Hence, these outflows are not purely magnetically driven.

This conclusion can also be reached by analysing directly the Bernoulli function of the outflow (e.g. figure 45). Such an analysis shows that thermal effects (enthalpy and heating terms) both contribute significantly to the energetics of the outflow (Béthune et al. 2017). Still, magnetic effects are clearly not negligible at the base of the outflow, where magnetic pressure helps pushing the flow upward.

Figure 45. Bernoulli invariant $\mathcal {B}$ measured in an outflow driven from a non-ideal PPD. Magnetic contribution are $B_\perp ^2v$ and $w$, thermal contributions comes from the enthalpy $\mathcal {H}$ and external heating $\mathcal {Q}$ whereas kinetic energy terms are represented by $v_\phi ^2$ and $v_p^2$. Magnetic terms contribute significantly close to the launching point, whereas thermal energy becomes important higher up in the outflow. The ideal MHD region starts for $z\gtrsim 5h$. Figure from Béthune et al. (2017).

For these reasons, these outflows have been labelled ‘magneto-thermal’. This kind of outflow has already been identified in self-similar solutions by Casse & Ferreira (2000). Compared with historical cold wind solutions, they are (obviously) warmer, denser and slower. They reach high $\xi$ values (typically $\xi > 0.1$) and have moderate $\lambda$. Of course, the fact that they extract angular momentum from the disc and that the initial acceleration is a result of magnetic effects implies that they are not purely thermal.

13.4. Self-organisation

Self-organisation was unambiguously identified in global simulations by Béthune et al. (2017) for simulations with $\beta _{\mathrm {mid}}\lesssim 10^3$ (e.g. figure 46), but it is absent from the models of Gressel et al. (2015) and Bai (2017), who only considered $\beta \gtrsim 10^4$. It is most of the time found in simulations exhibiting dissymmetric wind configurations, but does not seem to prefer one given field polarity, as would be expected from Hall-driven self-organisation. A careful examination of the flow shows that the poloidal field lines are concentrated in low-density regions. Hence, magnetic effects are playing a very important role in the mechanism.

Figure 46. (a) Self-organisation in a simulation with $\beta _{\mathrm {mid}}=10^2$ computed in 2.5 dimensions. The density is represented in colormap whereas magnetic field lines are in white lines and velocity field is shown in green arrows. Note that field lines are accumulated in regions of reduced density in the midplane. Figure from Béthune et al. (2017). (b) Volume rendering of a similar model, this time computed in the full three dimensions. Note that the flow remains axisymmetric.

It is possible to identify which process is responsible for self-organisation by looking closely at the non-ideal induction equation. It is then found that ambipolar diffusion is the only term responsible for the accumulation of magnetic flux in narrow regions, whereas Ohm and Hall effects are both diffusing the field away (Béthune et al. 2017). Hence, despite the presence of a rather strong Hall effect in these simulations, it is not the Hall-driven self-organisation that is at work in these models, but ambipolar-driven self-organisation. In essence, the mechanism seems to be similar to that driving self-organisation in stratified shearing box models subject to ambipolar diffusion only (Bai 2015). The local configuration found in the global simulations is indeed identical to the configuration found in shearing boxes (see figure 40), making the shearing box model a valuable tool to understand self-organisation in this regime.

Unfortunately, there is today no general theory predicting in which situation self-organisation is occurring nor what are the general properties of the structures that are formed.

14. Conclusions

In conclusion, the global modelling of the weakly ionised part of PPDs ($R\gtrsim 1\ \mathrm {AU}$) is still in its infancy. The inclusion of non-ideal MHD effects confirms that the disc midplane is mostly laminar, whereas the disc is still accreting thanks to magnetised outflows. Owing to the low magnetisations used in these models ($\beta _{\mathrm {mid}}\gtrsim 10^2\text {--}10^3$), the outflow is not purely magnetically driven as external heating and thermal pressure are found to contribute to the global energy budget. This implies that it is difficult to draw systematic conclusions regarding the accretion rate and mass loss rate, as these can depend both on the ionisation structure but also on the thermodynamics of the disc wind, both of which are plagued by huge uncertainties.

It would be tempting to consider stronger field models with $\beta =O(1)$ which are known to lead to historical ‘cold wind’ models (Ferreira 1997). However, such models would also have larger average accretion velocities which are typically sonic. If one wants to keep an accretion rate approximately compatible with observation, such a model would imply a disc surface density reduced by at least one order of magnitude (most likely two), which would be incompatible with the disc masses inferred from observations. Hence, the fact that the disc is massive, with an average accretion velocity much smaller than the speed of sound implies that if outflows exist in these regions, they must be due to a weak field: $\beta _{\mathrm {mean}} \gg 1$.

Self-organisation is also a very promising mechanism to explain some of the observables in PPDs. However, the theoretical background for these features remains limited. The fact that they are seen both in global and local simulations is encouraging, but, clearly, a detailed theoretical work is needed before any satisfactory prediction can be made and tested against observational data.

PART FIVE: Summary and future directions

15. Summary

These are exciting times for PPD modelling and planet formation theory in general. Indeed, we now start to have direct observational constraints for these astrophysical objects, ranging from large resolved structures, such as rings and non-axisymmetric bumps, to turbulent velocity dispersion measurements. Even magnetic field strength and topology are now beginning to be probed at the disc scale, giving more constraints to models.

The plasma in these discs is, however, relatively cold and therefore weakly ionised, reaching an ionisation fraction of $\xi \sim 10^{-13}$ in the disc midplane around 1 AU. This implies that non-ideal MHD effects (Ohmic, ambipolar and Hall effect) are essential to obtaining a proper description of the plasma. The amplitude of these effects, however, is poorly constraint. Because ionisation is mostly non-thermal, the amplitude of these non-ideal effects depends on the details of the ionisation sources, the disc structure and the plasma composition (especially the abundance of tiny dust grains). Overall, there is an uncertainty of several orders of magnitude on these effects, implying that very detailed models including complex reaction networks are likely unnecessary at this stage because the input parameters (disc composition and environment) are largely unknown.

In order to explain the observed accretion rates in these discs, and because angular momentum is a conserved quantity, one needs to find a way to remove the disc angular momentum, either by transporting it radially outwards in the disc bulk, or by transporting it vertically away in a magnetised wind. Although radial transport has been historically favoured thanks to its simplicity and elegance (the well-known $\alpha$ disc model), vertical transport is now believed to be key in several astrophysical objects because of its high efficiency.

The most favoured mechanism to explain angular momentum transport in discs is the MRI, a linear, ideal MHD instability found in rotating sheared flows. This instability has been the subject of intense studies since the early 1990s and it is known that it is strongly affected by the non-ideal effects present in PPDs. Most notably, it is suppressed in the regions where Ohmic and ambipolar diffusion are strong, and it gives a new branch, known as the HSI, when the Hall effect is dominant, in the case where the poloidal field is aligned with the rotation axis.

In the non-linear regime, the MRI behaviour strongly depends on the presence and strength of a mean vertical field threading the disc. Historically, most of the simulations published until early 2010 were in a regime without a mean field, commonly known as the ‘MRI dynamo’. In ideal MHD, this regime is known to produced vigorous 3D turbulence and radial transport of angular momentum, but in non-ideal MHD, turbulence is suppressed, leading to a laminar flow, no angular momentum transport and no accretion. This ‘dead-zone’ problem is circumvented by considering a mean vertical field threading the disc. Doing so, the MRI can indeed be revived in the upper layers, as expected from the linear analysis, but it then saturates into magnetised outflows, and the flow remains mostly laminar. In this situation, angular momentum is transported, mostly in the vertical direction, so accretion is saved, but its physical description then becomes fairly different from that of an $\alpha$ disc.

This connection between the MRI and magnetised outflow in discs threaded by a mean field was only realised during the past 10 years. In the presence of outflows, local models are insufficient since the dynamics of outflow is dictated by the global geometry of the system. In the case of PPDs, it is found that accretion driven by magnetised winds are compatible with observed accretion rates for relatively weak mean fields, $\beta _{\mathrm {mean}}=O(10^4)$, which are compatible with the upper bounds on the field strength from observations. Because the outflow is in a weakly magnetised regime, it is usually found that the mass loss rate can be of the order of, or even larger than, the mass accretion rate measured at the inner radius of the disc (this result being perfectly consistent with mass, angular momentum and energy conservation). In addition, the outflow is highly sensitive to thermal effects, which contribute significantly to the flow energetics, as is found in many models. Hence, these outflows have been called ‘magneto-thermal’.

16. Perspectives

Research on this topic is now following several paths. First, the fact that thermal effects can play a significant role implies that they should be modelled accurately. Several groups are now working actively on this problem (e.g. Wang, Bai & Goodman 2019; Gressel et al. 2020). It should be realised that thermal driving is a very complicated problem, as it involves heating by X-rays and UV photons, in addition to cooling, mostly by molecular and atomic lines. The computation of thermal processes therefore rely on complex chemical networks, coupled to radiative transfer codes, which are all computationally very intensive. Eventually, one hope is to find a way to simplify this physics using prescribed heating and cooling functions tested on complete models. This would allow a more systematic exploration of the long-term impact of thermodynamics on these systems, and make a connection to the winds observed in the sub-millimetric range.

A second question is the dynamical evolution of the mean magnetic field threading the disc. As shown previously, this mean field is key for magnetised outflows. It is strongly suspected that such a field should be present, as a direct result of the disc formation process, which relies on the collapse of a magnetised molecular cloud. During the collapse, a fraction of the magnetic field is trapped in the forming disc, and then plays the role of the mean field for outflows. However, once the disc is formed, it would be desirable to describe how this mean field evolves with time. It could be advected inwards by the accretion flow, leading to a strongly magnetised inner region, or it could inversely diffuse outwards because of non-ideal MHD effects. At the time of writing, this question is not settled, even qualitatively. Numerical models suggest that the field is diffusing outwards (Bai & Stone 2017; Gressel et al. 2020) whereas analytical models suggest inwards transport (Leung & Ogilvie 2019). If magnetised outflows are the dominant mechanism of accretion, then it is essential to address quantitatively this question in order to be able to model the long-term dynamics of these discs, because mass and magnetic flux are tightly linked. This flux transport can be at the origin of complex dynamics in the disc, such as time variability or even eruptions, which are also observed in these systems.

A third axis of research is the effect of this dynamics on planet formation, from the dynamics and growth of dust grains to the migration of giant planets. Most of the literature published to date rely on the $\alpha$ disc paradigm, which itself assumes that the disc is turbulent. However, if accretion is driven by magnetised winds in a mostly laminar flow, this framework has to be revised. Indeed, the lack of turbulence affects the dynamics of large grains (${\gtrsim }10\ \mathrm {\mu }\mathrm {m}$): vertical and radial settling, coagulation and disruption efficiency, etc., are all strongly modified. In addition, the fact that accretion is driven by surface stress, and not by turbulence is also going to reshape planet migration. It is not clear yet how type I and type II migration processes react to this shift in accretion paradigm, but the first attempts at including the wind stress in 2D planet migration numerical models already show a very significant effect (Kimmig, Dullemond & Kley 2020). Even the long-term evolution of the disc is quite different from that of a viscous disc, as viscous spreading is absent for a wind-driven disc. These are only a few example, but it shows that many things which were thought to be well established are now standing on wobbly foundations.