{'fraction_non_alphanumeric': 0.03481592017639591, 'fraction_numerical': 0.015425733312545936, 'mean_word_length': 4.385921980711099, 'pattern_counts': {'":': 1, '<': 57, '<?xml version=': 0, '>': 48, 'https://': 0, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 18, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Magnetohydrodynamics of protoplanetary discs

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>1 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 K) and dense, which implies a very low ionisation fraction close to the disc midplane. In these notes, I deliberately ignore the innermost R~0.1 AU region which is influenced by the star-disk 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 MHD, 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 which 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 since magnetised winds are also emitted from the surface of these objects. After a short introduction on winds physics, I present global models of protoplanetary discs, including both a large-scale wind and the non-ideal dynamics of the disc.

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 1 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 K) and dense, which implies a very low ionisation fraction close to the disc midplane. In these notes, I deliberately ignore the innermost (R ∼ 0.1 A.U.) region which is influenced by the star-disk 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 MHD, 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 which 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 since magnetised winds are also emitted from the surface of these objects. After a short introduction on winds 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 & physical description

Observational context
Recent years have seen a dramatic change in our understanding of protoplanetary discs, 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 [a.u.]) 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 protoplanetary discs is much more complex and rich than initially anticipated. In this part, we review the most recent evidence for proto-planetary disc structure and evolution, which can be used to constrain the most recent theoretical models.

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 protoplanetary disc 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 (Fig. 1). The disc outer edge can range from 100 AU to more than 1000 AU depending on the object under consideration.
Observations typically probe the following regions: • 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.
• The near and mid-infrared continuum (also known as infrared excess) is due to stellar photons scattered by small dust grains (typically less than a µ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.
• The (sub)-millimetre continuum probes the thermal emission of bigger dust grains (typically with a size of the order of a 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.
• 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.

Accretion
Because the thermal equilibrium of protoplanetary discs for R 1 A.U. 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. respectively to 0.65 dex theṀacc variability of tored on timescales of 1, nt of accretion from the y, the authors inferred a theṀacc variability in sing shorter ( hour to Costigan et al. (2014). es and their variability le of 14 objects; targets , in some cases iterated few nights. Combining monitoring with multiell as with results from at the dominant contriacc arises from the days ates ranging from 0.04 t variability on day-tonificantly on the large verageṀacc-M * trend; y linked to the differobjects (see Sect. 5.3). ge distributions for the non-accreting (WTTS) sue of an accurate age ests that an intrinsically ss and age, which may .g., the properties of the of the disk). e estimates for theṀacc timescales (this study; ear timescales (Nguyen consistent; this suggests iability over these basefrom the mid-term variorted by the results of een the L exc UV measuret al. (2002)'s survey of rve in this case is that, ed at the two epochs are about the equality line, average amount of vari-CFHT monitoring.
ter spread inṀacc ? provides a debated indication on nominal ages. In order to probe the evolution of accretion properties with stellar age, we inferred an age estimate for each source from isochrone fitting, adopting the models of Siess et al. (2000). The results for accreting sources are displayed vs. theṀacc estimates in Fig. 19. Only objects which are actual detections inṀacc and whose age estimate from models is between ∼10 6 yr and 10 7 yr have been retained for the comparison.
An overall decreasing trend of the averageṀacc with increasing stellar ages is traced in Fig. 19, consistently with previous results regarding the average evolution of accretion with stellar age (e.g. Hartmann et al. 1998;Sicilia-Aguilar et al. 2010;Manara et al. 2012). The anticorrelation trend is well described in the M * < 1 M⊙ regime, reaching a significance of over 5σ according to a Spearman's rank order test (Press et al. 1992). A linear regression to this group of objects yieldṡ D. Fedele et al.: Accretion Timescale in PMS stars   He I 5876 Å in emission (EW = −0.5 Å, -0.6 Å respectively). The evidence of large Hα10% together with the He I emission is most likely due to ongoing mass accretion, and these two stars are classified as accreting stars. We estimate a fraction of accreting stars in NGC 6231 of 11/75 or 15 (±5%). We warn the reader that this might be a lower limit to the actual fraction of accreting stars; further investigation is needed to disentagle the nature (accretion vs binarity/rapid rotation) of the systems with large Hα10% (>300 km s −1 ) but low EW [Hα].

NGC 6531
We identified 26 cluster members in NGC 6531 based on the presence of Hα emission and presence of Li. 13 other sources show presence of Li 6708 Å, but have Hα in absorption. As in the case of NGC 6231, these might be cluster members with no or a reduced chromospheric activity. We measured the EW of Li 6708 Å of these 13 sources and compared them with the typical EW of the 26 stars in NGC 6531 showing also Hα emission Page 5 of 7 Figure 2. (left) 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).
(right) 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 first is the observational signature of accretion columns at the stellar surface. These accretion columns are formed when the disc material gets lifted up 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 Atomic Jet

Molecular outflow
Disc seen edge on 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 protoplanetary disc seen edge-on. Courtesy of F. Louvet. 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 /yr with uncertainties of the order of an order of magnitude depending on the object under consideration † (e.g. Fig. 2-left). These accretion rates tend to decrease over timescales of a few million years.
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 .
By combining this information, and assuming that accretion is approximately constant during the lifetime of these objects, one deduces that typical protoplanetary disc masses range from 10 −3 M to 10 −1 M , which is consistent with mass inferred from the total dust content of the disc (Andrews et al. 2013).

Ejection: winds and jets
Protoplanetary discs 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 km.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 ∼ 1 − 10 km.s −1 (Fig. 3). They could are indicated by arrows in Fig. 1: (1) a spiral arm located in the south and east (hereafter, the SE spiral) at deprojected distances r ⇠ 0.26 00 0.53 00 ; (2) an arc in the west at r ⇠ 0.2 00 -0.3 00 ; (3) a spiral arm in the northwest (hereafter, the NW spiral) at r ⇠ 0.22 00 -0.34 00 ; (4) an arc located north of the NW spiral at r ⇠ 0.32 00 -0.53 00 ; (5) substructures inside the spirals at r ⇠ 0.1 00 -0.2 00 and (6) a substructure branching o↵ the SE spiral at r ⇠ 0.21 00 -0.27 00 . We note that although there is some clear scattered light signal very close to the coronagraph, the exact morphology of feature (5) may be a↵ected by the e ciency of the AO correction and the coronagraph centering. Appendix A provides the U images ( Fig. A.1) and S/N maps (Fig. A.2). The features are detected above the disk background at more than 3 and do not appear in the U images. This suggests that they are real and not di↵raction residuals that would a↵ect the Q and U images alike. Although one could naturally see a single spiral arm in the west, extending from PA ⇠ 180 to 350 , the radial mapping ( Fig. 1, right) shows that it consists of different features that have a very distinct dependence on the PA, with a sharp transition at PA ⇠ 240 . In the December 2014 image, feature (2) hosts an area that is ⇠1.5 times brighter than the NW spiral, at PA ⇠ 240 260 . In the March dataset, the area at PA ⇠ 260 in contrast seems fainter than the background spiral. Since it is very close to the residual spiders, it is di cult to conclude whether it is an instrumental e↵ect, or if it has a physical origin which, considering the timescale of the variations, would point to a variability in the inner disk. While the detected structures agree with the HiCIAO H-band image (Grady et al. 2013), we did not detect a significant di↵erence in the extent of the scat-only fit the corresponding dataset. scattered light traces small dust gr the gas, and thus, that it indirectly prescription of Rafikov (2002) and wake in polar coordinates (r, ) fol where ↵ and are the exponen sound speed profiles, respectively c s (r c )/r c ⌦(r c ) is the ratio of the p radius (also called disk aspect ratio (r c , c ). The pitch angle of the spir perature (hence on the aspect rati launching planet. The flaring inde variation of H/r with radius. Equ ear or weakly nonlinear regimes, i is launched by an embedded plane shape of the density wave given b 2002).
We attempted to fit both spira that they are launched by two plane radii in the disk. We assumed that tion and fixed ↵ = 3/2. Varying and we fixed = 0.45, following A be due to the interaction of the jet with its environment, or they could be a genuine outflow component, emitted from the disc at R 1 A.U..

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 (Fig. 4).
The second class of observation is based on interferometry at millimetric and submillimetric wavelengths. The ALMA observatory has been very successful at probing the very structure of protoplanetary discs with incredible resolution and unexpected results (Fig. 5, Andrews et al. 2018).
Although these observations probe the dust distribution in the disc, they also tell us about the gas distribution and dynamics, because the grains which 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 et al. 2019).
All these observations indicate that discs are not smooth and symmetrical. They are instead structured on length scales comparable to our solar system. Structures † These structures trace the disc surface and not the column density 8 G. R. J. Lesur executions. The AO staff. After vapor radiometer into 2 s integraand times. The w was calibrated e amplitude scale 1107-4449. The librated using the lthough images e of artifacts and t (1.5 Jy), folding er density of short truction. val ALMA data , 2012 November and 34 antennas first two observadows centered at latter had two z), one 469 MHz Hz window (at a gain calibrator, and J0522-364 33 (December 31) weather for these of 0.5-1 mm. The inutes. The basic ck, we compared et on overlapping nsistency. ation were shifted target and then tial emission from ere made with a uum images at a ere generated by ing with a multiof the CLEAN esized beam. All with the CASA an analysis of two et. The first used a f 0) to provide a 78°). While this he cost of a dirty coverage at long second image was elliptical taper to e sidelobes. Both differences) and .
of the 870 μm revealing a series radial distance of pole-on viewing this substructure, polar coordinates and azimuthally averaged into a radial surface brightness profile.
The inner disk includes an unresolved (<0.5 au in radius) 0.93±0.04 mJy source coincident with the stellar position and a bright ring that peaks at 2.4 au; between them is a dark annulus centered at 1 au. The bright ring and dark annulus are unresolved (<1 au across). Because it is unresolved, the depth of the dark annulus is difficult to determine unambiguously: we find at least a 30% brightness reduction.   optically thick and thin regions, driving any derived values of dust opacity spectral index to be lower than reality.

Compact Spectral Lines
At 1″ . 1 resolution, the HCO + (1-0) emission shows a morphology similar to the 12 CO(1-0) shown in Figure 1( At 2.9 mm, we h time in the millim with an inclination Figure 1(c) and Ta  Figure 2(e)). In panel (c) the gray region delineates the statistical 1σ spectral index uncertainty; it does not account for the absolute flux uncertainty. For the fully fitted rings, panels (b) and (c) show dashed lines for the locations of dark rings and solid lines for bright rings. The fitted offsets from the continuum peak (Table 2) have been taken into account. Panel (d) shows the same cross-cut as panel (b) but on a Planck brightness temperature (TB) scale. Panels (b) and (d) are shown on a log scale on the yaxis. The dashed curve on panel (d) shows a representative power law for TB as a function of radius with an exponent of −0.65 extending from the B1 peak; it is not a fit.

