Dust settling and rings in the outer regions of protoplanetary discs subject to ambipolar diffusion

MHD turbulence plays a crucial role in the dust dynamics of protoplanetary discs. It affects planet formation, vertical settling and is one possible origin of the large scale axisymmetric structures, such as rings, recently imaged by ALMA and SPHERE. Among the variety of MHD processes, the magnetorotational instability (MRI) has raised particular interest since it provides a source of turbulence and potentially organizes the flow into radial structures. However, the weak ionization of discs prevents the MRI from being excited beyond 1 AU. The strong sedimentation of millimetre dust measured in T-Tauri discs is also in contradiction with predictions based on ideal MRI turbulence. In this paper, we study the effects of non-ideal MHD and winds on the dynamics and sedimentation of dust grains. We consider a weakly ionized plasma subject to ambipolar diffusion characterizing the disc outer regions (>1 AU). For that, we perform numerical MHD simulations in the stratified shearing box, using the PLUTO code. Our simulations show that the mm-cm dust is contained vertically in a very thin layer, with typical heightscale ~0.4 AU at 30 AU, compatible with recent ALMA observations. Horizontally, the grains are trapped within pressure maxima induced by ambipolar diffusion, leading to the formation of dust rings. For micrometer grains, dust and gas scaleheights are similar. In this regime, the settling cannot be explained by a simple 1D diffusion theory but results from a large scale 2D circulation induced by both MHD winds and zonal flows. Overall, our results show that non-ideal MHD effects and their related winds play a major role in shaping the radial and vertical distribution of dust in protoplanetary discs. Leading to substantial accretion efficiency, non-ideal effects also a promising avenue to reconcile the low turbulent activity measured in discs with their relatively high accretion rates.