7
The Astrophysical Journal Letters, 808:L3 (10pp), 2015 July 20 eak emission of~390, and n the opposite r of 130. The in the north of ast make the me than earlier spectral slope .44 mm comencies n (14)] ggesting that ominate the However, the rom the same rian disk prowith an inner 1B). 12 CO 6-5 0 AU radius, ring at 30 AU eed still some gnificant drop f~25 T 5 AU. nd continuum ut the relative ate that large ed in an asymll micrometert the ring. To echanism that he distribution long-lived gas azimuthal diles get trapped respect to the ssure gradient: re v → dust and ities and p is disks without oints inward, ve-mentioned r, there exists ximum of the ere ∇ → p gas ¼ 0 would conmain trapped ard drift and se small dust the gas, they trated toward azimuthal dis mechanisms oduce a local nstance, when ubstellar comisk, hindering , the presence ly speculative, B C Fig. 1. IRS 48 dust and gas observations. The inclined disk around IRS 48 as observed with ALMA Band 9 observations, centered on the star (white star symbol). The ALMA beam during the observations is 0.32′′ × 0.21′′ and is indicated with a white ellipse in the lower left corner. (A) The 0.44-mm (685 GHz) continuum emission expressed both in flux density and relative to the root mean square (rms) level (s = 0.82 mJy per beam). The 63 AU radius is indicated by a dashed ellipse. (B) The integrated CO 6-5 emission over the highest velocities in contours (6,12,...,60s CO levels, s CO = 0.34 Jy km s −1 ): integrated over -3 to 0.8 km s −1 (blue) and 8.3 to 12 km s −1 (red), showing a symmetric gas disk with Keplerian rotation at an inclination i = 50°. The green background shows the 0.44-mm continuum. The position angle is indicated in the upper right corner. (C) The Very Large Telescope Imager and Spectrometer for the midinfrared (VISIR) 18.7-mm emission in orange contours (36 to 120s VISIR levels in steps of 12s VISIR , s VISIR = 0.2 Jy arc sec −2 ) and orange colors, overlayed on the 0.44-mm continuum in green colors and the 5s contour line in green. The VISIR beam size is 0.48′′ in diameter and is indicated with an orange circle in the bottom right corner.
www.sciencemag.org on July 5, 2013 www.sciencemag.org Downloaded from parent molecular cloud (22), accretes material from its surrounding disk with a mass accretion rate of 8 × 10 −8 M ⊙ year −1 (25). Previous observations at 0.6″ to 1.1″ angular resolution were well described by a smooth and axisymmetric distribution of material in the disk that extends near the stellar photosphere and decreases monotonically with distance from the star (20, 23).
To estimate the optical depth of the observed dust continuum emission, we performed radiative transfer calculations using RADMC-3D (26) at 1.3 mm (17), using the previous surface density constraint found for the Elias 2-27 disk (20). This model reproduces the azimuthally averaged radial profile of the observed ALMA 1.3-mm continuum emission (fig. S1) (17). At a radial distance from the star (hereafter referred to as radius and denoted by R) larger than R ≈ 10 AU, the emission is optically thin and thus traces the density of solid material down to the midplane of the disk. At the location of the spiral structures (from R = 100 to 300 AU), the azimuthally averaged optical depth t of the dust continuum emission is t = 0.1 at R = 100 AU, decreasing to t = 0.02 at R = 300 AU (fig. S1B), which is consistent with the measured peak brightness temperature on the spirals of 1.2 K at R = 150 AU.
The spiral structures are even more evident in Fig. 2A, in which the data has been projected into a polar coordinate grid that accounts for the viewing geometry of the disk. In polar coordinates, a ring with zero eccentricity would have a constant radius for all polar angles. However, shown in Fig. 2A are two bright structures that grow in radius from~100 to 300 AU as the polar angle increases. The brightest of these two structures lies northwest of the star, labeled "NW"; the spiral structure southeast of the star is labeled "SE." In Fig. 2B, we present the surface brightness contrast of the NW and SE arms, defined as the ratio between the peak of emission at the arm and the background surface brightness (17). We found that both arms have similar contrasts ranging between values of 1.3 and 2.5. The spiral arms reach their highest contrast at R = 150 AU, coinciding with the location in the disk where gravity has the most influence over thermal pres-   Partnership et al. (2015). Bottom left: horsehoe-like structure observed in Oph IRS 48 at sub-mm wavelengths (green, tracing mm-sized dust) and corresponding scattered light infrared emission (yellow, tracing µm size dust) from van der Marel et al. (2013). Bottom right: spiral structures seen at sub-mm wavelengths in the young and massive disc of Elias 2-27, from Pérez et al. (2016) 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 which are presumably in the final evolution stage of protoplanetary discs. All of these structures could be the signature of embedded planets perturbing the disc structure by gravitational interaction. However, other processes have been proposed which do not assume planets. One of the key questions is therefore whether or not these structures are necessarily a signature of embedded planets.

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 since the turbulent broadening of spectral lines is comparable to 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, since their thermal velocity is lower compared to lighter molecules at a given equilibrium temperature. High-resolution spectra obtained from ALMA for CO lines indicates that turbulence is very weak, or inexistent (Flaherty et al. 2015;Flaherty et al. 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 MHD turbulence which typically predicts δv 0.1c 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 ). The result is that 100 µ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 δv ∼ 10 −2 c s (α ∼ 10 −4 , see §4.4).

Magnetic fields
Evidence for magnetic fields in protoplanetary discs is scarce. Typical values are expected to be of the order of a Gauss at 1 au down to a few mG 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 kG 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 (10s of AU), attempts at measuring the field strength through Zeeman splitting in molecular lines have only led to upper limits, with B z < 0.8mG and B < 30mG (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 submillimetric 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 evidences 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.1G around 1 AU are inferred from remnant magnetisation in meteorites following this idea (Fu et al. 2014), while upper limits with B < 30mG in the region around 15-45 AU is deduced from the magnetisation of Comet 67P/Churyumov-Gerasimenko (Biersteker et al. 2019).

Fluid properties
Protoplanetary discs (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): Here, we have defined the main physical properties of a disc: its surface density Σ, which correspond to the usual mass density integrated in the direction perpendicular to the disc plane, its temperature T , and its angular velocity Ω around the central object. We also define for convenience a dimensionless distance from the central object, in astronomical unit: R AU ≡ R/1 AU. This simple model leads to a 0.04 M 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 ≡ P/ρ and using the vertical hydrostatic equilibrium to define the disc vertical scale height ( §4.2) H = c s /Ω, one gets:

Magnetic fields
Magnetic fields in PPDs are poorly constrained ( §1.6). It is widely believed that fields are largely subthermal: 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 β parameter In practice, β 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 β cannot reach a value lower than β ∼ 1, hence this value is actually a lower limit for the typical plasma β expected in these discs. It is possible to connect the field strength to the plasma β using the properties defined above and obtain: The upper bound B 10 mG for R ∼ 10 AU mentioned in §1.6 tend to suggest β 10 4 in these regions, which confirms that the field strength is expected to be strongly subthermal.

Fluid approximation
Protoplanetary discs 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 cm, we can estimate the cross section of neutral molecules as σ nn = 3 × 10 −16 cm 2 . This gives us an approximate mean free path mfp and collision frequency ω coll We therefore have mfp R and ω coll Ω, which validate the fluid approximation to describe the dynamics of protoplanetary discs 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, mfp increases significantly. One finds that mfp H when n 10 4 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. Fig. 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.

Grain population
The question of grains is of importance in PPDs. As it 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 ρ o = 3 g.cm −3 . The density of grains is therefore ρ grain = 6 × 10 −12 R −9/4 AU g.cm −3 n grain = 1.4 R −9/4 AU a −3 µm cm −3 In this last estimate, we have assumed that all the grains had the same size. This is an over-estimation since 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 protoplanetary discs. However, this order of magnitude estimate points to an important fact: the abundance of grains n grain /n ∼ 10 −14 a −3 µm . Hence, if grains are smaller than 1 µm, the typical ionisation fraction of protoplanetary discs ( 10 −14 ) suggest that grains are more abundant than free charge carriers. As we will see in §3.4.3, this has a huge impact on the plasma conductivity tensor as grains can become the main charge carriers.

Ionisation fraction
The ionisation fraction ξ ≡ 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 -10 −13 at 1 au to 10 −13 -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-Ray, UVs, cosmic rays and radioactive decay).

Sources of ionisation
Since we focus on the outer part of protoplanetary discs (R > 1 AU), the gas is mostly cold with T < 300 K. This implies that thermal ionisation (due 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 ζ: • X-ray ionisation due to bremsstrahlung emission from an isothermal T = 5 keV corona localised around the central protostar (Igea & Glassgold 1999;Bai & Goodman 2009, see their Eq. 21); • cosmic-ray ionisation with ζ CR = ζ CR,0 exp(−Σ/96 g cm −2 ) s −1 (e.g. Umebayashi & Nakano 1981) and ζ CR,0 = 10 −17 s −1 , corresponding to the interstellar value ; • radioactive decay with ζ rad = 10 −19 s −1 (Umebayashi & Nakano 2009). The amount of ionisation due to CRs is highly disputed. Some authors have proposed that due to the wind coming from the young star, CRs are magnetically mirrored from the protoplanetary disc, resulting in a significantly reduced ionisation rate due to CRs (ζ CR,0 ∼ 10 −20 s −1 , Cleeves et al. 2013). On the opposite, 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 ζ CR,0 ∼ 10 −13 s −1 (Padovani et al. 2018). Observations of TW Hya tend to suggest a low ionisation rate due to CR (ζ CR,0 10 −19 s −1 , Cleeves et al. 2015), though this is still highly model dependent. Because of 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 Fig. 6 the resulting ionisation rate following the disc structure presented in §2.1. We find that cosmic rays are shielded only in the innermost parts of the disc, where the column density goes above 100 g cm −2 . Most of the disc midplane up to z ∼ h has ζ ζ CR,0 , indicating that CRs are indeed the dominant source of ionisation in this region. Above z ∼ h, X-rays start to penetrate the disc and the ionisation rate rises.

A simple chemical model
To illustrate the typical ionisation fractions expected in protoplanetary discs, we follow Oppenheimer & Dalgarno (1974), Fromang et al. (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 : where ζ is the ionisation rate, δ is the dissociative recombination rate for molecular ions, δ r the radiative recombination rate for metal atoms, and δ t the rate of charge transfer from molecular ions to metal atoms. Following Fromang et al. (2002), we take δ r = 3 × 10 −11 T −1/2 cm 3 s −1 δ = 3 × 10 −6 T −1/2 cm 3 s −1 δ t = 3 × 10 −9 cm 3 s −1 In the absence of metals and grains, the rate equations admit a simple solution in steady state ξ = ζ δn n . (2.1) In the opposite metal-dominated limit, still without grains, one gets ξ = ζ δ r n n .
Since δ r δ, one clearly sees that the absence of metals leads to a dramatic decrease in the ionisation fraction (Fromang et al. 2002).

Typical ionisation fraction profile
Grain-free case, Metal-free case: combining (2.1) to 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 Σ = 1700 R −3/2 AU g cm −2 (Wardle 2007;Lesur et al. 2014). This makes the disc much denser in the inner part, resulting in a stronger shielding of cosmic rays and a lower ionisation fraction than less dense discs. As an illustration, we show in Fig. 7 the resulting ionisation fraction with the disc structure presented in §2.1 and with a MMSN disc model. We observe that the lowest ionisation fraction reach 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 CRs+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 AU since 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 protoplanetary discs (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 impact of grains, let us add the following reactions to our simplified reaction network: grain + m + → grain + + m grain − + m + → grain + m grain + e − → grain − grain + + e − → grain grain + M + → grain + + M grain − + M + → grain + M grain + + grain − → grain + grain The above reaction network only considers singly charged grains, while it is well known that grains can have many charges (Ilgner 2012). We chose this approach to illustrate in the simplest model the impact of grains on the ionisation fraction, and later on the diffusivities since 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 due to electrons and charged grains are presented in Fig. 8. We observe essentially two trends. First, the smallest ionisation fraction is found when grains are present, while the highest ionisation fractions correspond to grain-free metal-rich cases, with variations due to this composition effect of the order of 3 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 due 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 AU . In the following, we will use the value ξ = 10 −13 to evaluate several plasma parameters, keeping in mind this corresponds to a lower bound in our disc model.

Plasma description in protoplanetary discs
In this section we explore the properties of the plasma constituting protoplanetary discs and ask whether they can be described using non-ideal MHD. For this limit to be valid, we have to satisfy 3 criteria which we discuss below (i) Binary Coulomb interactions should be negligible. This implies that the plasma parameter (defined below) is much larger than 1.
(ii) Electro-neutrality is satisfied on timescales of interest, i.e. any charge separation is quickly eliminated by electrostatic interactions.
(iii) The behaviour of each fluid component (electrons, ions, neutrals, charged grains) can be described using a single fluid approximation.

Plasma parameter
Several quantities allow one to characterise a plasma, the first one being the nature of the electromagnetic interaction. The most fundamental quantity characterising a plasma is the Debye length which may be written in an electron-ion plasma where Z is the averaged number of charges on the ions and we have assumed electroneutrality so that n i = n e /Z. The Debye length is clearly below the scales of interest in protoplanetary discs. Even if one considers charged grains, the same Debye length can be derived since 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 that screening of short-range Coulomb interaction. This is quantified by the plasma parameter 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 10 3 in protoplanetary discs, a value which is never encountered, even in chemical models including grains (e.g. Wardle 2007).

Electro-neutrality and drag
Protoplanetary discs 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 multifluid 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 since the collision timescales are short, as demonstrated above. We therefore start from the following dynamical equations: where n j , m j , v j , P j , q j , 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 R j between this species and all of the other species. This force is due to inter-species collisions and can be written as since each fluid component is collisional and therefore has a Maxwellian velocity distribution. Here, γ = σv jk /(m j + m k ) and σv jk is the momentum exchange rate between species j and k. As expected from momentum conservation, we have j R j = 0. It is usually assumed that electro-neutrality follows from the fact that the plasma frequency ω p is much larger than any frequency of interest. While 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 other ones unperturbed. We moreover assume that the fluid pressure and other external forces are negligible compared to electromagnetic forces. The linearised equation of motion reads Solving these equations requires an equation for the electric field, which is obtained from one of Maxwell's equation We can combine these equations to obtain a second order relation on the density fluctuation where we have introduced the plasma frequency ω p and the stopping time τ s as Dynamically, this equation describes damped plasma oscillations with frequencies for which we can distinguish two physical limits.
(i) ω p τ −1 s in which case the plasma is subject to plasma oscillations at frequency ω p with a damping timescale equal to τ s . If we consider phenomena on frequencies much lower than ω p , we can average out the highest order time derivative and get a simple closure relation between v x and E x : v x = qE x /γmρ which constitutes the base of Ohm's law. Once these oscillations are time-averaged, the plasma can be assumed to be electrically neutral.
(ii) ω p τ −1 s in which case the plasma is subject to over-damped oscillations with two imaginary frequencies ω + = iτ −1 s and ω − = iω 2 p τ s ω + associated to two damping timescales τ ± = (ω ± ) −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. Because of the drag, however, it very rapidly reaches an asymptotic velocity, given by v x = qE x /γmρ. τ + 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 τ − to actually fill the charge deficit. This implies that Ohm's law, given by the asymptotic velocity, is valid on timescales longer than τ + , and that charge inertia can be neglected in that limit. However, charge neutrality is restored on the much longer timescale τ − . 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 max(τ s , ω −1 p ), and recover Ohm's law without time derivative. Note that this condition is different from electroneutrality, which requires timescales longer than max ω −1 p , (ω 2 p τ s ) −1 , which are significantly longer than τ s when ω p τ s 1. It should be pointed out that this analysis was done for a single species, while plasmas in protoplanetary discs can be made of many different species. Hence, the condition for electroneutrality need to be satisfied only by the most mobile specie of the plasma, which can then compensate for charge fluctuations, and not necessarily by all of the species present.
In protoplanetary discs, we get the following values for the plasma frequency, depending on the type of charge carrier ω p,e = 2.2 × 10 5 ξ 10 −13 where e, i, g stands for electrons, ions and grains. As it can be seen this frequency is always short compared to the timescales of interest, but grains tend to have significantly lower frequencies due to their higher inertia. The stopping times can be estimated starting from the momentum exchange rates σv ij . Since we are interested only in weakly ionised plasmas, collisions between charged species will be extremely rare. We will 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 σv e = 8.3 × 10 −9 × max 1, where σv e is deduced from Draine et al. (1983) and σv i is obtained from Draine (2011), following Bai (2011a) †. For grains above a size of a few 10 −2 µm, collisions mainly behave as billiard balls. In other words, σv is roughly equal to the velocity of the incident neutral times the cross-section of the grain. For spherical grains, this leads to ‡ σv g = πa 2 2k B T m n = 2.6 × 10 −3 a µm T 100 K 1/2 cm 3 s −1 .
These rates allow us to compute stopping times for each species following the definition above. τ s,e = 6.7 × 10 −7 R 9/4 AU s, τ s,i = 4.9 × 10 −5 R 9/4 AU s, which shows that because of the low ionisation fraction and the neutral drag, ω p,j τ s,j < 1 for ions and electrons, while it is > 1 for grains. This means that plasma oscillations are over-damped for ions and electrons (case ii above) and are not directly relevant for quasineutrality. Nevertheless, ω p τ s > 10 −2 , so even in this case, electroneutrality is recovered † The momentum exchange rate σv i estimated by Bai (2011a) is actually the collision rate. The momentum exchange rate quoted here is larger by a factor approximately 1.21 than the collision rate quoted by Bai (2011a) (see Draine (2011), eq. 2.39) ‡ Note that the expression provided by Bai (2011a) for this rate is incorrect by more than 3 orders of magnitude 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µ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 will not be explored here, and we will only consider situations where ions and electrons are still present in the system. ). However, it is numerically expensive since the numerical time steps are usually limited by τ s which is much smaller than the timescales of interest (see above). 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 τ 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 UVs 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 et al. 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 center of mass of the fluid, defining the total mass density ρ = j n j m j , the flow velocity v = j n j m j v j /ρ and the drift speed for each species w j = v j − v we sum equations (3.1) and (3.2) to obtain where we have introduced the total pressure and force P and f as well as the total current J = j n j q j v j . These equations are exact. However, they don't correspond to the usual dynamical equations one is used to, and it's important to understand why each extra term can be neglected. The first term on the righthand side corresponds to transport of momentum by the drift velocity. Physically, it can be interpreted as a diffusion of momentum due to drifting particles. It can be neglected, provided that drift velocities are small, i.e. that w j < LΩ ρ/ρ j where L is the typical length scale of interest and Ω is the typical frequency †. The presence of the density ratio ensures that even for drift velocities comparable to LΩ, this term is negligible.
We also have a term involving the total charge of the flow j n j q j . As shown above, 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

Ohm's law
In the equation of motion for the center of mass, we have left aside the fact that additional equations were required to obtain B and J . Indeed, Maxwell's equations give us The remaining unknown is therefore the electric field. Because of 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.1), and let us separate the velocity into a velocity for the center of mass, and the drift velocity for species j.
where we have defined the electric field in the center of mass frame E cm ≡ E + v × B/c and the forces on the center of mass F cm = −∇P + f + J × 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.
• d t w j can be neglected provided that Ω τ −1 s (i.e. this assumption is identical to the quasi-neutrality assumption discussed above). In other words, the inertia of charged particles is negligible and they instantaneously reach their asymptotic velocity.
• Similarly, the inertial term (second line) and external forces f j can be neglected since they modify the impulsion on timescales long compared to τ s .
• ∇P j ∼ ρ j c 2 s,j /Λ is negligible provided that c s,j ΩΛ. The equations of motion for charged particles in the frame of the center of mass therefore read where we have assumed that dominant collisions were due to neutrals. This is usually recast as whereb is a unit vector parallel to B and is the Hall parameter (Wardle & Ng 1999). Equation (3.3) can be solved for w j , which gives the asymptotic velocity We eventually obtain an expression closing our set of equations by relating the drift velocities to the current in the flow J = j n j q j w j and assuming quasi-neutrality j n j q j = 0: These expression constitute the base of Ohm's law. We can identify 3 conductivity tensors, the Ohmic, Hall and Petersen conductivity tensors: defined so that Ohm's law can be written in the more familiar form: This relation can be inverted one final time to obtain the electric field in the observer frame and write the induction equation as where the magnetic diffusivities are defined as where the subscript O,H,A stands for Ohmic, Hall and ambipolar.

Simplified case of two charged species
In the simplest case of an plasma made of 2 singly charged species (+) and (-), we get the following simplified expressions from (3.4), (3.5) and (3.6) First, all of these coefficients are proportional to n −1 + , i.e. are inversely proportional to the ionisation fraction. Secondly, since µ j ∝ B, we find that η O does not depend on B, while η H ∝ B and η A ∝ B 2 . Finally, we find that η H may have either sign. If |µ − | > |µ + |, we find η H > 0, and η H < 0 otherwise. Since µ 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 γ so that µ − = −µ + , the Hall effect vanishes.
In the case of an electron-ion plasma, we have |µ e | m i /m e |µ i | |µ i |. Hence η O ∝ |µ e | −1 and η A ∝ |µ i |, which justifies the usual statement that Ohmic diffusion is due to electron-neutral collisions and Ambipolar diffusion to ion-neutral collisions. We also have η H = |µ e |η O and η A = |µ e µ i |η O . Hence, we can distinguish three regimes depending on the Hall parameter of the ion and electrons: This allows us to delimit the Ohmic, Hall and Ambipolar regime as a function of the neutral density and the field intensity ( Fig. 9). As it can be seen, the midplane of protoplanetary discs is expected to lie mostly in the Hall regime and possibly in the ambipolar regime in the outer most parts of the disc.
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 since J × B · B = 0). Being in the Ohmic or ambipolar dominated regime does not automatically imply that the Hall effect is dynamically unimportant. . Non-ideal regimes as a function of the neutral density and magnetic field intensity, computed for an electron-ion plasma and assuming T < 100 K. The blue and green lines correspond to the typical values of a protoplanetary disc midplane, for various plasma β parameters.
Finally, we obtain the usual expressions for the diffusivities in the electron-ion case: These values can be compared to diffusivities of everyday material such as iron (η = 8 × 10 2 cm 2 /s), demineralised water (η = 1.4 × 10 15 cm 2 /s) and dry air (η = 1.6 × 10 24 cm 2 /s). Even though one might wrongfully conclude from this that MHD effects are irrelevant, the time scales (∼ year) and length scales (∼ 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 will see, one obtains magnetic Reynolds numbers of O(1), which put these flows in a regime comparable to liquid sodium experiments on Earth.

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 where V A is the Alfvén speed and Ω is the rotation rate of the system. Note however that Λ O ∝ B 2 and Λ H ∝ B, which makes these numbers less useful when it comes to predicting the saturation of MHD instabilities since B is a priori unknown. It is therefore useful to define two additional dimensionless numbers, the magnetic Reynolds number and the Hall Lundquist number 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 Fig. 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" . In addition, we find 10 −1 < L H < 10 in most of the disc midplane with a sharp increase at the disc surface while Λ A 1 in most of the disc.

The role of grains
When it comes to the conductivity of protoplanetary discs, grains play essentially two roles.
(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).
(ii) Due 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 Fig. 11 an example of diffusivity computation in plasmas of different compositions (these compositions are identical to the ones discussed in §2.5.3). As it can be seen, the addition of 0.1 µm size grains to the system has a dramatic impact on the diffusivities: all of the dimensionless numbers decrease by several orders of magnitude close to the midplane, while ambipolar diffusion becomes stronger than Hall in the case with grains. Overall, ambipolar diffusion increases by 10 4 while Hall and Ohmic increase by 10 2 compared to 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 (electrons+grains-). In the end, the Hall conductivity becomes dominated by ions when they become sufficiently abundant: at the disc surface.
In the grain-free case, the presence of metal atoms tends to decrease the diffusivities by typically 1-2 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 impact 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 a very low grain size 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).

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 AU) of protoplanetary discs. 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 cosmic rays can be both shielded by the stellar wind, or"locally" produced in shocks surrounding the forming star, there is an tremendous uncertainty of 6 orders of magnitude on the ionisation rate due to cosmic rays (see the discussion in §2.5.1). Since this mechanism is the main source of ionisation in the disc below two scale heights, this implies a 3 orders of magnitude uncertainty in the ionisation fraction ξ and similarly on the diffusivity coefficients. Ionisation due to X-rays is also subject to caution since 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 which determines the penetration depth of ionising radiations but also the recombination rate. Denser disc models, such as the minimum mass solar nebula, tend to have lower ionisation fractions and larger diffusion. We have shown that by comparing a theoretical MMSN disc model to a model favoured by observations ( §2.1), one can change the ionisation fraction by 2 G. R. J. Lesur orders of magnitude (see §2.5.3). Since the gas column density profile is largely unknown for R ∼ 10 AU, one is forced to use the gas column density as a free parameter.
Grains: as shown above, grains affect both the ionisation fraction and the dependence of the diffusivities on ξ. Overall, grains tend to reduce the ionisation fraction by several orders of magnitude (typically 2-3). Because of the change in composition (grains become the dominant charge carrier close to the midplane), diffusivities increase by 2-4 orders of magnitude in our worst-case scenario, compared to the metal-free case. In addition, the Hall diffusivity η 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 µm do not impact 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 impact 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 which are within an order of magnitude of the diffusivities discussed in our grain-free metal-free scenarios. Last, the physics at 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), while others include dependencies 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 impact on the resulting diffusivities.
Overall, one is forced to conclude that the conductivity tensor of protoplanetary discs is plagued by uncertainties and that no chemical/ionisation/grain model is better than the other. Given the discussion above, the uncertainty on the diffusion coefficient is at least ±3 orders of magnitudes, 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.
Because of these uncertainties, we will focus in the following on the "intermediate" case of a metal-free grain-free case which we discussed in §3.4.2 and for which diffusivities are given in Fig. 10.

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 by unforeseen difficulties. The main idea which 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 largescale quantities such as the disc surface density Σ and thickness H, the rotation rate Ω, the diffusivities, the magnetic field strength, etc. This approach is, in essence very similar to the mixing length theory of convection, except that in the disc case, one doesn't transport heat but angular momentum. In discs, it is called the α-disc theory.
Here, we present the basic concepts behind accretion and the α-disc theory. Then, we introduce the magnetorotational 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 protoplanetary discs, taking into account non-ideal MHD effects.

Disc equilibrium
A protoplanetary disc 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 where (R, z) are cylindrical coordinates and Ω is the angular velocity of the flow, which we assume only depends on R and ψ = −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: 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 where k is Boltzmann's constant and µ is the mean molecular mass. Since 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 R) since we assume the disc is thin: We deduce from this the vertical density profile (4.1) The thin disc approximation H R implies that the disc is cold, or in other words that c s RΩ K . In the radial direction, we first have to compare the radial pressure gradient to the gravitational potential 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 Ω = Ω 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.

Accretion theory
The energetics of MHD-driven discs has been extensively discussed by Balbus et al. (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: where E NI are electromotive forces due 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 to viscosity are always negligible for practical applications.
In order to capture the dynamics of the disc, we separate the gravitational potential ψ as a midplane potential Ψ and a deviation as one moves away from the disc midplane Φ: We also separate the mean rotational motion of the disc from its deviations (not necessarily small): where we only assume that Ω satisfies the radial equilibrium in the disc midplane Under these assumptions, it is possible to rewrite the angular momentum conservation as where we have used the continuity equation to eliminate the terms proportional to ΩR 2 . A similar procedure can be followed for the energy equation, which can be written as We recognise the mass conservation equation in the first line, and the angular momentum conservation equation in the second line. Substituting (4.2) and (4.6) in the equation above allows us to recast energy conservation as where we recognise the radial equilibrium in the first term, which can be cancelled out.
Hence, we get an energy equation for the velocity fluctuations which reads In order to compute averaged conservation equations, we define an azimuthal average as and a vertical integration of the azimuthal average so that the continuity equation (4.2) reads where Σ ≡ ρ is the gas surface density. The equation of angular momentum conservation (4.6) can be recast using the same averaging procedure (4.8) defined above to get an equation relating the mass accretion rate ρv r as to the radial and surface stresses Radial stress where we have assumed that v Ω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.7): where we have the mechanical energy of the fluctuations E m = 1 2 ρv 2 + ρΦ + B 2 8π and its associated energy flux 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 which 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 Since 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.

α 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 as a function of vertically averaged quantities such as Σ, P , etc. Historically, and based on a purely dimensional argument (Shakura & Sunyaev 1973), it is usually assumed that where α 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 Ω K (these two quantities are the only length and frequency accessible to an ideal system). Hence, we expect v = θHΩ K where θ is a dimensionless constant, of order unity. Therefore W rφ = θ 2 ρv 2 = θ 2 ρH 2 Ω 2 K . Using (4.1), one gets W rφ = θ 2 ρc 2 s = θ 2 P . Hence, thanks to the vertical equilibrium of a thin disc, Shakura & Sunyaev (1973) prescription shows up as a mixing length theory with a length H, a frequency Ω K and α = θ 2 .
Interestingly, since H = c s /Ω K , θ is actually a measure of the Mach number of the flow θ = 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 θ 1 and therefore α < 1.
This prescription may be seen as a viscous theory. Indeed, the α-disc prescription leads to W rφ = αP = αΣc s HΩ K . Since RdΩ K /dR = −3/2Ω K , the stress can be recast as where we have defined an effective viscosity ν t = αc s H. Here, we clearly recognise the usual R − φ component of the viscous stress in the Navier-Stokes equations.
Plugging the α prescription in (4.10) and neglecting surface (wind) contribution leads This allows us to express the mass accretion rateṀ ≡ −2πRρv r aṡ We can then use mass conservation (4.9) to get an equation for Σ which essentially constitutes a diffusion equation for the surface density. The diffusion timescale associated to accretion can be estimated using c s = Ω K H. One finds 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 inexistent on this timescale.