Introduction
One of the most challenging problem in astrophysics is to understand how protoplanetary (PP) discs accrete and form planets. For more than forty years, accretion has been thought to originate from turbulent motions, acting like an effective viscosity (Shakura & Sunyaev 1973). The magneto-rotational instability (MRI, Balbus & Hawley 1991;Hawley et al. 1995) has long been considered as the main source of turbulence, among various other instabilities (of gravitational or thermodynamics nature) whose role is still debated. Local and global MRI simulations in the ideal limit suggest in particular that the rate of angular momentum transport is compatible with that inferred from observations (Flock et al. 2011(Flock et al. , 2013. However the MRI has been shown to survive only in the very inner part of the disc, below 0.1 -1 AU where the ionization sources are strong enough to couple the gas with the ambient magnetic field (Gammie 1996). Further away, non-ideal MHD effects, such as ohmic diffusion, Hall effect and ambipolar diffusion prevail and tend to suppress any form of MRI-driven turbulence (Fleming et al. 2000;Sano & Stone 2002;Wardle & Salmeron 2012;Bai 2013;Lesur et al. 2014;Bai 2015).
Observationally, the prevalence of a sustained and vigorous MRI-turbulence in the regions 1 AU is also largely questioned. The high resolution imaging at optical, near-infrared (VLT/SPHERE) and sub-millimetric wavelenght (ALMA) have provided precious information about the disc structure and turbulence, essentially by probing the spatial distribution of micrometer and millimetre dust grains. Few observations have been able to map the dust heightscale, by measuring the rings/gaps contrast (Pinte et al. 2016) or directly by imaging edge-on discs (Wolff et al. 2017). In particular, the survey of Pinte et al. (2016) in HL Tau suggests that the millimeter dust forms a very thin layer around the midplane with typical scaleheight 2 AU at a radius of 100 AU. This result, together with the measurements of gas velocity dispersion in CO (e. g. Flaherty et al. 2017) indicate that the turbulence in protoplanetary discs is weak, contrary to the current theoretical predictions based on the ideal MRI. In particular the accretion efficiency that one would deduce from such low turbulent levels is smaller by one or two orders of magnitude than those inferred from direct observations (see Table 4 of Venuti et al. 2014). These elements suggest that accretion is happening through an essentially laminar process which is yet to be identified Article number, page 1 of 22 arXiv:1805.00458v1 [astro-ph.EP] 1 May 2018 A&A proofs: manuscript no. article_aa In parallel to that, ALMA and SPHERE have unveiled a lot of structures in the scattered light emission or the dust continuum, such as rings, spiral arms, gaps or vortices (Garufi et al. 2017). Numerous effort has been deployed to explain the formation of concentric rings, a feature observed in various protoplanetary discs such as HL tau (ALMA Partnership et al. 2015), TW Hydra (Andrews et al. 2016) or the disc around Herbig Ae star HD 163296 (Isella et al. 2016). If the presence of planets is often invoked to explain these structures, there is today no consensus on their origin. Many other scenarii have been proposed such as dust-drift-driven viscous ring instability (Wünsch et al. 2005;Dullemond & Penzlin 2018), snow lines (Okuzumi et al. 2016), MHD turbulence with dead zones (Flock et al. 2015), zonal flows (Johansen et al. 2009) or secular gravitational instabilities in the dust (Takahashi & Inutsuka 2014).
A possible solution to account for all these observations is to consider the atypical dynamics induced by non-ideal MHD effects in protoplanetary discs, in particular the Hall effect and ambipolar diffusion. Indeed, if the latter tend to weaken considerably the disc turbulence, both are able to produce large scale horizontal magnetic fields within the disc, producing strong currents. These currents allow the development of vigorous outflows at the disc surface and a laminar Maxwell stress in the midplane (Lesur et al. 2014), both transporting significant angular momentum (Bai 2013;Simon et al. 2015;Gressel & Pessah 2015). Moreover, these non-ideal effects lead in the nonlinear regime to the formation of self-organized axisymmetric structures, such as zonal flows (Kunz & Lesur 2013;Bai 2015), corresponding to concentric rings in a global disc view (Béthune et al. 2016). These rings, which are very stable and coincide with gas density maxima, are likely to be ideal locations for dust accumulation and planet formation.
Whether the Hall effect and ambipolar diffusion are key physical processes accounting for the observed disc quiescence (in terms of turbulent motions) and ring-like structures remains an open question. Since current sub-millimeter or near-infrared observations are sensitive to dust, direct comparison with them requires to go one step forward and determine the dust dynamics in such non-ideal MHD flows. Ultimately it requires to compute its scattered and direct emission, but this remains out of the scope of the paper. Although the dust behaviour has been widely studied in local and global MRI disc simulations (Johansen & Klahr 2005;Carballido et al. 2006;Balsara et al. 2009;Zhu et al. 2015), it remains relatively unexplored in non-ideal MHD simulations characterizing the disc outer regions (r 0.1−1AU). Ruge et al. (2016) studied the case of zero net flux MRI with ohmic resistivity but did not include other relevant non-ideal effects. Several simulations by Zhu et al. (2015) have combined the dust dynamics with ambipolar diffusion but focused into a regime of large particle size ( mm), in relatively small domains compared to the characteristic length of the zonal MHD flows.
Our aim is to generalize the study to smaller particles sizes and more extended domains. A first step is to characterize the degree of dust sedimentation, and the ability of non-ideal MHD effect to concentrate the grains into rings. The goal is to make comparison with recent ALMA observations for millimeter size particles. Our ambition is also to provide a theoretical model that predicts the dust scaleheight and its vertical distribution, as a function of the disc magnetization and particle size. Current models of dust settling in turbulent discs are based on a 1D diffusion theory (Dubrulle et al. 1995;Dullemond & Dominik 2004) which is not necessarily appropriate to describe non-ideal MHD flows characterised by nearly laminar midplane and windy corona.
To address the problem, we performed 3D MHD shearing box simulations of stratified discs with an imposed vertical magnetic field, using a modified version of the PLUTO code. To make the interpretation easier, we restrict our study to a large radius where the Hall effect can be neglected. The dust population is approximated as a pressure-less fluid made of different particle sizes. We explore different grain sizes from few micrometers to few centimetres and different disc magnetization with plasma beta ranging from 10 3 to 10 5 . We enable the back reaction of the dust on the gas. As a preliminary step, we do not include radiative transfer and neglect coagulation or fragmentation processes. The paper is organized as follows: in Section 2, we describe the model and review the main characteristics of non-ideal physics and dust/gas interaction. We also present the numerical methods used to simulate the dust/gas dynamics. In Section 3, we perform numerical simulations without dust in order to characterise the properties of the flow subject to ambipolar diffusion (transport, turbulent levels, winds). In Section 4, we add the dust in the simulations and determine its vertical and radial distribution. We show in particular that the sub-millimetre dust is highly sedimented and form ring structures with high density contrast. We also explore the effect of a mean radial pressure gradient in the disc. In Section 5, we compare different settling models and show that the classical 1D model fails to predict the vertical dust distribution for small particle sizes or large magnetization. We then propose a 2D model including the effect of coronal winds and radial gas density structures (zonal flows). Combined each other, these two features induce a large scale circulation of the dust in the corona, which affects the sedimentation process. Finally, in Section 6 we compare the vertical and radial dust distributions with those inferred from ALMA observations (Pinte et al. 2016). We conclude in Section 7 by discussing the applications of our work on protoplanetary discs evolution and planet formation.

MHD and dust equations in the shearing sheet
To study the gas and dust dynamics in the outer part of protoplanetary discs, we use the local shearing sheet model (Goldreich & Lynden-Bell 1965). This corresponds to a Cartesian patch of the disc, centred at r = R 0 , where the Keplerian rotation is approximated locally by a linear shear flow plus a uniform rotation rate Ω = Ω e z . We note (x, y, z) respectively the radial, azimuthal and vertical directions. We adopt a multi-fluid approximation in which the ionized gas and the dust interact and exchange momentum through drag forces.
The gas is coupled to a magnetic field B and is assumed to be inviscid and isothermal, its pressure P and density ρ related by P = ρc 2 s , with c s the uniform sound speed. The evolution of its density ρ, and total velocity field v follows: where the total velocity field can be decomposed into a mean shear and a perturbation u: v = −qΩx e y + u.
(3) with q = 3/2. The term g = −Ω 2 z e z + 2qΩ 2 x e x corresponds to the tidal gravity field in the local frame. The last term in the momentum equation (2) contains the acceleration γ d−>g exerted by the dust drag force on a gas parcel (detailed below). The magnetic field B, which appears in the Lorentz force (∇ × B) × B, is governed by the non-ideal induction equation, The term ∇ × η A (J × e b ) × e b corresponds to ambipolar diffusion, with J = ∇ × B the current density, e b = B/ B the unit vector parallel to the field line and η A the ambipolar diffusivity. This diffusion is due to ions-neutral collisions and is assumed to be the dominant non-ideal effect in the regions considered in this paper, i.e R 0 = 30 AU (see justifications in Simon et al. 2015). Therefore, the Hall effect and ohmic diffusion are neglected. Due to our assumption of ionization, η A is not uniform and its vertical profile is detailed in Section 2.4.
Finally, the dust is composed of a mixture of different species, characterizing different grain sizes. Each specie is described by a pressure-less fluid, with a given density ρ d and velocity v d . We assume in this paper that dust grains remain electrically neutral so that they do not feel the magnetic field. The equations of motion for each specie are: with γ g−>d the drag acceleration induced by the the gas on a dust particle. If we label each dust specie by a subscript k, we obtain by conservation of total angular momentum: In the next section, we specify the form of the drag term γ g−>d and its dependence on the particle size.

Drag and stopping times
In this study we consider that dust particles are spherical and small enough that they are in the Epstein regime (Weidenschilling 1977). The acceleration associated with the drag and acting on a single particle is given by: For a spherical particle of size a and internal density ρ s (should not be confused with the fluid density ρ d ), the stopping time τ s is This corresponds to the timescale on which frictional drag will cause an order-of-unity change in the momentum of the dust grain. This is a direct measure of the coupling between dust particles and gas. A useful dimensionless quantity to parametrize this coupling is the Stokes number If not precised, St denotes the Stokes number in the midplane. Note that the effective Stokes number in the disc atmosphere is larger than St, since it is inversely proportional to the density.

Converting Stokes to particle size
In this paper, we preferentially use the Stokes number rather than particle size to describe the dust dynamics. Indeed St is a dimensionless quantity and does not depend on the disc properties and geometry. However, to allow comparison between our simulations and real systems, it is helpful to associate the Stokes number to a grain size.
In the limit of small magnetization, the disc is in hydrostatic equilibrium and its column (or surface) density is Σ = ρ 0 H √ 2π where ρ 0 is the midplane density and H = c s /Ω is the disc scaleheight. Thus, combining these different relations, we obtain: To evaluate Σ, a first possibility is to assume that the disc can be described through the standard MMSN model (Hayashi 1981), with surface density This gives Σ = 10.34 g/cm 2 at R 0 = 30 AU. However, this model is probably not realistic since recent observations suggest flatter profiles with a power law between -0.5 and -0.2 (Williams & McPartland 2016), and in some cases Σ varying up to 20 g cm −2 at R 0 = 30 AU (see Pinte et al. 2016, in the HL Tau disc for instance). By assuming ρ s = 2.5 g/cm −3 , we obtain the following conversion 2.4. Ionization profile and ambipolar diffusivity η A If the only charged particles are the ions and electrons, the ambipolar diffusivity η A is given by: (Wardle 2007), where γ i = σv i /(m n + m i ) and σv i = 1.3 × 10 −9 cm 3 s −1 is the ion-neutral collision rate. ρ n and ρ i are respectively the neutral and ion mass density; m n and m i their individual mass. By introducing the Alfvén speed v A = B/ √ ρ, we also define the dimensionless ambipolar Elsasser number To calculate ρ i , one generally needs to model the various ionization sources of the disc (X rays, cosmic rays, radioactive decay and FUV) and compare them to the dissociative recombination mechanisms. However, this requires to compute the complex chemistry occurring in the gas phase and on grain surfaces.
To simplify the problem, we use the same model as Simon et al. (2015) which captures the essential physics of disc ionization, i.e the effect of FUV in the disc corona. In this model, Am is constant and of order unity in the midplane, and increases abruptly up to a certain height z io . The height z io is free to vary during the simulation and corresponds to the base of the ionization layer, above which the FUV can penetrate. If we note Σ i (z) the horizontally-averaged density, integrated from the vertical boundary, we assume that FUV are blocked for Σ i (z) > Σ ic = 0.1 g/cm 2 . Given that carbon and sulfur atoms are fully ionized by FUV and immune to recombination on grains, the ionization fraction x i = ρ i m n /(ρm i ) above z io is 10 −5 (Perez-Becker & Chiang 2011). To avoid any discontinuity of this quantity at z = z io , we use a smooth function to make the transition between the midplane regions and the ionized layer, so that the ionization profile is The profile of Am is then given by: where Am 0 is the constant midplane value. Note that we have supposed an MMSN disc model (see Eq. 12) to convert numerical Σ i into g/cm 2 .

Numerical methods
We use a modified version of the Godunov-based PLUTO code (Mignone et al. 2007) to simulate the gas and dust motions. This code provides a conservative finite-volume scheme that solves the approximate Riemann problem at each inter-cell boundary. Physical quantities are evolved in time with a second order Runge-Kutta scheme, and the orbital-advection algorithm FARGO is used to reduce numerical dissipation. The implementation of ambipolar diffusion is described in Lesur et al. (2014) (Appendix B) and has been used in various studies (Lesur et al. 2014;Béthune et al. 2016Béthune et al. , 2017. For the dust, we implemented our own version which is described and tested in Appendix A. Our simulations are computed in a shearing box of size L x = 8H, L y = 8H and L z = 12H with resolution N X = 128, N Y = 64, N Z = 192. The inverse of the orbital frequency Ω −1 = 1 defines our unit of time and H = c s /Ω the disc scaleheight, our unit of length. Boundary conditions are periodic in y and shear-periodic in x. In the vertical direction, we use a standard outflow condition for both the gas and dust velocities but compute a hydrostatic balance in the ghost cells for the gas density. For the magnetic field, we adopt the so called "vertical field" boundary with B x = B y = 0 at z = ±L z /2. The presence of magnetic fields produces outflows, which depletes the mass contained in the box. Then, at each time step, we artificially multiply ρ in each cell by a uniform factor such that the total mass in the box is maintained constant. We checked that the mass injected at each orbital period is negligible compared to the total mass (≈ 1%) and does not affect the results. Note that the mean B z is conserved in the box but not the mean horizontal magnetic field, whose evolution depends on the electromotive forces (EMFs) at the boundary.

Initialization and main parameters
We fix R 0 = 30 AU in all runs. Non ideal MHD simulations without dust are initialized with hydrostatic equilibrium and weak random motions. The box is threaded by a net vertical field B z , which allows us to define the beta-plasma parameter with ρ 0 = 1 the midplane density. The midplane ambipolar Elsasser number is fixed to Am 0 = 1 in all cases. Simulations with dust are systematically initialized from developed non-ideal MHD states. Their detailed setup and the range of Stokes numbers considered are given in Section 4.1.

Diagnostics
To analyse the statistical properties of the MHD flow, we define the box average: where L x , L y and L z are the dimensions of a Cartesian portion of volume V. We also define the horizontally averaged vertical profile of a physical quantity: and its mean radial profile, averaged in y and z (over a size l z ≤ L z ): Finally we note X the time average of a quantity over the simulation time.
An important quantity that characterizes the turbulent dynamics is the coefficient α which measures the radial angular momentum transport efficiency. This quantity is the sum of the total stress, which includes Reynolds H xy = ρu x u y and Maxwell stresses M xy = −B x B y , divided by the box average pressure.

Ambipolar-MHD simulations without dust
A first and necessary step before computing the dust/gas interaction is to characterize the properties of the flow in outer regions of PP discs. For that purpose, we start with non-ideal MHD simulations without dust that include ambipolar diffusion. We examine the turbulent transport, velocities dispersion and flow structures at R = 30 AU for different vertical magnetization. The turbulent states obtained will then be used in Section 4 to initialize simulations with dust.

Transport and turbulence properties
To characterise the properties of non-ideal MHD flows with ambipolar diffusion, we perform three different simulations (assuming Am 0 = 1) with β = 10 5 , 10 4 and 10 3 respectively, which corresponds to B z = 0.0045, 0.014 and 0.045 in our code units. Figure 1 shows the evolution of the box averaged radial transport efficiency α for the three different cases. For all runs, a quasisteady state that transports a significant amount of angular momentum is obtained. The transport is an increasing function of B z , with values in agreement with those found in the literature. For instance, we obtain α 0.0034 for β 0 = 10 5 and α 0.02 for β 0 = 10 4 , which is identical to what Simon et al. (2013) found for a similar setup with the ATHENA code (see their Table 1). Figure 2 (top panels) shows the vertical profiles of horizontally-averaged Maxwell and Reynolds stress. Each of them contributes positively to the radial transport efficiency α (except H xy in the midplane for β = 10 3 ), although the Maxwell stress always dominates over the Reynolds stress. For all β, both quantities peak around z 2H where the ionization is sufficiently large to excite the MRI. Their profiles show a dip in the midplane where the ionization is weak and the MRI absent.
However, as the vertical field increases, a laminar stress due to the large scale B x and B y starts to be important in the midplane and induces a shallower dip in Maxwell stress. In particular, the ratio of Maxwell to Reynolds stress sharply increases with magnetization, from 2.5 at β = 10 5 to ∼ 100 at β = 10 3 , and thus deviates significantly from that of ideal MRI. We will see in Section 6 that such result has consequences on the so-called Schmidt number, ratio between the gas effective viscosity and dust turbulent mixing. Note that a transition to a magnetic and quasilaminar stress occuring around β 10 3 has been already pointed out by Simon et al. (2013) and also exists when the Hall effect is included (Lesur et al. 2014). We emphasize also that the radial stress W xy is not the only source of angular momentum transport (and thus accretion) in these type of flows. The vertical magnetic field also removes angular momentum from the disc via strong winds powered by the MRI in the active layers. This leads to a zy component of the stress tensor which is smaller than W xy but can influence the accretion rates (due to a factor R/H appearing in the equation forṀ , see Fromang et al. 2013;Simon et al. 2013) .
Finally, the bottom panels of Fig. 2 show the vertical profiles of the horizontally averaged r.m.s velocity u rms (z)=(u 2 ) 1/2 in all directions. These quantities are important to quantify the diffusion of dust particles (see Section 5). We checked that their dependence in z is similar to that of Bai (2015). For weak magnetization β = 10 5 and 10 4 , radial and vertical velocity dispersions in the midplane are comparable in both cases (≈ 2 × 10 −2 ) while they reach unity near the vertical boundaries. In the midplane region, they are mainly associated with the fluctuating part of the velocity field u − u(z). The mean component u(z) (orange dashed lines in Fig. 2) associated with a large scale wind remains negligible compared to the fluctuations up to z 4H. However, for higher magnetization (β = 10 3 ) , this wind component is stronger and becomes similar in amplitude to the fluctuating part for z 1.5H. Note that the azimuthal velocity dispersion increases significantly with B z in the midplane. This is a consequence of the emergence of large scale radial zonal flows whose amplitude is particularly enhanced at large B z (see next section).
3.2. Zonal flows, pressure maximas and sound waves Self-organized axisymmetric structures, or commonly called "zonal flows" are found in a variety of non-ideal MHD flows threaded by a net vertical field: Hall-dominated turbulence (Kunz & Lesur 2013), ambipolar-dominated flows (with or without ohmic diffusion, see Simon & Armitage 2014;Bai 2015) and flows combining all non-ideal effects (Lesur et al. 2014;Bai 2015). It has been pointed out that ambipolar diffusion enhances these features (Bai 2015) and can lead to zonal flows with density amplitude of 10-20% the background disc density (Simon & Armitage 2014). These zonal flows become large scale rings in the global configuration, as revealed by the simulations of Béthune et al. (2017).
To check whether these axisymmetric structures exist in our stratified simulations, we show in Fig. 3 the spacetime-diagram (t, x) of the mean gas density (top) and B z (bottom), averaged in y and z, for our three different β. Note that the average is done between z = −1.5H 0 and z = 1.5H 0 to exclude the ionized layers where coherent structures are mixed up by the strong MHD turbulence. The top panels show that the disc develops radial density variations associated with zonal flows whose amplitude increases with mean B z (decreases with β). In case β = 10 5 , the density structures are very faint and do not exceed 5% of the mean background density. For β = 10 4 and 10 3 , the zonal flows are much more developed and reach amplitudes respectively 30 − 40% and 170% of the mean background density. Within the simulations time (∼ 100 orbits), these structures remain almost steady, although variations on longer timescale occur (see Bai 2015). Note that for β = 10 4 and 10 5 , the density plots of Fig. 3 reveal the presence of large scale acoustic waves propagating radially in the disc. By computing the Fourier transform of ρ in time and space (along the x direction), we checked that the signal frequencies match the dispersion relation of p-modes (sound waves). These waves might be generated via the turbulent stress, through a process called "aerodynamic noise" (Lighthill 1952;Heinemann & Papaloizou 2009).
The bottom panels show that for β = 10 4 and 10 3 , the vertical magnetic flux is concentrated into very thin shells (or filaments), while outside these shells, the net vertical flux is close to zero. These filaments are located within the gas density minima. Their existence have already been pointed out by Bai (2015) but their origin is still unclear. One possibility suggested by Bai (2014) is that the flux concentration is driven by the anisotropic turbulent diffusivity associated with the MRI turbulence. However ambipolar diffusion seems to reinforce the magnetic shells and might also play an important role. It is likely that the disc stratification is also essential to the formation of the structures. We explore in next section the effect of zonal density and B z structures on the dust sedimentation and radial concentration.

Simulations with dust and ambipolar diffusion
We now study the dust dynamics in the MHD flows described in Section 3. For each β, we simulate the evolution of different grain sizes, corresponding to different midplane Stokes numbers St. We first study the settling of the dust and its evolution toward a quasi-steady state. We then characterize their vertical distributions and typical scaleheight with respect to St and β. We analyse the effect of zonal flows (pressure bumps) on their radial struc-tures and vertical dynamics. Finally, the presence of a mean radial pressure gradient in the disc is investigated.

Initial condition, settling times and convergence
The initial conditions for the dust are similar for all β and Stokes number. The distribution at t = 0 is with H d 0 = 0.35 H and ρ 0 a constant evaluated so that the ratio of surface densities Σ d /Σ is 0.014. The dust velocity is initially Keplerian with zero perturbation. For β = 10 5 , we integrate simultaneously, in a same simulation, the motion of four different grain sizes with Stokes numbers 0.1, 0.025, 0.01 and 0.001. For β = 10 4 , we compute simultaneously six different grain sizes with Stokes numbers 0.1, 0.025, 0.01, 0.0025, 0.001 and 0.00025 ( 5 mm to 10 µm) . Finally, for β = 10 3 , we ran two different simulations, one with four dust species (St = 0.1, 0.025, 0.01 and 0.001) and the second with two dust species (St = 0.0025 and 0.00025). Note that for simplicity the dust mass distribution is initially independent of the particle size, which is far from being the case in real protoplanetary discs. However, since the gas to dust ratio remains small when summed over all sizes, this approximation has little effect on the results. We checked in particular that the back reaction of the dust onto the gas does not produce substantial effects (at least within the simulation time). Moreover, we found that the dust evolution and distribution are independent of the initialization (provided that k ρ d k /ρ 1), and whether the dust components are simulated altogether or separately. Figure 4 shows the space-time diagrams of ρ d (z, t) (horizontally averaged vertical density distribution) obtained from the simulations, for β = 10 4 (left) and β = 10 3 (right), and for different Stokes numbers. The density profiles seem to converge to a steady state depending on the Stokes number after a certain period of time. To quantify more precisely this timescale, we plot in Fig. 5 the time-evolution of ρ d at z = 0 (midplane) and z = 2H, for β = 10 4 . If St 0.001, the initial evolution is purely exponential and dominated by gravitational settling. The dust particles fall toward the midplane within a time proportional to Ω −1 /St, which is characteristic of such process (Dullemond & Dominik 2004). Ultimately, when turbulent diffusion and mixing start to be important, dust particles stop settling toward the midplane and their mean vertical distribution stabilizes. Note that at z = 2H, the gas density is smaller and dust particles have an effective Stokes number greater than in the midplane. Therefore, the settling and convergence toward a steady state are more rapid at this location. Nevertheless, Fig. 5 (bottom) shows that the dust density at z = 2H undergoes brutal variations (especially for St = 0.0025 between t = 400 and 500 Ω −1 ) probably associated with stochastic events, like jets or wind plumes. The case St = 0.00025 (our smallest particle size) is a bit more delicate since the typical settling time can be, in theory, much larger than the simulation time. In this case, an equilibrium is not necessarily reached and the vertical dust distribution keeps evolving slowly in time. Note that for β = 10 3 , the settling rates Article number, page 7 of 22  are similar but equilibria are reached sooner. The reason is that turbulence is more vigorous and the winds more efficient when the magnetization is higher (see Fig.2).

Vertical profiles and dust scale height with respect to St
We now characterise in more detail the shape of the vertical equilibrium and estimate the typical dust scaleheight as a func- The dotted curves are derived from a simple diffusion model assuming that gravitational settling is balanced by homogenous turbulent diffusion (see Dubrulle et al. 1995;, and section 5). The dashed/red horizontal line denotes the size of the numerical cells.
tion of the Stokes number. A quick inspection of Fig. 4 shows that the size of the dust layer increases with decreasing Stokes number. This is physically expected, as small dust particles are less sensitive to gravitational settling and tend to follow the turbulent gas motion. To be more quantitative, we show in Fig. 6 the vertical density profiles, averaged in time, for β = 10 4 and different dust sizes. For comparison we superimposed the gas vertical density profile (dashed blue line) fitting perfectly a gaussian of width H = 1 in the domain |z| < 2H. For large Stokes numbers 0.0025 St 0.1, the dust density profiles fit also quite well with gaussians but with a width much smaller than H. However, for St 0.0025, the profile turns out to be more flared and follows an exponential distribution of the form ρ d = ρ 0 exp(−z/λ d ) with λ d H for the minimum grain size.
To be more quantitative, we define the dust scaleheight H d (St, β) as the altitude z such that In case of large St, this corresponds to the standard height of the Gaussian distribution. For each β and St, we calculate H d from the simulations data and plot their values in Fig. 7. For the largest St = 0.1, almost all dust material is contained within two grid cells around the midplane, so that our estimation of H d is probably incorrect. We remind also that for St = 2.5 × 10 −4 , the distribution is probably not fully converged and one might expect that H d is slightly underestimated.
By examining the scaling relations of Fig. 7 Zhu et al. 2015) and can be understood, in a first and crude approach, within the framework of the classical diffusion theory (Morfill 1985;Dubrulle et al. 1995). In this theory, the equilibrium distribution in the vertical direction is simply achieved by the balance between the gravitational settling and the diffusion due to the turbulence. The latter is assumed to be homogeneous and isotropic, and can be merely described through a uniform anomalous diffusion coefficient. A simple 1D advection-diffusion equation is solved in the vertical direction which gives the relation (25) in the limit of H d H and a flatter dependency in the limit of small St (see Sect 5.1.1). For a comparison with the numerical data, we superimposed in Fig. 7 the predicted scaleheight from this model for the three different magnetizations (see section 5 for details about their calculation).
Although the turbulent diffusion theory seems to reproduce quite well the numerical dependency for β = 10 5 and for β = 10 4 , it does not explain the flared shape of the density profiles at small St (Fig. 6). In the case β = 10 3 , it is even worst; the dust scaleheight is 2 or 3 times larger than the predicted value in the regime of small St ( 0.01). There are actually many caveats to this simple approach. First it assumes that the turbulence is homogeneous and isotropic and neglects the presence of a wind (mean u z ). Indeed, for small St, the typical timescale associated with MHD outflows in the simulations becomes smaller than the stopping time, which suggests that the wind might lift small dust particles and modify their vertical equilibria. For β = 10 3 , the association of a powerful wind with a zonal flow drastically changes the settling process (see next section) and breaks the assumptions usually made in the classical model. All these effects are studied quantitatively and in detail in Section 5.

Effect of zonal flows
As we showed in Section 3, MHD flows dominated by ambipolar diffusion develop quasi-steady zonal density structures (see Fig. 3). What are their impact on the radial distribution of grains and, more indirectly, on their vertical settling?

Dust trapped in rings
First, zonal flows or annular structures are expected to have an effect on the dust radial distribution. Indeed, they correspond to gas pressure maxima in isothermal discs, which are known to collect dust particles and play a fundamental role in the formation of planetesimals. To concentrate dust in the radial direction, the drift associated to the pressure maxima has to be stronger than the radial turbulent diffusion of grains. Since the drift velocity is directly proportional to St (in the limit St 1), we expect that dust will be efficiently concentrated into radial bands for sufficiently large values of the Stokes number. The dependence of the process on β is however less clear. If the radial zonal density structures are less prominent in the low magnetization regime, the turbulent diffusion might be also weaker. Assuming that the dust rings form, another unknown is the stability of these structures. One might ask whether they reach a steady configuration before the dust to gas ratio becomes excessively high and leads to instabilities (e.g. the streaming instability).
To address these questions, we show in Fig. 8 (second, third and fourth rows) the radial density distribution of the dust as a function of time, for St = 0.1, 0.01 and 0.001 and different β. For β = 10 5 (left column), the grains concentrate into very thin radial bands, spaced out by less than a scaleheight H in the radial direction. As expected, the concentration is stronger for the largest particles (high St). For St = 0.1, the density contrast can be quite high, of order 1000, while for St = 0.001, it is never greater than a factor 2. Although the density contrasts are pretty high, the dust to gas ratio remains bounded and smaller than 0.06. For higher magnetization, (β = 10 4 and β = 10 3 ), dust rings also form, with density contrast and concentration much higher than in the case β = 10 5 but with similar dependence on the Stokes number. Moreover, the spacing between the rings is larger than H and comparable to the box size. In particular for the largest magnetization (β = 10 3 ) and St = 0.1, all the dust material accumulates into one single narrow band with averaged dust to gas ratio ∼ 0.3. In this particular case, the stability of the ring is not guaranteed as its density keeps growing in time, even after 300 Ω −1 . The evolution of such structure with high dust concentration is hard to predict from our simulations, since the bi-fluid assumption is not valid anymore when ρ d /ρ 1.
We check that in most cases, the location where the dust is trapped corresponds to a gas pressure maximum (for comparison, we plotted the radial density profiles of the gas in the top panels of Fig. 8). However, for small St, some rings develop outside gas pressure maxima, like the upper one in the right panels of Fig. 8, for St = 0.01 (β = 10 3 ). Actually, they correspond to locations nearby magnetic shells where B z is strong. At these radii, the radial velocity associated with the steady zonal MHD structure is strong enough to advect and concentrate the small dust particles. This occurs because the gas velocity exceeds the radial drift velocity due to the local pressure gradients.

Dust vertical circulation
In the previous paragraph, we showed that zonal flows directly impact the radial dust distribution by efficiently concentrating the grains. But do they also impact vertical settling? First, the stopping time in the density (or pressure) maxima is reduced so that the effect of gravitational settling is weaker. We checked that the dust layer is puffed up in these regions (and shrinked in the density minima by using similar arguments). Note however that this effect, when averaged in x, should not change significantly the mean vertical dust profile. Another and more subtle effect originates from the wind topology and vertical structure of the zonal MHD flow. The top panel of Fig. 9 shows the gas density and poloidal streamlines averaged in time and y, for β = 10 3 . We see that the wind flow is not distributed uniformly along x. Instead, a strong and coherent windy plume emanates from the right part of the box where ρ is minimum and B z is maximum (see Fig. 3 for comparison). This localized structure is maintained during the entire simulation time 300 Ω −1 . It is inclined and connected to a large scale roll in the midplane. Such asymmetric outflow with odd parity is actually reminiscent of global simulations (Fig. 25 of Béthune et al. 2017; and does not appear as an artefact of the shearing box. There is locally, near the midplane, a spontaneous breaking of vertical symmetry of the flow and magnetic fields. Our local model however does not capture the outflow behaviour at large z and its connection to the global disc geometry. Global models predict that at a certain height, the stream that flows toward the star (on one side of the disc) strongly bends and ends up flowing away from the star (see Fig. 25 of Béthune et al. (2017)). Note finally that the rest of the disc is characterized by fainter and smaller scale axisymmetric rolls (or eddies) with smaller correlation time. Some of these features are probably residuals of our averaging procedure, which contains a finite number of datasets (typically 50 in time).
The second and third panels of Fig. 9 show the dust density and mass flux streamlines in the poloidal plane. For intermediate Stokes number (St = 0.01), the dust remains confined in the pressure maxima, out of the windy plume structure. A tiny fraction of the material is advected in the large scale roll but the gravitational settling is too strong for the dust to be lifted away by the plume. The effect of the wind is then partially impeded. However for St = 0.001, more material enters the large scale roll (which also corresponds to gas density minimum) and get lifted by the plume. The flux streamlines of Fig. 9 (bottom) shows that the dust is then dragged horizontally toward the pressure maxima, because of the plume inclination, and finally falls back onto the disc. Therefore a large scale circulation of small dust grains develops between the midplane regions and the disc corona. This circulation, driven by the wind plume associated with the zonal flow, allows grains to reach altitudes higher than those allowed by a simple turbulent diffusion (see Section 5 for a more quantitative evidence of this result).

Simulations with an imposed mean pressure gradient
In protoplanetary discs, gas pressure decreases with radius. In our previous simulations, we neglected the effect of a mean radial pressure gradient, which is known to cause a mean drift of dust particles (Birnstiel et al. 2016), and potentially lead to streaming instabilities (Youdin & Goodman 2005;Youdin & Johansen 2007;Jacquet et al. 2011). The aim of this paragraph is to understand whether the mean drift changes the radial and vertical distribution of solids.
We introduce a dimensionless parameterĝ e that characterises the acceleration induced by the pressure gradient: with R the radius. In the Minimum Mass Solar nebula model (MMSN) of Chiang & Youdin (2010), the aspect ratio is H/R ≈ 2.8 × 10 −2 (R/AU) 2/7 , which gives g e 0.07 at 30 AU. Observations of T-Tauri discs suggest a disc aspect ratio H/R = 0.08 at R = 30 A.U (see Fig. 8 of Pinte et al. 2016). We recall that in this configuration, the mean drift velocity ∆v = v d − v is: where f g = ρ/(ρ + ρ d ) (Youdin & Goodman 2005). The above relation shows that the drift increases with St for St 1 and reaches a maximum for St ≈ 1. Therefore, to maximise the effect of the radial pressure gradient, it is suitable to consider the largest size simulated St = 0.1. We ran two different simulations, one withĝ e = 0 (no pressure gradient) and the other witĥ g e = 0.05. The latter corresponds to a radial drift of 0.005 c s in the midplane. To allow comparison, we integrate a single size (St = 0.1) and use identical initial conditions in both cases. Note that the radial pressure gradient in PLUTO is implemented by adding a force term S drift = ρHΩ 2ĝ e per unit volume in the gas momentum equation, with uniform and fixedĝ e . Figure 10 shows the vertical and radial density profiles for both casesĝ e = 0 andĝ e = 0.05, at a given time t = 220Ω −1 . When a radial drift is present, the dust layer is slightly wider in z, although the difference remains very small. Axisymmetric rings are still formed but the dust concentration inside is weaker. We checked that the gas radial density structures as well as the dust rings remain at the same location as in the the caseĝ e = 0.
In the limit St 1 and for perturbations wavelength λ 2πĝ e St 2 , Jacquet et al. (2011) showed that the pressure gradient induces a streaming instability with growth rate σ: where k x and k z are the radial and vertical perturbations wavenumbers and f p = ρ d /(ρ + ρ d ). This instability is only relevant for Stokes numbers 0.1 as σ depends on the cube of St. However it is not detected during the course of the simulations even for our largest St; the reason is that our resolution is probably too low. Indeed, the dust settles in a very thin layer, with a size not exceeding 2 or 3 grid points. We do not exclude that the streaming instability might occur on very short wavelength ( H) and destabilize the rings, but this requires to increase the number of grid points in z and probably in x as well.

Settling models
The aim of this section is to elaborate a settling model that accounts for the vertical density profiles obtained in Section 4.2.
We start with the historical 1D diffusion model of Dubrulle et al. (1995) and attempt to apply it to our non-ideal MHD flows. We examine in particular the influence of a mean wind component in the model. In a second part, we propose a 2D model particularly adapted to the case β = 10 3 , including the effect of zonal flows and vertical circulation induced by the MHD wind plumes.

1D models
5.1.1. The turbulent diffusion model of Dubrulle et al. (1995) without winds The classical 1D diffusion theory (Dubrulle et al. 1995) assumes that all quantities can be decomposed into a mean component (that depends only on z) and a small fluctuating part: By using such decomposition, the mass conservation equation (5), averaged in x and y, becomes: where ∆v = v d − v is the drift velocity between dust an gas. The first term in the z-derivative corresponds to the advection/stretching of dust by the mean vertical gas flow (wind) that we leave behind here. The second term is due to the correlation of turbulent fluctuations which is reduced to a diffusion operator in Dubrulle's theory. The third term accounts for the mean vertical drift motion of dust due to gravitational settling. By making some assumptions, listed in Section 5.1.2 , it is possible to show that Eq. (31) takes the form of an advection-diffusion equation: (Morfill 1985;Dubrulle et al. 1995;Schräpler & Henning 2004), with D z v 2 z τ corr > 0 the diffusion coefficient and τ corr the correlation time of the turbulent eddies. This equation and the relation liking D z to the vertical r.m.s fluctuations are derived in appendix B. To find equilibrium solutions of this equation, we assume the ergodicity of the system and ∂ρ d /∂t 0. If the gas is in hydrostatic equilibrium ρ = ρ 0 e − z 2 2H , an analytical solution is For uniform diffusion coefficient D z and z H, this gives: The distribution is Gaussian and the dust heightscale has a dependence on St −1/2 for St D z /(ΩH 2 ). It tends toward unity in the limit of small St.

Check of the assumptions
To obtain Eq. (32) and the relation (34), one needs to make several assumptions (see appendix B): 1. Gas density fluctuations small compared to the background density. 2. Vertical drift ∆u z = v z d − v z proportional to z St/ρ. This is called the "terminal velocity approximation". 3. Wind advection u z ρ d negligible compared to the turbulent transport δρ d δv z d . 4. Turbulent transport in the form of a diffusion operator −D z ρ ∂ ∂z ρ d ρ with D z uniform diffusion coefficient (homogeneous and isotropic turbulence) The aim of this paragraph is to discuss the validity of these assumptions in the different regimes probed by our simulations.
The first assumption (1) is verified for β = 10 5 and β = 10 4 . Indeed, in both cases, the density fluctuations never reach more than 30% of the background density. However it is not valid for β = 10 3 for which the density contrast can be higher than 3. To check assumptions (2) and (3), we plot in Fig. 11 the three different flux terms appearing in Eq. (31) for β = 10 4 . The turbulent flux δρ d δv z d is computed directly from the 3D numerical data, with perturbations δρ, δu obtained by subtracting to each quantity its mean value. We superimposed in the same figures the approximated settling term zStρ d /ρ. We show that there is an excellent agreement between the actual vertical drift and its value in the terminal velocity approximation. A tiny difference is however obtained at z 2H for St = 2.5 × 10 −4 . Regarding hypothesis (3), Fig. 11 shows that the flux associated with wind advection (blue/plain lines) remains small compared to gravitational settling and turbulent transport, provided that St 0.0025,. However for St < 0.001, it is comparable to other terms and carries a large fraction of the dust vertically. By neglecting the wind component, the 1D model of Dubrulle et al. (1995) is probably incomplete for small particles with St 0.001. Finally, to check assumption (4), we plot in Fig. 12 the ratio for different β and St. The diffusion coefficient D z is computed by averaging the numerical 3D datasets in time with a sampling period of 5Ω −1 . For β = 10 5 and β = 10 4 (top panel), D z is positive and well-defined. In the midplane, within z ± 0.5H, we found that D z 0 1.9 × 10 −4 for β = 10 5 and D z 0 1.6 × 10 −4 for β = 10 4 . However D z is not uniform in z and is multiplied almost by a factor 10 at z = 1.5 H. The diffusion coefficient increases rapidly with altitude, fitting quite-well with an exponential function D z (z) D z 0 exp(1.6 |z|/H). These results are consistent with the theoretical prediction that D z v 2 z τ corr with τ corr 0.5Ω −1 . In particular, we showed in Section 3.1 that the vertical r.m.s fluctuations increase quasi-exponentially with altitude and have similar values for both β = 10 5 and β = 10 4 . Note that the diffusion coefficient does not depend too much on the particle size.
To summarize, as long as the size of the dust layer remains much smaller than H, the approximation of uniform diffusion coefficient holds and leads to the simple solution given in (34). However for St 0.001, this assumption does not hold; winds and non-uniform D z have to be included in the model to obtain acceptable solutions. For β = 10 3 (bottom panel of Fig. 12), we found however that the 1D diffusion coefficient is not welldefined, as it can be either positive or negative, and undergoes brutal variation in z. We checked in particular that the flux associated with the turbulent fluctuations is negative for z > 0, carrying mass downward. This result suggests that the simple 1D diffusion theory is inappropriate to describe dust settling in this regime.

Numerical vs theoretical density profiles
We focus here particularly on simulations with β = 10 4 (the case β = 10 3 requires to go beyond the 1D model and is treated further in Section 5.2). Figure 13 shows a comparison between the numerical density profiles in z and the theoretical prediction of the 1D model, in case of a uniform diffusion coefficient D z = D z 0 = 1.6 × 10 −4 (blue lines) and a varying D z = D z 0 exp(|z|/h i ) (orange lines) with h i = 0.6 andz = z/H. The theoretical profile in case of a non-uniform D z is obtained analytically by integrating (33): with Article number, page 13 of 22 The solution with uniform diffusion coefficient (33) is recovered by setting h i = ∞ and expanding the exponential in ξ(z) to first order.
For St = 0.01 D z 0 /(ΩH 2 ), Fig. 13 (top) shows that the analytical solution fits closely the numerical profile. As expected, the difference between uniform D z and varying D z solutions is very small. However, for St = 2.5 × 10 −4 D z 0 /(ΩH 2 ), both solutions fail to reproduce the numerical vertical distributions. If the heightscale calculated from the theoretical solution is not too far from the simulated one, the shape of the profile differs significantly. In particular, the 1D model does not account for the flaring at large z.
In conclusion, we showed that the 1D diffusion model of Dubrulle et al. (1995) applied to a non-ideal MHD turbulent flow works pretty well for β = 10 5 and β = 10 4 , provided that the dust particles are not too small. For St 0.001, the model does not reproduce very well the shape of the vertical profile, even when the assumption of uniform diffusion in z is relaxed. The wind starts to be important in this regime and has to be included in the model. This is investigated in the next sections.

1D models with a uniform and mean wind
We showed that for small particles (typically St < 0.001), the advection of the mean wind ρ d v z in Eq. (31) becomes comparable to the gravitational settling. Then, assumption (3) needs to be relaxed. A first and naive approach is to model the wind by a linear function of z but uniform in x: with a characteristic rate Ω w . We checked that this simple dependence matches relatively well with the simulations data, at least for z within ±H. We found respectively Ω w = 1.6 × 10 −4 , 7 × 10 −4 and 6 × 10 −3 Ω for β = 10 5 , 10 4 and 10 3 . The modified equilibrium, including the effect of the wind, is: Assuming a uniform D z = D z 0 , it is straightforward to show that for the dust density increases with altitude in the vicinity of z = 0. In case β = 10 4 with Ω w = 7 × 10 −4 , this occurs for a critical St = 5.4 × 10 −4 . We plot in Fig. 14 the solutions of Eq. (39) for different Stokes numbers and β = 10 4 . For St = 2.5 × 10 −4 , the density profile exhibits two bumps at z 1.5H where the dust settles. This configuration corresponds to "floating" or "levitating" dust layers. Particles are trapped between z H where the wind dominates and regions of higher altitudes where gravitational settling dominates (the latter is an exponential function of z while the wind advection is only proportional to z). We checked that assigning a z-dependence of D z leads to similar results. These floating layers are actually reminiscent of the simulations by Miyake et al. (2016), who found that magnetorotational driven winds can concentrate dust grains with size 20 − 40 µm around four disc scaleheights by a similar mechanism.
In conclusion, the simple 1D wind model does not reproduce at all the profiles computed numerically. Indeed, none of them exhibit such distribution with floating particles. The reason behind that is intimately related to the baseline assumptions of the 1D model and particularly to the decomposition (41). While the gas density fluctuations remain small compared to the background density, it is not necessarily true for the gas vertical velocity (see Fig. 2 for example). The term "fluctuations" here is actually misleading since it does include both random turbulent motions and the x-dependent self-organized and coherent axisymmetric flow induced by ambipolar diffusion. In particular, we showed in Section 4.3.2 that the wind is distributed non-uniformly in x and rather active in the gas density minima, where the dust concentration is low. By averaging quantities in x, we thus overestimate the wind effect. Moreover, the presence of zonal flows in x leads to a large scale meridional circulation (see Fig. 9) which is not taken into account in the 1D model. To better understand the dust distribution in magnetized discs with zonal flows (typically for β = 10 3 and marginally for β = 10 4 ), one needs to treat the advection/diffusion equation in 2D. This is the object of the next section.

2D models with MHD zonal flows and wind plumes
In this section, we develop a 2D model that takes into account the wind geometry (in radius) and the zonal flow structures. We focus particularly on the case β = 10 3 for which the settling process cannot be described by a 1D turbulent diffusion model (see 5.1.2), Note that the 2D model described below applies also for the case β = 10 4 in the regime of small dust particles.

Formalism and toy model
To begin, we assume that the dust and gas motion is the sum of a steady axisymmetric background flow (physically representing the zonal structure) and turbulent fluctuations δv: where X a (x, z) = X dy denotes the average in time and y of a quantity X. Like in the 1D theory, we model the turbulent Fig. 15. Gas density (colormap) and velocity streamlines (cyan arrows) of the 2D toy model. The lines thickness is proportional to the velocity v . The outflow inclination, width and averaged intensity are fixed and equal respectively to i = 40 • , δ = 0.7H, and Ω w = 6 × 10 −3 Ω. The density is Gaussian in z and has a cosine perturbation in x with amplitude A x = 0.5. The outflow configuration can be directly compared with that of Fig. 9 (here the plume has been shifted toward the center of the 2D domain) transport (cross-correlation of density and velocity fluctuations) through a diffusion operator: with D x and D z the radial and vertical diffusion coefficients. To simplify notation we drop the brackets and the overline denoting the time and y averages. The 2D steady axisymmetric dust distribution is then governed by a second order partial differential equation: Before solving this equation for ρ d , it is necessary to prescribe the background axisymmetric gas density and velocities. To model the zonal flow, we assume the following form for the gas density: with A x the amplitude of the radial density pattern. Instead of a uniform wind v z = Ω w z prescribed in 5.1.4, we consider here a x-dependent steady axisymmetric outflow modelled by: where x i = x cos i + z sin i. This simple prescription mimics the inclined outflow obtained for β = 10 3 and illustrated in the top panel of Fig 9. The inclination i refers to the angle relative to the vertical axis. The plume is confined within the gas density minimum, with a Guaussian profile in x and a given width δ.
The infinite sum ensures the strict periodicity of the flow in the radial direction (with period L x ). Figure 15 shows the gas density and the streamlines associated with the 2D outflow for a given set of parameters i, δ, Ω w and A x . Given such prescription for the gas motion, it is possible to calculate the dust velocities. To simplify the problem as much as possible we neglect the radial drift induced by the pressure maximum, so that v x d = v x . This is correct in the limit of small Stokes numbers. In z, we use the terminal velocity approximation (see Appendix B) to estimate the vertical drift:

Parameters of the 2D model
For a given Stokes number St, the 2D model described in 5.2.1 contains exactly six free parameters, three associated with the geometry of the outflow (i, δ and Ω w ), two associated with the turbulent diffusivities (D x and D z ) and one related to the amplitude of the zonal flow. To match with the simulations, we choose i = 40 • , δ = 0.7H and Ω w = 6 × 10 −3 L x /(δ √ 2π). The latter is calculated such that the horizontal average v z (z) matches with the 1D prescription of 5.1.4. The amplitude of the zonal flow A x is fixed to 0.5. The two diffusion coefficients remain to be estimated from simulations. We remind that their definition in the 2D formalism are:  (Dubrulle et al. 1995) respectively without wind and with uniform wind. All profiles are normalized to unity in the midplane to allow comparison.
The numerical calculation for β = 10 3 gives D x 2.5 × 10 −4 and D z 3×10 −4 in the midplane regions within z±1H. However we point out that these coefficients are difficult to evaluate further in the disc atmosphere.

2D solutions and comparison with simulations
Once the dust velocities and parameters of the model are fixed, we compute numerically the 2D equilibrium solutions of (44).
The numerical methods to obtain these solutions are detailed in Appendix C. Figure 16 shows two examples of solutions computed for St = 0.001 and St = 0.00025. The top panels correspond to a uniform wind in x with δ = ∞ while the bottom panels correspond to a wind distributed non-uniformly in x with δ = 0.7H. Both have non-zero inclination i = 40 • . In the first case, the dust forms two layers off the midplane, located around z = 2H, exactly as in the 1D case (see red line in Fig. 14 for comparison). Note that the dust distribution have a radial dependence due to the gas density modulation in x. In the second case however, the dust remains confined in the midplane. The strong collimated wind associated with the plume locally depletes the dust and forces some material to float around z = 2H, but most of the grains fall back into the disc once they leave the plume region. A large scale poloidal circulation, similar to that described in Section 4.3.2 is then sustained in the disc atmopshere. One critical parameter that controls the rate of mass entering the plume and therefore the amount of dust accumulating at large z is the radial diffusion coefficient. Indeed strong radial diffusion will allow a rapid replenishment of the depleted region (plume) and thus massive dust outflow. On contrary, small D x will tend to reduce the mass loading. The simulations are rather in a regime of inefficient wind since the diffusion timescale H 2 /D x ≈ 4000 Ω −1 is much longer than the outflow timescale Ω −1 w ≈ 200 Ω −1 . Note that in the simulations (Fig. 9) the outflow region was not that depleted because a large scale poloidal roll was trapping the dust inside. This roll has not been included in our model as its effect is probably not dominant in the vertical settling process.
To check that our 2D model is compatible with simulations, we compare in Fig. 17 the horizontally averaged density profiles ρ d of the 2D solutions (red dashed lines and orange lines) with those obtained in the shearing box (green lines). For both Stokes numbers St = 0.001 and St = 0.00025 our model is able to accurately reproduce the shape of the vertical profiles. For St = 0.00025, the result is independent on whether the diffusion coefficients are uniform or not. However this is not necessarily true for St = 0.001 for which the 2D model with uniform D z fails to reproduce the simulations data. We also superimpose in Fig. 17 the profiles obtained with the 1D Dubrulle's model (blue lines) and the 1D model with wind (purple dotted line). Clearly the 2D models are better to explain the simulations. This indicates that the dust poloidal circulation induced by the zonal flows is crucial to model the dust vertical equilibrium.

Conclusions on settling models
In the regime of β 10 4 and for St 0.001, we found that the 1D model of Dubrulle et al. (1995) correctly predicts the shape of the density profile and the ratio H d /H. However when the magnetization is stronger or the particles size smaller, the 1D model fails to reproduce the simulations profiles. We found that the wind dynamics and its dependence on x has a crucial effect on dust distribution and needs to be modelled in a 2D framework. The dependence of diffusion coefficients on z seems also to matter but does not fundamentally change the physical processes involved.

Comparison with other MHD simulations
Our results suggest that the dust is highly sedimented in presence of ambipolar diffusion, regardless the magnetization. In particular for particles size larger than the millimetre, we found that the ratio H d /H is smaller than 0.1 and the vertical diffusion coefficient less than 3 × 10 −4 . This corresponds to values of the Schmidt number S c = αΩH 2 /D z between 10 for β = 10 5 and .01, non-ideal simulation β = 10 3 St=0.01, ideal MRI simulation  ALMA constraints by Pinte et al. 2016 (0.87-2.9 mm) Fig. 18. Vertical density profile of the millimetre dust simulated at 30 AU (St 0.025) and comparison with ALMA observations. The red profile corresponds to our simulation with β = 10 3 and ambipolar diffusion (Am = 1). The purple profile is that obtained in the ideal (zero net flux) MRI simulation of  for the same Stokes number. The orange area corresponds to the range of profiles compatible with the 0.87-2.9 mm dust continuum emission measured by ALMA (Pinte et al. 2016). To help the comparison, the density profiles are all normalized with ρ 0 = 1; the gas profile is represented by a blue line. The x-axis unit is in AU, the conversion is about 1H 2.5 AU at R = 30 AU.
200 for β = 10 3 . Hence, the gas transport is always much larger than the vertical dust transport. This in contrast with ideal MRI turbulence simulations (either zero or non-zero net vertical flux) for which typical Schmidt number is never larger than 10 (Johansen et al. 2006;. We can interpret this result as a consequence of enhanced Maxwell to Reynolds stress ratio at large B z (see Section 3.1). In fact, in non-ideal MHD flows, angular momentum is mainly extracted through coronal winds and laminar Mawell stress, while dust is more sensitive to turbulent gas motions. Note that the large S c found in our simulations, as well as the the small H d /H and τ corr Ω 0.5, seem incompatible with the recent non-ideal MHD simulations of Zhu et al. (2015) with ambipolar diffusion. Indeed in their case, they found S c 1 and diffusion coefficients 10 to 20 times larger than ours. The main difference between their simulations and ours is the vertical and radial box size (they used L x = 4H, L z = 6H. while we use a larger domain with L x = 8H, L z = 12H). We ran few simulations by varying both L x and L z and found that the discrepancy is mainly due to the vertical extent L z . The reason is that the wind properties, in particular its strength and mass loss rate, depend strongly on the the vertical extent of the box . We checked that for a box of size L z = 6H, the vertical r.m.s turbulent fluctuations are stronger by a factor almost 2 and the wind by a factor 3-4, compared to the case L z = 12H. This leads to a much higher diffusion coefficient and a ratio H d /H three time larger. We emphasize that the convergence of the outflow properties with L z is probably not reached even for L z = 12H. The ratio H d /H may therefore be slightly over-estimated by our simulations. To avoid such bias, it is necessary to model the wind in a global way and characterize the dust settling in global simulations, which will be the object of a future paper.

Vertical settling and comparison with ALMA observations
We compare here the dust scaleheight H d measured in our simulations with that inferred from observations. We discuss also about the limitations of such comparison. Most of the constraints available today on H d come from the sub-millimetre ALMA observations of T-Tauri discs like HL-Tau (Pinte et al. 2016). These constraints, based on the measure of the rings contrast, suggest that sub-millimetre particles are contained in a geometrically thin disc with scale height between 0.7 and 2 AU at 100 AU. If we assume that the ratio H d /H does not vary strongly with R and consider disc aspect ratio H/R 0.1 at R = 100AU (suggested by Pinte et al. 2016, for HL Tau), we find H d /H 0.07 − 0.2.
Using either the MMSN model or current estimation of disc surface density profiles, it is possible to show that sub-millimetre dust particles measured by ALMA lies into the range St = 0.005 − 0.01 at 30 AU (see Section 2.3). To obtain this range, we have considered that most of the emission in the band 6-7 of ALMA (with wavelength λ 1 mm) is due to particles of size a λ/(2π) 160 µm (Kataoka et al. 2017). For such St, our simulations with ambipolar diffusion indicate that the dust is contained within H d 0.11 − 0.19H, depending on the magnetization of the disc. This scaleheight is then compatible with observation. To facilitate the comparison, we superimpose in Fig. 18 the normalized density profile in z obtained numerically at β = 10 3 , St = 0.01 (red line) and a surface (in orange) covering the range of Gaussian profiles estimated by the observations of Pinte et al. (2016). Note that the profiles obtained for smaller β lies only marginally within the error bar. We stress that the scaleheights considered here are 4 to 10 times smaller than those found in ideal MRI simulations . In particular, the typical vertical diffusion coefficient in this case ( 5.5 × 10 −3 ) is 20 times larger than those obtained in our simulations with ambipolar diffusion. All these results are indications that the dynamics of discs like HL-Tau, in regions beyond 1 AU, is strongly affected by non-ideal MHD processes such as ambipolar diffusion.
There are however many caveats associated with such comparison. First the uncertainties on the gas surface density in the outer regions can be more important than those estimated. The current MMSN models and the observations can differ by one order of magnitude at most. We have also assumed that the dust to gas scaleheight ratio H d /H does not vary with radius, which is a crude approximation. In reality, it may encounter large variations between 30 AU and 100 AU, either due to radial variation of turbulence strength or local change in the gas to dust ratio (Pinte et al. 2016). The estimation of H d from the rings contrast is also largely arguable, since it depends on the dust opacities and radiative transfer used in the synthesized model of the disc. High resolution imaging of edge-on discs is probably a better and promising avenue to measure directly the dust scaleheight. Finally, although disc emission at λ = 1mm peaks around a 150 µm, there is a wide range of particles size from 50 µm to 500 µm that also contributes, to a lesser extent, to such emission. Since particles are not distributed uniformly in size in real discs, there is potentially a bias in our estimation of the Stokes number at 30 AU. This last issue however should not change drastically the conclusions of the present work.
Note that ongoing ALMA observations of edge-on discs (in particular in HH30,Menard et al.,in prep), combined with previous HST surveys, are going to provide accurate measure of both millimetre and micrometer dust thickness. We achieved a preliminary confrontation between our model and these measures and found that the vertical extent of both dust components fit with the observations when using β = 10 3 . A forthcoming paper will be dedicated to the details of such comparison.

Rings and comparison with ALMA observations
The Table 1 of Pinte et al. (2016) summarizes the properties of the rings/gaps in HL Tau estimated from the millimetre emission. In particular information about the gap width and depth is given. The width associated with the gap at 30 AU is about 10 AU, although it varies moderately from one gap to another. In our simulations, the typical width of the sub-millimter dust rings (St 0.01), averaged over the different structures, is respectively 0.4H, 1.7H and 3.5H for β = 10 5 , 10 4 and 10 3 . Given that H 2.5 AU in our case, there is again a good agreement between the numerical result and the observation for β = 10 3 . We stress however that the boundary conditions and the lack of global radial structure in shearing-boxes might influence the distance between the rings. Note also that other discs than HL-Tau, like Elias 24, exhibit much wider rings or gaps with size of few tenth of AU. To explain such systems, our simulations are probably not relevant anymore ; the presence of an embedded planet with a mass comparable to that of Jupiter might be required to form such gaps (Dipierro et al. 2018).
As for the gap depth, which measures the ratio between the average dust density in the ring and the density at the minimum, direct comparison is much less obvious and has to be taken with caution. Indeed, due to the instrumental resolution ( 3 − 5 AU for HL-Tau) the contrast between the rings cannot be measured accurately. At best, a lower bound for the gaps depth can be estimated, which is around 20 for HL-Tau. In the simulations, this quantity is also extremely difficult to estimate since it drastically changes between the different gaps. In the case β = 10 3 , the density contrast in the region lying between the two rings (see Fig. 8 and St = 0.01) is 35. However in the wider region across the radial boundary, this ratio is huge, of order 10 6 , since the dust is completely depleted there. Morevoer, we remind that the dust density within the rings has not necessarily reached a steady value in simulations.

Conclusions
In summary, we characterised the dust dynamics in the outer regions of magnetised accretion disks with ambipolar diffusion. By using 3D MHD stratified shearing box simulations, we showed that centimetre to millimetre dust particles are highly sedimented in the midplane with scaleheight H d /H 0.1 − 0.2. This result fits particularly well the observational constraint of Pinte et al. (2016) (ALMA) in case of a large enough magnetization β 10 3 . The vertical settling associated with particles of intermediate size (mm-cm) can be described through the 1D diffusion model of Dubrulle et al. (1995) which nicely reproduces the density profiles and the ratio H d /H. We also found that the dust particles in this range accumulate within the pressure maxima and form large scale axisymmetric structures (rings in a global view) whose properties depend strongly on the magnetization. For β 10 3 , the rings separation (or gaps width) is of order 3-4 H, which is in agreement with measures in HL-Tau (Pinte et al. 2016). We checked that the presence of a mean pressure gradient in the box does not change drastically these results, although our simulations are not enough resolved to trigger the streaming instability and thus determine its impact on the dust distribution.
For small Stokes number (sub-millimetre to micron size particles), we showed on contrary that the dust is less sedimented in the vertical direction, with H d /H ≈ 1 and forms much fainter rings. In case of a strong magnetization (β = 10 3 ), particles are subject to a large scale poloidal circulation driven by the zonal flows and the winds associated with the ambipolar MHD flow. Given their small size, particles in the midplane can diffuse toward density minima where the magnetic activity is high. Once they enter such regions, they are lifted by powerful gas plumes and reach the base of the ionized layer. Due to gravitational settling and inclination of the plumes, they fall back onto the disc. Such circulation makes the 1D diffusion model of Dubrulle et al. (1995) inappropriate to explain the dust settling in this regime. A 2D model taking into account the radial structure of the gas flow is necessary to reproduce the vertical dust profiles of the simulations. One crucial element of this model is the non-uniformity of the wind in the radial direction. Indeed, models assuming a uniform outflow (inclined or not) lead to "floating" dust layers which are not seen in the present simulations.
We discuss finally the implications of this work on accretion and planet formation scenario in protoplanetary discs. Since they are cold and dense, such objects are inevitably subject to nonideal MHD effects when threaded by a magnetic field. Our simulations suggest that both accretion efficiency (up to α 10 −2 ) and dust scaleheight in discs dominated by ambipolar diffusion are compatible with observations. The turbulent r.m.s fluctuations obtained in these simulations are also in agreement with the dispersion velocity in CO measured recently by Flaherty et al. (2017) who predict an upper limit v turb < 0.05c s in the midplane. Non-ideal MHD effects are therefore a possible avenue to reconcile the accretion theory with the observations. A next step is to perform global simulations, including a more realistic dust size distribution and the three non-ideal effects. In particular the role of the Hall effect (Lesur et al. 2014;Bai 2014) on the dust is still unknown and has to be investigated.
This work is also directly relevant for planet formation theory and may be important to understand the growth of grains to pebbles ("radial-drift" barrier) and from planetesimals to planets ("fragmentation" barrier) (Birnstiel et al. 2016). Indeed, the confinement of solids in pressure maxima (here the rings) is known to be a very efficient way to overcome these barriers. In addition to stop their radial migration and accelerate their growth, such trapping regions also reduce their relative velocities, which prevents grain fragmentation (Gonzalez et al. 2017). The high dust to gas ratio in these regions can ultimately lead to gravitational runaway of the solids. Non-ideal MHD effects appear to be an alternative to other trapping mechanisms such as vortices, whose stability and origin remain unclear (Lesur & Papaloizou 2009). It may be also more generic and efficient at collecting dust particles than the streaming instability (and related self-induced traps) whose excitation requires very specific conditions. The high dust to gas ratio in the midplane due to the strong sedimentation induced by MHD zonal flows could also enhance the so-called "pebble accretion" (Ormel & Klahr 2010;Lambrechts & Johansen 2012), and therefore reduce the growth timescale of massive cores.
Decay rates of 1D sound waves in a gas/dust mixture as a function of the stopping time. The dashed/red curve is the analytical solution whereas the diamond markers are the rates measured from PLUTO simulations for different resolution in x. The background dust to gas ratio is χ = 1 and the wave number is k x = 2π/L x .

Appendix A: Tests of dust implementation in PLUTO
Appendix A.1: Numerical methods The numerical integration of a given dust specie is similar to that of the gas, except that the pressure is not computed. We use a HLL Riemann solver to compute the density and momentum flux at cells interfaces. A standard monotonized central flux limiter is used to bound the spatial derivatives. The drag force is treated as a source term in the right hand side of the Runge Kutta solver. The time step is adapted to take into account the large drag force encountered by small particles and ensure the stability of the scheme. The upper bound for the time step is taken from Laibe & Price (2012) (cf their section 3.2): where the subscript k labels the different particle species. This prescription, derived from a Von Neumann analysis, guarantees the stability of the simple explicit Euler-scheme. It is then robust but perhaps not optimal for Runge-Kutta solvers. Finally, in order to avoid strong numerical diffusion, we implement a version of the FARGO algorithm for the dust components, which treats the shear advection separetely from the other flux and source terms.

Appendix A.2: Sound waves
To check our implementation, we first perform a series of tests involving the propagation of simple axisymmetric sound waves in a dust-gas mixture along the radial direction. As a first step, we neglect the Coriolis force (no rotation) and the shear. Assuming a single dust specie of size a and perturbations δρ, δρ d , δv x and δv x d ∝ exp(ik x − ωt), the equations governing the linear sound where ρ, ρ d and χ are respectively the gas density, dust density and dust to gas ratio related to the background. By calculating the determinant of this system, and eliminating the trivial mode ω = 0 (gas components reduced to 0), we obtain the following dispersion relation: The solution is composed of one standing wave and two oscillatory modes with frequencies close to k x c s (in the limit of small τ s ). Because the dust and gas exert a mutual drag, the oscillatory modes decay in time with a rate that depends on the dust to gas ratio and τ s . In Fig. A.1, we plotted in red/dashed line the theoretical decay rate as function of τ s , calculated analytically from the dispersion relation, for a fixed dust to gas ratio χ = 1.
To check that our PLUTO implementation is correct, we simulate a linear sound wave for different τ s between 0.1 and 1000, and for χ = 1. The wave corresponds initially to a pure eigenmode of the linear system with an amplitude of 10 −3 . We choose k x = 2π so that the wave fills the box of size L x = 1. The diamond points in Fig. A.1 are the decay rates obtained in the simulations for different τ s and numerical resolution (per wavelength). Clearly, there is a perfect agreement with the teoretical prediction, provided that N x 64. However, the difference remains small even for N x 32. Below N x 16, the decay rate is not reproduced accurately and numerical dissipation starts to be prominent. Note that the wave never dissipates entirely and its kinetic energy saturates around a very low value, regardless of the resolution. This behaviour results probably from the initial numerical noise that produces spurious modes and pollutes the initial wave. However we emphasize that such noise is completely negligible compared to the initial wave amplitude.
Note that this simple problem has been already tested in SPH codes (Laibe & Price 2011. Unlike SPH implementation, we did not experience any convergence problem when consider-ing stopping times much smaller than the wave period. In particular the rate at which it converges with resolution seems independent on the stopping time. Laibe & Price (2012) argued that at small τ s the spatial de-phasing between gas and dust is small, of order c s τ s , and must be resolved numerically in order to capture properly the physics of the wave dissipation. Our tests in PLUTO suggest however that the dephasing length does not need to be resolved to obtain correct dissipation. Such length is actually purely geometrical and is not associated with any energy or momentum transfer.

Appendix A.3: Shearing waves with rotation
We performed similar tests for 2D non-axisymmetric waves of the form exp(ik x + ik y − ωt) with k y 0, in the presence of rotation and shear. Since these waves have a wavenumber k x = k x 0 + S k y t that increases linearly in time, analytical solutions are not straightforward to obtain. Waves are rapidly sheared out and their evolution on long time scales cannot be described through simple decaying trigonometric functions. To deal with this issue, we solve numerically the linearised problem with a simple Runge Kutta integrator and compare the solutions with those obtained numerically with PLUTO.
We consider an initial wave with k x = −2π/L x and k y = 2π/L y in a box of size L x = 1 and L y = 3. The amplitudes of the density and velocity fields of the mode are initialized randomly and uniformly between 0 and 10 −3 . We tested three different stopping time by fixing the dust to gas ratio χ = 0.4. The equation of state is isothermal and the gas is unmagnetized. The results are shown in Fig. A.2, where the blue curves represent the time-evolution of kinetic energy integrated with our linear solver. The other curves correspond to PLUTO simulations with different resolutions. Our code reproduces quite well the desired solution during the first shearing times but the gas component dissipates at longer times (t > 10 Ω −1 ) as the wave is strongly sheared out and damped by numerical diffusion. We checked in particular that increasing the numerical resolution makes the evolution of the waves closer to the theoretical prediction. Note that for τ s Ω = 0.2, the dust component tends to follow the gas oscillations within the first orbits, since τ s is smaller than the oscillation period. Dust and gas are well-coupled and the wave is rapidly damped due to the strong drag between the two components. However for τ s = 20Ω −1 , there is only a weak coupling between gas and dust, since τ s is much larger than the gas oscillation period. The dust motion follows its own oscillation with larger period and slowly decays due to the small drag.

Appendix A.4: Streaming instability
Finally we tested the streaming instability in the 3D unstratified box with a mean pressure gradient. The properties of this instability and its linear derivation can be found in Youdin & Goodman (2005); Jacquet et al. (2011). The setup is similar to that described in Section 4.4, but withĝ e = 0.02. The background flow has a dust to gas ratio χ = 0.2 and a mean radial drift given by Eq. (28). We first checked the 2D linear unstable axisymmetric modes in x and z. The theoretical growth rates and eigenfunctions of these modes are computed semi-analytically by solving the linearised system given in Youdin & Goodman (2005).   We ran PLUTO simulations for different stopping times (0.025, 0.1, 0.2, 0.5, 2, 10, 50, 200), with initial states composed of a background drift flow plus an unstable eigenmode of amplitude 2 × 10 −3 . The box size is L x = 1H, L z = 4πH and the numerical resolution is N X = 128, N Z = 64. The initial mode has k x = 2π/L x and k z = 2π/L z . The red diamond markers in Fig. A.3 correspond to the numerical growth rates obtained with PLUTO. We checked that for τ s 100 there is a very good agreement between the theoretical and numerical growth rates, with a relative error less than 2% in all cases. For τ s = 200 Ω −1 , however we do not reproduce the correct growth rate. Actually we found that this discrepancy is not due to our implementation but to the presence of unstable harmonics with much higher growth rates, probably excited by the initial numerical noise. To make sure that our implementation is robust, we also simulate the nonlinear evolution of such mode in a bigger box (L x = L z = 6H) with higher resolution (N X = 512, N Z = 256) during 1000 orbits. We recover the classical result of Youdin & Johansen (2007) that dust concentrates into small-scale clumps with local dust to gas ratio greater than 1.