α-υ 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, due to a hypothetical wind. To do so, let us define and in a way similar to the α prescription, we assume where P mid is the midplane pressure of the disc. Using the same procedure as for the α disc, we can express the mass accretion rate as a function of α and υ The comparison between the α term and the υ term is revealing as it compares the role played by the radial and vertical stresses. One can assume that in first approximation P P mid H so that the vertical contribution is R/H(υ/α) times larger than the radial one. This implies, in particular in thin discs where R/H 1, that magnetised winds can easily be the dominant source of accretion.
In addition, using (4.14) in the continuity equation, the vertical stress term shows up as a first order radial derivative of Σ (=advection) while the radial term appears as a second order derivative as in the usual alpha disc theory. For this reason, winddriven discs cannot be treated as viscous discs, since the wind component appears as an advective term in the surface density evolution. The α 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 directly the turbulent stress W rφ from the equations of motion for the gas. This is however a rather complicated task which often implies using numerical tools, as the equations of motion cannot, in general, be solved analytically. Since the disc is thin, and turbulence, in the α disc theory, is supposed to be confined by the scale height H 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 et al. (1995).

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's been used more recently as an efficient tool to model 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 Ω O ≡ Ω 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 (Fig. 12).
In this frame, the system follows the usual single-fluid equations of motions (MHD). Since it is rotating, we have in addition a Coriolis force and a centrifugal force, so that the equations of motion read where ψ is the gravitational potential, E NI is the non-ideal electromotive force and R ≡ (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 γ. As is well known, the centrifugal force derives from a potential of the form ψ c = −Ω 2 0 R 2 /2. The effective potential (gravitational plus centrifugal) in the corotating frame therefore reads The 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 ∼ y ∼ z H to get the Hill's potential This effective potential has been truncated at the first non-trivial order. It is however interesting to note that it does not depends 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 dependencies, making calculations much more tedious. Let us emphasise already at this stage that the Hill's potential is not adapted to global phenomenon. This can be seen by comparing the iso-potentials of ψ eff and ψ eff,Hill (Fig. 13). The Hill's approximation is found to be symmetrical in x → −x, implying that one doesn't know where the center of attraction is located (both x → −∞ and x → +∞ are technically valid). Moreover, the neutral isopotential ψ(x, z) = −3Ω 2 0 R 2 0 /2 has an asymptote for z → +∞ at x = R 0 ( √ 3 − 1) which is absent in the Hill's approximation. This asymptote is key for outflows to be ejected to z → ∞, and is the main reason why local models always produce outflows which depend on the location of the z boundary conditions (see §9.2).
Finally, it can be shown that in the more general case of a central gravity of the form ψ grav = r −2(q−1) , the equilibrium rotation profile is Ω ∝ R −q and the Hill's potential reads ψ eff,Hill = Ω 2 0 R 2 0 (qx 2 + z 2 )/2 + constant. In the case of a gravitational potential from a central point mass, we simply have q = 3/2.
The equations of motion in the Hill's approximation finally reads where we have differentiated the Hills potential to obtain the tidal acceleration Ω 2 0 (qxe x − ze z ). This set of equation admits a simple steady solution V 0 for the velocity field by simply balancing the Coriolis and tidal forces in the x direction: 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. Since the Hill's approximation only retains the first terms of the effective potential, the shear in the 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 w ≡ v − V 0 (where w is not necessarily small compared to v), for which the equations of motion read where we have defined the comoving derivative D t ≡ ∂ t − qΩx∂ y . It is worth noting that this last set of equation 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 will present later. In exchange, new source terms have appeared when going from v to w: the "lift up" effect qΩ 0 w x e y , which is essentially the advection of the mean keplerian flow by radial motions, and the "Ω effect" of dynamo theory (actually due to shear, not rotation) −qΩ 0 B x e y . As we will 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.3). By defining L = w y + (2 − q)Ω 0 x, we get the conservation equation where the Coriolis and tidal terms have vanished in the definition of the local angular momentum.

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 the 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 the 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 H. In this case, vertical gravity can be neglected in (5.3) and since one expects w ∼ ΩL, the flow is strongly subsonic so that an incompressible approximation † can be made in place of (5.2). The asymptotic of these two models are discussed in details by Umurhan & Regev (2004). Here we will mostly use the LSB, keeping in mind that an incompressible approximation is possible to study the basic effects of the shear.

Boundary conditions
The shearing box model makes use of the Hill's approximation in the simplest possible numerical setup: a periodic box. Because the 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 ∼ H R 0 . However, the presence of the comoving derivative makes things a bit more difficult. The fields ρ, w, etc. are advected everywhere at the azimuthal velocity −qΩxe 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ΩL x t, z), which can be represented graphically as in Fig. 14.
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.

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.10). One therefore needs to quantify which is the Hill's equivalent of the radial stress term in (4.10). Despite a non-zero radial stress, the shearing box model does not exhibit any accretion due 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 the Hill's potential.

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 · ≡ dx dy dz (5.7) Mass, momentum and flux conservation equation eventually read 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 while 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 Ω-effect which shows up as a source term of toroidal field. These simple conservation laws allows one to control very carefully the box physics. Moreover, since 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 w and the induction equation with B. One eventually gets an equation for the mechanical energy (kinetic plus magnetic) in the box This set of equations if the local equivalent of (4.7). 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: • P dV work. This term is usually small when thermal effects are negligible, but can become important in thermally-driven winds for example.
• Nonideal effects. For Ohmic and ambipolar diffusion, E NI ∼ −J , hence for these two terms, E Ohm/Ambipolar · J < 0. Meanwhile, the Hall effect has no contribution to this term since J × B · J = 0. Overall, the non-ideal term therefore appears as a sink of mechanical energy, as expected. Energy is dissipated into heat.
• Radial stress. We recover the radial stress term of the global angular momentum equation (4.10) which appears as a source term here.
• 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.
• Viscous terms (not shown). When the flow is turbulent, viscous terms can become important at 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 F Mech is periodic/shear periodic. If we average the energy equation, we find that the only relevant flux component is the one 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 the 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 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 Ω 0 . In the horizontal direction, the equations describes epicycles. To show it, let us first integrate (6.2): where L is the local angular momentum of the particle, as defined in (5.6). Our particle being initially in equilibrium at (x, y) = 0, it has L = 0 and we can write the radial equation of motion as hence, our effective gravitational potential which was initially unstable (∂ x ψ eff,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 This characteristic frequency is named epicyclic frequency. In the particular case of a Keplerian disc (q = 3/2) we find κ 2 = Ω 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. Fig. 15).
This shows that at the linear level, pure Keplerian flows are stable. This doesn't mean that nonlinear (or subcritical) instabilities cannot exist in these flows, given that the shear is a natural reservoir of free energy to trigger instabilities. Such nonlinear 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 number up to a few millions. However, thermal effects, such as heating, vertical stratification, etc. are known to affect this picture, leading to thermally-driven linear and nonlinear instabilities. This is a whole field of research, which will not be covered here. Interested reader may consult Fromang & Lesur (2019) for a more detailed overview of this topic.

Linear Magnetohydrodynamic stability
As we have seen, protoplanetary discs are hydrodynamicaly stable at the linear level. In MHD, however, things start to be 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 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, 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 the resulting magnetic perturbation can be obtained by integrating the induction equation with respect to time: In order to model how the field impacts the dynamics, we have to include the Lorentz force F L in the equation of motion. In the horizontal direction, only the magnetic tension term B · ∇B appears, so we have The horizontal equation of motion are therefore reduced to 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 which leads to an instability. To show this, let us assume x = x exp(σt). The equations of motion lead to the following eigenvalue problem A k 2 )y = −2Ω 0 σx, which allows us to get the dispersion relation where we recover epicyclic oscillations when V A = 0 with σ 2 = −2Ω 2 0 (2 − q) = −κ 2 and pure Alfvénic oscillations when Ω 0 = 0 with σ 2 = −V 2 A k 2 . Expanding this dispersion relation leads to This dispersion relation describes a linear instability when σ 2 is real, i.e. when This instability is the magneto-rotational instability (or MRI). It appears when the magnetic tension force is not too strong, as suggested by (6.6). It is possible to solve the full dispersion analytically to get the eigenvalues (see Fig. 16). When V A k < √ 2qΩ 0 , positive eigenvalues are found which are the signature of the MRI. The maximum growth rates are obtained for V A k = √ 2qΩ 0 /2 with σ max = qΩ/2. Above the limit (6.6), 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 (Fig. 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). The physical interpretation of the MRI is straightforward: consider 2 fluid particles attached to a vertical field line and assume we slightly move these particles radially. At first, they will start an epicyclic motion and drift azimuthally (Fig.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 reversely for the outer particle. This mechanism can only work if the radial magnetic tension is sufficiently weak otherwise, the particles come back to their initial point resulting in an Alfvénic oscillation. It is this radial component of the Lorentz force which is the stabilising agent of the MRI.

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,min 1/H which implies that any disc model threaded by a vertical field has a minimal Alfvén frequency ω A ≡ k z V Az . 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,min and ω A is also quantised. As an example, we show in Fig. 18 the modes accessible to a disc threaded by a vertical field with V Az = 0.2 ΩH. In this particular example, only the four smallest n are MRI-unstable, the most unstable mode being n = 2.
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 Fig. 18 since k z,min increases. The MRI is entirely stabilised when the lowest n gets into the stable regime, i.e. when V Az > √ 2qΩ/k z,min . With k z,min 1/H this implies in Keplerian discs 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 it is sometimes thought. Indeed, since ΩH = c s , the criterion above 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).

Non-axisymmetric MRI
We start from equations (5.2)-(5.5) in which we assume the disc is threaded by a mean field having a mean vertical and azimuthal components: B = B 0,y e y + B 0,z e z . In this subsection, we relax the axisymmetry hypothesis used above, but we still assume the flow is incompressible, and neglect vertical gravity and stratification. We then consider the following set of equations: We are going to linearize this system, assuming the field can be decomposed as B = B 0 + b. The presence of the advection term qΩx∂ 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 quantityX assumes as spatial decompositioñ Using this decomposition, it is easy to show that Since it is the only term which shows this explicit x dependency in the equations of motion, and since these equations are assumed to be valid for all x, we are forced to conclude that the term in parenthesis cancels out: Without loss of generality, we can therefore assume a decomposition into "shearing waves", defined as k(t) = k 0 + qΩtk y e x † Although these solutions are often called "waves", these are not normal modes of the physical system, but they are just a convenient way to decompose solutions in a sheared flow 46 G. R. J. Lesur which is solution to the set of equations above. We using this shearing wave decomposition we obtainẇ where we have dropped the explicit time dependency of k and the˜symbols for simplicity. We next take the time derivative of (6.12): This allows us to express the generalised pressure So that in the end, the equations of motion read: − qΩw x g yy e y − (2 − q)Ωw x (1 − g yy )e y − 2Ωw y g xy e y + 2(1 − q)Ωw x g yz e z − 2Ωw y g xz e z where we have introduced g ij = k i k j /k 2 . It is not possible to go any further without making any approximation. Indeed, although by construction k · 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 which only grow for a finite time. To understand why this is always the case, let us continue our analysis using a first order WKB approximation.
To get a dispersion relation, one needs to assume that k is "almost" steady, i.e. k y k x so that d log k/dt Ω. 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 above expansion. Assuming X =X exp[σt + ik(t) · x)], we then get These equations clearly exhibit an Alfvén mode in the z direction, with σ = ±i(k · V A,0 ) and V A,0 ≡ B 0 / √ 4πρ as one would expect. We can then solve independently the horizontal problem to obtain the dispersion relation: where we have defined the Alfven frequency ω A = k · V A,0 . This equation describes both the traditional MRI mode and the non axisymmetric MRI, as its close ressemblance with (6.5) suggests. It should however be noted that in the k z /k → 0 limit, the instability is lost since g zz → 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 k is slowly evolving with time. As t → ∞, one then expects |k x | → ∞ and therefore g zz → 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, since one needs k · B 0 = 0. Therefore, there is no "toroidal field MRI" as it is sometimes seen 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 shown above.

Historical background
The linear MRI in the non-ideal MHD regime has been explored by many authors. After the discovery of the MRI in the disc context by Balbus & Hawley (1991), it was soon realised that protoplanetary discs, 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 protoplanetary discs . Later, Wardle (1999) considered the 3 non-ideal MHD effects in simple axial geometry while 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).

Linearised equations
As in the cases above, we start from equations (5.2)-(5.5) in which we assume the disc is threaded by a mean field having a mean vertical and azimuthal components: B = B 0,y e y + B 0,z e z . In the following, we consider small axisymmetric † perturbations of the equilibrium. We seek solutions of the form w = u exp[k · x + σt] and B = B 0 + b exp[k · x + σt]. Where k = k x e x + k z e z is the wave number, and σ 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 small shearing box approximation. As we will show later, this is enough to capture most of † in the Hill's approximation, non-axisymmetric linear perturbations only lead to transiently growing solutions which ultimately decay (Balbus & Hawley 1992). These transient "modes" are however important for the MRI dynamo once nonlinear feedback is taken into account (Lesur & Ogilvie 2008a).
the physics relevant to the problem. We will explore in §6.4.6 the impact of the vertical stratification on linear modes.
Under these assumptions, the linearised equations read 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 k × b, which means it rotates the magnetic perturbation around the k direction, keeping its norm constant. Note that the direction of rotation is given by the sign of η H , which shows that the handedness given by the Hall effect is directly connected to the microphysics of the plasma (see §3.4.1). Finally, ambipolar diffusion involves an anisotropic diffusion term which we will 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 4 th order problem The role played by ambipolar diffusion is here a bit more self-explanatory. We observe that the diagonal terms in the 2 nd 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 B 0 . However, there are also off-diagonal terms proportional to k x B 0,y . As we will see, these terms can lead to oblique unstable modes (see also Kunz & Balbus 2004). We next solve the above set of equations for σ and look for unstable eigenvalues. We follow Pandey & Wardle (2012) and first obtain an equation for the velocity fluctuations So that the induction equation becomes in matrix form: where we have introduced the Alvén speed V A = B 0 /(4πρ) 1/2 , the Alfvén frequency ω A ≡ k · V A , the Hall length H ≡ η H /V A , the ambipolar time τ A ≡ η A /V 2 A and the epicyclic frequency κ 2 ≡ 2Ω 2 (2 − q).
After a long but straightforward calculation, one eventually gets the dispersion relation which can be written The stability of this linear system can be analysed in the vicinity of σ = 0. In this case, a necessary and sufficient condition for instability is 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 = 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 don't 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 will see later.

Ohmic diffusion
The impact of Ohmic diffusion on the stability of the MRI is physically very intuitive: it stabilises MRI modes, starting from the largest ω A of the system (see Fig. 19). The stability condition in the presence of Ohmic diffusion is deduced from the condition C 0 = 0 and reads . MRI growth rate as a function of the Ohmic Elsasser number ΛO for a Keplerian disc (q = 3/2) with kx = 0. Note the damping of the most unstable mode for ΛO = 1/ √ 3 and the survival of low growth rate modes in the limit ωA → 0, ΛO → 0.
Clearly, the modes with k x = 0 are the last ones to be stabilised when one increases Ohmic diffusion. Let us focus on this case and assume the flow is Keplerian so that κ = Ω 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 ω A = √ 3Ω/2 so that the above inequality reduces to Λ 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 Ω → 0 and neglecting Ohmic and ambipolar diffusion. In this case the dispersion relation (6.19) gives We recognise two waves with frequency ω ≡ iσ given by where "+" waves are known as whistlers or electron-cyclotron modes, while "-" waves are ion-cyclotron modes. The whistler frequency ω H = ω A H k increases as k 2 in the limit k → ∞ while the ion-cyclotron frequency tends to a constant ω IC = ω A /( H k). Hence these two waves behave very differently at small scale. For positive H , whistlers are right-handed polarised wave while 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 → 0.
Hall-shear instability: The Hall-shear instability † (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 ω A → 0 while keeping H ω A > 0. This "low magnetisation" limit allows one to decouple the ions from the electrons, as is evident from (6.17). Neglecting Ohmic and ambipolar diffusion, one gets the following dispersion relation which exhibit a linear instability when Interestingly, the HSI shows up only when q H V Az > 0, or in other words, when the vertical field points in the same direction as the rotation axis in Keplerian discs ‡, assuming H > 0. When the whistler frequency becomes too large (k 2 H V Az > qΩ), 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 way similar to the MRI. Last, the maximum growth rate is identical to the MRI σ max = qΩ/2 and is obtained for k 2 H V Az = qΩ/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 (Fig. 20). If the vertical field is positive (i.e. aligned with the rotation axis), the Keplerian shear is going to shear the perturbation in the opposite direction, amplifying the toroidal field from the radial field and feeding the instability. On the contrary, 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. Ion-cyclotron instability: The ion-cyclotron instability ¶ is more difficult to isolate compared to the HSI since ion inertia cannot be neglected in this case. However, it is still possible to filter out whistler waves by letting k z → ∞ and keeping constant k x k z . In this case, the whistler wave frequency ω H becomes infinite while the cyclotron frequency stays finite. Since we're looking for a finite growth rate, we will assume σ is finite when k → ∞. Neglecting ambipolar and Ohmic diffusion and keeping only O(k 4 ) terms, the dispersion relation (6.19) becomes (see also Simon et al. 2015): This relation clearly describes ion-cyclotron modes since in the limit Ω → 0, we recover the ion-cyclotron frequency σ = ±iω CI . It describes unstable modes, resulting from the interaction between ion-cyclotron waves and epicyclic motions, provided that → ion-cyclotron unstable (6.24) This inequality brings up two remarks: 1-the ion cyclotron instability appears for antialigned field configurations, i.e. field configuration opposite to the HSI. 2-this instability does not vanish in the limit k z → ∞. 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 Fig. 21 the growth resulting from (6.19) in the Hall only case. To quantify the intensity of the Hall effect, we have defined a Lundquist number based on the vertical wavenumber We have chosen the most unstable modes k x = 0 and assumed k z > 0 so that ω A < 0 corresponds to B 0z anti-aligned with Ω. As can be seen in Fig. 21, we recover the MRI in the limit L H → ∞ which gives identical growth rates under the symmetry ω A → −ω A . As L * H decreases and the Hall effect increases, this symmetry is broken. For ω A > 0 (aligned field configuration), the MRI becomes the HSI, and the optimum growth rate starts to move to lower ω A . This is expected from our analysis since the Figure 21. Hall-MRI growth rate as a function of the modified Hall Lundquist number L * H = ( Hkz ) −1 for a Keplerian disc (q = 3/2) with kx = 0. We have assumed kz > 0 so that ωA < 0 corresponds to the anti-aligned case VAz < 0. The white plain line corresponds to the HSI stability limit (6.23) and the white dashed lines to the ion-cyclotron stability limits (6.24). maximum growth rate of the HSI is found for ω A /Ω L * H /2. We show in white filled contours the stability limit of the HSI obtained from (6.23). This limit matches the full dispersion relation for L * H 0.3. For ω A < 0 (anti-aligned field configuration), the MRI becomes the ion-cyclotron instability and the optimum growth rate moves to higher |ω A |. The stability contours of the ion-cyclotron instability (6.24) are shown in dashed white line. As for the HSI, they match the full dispersion relation for L * H 0.3. Impact of Ohmic diffusion on the Hall-MRI: As shown above, Ohmic diffusion tends to damp MRI modes, starting from the largest ω A and moving to lower ω A as Λ O decreases. On the other hand, the Hall effect creates two distinct branches from the MRI, depending on the field alignment configuration. Since Ohmic diffusion suppresses first high |ω A | modes, the ion-cyclotron instability will be first to be stabilised. On the other hand, since the HSI is living at lower |ω A | compared to the MRI, it will be less affected by Ohmic diffusion than the MRI. This effect is illustrated in Fig. 22, demonstrating that HSI modes survive more easily to Ohmic diffusion, even when ideal MRI modes are suppressed. We can go further and estimate the stability limit of the HSI in the limit ω A → 0 while keeping H ω A > 0. In this limit, the dispersion relation of the HSI including Ohmic diffusion reads from which we can deduce a general instability criterion. For the most unstable HSI where Λ H is the Hall Elsasser number (see §3.4.2). This expression can be compared directly to the Ohmic-only criterion (6.21), and shows that as expected, the HSI can make the system unstable even when Λ O 1, provided that Λ H 1 as well. For this reason, some authors 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 ω A ∼ Ω, in other words, the Alfvenic and rotation frequencies match. In the HSI case, it is the whistler and rotation frequencies ω H ∼ Ω which have to match. In the limit of strong Hall effect, ω H = H k z ω A . Therefore when H kz 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, since diffusion first damps small scales 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 will see in §8.4, the HSI in the non-linear regime is indeed full of surprises. . .

Ambipolar diffusion
The linear ambipolar diffusion is a non-diagonal operator as shown in (6.17). The nondiagonal terms are non-zero only when k x V Ay = 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 Ay = 0.
Diagonal case: (k x V Ay = 0) In this case, the stability criterion is very similar to that of Ohmic diffusion. By requiring C 0 > 0, one gets 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 Ay 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 → Stability of the most unstable ideal MRI mode, (6.26) and a criterion for the stability of all the modes available smaller than the minimum vertical wavenumber k z,min : This shows that the MRI can be stabilised even for Λ A > 1 provided that the toroidal field dominates over the poloidal one. 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. On the contrary, ambipolar diffusion increases with the total field strength. Hence, an increase in V Ay results in an increase of the effective diffusion η A , while the MRI feedback loop is left unperturbed.
Non-diagonal case: (k x V Ay = 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.19): qΩk x k z V Ay V Az < 0. An illustration of the impact of this term on the general stability property is shown in Fig. 23 where we have assumed V Ay = V Az , q = 3/2 and Λ A = 0.4. We observe that oblique modes (k x = 0) tend to be more unstable than axial modes when Λ A < 1. It can be shown that oblique modes are always unstable, albeit with a vanishing growth rate, in the limit Λ A → 0 provided that k x /k z is sufficiently large (Kunz & Balbus 2004).

Vertical stratification
All of the above was computed ignoring vertical stratification. It should however be pointed out that since discs are vertically stratified, unstratified results should be taken with care. Let us emphasize here a few important results regarding the effect of stratification on MRI modes. In this section, we will 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 . MRI growth rate with ambipolar diffusion and ΛA = 0.4. We have chosen VAy = VAz and kx = 0 to illustrate the effect of non-diagonal ambipolar terms. The most unstable mode in this case is found for kx < 0 and a weakly unstable branch exists for kx/kz 3.5 in the limit ωA → ∞. We name these modes "oblique ambipolar modes". Note however their low relative growth rate midplane: . Ideal MRI: the ideal MRI case with an isothermal disc is presented in details by Latter et al. (2010). We will therefore just summarise here the main conclusions. First, stratified eigenvalues (growth rates) satisfy the same dispersion relation (6.5) as non stratified modes with V A = V A (z = 0). Vertical wave numbers 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 n H = 1.1584, 2.0796, 2.9829, 3.8798 . . . for n = 1 . . . 4. Taking k 1 as the lowest wavenumber accessible to the system, the MRI is stable for all possible eigenmodes in Keplerian discs This value can vary slightly depending on the boundary conditions for the perturbation as z → ∞. Overall, it is safe to assume that the MRI exists provided that β mid 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 Fig. 24). The MRI does not choose specifically any symmetry in these local stratified models: both even and odd n follow the same dispersion relation. Third, the fastest growing modes tend to be more oscillatory in the disc midplane at   Note the strongly oscillating behaviour close to the midplane.
higher β mid . This is because qualitatively, one expects kV Az ∼ Ω for the fastest growing mode, hence in the midplane, one gets k z H ∼ √ β mid . An example of such a mode is shown in Fig. 25 Ohmic and ambipolar diffusion: as shown in the non-stratified section, Ohmic and  Figure 26. Fastest growing MRI eigenmode in a stratified model for β mid = 10 3 using the diffusivity profiles from Fig. 10 at R = 1 AU. Left: including Ohmic diffusion only (σ = 0.72Ω), right including Ohmic and ambipolar diffusion (σ = 0.56Ω). 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. 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 Fig. 10). We show in Fig. 26 the most unstable MRI eigenmode with β = 10 3 resulting from our fiducial metal-free mode at R = 1 AU. As expected, the MRI is strongly suppressed where Λ O < 1. Ambipolar diffusion tends to reduce the overall growth rate but does not impact significantly the shape of the eigenmode. 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 1-2 scale-height essentially magnetically dead . Note that when the linear eigenmodes do not propagate in the disc, as it 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 above, 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 Fig. 27 where the MRI eigenmode now propagates in the midplane resulting in a fully active disc column at R = 1AU when V Az Ω > 0. In this case, however, the growing perturbation is mostly a magnetic one 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.

The Helicoidal MRI (HMRI)
The Helicoidal MRI 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 protoplanetary discs. For the sake of completness, let us show here the origin of this instability and discuss its application to Keplerian discs.  Figure 27. Fastest growing MRI eigenmode in a stratified model for β mid = 10 3 using the diffusivity profiles from Fig. 10 at R = 1 AU including Ohmic, Hall and ambipolar diffusion. Left: assuming VAzΩ > 0 (σ = 0.65Ω), right: assuming VAzΩ < 0 (σ = 0.51Ω). The "dead zone" is now subject to long wavelength perturbations when VAzΩ > 0, as a result of the HSI.
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 the Hills approximation (even though a WKB analysis can capture it in the global geometry, see Kirillov et al. 2014 and below). Therefore, let us start back from the evolution equations in global geometry.

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 u and the the magnetic field B. We consider the inviscid, incompressible and ideal equations of magnetohydrodynamics, which reads, component by components: for the equation of motion, which we combine to the induction equation: and the continuity equation We linearise the equations with respect to a background flow where Ω(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 B 0 = B φ,0 (R)e φ + B z,0 e z so the only spatial dependency is in the R direction for the toroidal component of the field. In addition, we're going to assume that the flow is axisymmetric, so that we can cancel ∂ φ derivatives. The velocity and magnetic fields are then expanded as where deviations are assumed to be infinitely small compared to the means. The linearised equations of motion eventually read In order to make the notations more concise and consistant with the usual MRI derivation, we define the following coefficients 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 these notations the equations eventually reads

WKB approximation and pressure
We are going to look for a local solution, i.e. a solution with fast variation compared to R. We will therefore zoom on a tiny region around a fiducial radius R 0 defining x ≡ R − R 0 , and expand the solution as where we have defined the Alfvén frequency ω Az = k z B 0,z / √ 4πρ 0 and the toroidal Alfvén . 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 will 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 get an equation for Π. This gives This eventually leads to the following linear problem where g zz = k 2 z /k 2 . This has non trivial roots provided that the matrix determinant cancels out. This condition leads to a dispersion relation on σ: σ 4 + σ 2 2ω 2 Az + κ 2 g zz + 2(p + 1)ω 2 Aφ g zz − 8iσω Az ω Aφ Ωg zz +ω 2 Az ω 2 Az − 2qΩ 2 g zz + 2(p − 1)ω 2 Aφ g zz = 0 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 will see later, the HMRI driving term is actually the linear term in σ.

HMRI with resistivity
Adding resistivity η means that the linearized induction equations are modified as follows

Inductionless limit
In the inductionless limit, we neglect magnetic induction by considering roots σ ηk 2 . We therefore introduce a small parameter ε = 1/ηk 2 , and expand the above dispersion relation at first order in ε. This suppresses the strongly damped magnetic modes which have σ ∼ −ηk 2 and 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: Az + 2(p + 1)ω 2 Aφ g zz + κ 2 g zz − 2iεω Az ω Aφ Ωg zz (4 − q) = 0 The roots are simple to obtain and give at first order in ε: As it 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 Since this condition is on a second order polynomial on ω Aφ , it is strictly equivalent to requiring that the discriminant of the polynomial is positive, i.e that Assuming κ 2 > 0, this is again equivalent to asking that which is (again!) a second order polynomial in q, which we ought to resolve. The stability condition is then simply In the current-free configuration (p = 1), we recover the so called Liu limit (Liu et al. 2006) q > 4( √ 2 − 1) 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 = 1), it is possible to recover the instability for Keplerian rotation profiles. This can be deduced from (7.2): 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.1) which is simply The maximum growth rate is then given by 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 (σ max )/Ω Λ O . This implies relatively low growth rates in protoplanetary discs, unless they are hosting a dynamically strong vertical field (with V Az ∼ ΩH where H is the disc thickness). In addition, these "optimum" growth rates can only be reached for ω Aφ ∼ ω Az (from Eq. 7.3). Combining these constrains, we find that V Aφ ∼ RΩ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 not Keplerian anymore 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 protoplanetary discs 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 ω Aφ Ω.

Nonlinear saturation of the MRI
The saturation of the MRI has been mostly studied in local shearing box simulations. Since 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 small shearing box 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 will first focus on unstratified models before moving to stratified ones.
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's 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 due to the technical difficulties associated to simulations with a mean field.
In all of these models, one uses box averages, defined for a quantity Q by One of the key element is then to quantify the α parameter from the turbulent stress, which is defined as which is also sometimes averaged in time.

Unstratified models
Except in rare circumstances, unstratified models are periodic in the vertical direction. In this case, the energy conservation equation (5.8) reads In numerical experiments, since MRI turbulence is subsonic, the contribution of the P dV term is small while the x − y stress term is definite positive (it is the term driving accretion!). Therefore, a quasi-steady turbulent state for which E mech ∼ constant necessarily implies that non-ideal effects (or viscous effect, 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 will keep the terminology "ideal MRI" in cases where dissipation is sufficiently small to be negligible on large scales Λ O,H,A 1, or when no physical dissipation is explicitly considered 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 see. . . Configurations with a mean field are often characterised by the plasma β parameter of the mean field, as defined by This should be confused with the turbulent plasma β parameter which is not a control parameter of the physical system since usually B B . In principle, the vertical box extension L z does not necessarily match the pressure scale height H = c s /Ω since 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 get an alternative expression for β mean which can be useful to interpret numerical simulations It is this expression for β which is used in incompressible simulations.

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, breakup presumably thanks to secondary instabilities which are essentially Kelvin-Helmholtz type modes (Goodman & Xu 1994;Latter et al. 2009) and end up in fully-developed MHD turbulence.
In the mean vertical field case, Hawley et al. (1995) found that and α 6.7β −1/2 mean . 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 α and β are overestimated by a factor of 2 compared to wider boxes (Bodo et al. 2008). We therefore correct Hawley et al. (1995) for this and get the following scaling law in the range 400 < β mean < 5 × 10 4 : a scaling which has also been verified in incompressible simulations (Longaretti & Lesur 2010). The extrapolation down to β → 1 suggest that α 1 could be reached. However, it is very difficult to explore this limit numerically since the MRI becomes in this case 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 α 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 < β mean < 1200 gives lower α values: 8.1.2. Zero mean field: the "MRI dynamo" The existence of an MRI-driven dynamo was first demonstrated by Hawley et al. (1996). In this configuration, no external field is imposed and turbulence regenerates a field on which the MRI can grow, feeding back turbulent motions. Because of 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 α ∼ 0.01. However, it was later realised that the value of α in ideal MHD simulations (where no physical dissipation is introduced) depends on the numerical resolution , with α ∝ 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 α from zero net flux simulation is intrinsically flawed. Obviously, this is most probably due to numerical dissipation (the only dissipation channel of these simulations), which is not acting as a real physical dissipation operator. If one introduces explicit viscosity ν and resistivity η O in the system, α then converges to a finite value, which seems to depend only on Pm ≡ ν/η O (Fromang 2010) and not anymore on the resolution. This, however, leads to another complication known as the "Pm effect".

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 ν and resistivity η 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 (Lesur & Longaretti 2007;. In simulations with a mean field (Lesur & Longaretti 2007;Simon & Hawley 2009), the Pm effect results in an increase of α when Pm increases. In the zero mean field case, the effect is even stronger since the MRI dynamo, and therefore MHD turbulence, simply disappears when Pm 1 Walker & Boldyrev 2017), giving α = 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's been proposed that non-local energy transfers in spectral space ) could be responsible for this effect. If true, since non-locality is necessarily bounded (Aluie & Eyink 2010), the Pm effect should disappear as (η O , ν) → 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. 2015). This, however, is not fully satisfactory as it does not prove that the MRI dynamo is inexistent in the limit Pm → 0 and Rm → ∞.
Despite numerous efforts, recent results indicate that the Pm effect is still very much alive in simulations, up to Rm = 5 × 10 4 and Λ 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 protoplanetary discs where Rm 1 is expected. For the regions above 1 AU, large-scale physics is dominated by magnetic diffusion and the P m effect is not relevant.

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.22). Using the expression for β mean , and assuming that k z,min = 2π/L z , the stability condition of a non-stratified shearing box reads The impact of a strong Ohmic diffusion on the non-linear saturation of the MRI was first explored by Sano et al. (1998) in 2D andFleming et al. (2000) in 3D. In the case with a mean field, it is found that when Λ O 1 (or equivalently Rm β mean ) and (8.2) is verified so that the system is MRI-unstable, MRI turbulence is affected: α becomes lower than in the ideal case, and turbulence becomes intermittent, with periods of linear growth followed by rapid decay due to reconnection events. As in the low-β ideal case, it is unclear whether this regime can be modelled with a constant α 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 ( Fig. 28). We have assumed that in the limit β mean → ∞, the zero net flux subcritical threshold was to be considered instead of the linear stability limit. As a result, the region β mean > 10 5 follows the zero mean field criterion, while β mean < 10 5 shows a transition region between fully developed turbulence with Λ O > 1 and a linearly stable flow, which we have named intermittent turbulence, in reference to Fleming et al. (2000). Since the data is sparse, we do not have any clear estimate of the values α in the intermittent turbulence region.

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 Λ A 100. Because protoplanetary discs are strongly collisional, the two-fluid approach is not very efficient since 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 Λ A decreases.
Starting from Bai & Stone (2011) results, it is possible to create a simple phenomenology for MRI saturation under the impact of ambipolar diffusion. First, let us recall from (6.27) that the linear stability criterion with k z,min = 2π/L z is given by The linear stability of a box with a pure vertical field is therefore very similar to the Ohmic case. However, as the instability grows, V A increases which eventually leads to the violation of the stability criterion above. This self-suppression effect of MRIturbulence allows us to deduce the saturation level by assuming that at saturation The ratio of Alfvén speeds can be evaluated by β : In principle, one would naively expect δ = 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 β computed from the fluctuations at all scales. Instead, we therefore choose δ = 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 β (2α) −1 , we get a simple estimate for α by combining the expression above: Of course, this estimates diverges as Λ A → ∞ since the saturation mechanism due to ambipolar diffusion becomes inexistent. In this case, the ideal-MHD estimate (8.1) should be used instead. A sample of predicted α from this saturation estimation is given in Fig. 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 ambipolardominated discs. This estimate also successfully recovers the ideal regime when Λ A 50, as reported by numerical simulations.
In the case of a pure azimuthal field, Bai & Stone (2011) have shown that no turbulence is sustained below Λ A 3, and α progressively drops to 0 from the ideal MHD value at Λ A = 100.
Finally, in the case of a mixed mean vertical and azimuthal field, the subsistance of oblique modes (see §6.4.5) even at Λ A 0.1 creates a weak transport with α ∼ 3 × 10 −4 for V Ay ∼ V Az (Bai & Stone 2011). Above Λ A ∼ 1, the scaling (8.3) is approximately recovered as oblique modes become unimportant.

Hall effect
The impact of the Hall effect on the saturation level of MRI turbulence was first explored by Sano & Stone (2002) with a relatively weak Hall effect (L H 20). They found that with a mean vertical field, α increases with decreasing L H in the aligned case, while α decreases in the anti-aligned case. However, protoplanetary discs are likely to have lower L H than the ones studied by Sano & Stone (2002) (typically L H ∼ 1 as in Fig. 10). For this reason, this problem was revisited by Kunz & Lesur (2013) with simulations in the L H ∼ 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 due to the HSI, the flow settles down into a quasi-laminar state for L H 5, with negligible turbulent transport (Fig. 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 In the zero net-flux case, Sano & Stone (2002) found evidence that Rm c decrease from around 10 4 in the ideal MHD regime possibly down to a few times 10 3 for L H 20. At the same time, α increases by a factor of a few compared to the ideal case. For stronger Hall effects (L H 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 (Fig. 31).
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.

Self-organisation
Self-organisation is a process by which a disorganised (i.e. turbulent) flow creates largescale and long-lived structures. There are several examples of self-organisation in nature, the most well-known one 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  10), α is larger than in the case without Hall effect by a factor ∼ 3. When Hall gets stronger, the system gets into the low transport state (LH 5). In contrary to the net field case, it periodically switches back to a high transport state, resulting in bursts of α. as Kelvin-Helmholtz or the Rossby wave instability (RWI, Lovelace et al. 1999) which 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 which have been proposed to explain this phenomenon. Hawley (2001) and Steinacker & Papaloizou (2002) were among the first to notice 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 while 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, α can reach values as high as 1, consistently with the fact that these gaps correspond to low β 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 β mean (assuming the mean field is kept at its initial value). Since α ∝ β −1/2 mean (see 8.1), α 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.

Ideal MHD
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 since several authors have measured the turbulent resistivity tensor of MRI turbulence (Guan  fig. 3). Hence, the turbulent resistivity prescription E = η turb ∇ × B likely breaks down altogether and should be replaced with a more elaborated 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 et al. (2009) also reported the formation of "pressure bumps" which are caused by large-scale fluctuations of α. 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 α which predicts long-lived axisymmetric structures in quasi-geostrophic equilibrium, as observed in their simulations.

Hall-MHD
The first mention of self-organisation in Hall-MHD appears in Kunz & Lesur (2013), where self-organisation has a dramatic impact 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 (Fig. 32). While in the ideal-MHD case, self-organisation appears as a "second order" effect on top of MRI turbulence, in the Hall-MHD case, selforganisation is the main saturation mechanism of the MRI. In other words, turbulence is mostly suppressed by self-organisation. Self-organisation shows up when L H 5, and α essentially vanishes as a result of the lack of turbulence (Fig. 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.
The origin of Hall-driven self-organisation can be tracked down to the induction equation in the presence of Hall effect. Indeed, When the Hall length is constant, the induction equation reads which highlights the role of the Maxwell stress in the induction equation. Guided by numerical simulations which show the appearance of a vertical magnetic field with variations in the x direction, let us define an average Q = dy dz Q so that the induction equation for the "mean" vertical field reads where we recognised the radial Maxwell stress term M xy = −B x B y /4π, also present in the angular momentum conservation equation (4.10). This demonstrates that in Hall-MHD, the transport of magnetic flux is tightly linked to the the transport of mass. Due to energetic constraints, M xy > 0 (see §5.2.4) in shear-driven instabilities/turbulence. Therefore, a concentration of magnetic field due 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 equation above. 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 B z < 0 and the systems settles down into a quasi-stationary state with very low-stress level (see also Fig. 33). Although this process was initially identified in unstratified shearing box simulations, it was later unambiguously identified in cylindrical unstratified simulations (Béthune et al. 2016). However, vertically stratified simulations do not seem to exhibit this process, for reasons not yet identified to date (Lesur et al. 2014).

Ambipolar diffusion
Self-organisation due 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 zmax /B zmin ) in the ideal-MHD case is larger than the ambipolar diffusion case by a factor 10 ( Fig. 2 in Bai & Stone 2014). In addition, the mechanism proposed for selforganisation 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 = 4H. However, a choice of box with L y > L x tends to break zonal flows (Fig. 34). Overall, it is very possible that self-organisation in ambipolar-dominated unstratified shearing boxes is a mere numerical artefact.

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). 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 β parameter of the mean vertical field threading the disc where averages are taken in the disc midplane. In the following, we use a definition of α using the box averaged pressure, i.e.
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 α from the vertically-averaged stress and the midplane pressure P 0 ≡ P (z = 0), which leads to predicted α values 2-3 times smaller than the ones obtained with the definition above (Davis et al. 2010). Note also that Stone's definition yields α which decreases as the box size increases. Indeed, the stress being concentrated in the region z 2H, the vertically averaged stress decreases as the box size increases above 2H. These differences in the definition of α 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 (Fig. 35, see also Davis et al. 2010 andSimon et al. 2012). This diagram shows quasiperiodic 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 (e.g. Hogg & Reynolds 2016 for an example). On a theoretical side, it is well reproduced by mean field models (Gressel 2010;Gressel & Pessah 2015), but we're still lacking a first principle theory for this dynamo.
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 α 0.02 ± 0.01 (stratified, isothermal, zero vertical flux) (9.1) (e.g. Stone et al. 1996;Davis et al. 2010). In the outer region of protoplanetary discs, 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 authors have explored non-isothermal models, such as Hirose and collaborators (Hirose et al. 2006(Hirose et al. , 2009(Hirose et al. , 2014, Flaig and collaborators (Flaig et al. 2010(Flaig et al. , 2012 and Bodo and collaborators (Bodo et al. 2013(Bodo et al. , 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  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 α 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 α ∝ 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.

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  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 (or not) 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 . It was quickly realised that ambipolar diffusion could dramatically reduce the strength of MRI turbulence in this layer, and even suppress it, since These results imply accretion ratesṀ 10 −10 M /yr in the region 1 − 30 au, which are too low by at least 2 orders of magnitude compared to observational constraints. However, the zero mean field configuration is rather artificial. Indeed, protoplanetary discs are expected to form from a collapsing cloud which 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 & 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 et al. (2010), thanks to more robust numerical techniques compared to the initial attempts of Stone et al. (1996) and carefully designed vertical boundary conditions. Their simulations have a relatively weak field β mean 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 β mean = 10 4 compared to the zero net field case (Suzuki & Inutsuka 2009). This problem was revisited by a large number of authors, both in the strong field limit β mean 1 (Ogilvie 2012;Moll 2012;Lesur et al. 2013, Fig. 36), and in the intermediate regime β mean = 10 3 − 10 4 . It was quickly demonstrated that these outflows were a variation of Blandford & Payne (1982) magneto-centrifugal paradigm Fromang 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.
First, the mass loss rate depends on the vertical box extension , 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 → ∞. Some authors have proposed to fix this issue by adding higher order terms to the vertical gravity force, making possible to escape at z = ±∞ with a finite amount of mechanical energy (Suzuki et al. 2010;McNally & Pessah 2015). This approach, however, violates the conservative nature of gravitation, since 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 the 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 (Fig. 37). At the linear level (see §6.4.6), the two symmetries show similar growth rates and properties. However, in the nonlinear regime, the odd symmetry is problematic for the physical interpretation of the outflow. By simply looking at the poloidal field topology (Fig. 37), it is clear that the odd symmetry connects the bottom side of the disc to the central star, while the top side is connected to R → +∞. Physically, it means that no angular momentum is actually extracted from the disc: angular momentum comes from an unknown source at z → −∞ and flows through the disc up to z → +∞. The angular momentum conservation equation (4.10) clearly describe this situation: in the odd symmetry case, one has B y B z (z) = B y B z (−z), so the surface stresses on both sides of the disc cancel out and the disc is not accreting any material! Naturally, this situation is rather unrealistic and the only physical configuration is Figure 36. Growth and saturation of the MRI in the presence of a strong vertical field (βmean = 10) in ideal MHD. Tubes are magnetic field lines while 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). Figure 37. Symmetries of the outflow configuration allowed in a shearing box. Poloidal field lines are represented in green while the toroidal field component By < 0 is shown in blue and By > 0 in red. Left: odd symmetry configuration, right: even symmetry configuration. The usual Blandford & Payne (1982) picture corresponds to the even symmetry in a shearing box model. 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 β mean 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 → ±∞. Some authors have enforced the even symmetry manually ) 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 β mean ∈ [10, ∞] (Salvesen et al. 2016) shows that α = 10.1 β −0.53 mean for β mean < 10 5 (9.2) α recovering the zero net flux value when β mean > 10 5 . This scaling leads to α 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, α increasing with increasing Pm , 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 above. 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 ζ and a torque parameter υ which respectively measure the mass and angular momentum evacuated by the magnetised wind.
It is customary to define the outflow rate as where "top" and "bottom" subscripts stand for the top and bottom of the disc (which in principle can be chosen freely) and 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. Since there is technically no accretion flow in shearing boxes, this replenishment is absent, leading to boxes which slowly loose mass. The typical timescale over which a box looses a significant fraction of it mass is given by τ loss ∼ (ζΩ) −1 . This implies that the shearing box model is really valid in the limit ζ 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 4H . In the vertical direction however, as the box gets higher, the outflow rate decreases. Fromang et al. (2013) and Lesur et al. (2013) show that doubling the vertical extension of the box divides ζ by a factor 3 for β mean = 10 4 and by a factor 1.6 for β 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 the Hill's approximation (see 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 10H, the combination of the data published by Suzuki & Inutsuka (2009),  and Fromang et al. (2013) leads to: ζ 5 × 10 −5 + 10 β mean for β 10 2 (9.4) This relation gets shallower for β mean < 10 2 Lesur et al. 2013), and eventually leads to a sharp decrease of ζ when β 1 (see Fig.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 . However, this dependence is complex since it likely depends on the magnetisation. It should also be pointed out that this relation is close to the one 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, υ, defined as: where we have defined the stress tensor T yz = ρw y w z − B y B z /4π. 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 υ is however notoriously difficult in shearing boxes because its sign is not well-defined (see the symmetry discussion above). A naive 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)  where the uncertainty is due to the different method of measurements found in the literature. As for ζ, this scaling probably gets shallower for β mean < 10 2 . However, this has never been properly evaluated in the literature.
It should be noted that for β mean > 1, υ < α. Hence, the vertical stress is always weaker than the radial stress, by a factor O(β 1/2 mean ). However, the respective contribution of these two terms to the mass accretion rate in the disc is proportional to (R/H)(υ/α) (see §4.5), implying that the vertical stress is the dominant mass accretion mechanism whenever β mean (R/H) 2 .

Non-ideal MHD
Simulations with a mean vertical field first focused on the impact of ohmic and ambipolar diffusion 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 (Lesur et al. 2014;Bai 2014;Simon et al. 2015;Bai 2015) showed similar outflows, with in addition the presence of a midplane laminar stress due to the Hall effect (Lesur et al. 2014).
The difference between Ohmic only, Ohmic+Ambipolar, and Ohmic+ambipolar+Hall effect is demonstrated in Fig. 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 since the linear analysis shows little difference between the Ohmic and Ohmic+Ambipolar cases regarding the localisation and growth rates of the eigenmodes (see Fig. 26). The difference between Ohmic and Ambipolar cases is due to the fact that η A ∝ B 2 while η O does not depend on B. As eigenmodes grow, the horizontal field grows as well and η 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.
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 which subsists. The midplane stress is tightly linked to the field polarity: a negative mean field (Ω · B < 0) leads to its disappearance. This laminar stress is due to the HSI (Lesur et al. 2014) which only shows up for positive polarities (see 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 ∼ 30 AU but disappears at 100 AU in MMSN models  since 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 β 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 , or that the HSI spontaneously saturates in an even configuration (Lesur et al. 2014), or a combination of these effects.
Using similar technics to the one 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 Fig. 39 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 α, while the vertical (wind) transport υ appears to be the less affected. This difference in behaviour is the main reason why winds are today favoured in protoplanetary discs: the efficiency of winddriven transport is less affected by non-ideal MHD effects than radial angular momentum transport. In addition to these remarks, one 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 (Fig. 10). As already guessed in the discussion above, aligned field tend to  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 βmean have been slightly shifted horizontally to improve readability.
have a larger α, ζ and υ. It is interesting to also notice that, for simulations at 30 AU, the scaling with β mean is similar to the one found in ideal MHD, apart from a constant offset.

Outflow induced self-organisation
There are several 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 β mean 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 intrication 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 (Fig. 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). 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 Fig. 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 fig. 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 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 above, 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.
! 2 Figure 41. Self-organisation feedback loop proposed by Riols & Lesur (2019). Consider a small density deficit (top left), this deficit induces a radially converging flow (top right), which drags the poloidal field line inwards (bottom right). The stronger field, leads to a more efficient local ejection index which empties even more the region (bottom left).

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 β mean (R/H) 2 . When non-ideal effects are introduced however, the radial transport of angular momentum is significantly reduced, while 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 due to the HSI and the existence of a self-organisation mechanism associated to 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 properly modelled.

Global models of protoplanetary discs
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 protoplanetary discs which exhibits accretion and outflows.

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, φ, z), the conservation of mass and momentum read ∇ · ρu p = 0 (11.1) where the index p denotes the poloidal components (R, z) of the vector fields, Ω ≡ u φ /R (Ω is not necessarily Keplerian) and ψ = GM/(R 2 + z 2 ) 1/2 is the gravitational potential, assumed to be due solely to the central star. (11.3) can be recast in conservative form to obtain an angular momentum conservation equation which shows that the magnetic field can carry angular momentum in the Maxwell stress B p B φ . In addition to the equation of motion, one needs to solve the induction equation in the ideal regime in the poloidal and azimuthal directions ∇ × u p × B p = 0, (11.6) Finally, the magnetic field has to satisfy the solenoidal condition ∇ · 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) 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 constant of motion which are useful to interpret outflow solutions.

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 from which one deduces that where β is an unknown scalar function (physically, β/R is the φ component of the electromotive force E φ ). Using the radial induction equation, we may see that ∂ z β = 0 so that β is a constant. To avoid any singularity for E φ at R = 0, we are then forced to have β = 0 and hence B p and u p are parallel, i.e. u p = µ(R, z)B p (11.9) We can inject this relation in the mass conservation equation to obtain B p · ∇(ρµ) = 0, which indicates that η ≡ ρµ is a constant along magnetic field lines, i.e. η = η(a). We can then recast (11.9) into which constitutes the frozen in condition where η describes the amount of mass loaded along a poloidal field line. η is a direct measure of the gas density ρ A at the Alfvén surface, which is the region defined by u p = V A,p = B p / √ 4πρ: (11.11) 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 φ can be non-zero in principle. This implies that in a stationary shearing box solution, the poloidal field and velocity are not necessarily parallel . 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.

Magnetic surface rotation
The rotation of magnetic surfaces is obtained by substituting the frozen-in condition (11.10) into the azimuthal component of the induction equation (11.7) : which allows us to define a new invariant along field lines 11.12) where the A indices denotes 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 Ω * (a), the total field is parallel to the total velocity. In other words, there is no induced electromotive force in the frame rotating at Ω * (a). Ω * can therefore be interpreted as the rotation speed of magnetic surfaces. By combining (11.10) and (11.12), we can express the total velocity as a function of the magnetic field, (11.13) This illustrates the fact that only the poloidal components of the field and the velocity are parallel in general. In the particular case where η = 0 (no motion along the poloidal field lines), we recover Ferraro's iso-rotation law Ω = Ω * (a).

Angular momentum
The angular momentum invariant is easily deduced from the angular momentum conservation equation (11.5): from which we deduce the angular momentum invariant L(a): L(a) ≡ ΩR 2 − RB φ η(a) (11.14) = Ω * (a)R 2 A (11.15) where the last equality has been obtained at the Alfvén surface using (11.11). 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 φ .

Bernoulli invariant
The Bernoulli invariant is obtained from the scalar product of the poloidal equations of motion (11.2-11.4) with u. One gets where ψ G = GM/(R 2 + z 2 ) 1/2 is the gravitational potential. We then use (11.13) to obtain the work of the Lorentz force: Additionally, we may express the work of pressure forces using the enthalpy per unit mass H u p · ∇P ρ = u p · ∇H + δQ in which we have also considered the effect of an additional heating term due to radiative heating/cooling δQ. Putting all the terms together and integrating along one particular streamline a, we find that the quantity δQ ds (11.16) is conserved along poloidal field lines and streamlines. Note that since the heating term is not a proper differential, it depends on the integral of the heating term along the chosen streamline S(a). Of course, for an outflow to escape up to z → ∞, on needs B > 0 on the streamlines (assuming H > 0). In the midplane of a Keplerian disc, one has u 2 /2 + ψ G = −v 2 K /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 : where 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 protoplanetary discs community).
Cold MHD winds :where thermal effects are neglected and the toroidal field acts as an energy reservoir to launch the outflow.

Dimensionless numbers characterising an outflow
Based on the MHD invariants, it is possible to define a series of dimensionless numbers which characterise an outflow streamline unambiguously. These numbers use the physical properties at the base of the outflow. Let us therefore write Ω 0 the rotation rate at the base of the outflow and R 0 its cylindrical radius (one typically has Ω 0 = Ω K (R 0 )), and B 0 the poloidal field strength threading the disc at the location where the outflow is launched. We then define The Bernoulli invariant can be easily expressed in terms of the other invariants, using the dimensionless rotation rate ω ≡ Ω/Ω 0 : where θ = 2 H + a=cste δQ ds /Ω 2 0 R 2 0 is the dimensionless thermal energy content of the flow. This expression clearly demonstrates the contribution of λ 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 Ω 0 R 0 , ω = 1 and R = R 0 : e 2ω * (λ − 1) + θ − 1. (11.18) Obviously, for the outflow to be able to propagate up to infinity, one needs e > 0. We recover here the two extreme limits discussed above: purely thermal winds, which have λ = 1 and require θ > 1, or cold MHD winds with θ = 0 which need For all practical applications of MHD outflows, one has ω * 1 to a very good approximation. This implies that an outflow can exist only if λ > 3 2 (11.19) 11.6. Ejection efficiency The existence of an outflow is tightly linked to the process of accretion happening inside the disc since 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 disċ 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 in and the radius under consideratioṅ These two quantities are connected by the continuity equation (4.9): At this stage, it is useful to introduce the ejection efficiency index ξ ≡ d logṀ acc d log R (11.20) (11.21) 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, ξ = 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.10) aṡ where we have define the stress tensor W iφ ≡ ρv i v φ − B i B φ /4π. This allows us to define the radial and vertical contribution to the accretion of the disc τ Rφ and τ zφ . It is then useful to introduce the ratio of these two quantities Λ ≡ τ z τ r so that the accretion rate is simplyṀ It is then possible to relate the vertical torque τ zφ 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 where the second equality assumes the outflow is top/down symmetric and that the poloidal field strength doesn't 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 . The vertical mass flux being directly to the radial derivative of the outflow rate, we can express the accretion rate aṡ From which we get an expression for the mass ejection index This relation reveals several key features of MHD outflows. First, it shows that large level-arms λ are associated to small ejection indices. Interestingly, this result does not depend on the radial contribution to the angular momentum budget Λ. Second, the energy constraint (11.19) imposes an upper bound on ξ in cold MHD winds: ξ Λ/(Λ + 1) 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 ξ 1 therefore necessarily require some thermal energy driving to escape the potential well, as one would expect.

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 constrains, 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 will make the assumption that the global outflow is top/down symmetric, so that the invariant are the same on both sides of the disc.
We first start with the outflow rate ζ (defined in 9.3), which can be easily connected to the mass loading parameter κ The magnetic level-arm can also be related to υ (defined in 9.5), 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 The energetic constrain (11.19) therefore leads to a new constrain in shearing boxes: 2ζ/υ < H/R. This constrain cannot be satisfied by the shearing box scalings (9.4) and (9.6) when H/R 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

Disc region
Wind 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 ΩK (R0) 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 Λ to α and υ. For this, let us assume that the disc surface density follows a power law: Σ = Σ 0 (R/R 0 ) −p , that α and H/R are constant with radius, and that the disc is vertically isothermal. Under these assumptions, the contributions to the mass accretion are: so that the ratio simply reads The scalings (9.4) and (9.6) then suggest that Λ → 0 when β mean → ∞. Combining these results with 11.22, this implies that the ejection index ξ tend to decrease as β mean → ∞, as one would expect.

Outflow phenomenology
The launching of an outflow is tightly linked to physics of the accreted material since 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 ( Fig. 42) and describe the physical processes in each region.

Disc region
In the disc region, the flow is mostly accreted, while 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: 1-assuming that the disc ionisation fraction is extremely low as in protoplanetary discs, so that ambipolar and ohmic diffusion are very large (Konigl 1989;Wardle & Konigl 1993), or 2-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 φ 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 ∂ t B φ B R R∂ R Ω. This effect is actually the main energy source for the outflow, which converts shear energy into magnetic energy stored in B φ 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 φ (z = 0) = 0 follows. The bending of poloidal field lines in the disc implies ∂ z B R > 0 for z > 0 and hence ∂ z B φ < 0 up to the disc surface. The actual value of B φ depends on the competition between the shearing of B R and the magnetic diffusion which damps the shear amplification. Overall, one can estimate B + φ at the disc surface by balancing shear and diffusion where h is the disc scale height, η is the magnetic diffusivity and θ is the inclination of the poloidal field at the disc surface (see Fig. 42). Accretion of the disc material naturally follows from the profile of B φ (z). The toroidal field is responsible for an azimuthal magnetic tension force B z ∂ z B φ /4π which slows down the rotating material and lead to accretion. We can deduce the accretion rate from (4.10): where ρ 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 due to a transfer of angular momentum from the accreted material to the toroidal field. In the end, the angular momentum is stored in B + φ 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 since B 2 R and B 2 φ are both increasing functions of z. Gravity is also directed downward. The magnetic tension term is usually small (it involves a radial derivative while the other terms involve vertical derivatives, which are larger by a factor 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 + φ cannot reach arbitrarily large values which would otherwise prevent thermal pressure from pushing materials to the wind base and vertically squeeze the disc. Ejection therefore requires (B + φ ) 2 /8π P + (Ferreira 1997). To summarise the physics of the disc region: • The region is non-ideal as the gas has to be allowed to stream through poloidal field lines which are stationary.
• the toroidal field at the disc surface stores the energy and angular momentum taken from the accreted material • the upward motion needed to replenish the outflow base is due to the vertical thermal pressure gradient. This sets a limit to the amount of toroidal field that can be stored at the disc surface since magnetic pressure prevents this upward motion.

Outflow region
In order to understand the dynamics of the outflow from its base, it is convenient to introduce the Alfvénic Mach number ξ defined as where the second equality is deduced from (11.10). In the outflow, we can expect a continuous acceleration of the flow, so that ξ is an increasing function as one moves along a poloidal streamline, with ξ ∼ 0 at the outflow base. We can eliminate B φ in favour of the angular velocity by combining the magnetic surface rotation invariant (11.12) and the angular momentum invariant (11.14) to obtain This expression looks singular for ξ = 1 which corresponds to the Alfvén point. However, at this particular point, R = R A so that Ω is actually smooth across this point. Second, in the limit of low Mach numbers ξ 1 which corresponds to the base of the outflow, we can expand this expression to get Hence, the outflow is rotating at a constant angular velocity up to the point where ξ 2 1/(R 2 A /R 2 − 1). Physically, the angular momentum stored in B + φ 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 ξ → ∞, we find Ω Ω * 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 above, angular momentum stored in B φ is converted into kinetic angular momentum. Hence B φ decreases along the streamline. This leads to a magnetic pressure force ∂ z B 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 ∂ z B 2 R accelerating the outflow vertically. Overall, the vertical acceleration is a magnetic pressure effect due to the decrease of B φ 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 + δR, δz). Since we focus on the launching region of the outflow, we assume that Ω Ω * . During this displacement, the Bernoulli invariant should be conserved, so we should have where the inequality comes from the assumption that u p (R 0 + δR, δz) > u p (R 0 , 0) (i.e. the flow accelerates). It can be converted into a criterion on θ: tan θ > 1/ √ 3 or θ > 30 • . 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 sufficiently inclined, 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: (i) Magnetic field lines do not behave as rigid poloidal wires. As shown above, the magnetic configuration has a strong toroidal field at the wind base B + φ which indicates that the field is wound up in this region.
(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 ∼ 1 will therefore experience a very small region of solid rotation (if any).
(iii) As discussed above, vertical acceleration is always a magnetic pressure effect which 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 • The flow is accelerated azimuthally by the magnetic tension due to B φ . This results in a transfer of angular momentum from B φ to the gas as it accelerates.
• Since B 2 φ (and possibly B 2 R ) decrease with z, the flow is accelerated vertically by magnetic pressure.
• Close to the launching point and far from the Alfvén radius, the outflow is approximately in solid rotation.
• In a cold wind, energy conservation at the wind base implies θ > 30 • .

Global numerical models
Many numerical models of protoplanetary discs have been published in the literature. Here, we will focus on models tackling the impact of non-ideal MHD effects and winds on the dynamics of a disc, since we focus on the outer parts of protoplanetary discs. We therefore exclude "ideal MHD" models and simulations without a large scale magnetic field.
This first model to investigate this regime was published by .
However, the limited vertical extension of this work (typically 4 scale heights) makes the interpretation of outflow properties difficult. Therefore, we focus on models with a larger vertical extension such as the ones published by Béthune et al. (2017) and Bai (2017).
13.1. Global topology Figure 43. Global simulation for a wind in a protoplanetary disc 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 φ . Close to the midplane, the configuration has an odd symmetry, as found in shearing box models and a current sheet is found at z ∼ 3H. Figure  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 since in this case, the gravitational potential is not symmetrical with respect to r. Many models seem to converge towards dissymmetric outflows, i.e. outflows which 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Ṁ w reaching a factor of a few). This kind of solution was found both with only Ohmic and ambipolar diffusion  and with all the three non-ideal effects (Béthune et al. 2017;Bai 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.  have proposed that the global field structure was playing a significant role in shaping the outflow. Indeed, by choosing initially a field configuration with ∂ R B z > 0, it is possible to produce outflows directed inwards (towards the star). It's been therefore proposed that the global radial magnetic pressure gradient was responsible for the symmetry breaking observed in these models. Also, the vertical rotation profile (and therefore the vertical temperature structure) might be playing a role by shearing the poloidal field . The question of the origin of the global outflow configuration therefore remains open.
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 (β mid 10 4 ) (Béthune et al. 2017;Bai 2017) and for anti-aligned cases (B z Ω < 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 Σ v r ρ τ r + τ z τ r τ z 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 due to the wind (τz term) while the radial stress does not seem to contribute to the angular momentum budget. Figure  a magnetic configuration as close as possible to the configuration expected from core collapse calculations (Bai 2017).

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 equatioṅ where we have defined the mass accretion rateṀ acc = −2πRρv r and the radial and vertical torques τ r,z . An example of such a measure is given in Fig. 44. We find that accretion is mostly due to the wind, the surface torque being the main contribution to angular momentum extraction in the disc. 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 α 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Ṁ acc ∼ 10 −8 − 10 −7 M /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 the one derived in the bulk of the disc. From the definition of the ejection efficiency (11.20), we havė Large ejection efficiencies (ξ = O(1)) as the one 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 the one 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.

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 to, or even larger than the mass accretion rate in the disc (Ṁ w 10 −8 − 10 −7 M /yr), which implies ξ ∼ 1.
The mass outflow rate is tightly connected to the amount of flux threading the disc, withṀ w ∝ β −1/2 mid (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 Lesur et al. 2014) while steeper dependencies are found in ideal shearing box models (Suzuki & Inutsuka 2009;. 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 λ < 2. This is coherent with the very high mass loss rates found in these simulations (high ξ, see Eq. 11.22). While it is in principle possible to obtain cold outflows with low λ and high ξ in discs threaded by a weak field (β mean 1, see for instance Jacquemin-Ide et al. 2019), some of the outflows published to date in protoplanetary discs have λ < 3/2 (e.g. Béthune et al. 2017), which violate the cold MHD wind constraint (Eq. 11.19). 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. fig. 45). Such an analysis shows that thermal effects (enthalpy and heating terms) both contribute significantly to the energetic 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.
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 to historical cold wind solutions, they are (obviously) warmer, denser and slower. They reach high ξ values (typically ξ > 0.1) and have moderate λ. Of course, the fact that they extract angular momentum from the disc and that the initial acceleration is due to 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 β mid 10 3 , but it is absent from the models of  and Bai (2017), who only considered β 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.
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, while 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 which is at work in these models but ambipolar driven self-organisation. In essence, the mechanism seems to be similar to the one driving selforganisation 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 Fig. 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 selforganisation is occurring nor what are the general properties of the structures which are formed.

Conclusions
In conclusions, the global modelling of the weakly ionised part of protoplanetary discs (R 1 AU) is still in its infancy. The inclusion of non-ideal MHD effects confirms that the disc midplane is mostly laminar, while the disc is still accreting thanks to magnetised outflows. Because of the low magnetisations used in these models (β mid 10 2 − 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 β = 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: β mean 1. Self-organisation is also a very promising mechanism to explain some of the observables in protoplanetary discs. 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.

Summary
These are exciting time for protoplanetary disc 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 nonaxisymmetric 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 weekly ionised, reaching ionisation fraction ξ ∼ 10 −13 in the disc midplane around one astronomical unit. This implies that non-ideal MHD effects (Ohmic, Ambipolar and Hall effect) are essential to get 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 nonideal 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 since the input parameters (disc composition & environment) are largely unknown.
In order to explain the observed accretion rates in these discs, and since 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. While radial transport has been historically favoured thanks to its simplicity and elegance (the well-known α 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 magnetorotational instability (MRI), a linear, ideal MHD instability found in rotating sheared flows. This instability has been the subject of intense studies since the early 90s and it is known that it is strongly affected by the non-ideal effects present in protoplanetary discs. Most notably, it is suppressed in the regions where Ohmic and ambipolar diffusion are strong, and it gives a new branch, known as the Hall shear instability (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 α 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 protoplanetary discs, it is found that accretion driven by magnetised winds are compatible with observed accretion rates for relatively weak mean fields, β 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 consistant 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".

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 et al. 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 winds observed in the sub-millimetric range.
A second question is the dynamical evolution of the mean magnetic field threading the disc. As shown above, 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) while 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, since 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 impact 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 α 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 ( 10µ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 impact (Kimmig et al. 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.