Issue 
A&A
Volume 617, September 2018



Article Number  A117  
Number of page(s)  22  
Section  Astrophysical processes  
DOI  https://doi.org/10.1051/00046361/201833212  
Published online  28 September 2018 
Dust settling and rings in the outer regions of protoplanetary discs subject to ambipolar diffusion
Univ. GrenobleAlpes, CNRS, Institut de Planétologie et d’Astrophysique de Grenoble (IPAG), 38000 Grenoble, France
email: antoine.riols@univgrenoblealpes.fr
Received:
11
April
2018
Accepted:
1
May
2018
Context. Magnetohydrodynamic (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 in discs, the magnetorotational instability (MRI) has raised particular interest since it provides a source of turbulence and potentially organizes the flow into large scale structures. However, the weak ionization of discs prevents the MRI from being excited beyond 1 AU. Moreover, the low velocity dispersion observed in CO and strong sedimentation of millimetre dust measured in TTauri discs are in contradiction with predictions based on ideal MRI turbulence.
Aims. In this paper, we study the effects of nonideal MHD and magnetized 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).
Methods. To compute the dust and gas motions, we performed numerical MHD simulations in the stratified shearing box, using a modified version of the PLUTO code. We explored different grain sizes from micrometre to few centimetres and different disc vertical magnetizations with plasma beta ranging from 10^{3} to 10^{5}.
Results. Our simulations show that the mmcm dust is contained vertically in a very thin layer, with typical heightscale ≲0.4 AU at R = 30 AU, compatible with recent ALMA observations. Horizontally, the grains are trapped within the pressure maxima (or zonal flows) induced by ambipolar diffusion, leading to the formation of dust rings. For micrometre grains and strong magnetization, we find that the dust layer has a size comparable to the disc heightscale H. In this regime, dust 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.
Conclusions. Our results suggest that nonideal MHD effects and MHD winds associated with zonal flows play a major role in shaping the radial and vertical distribution of dust in protoplanetary discs. Leading to effective accretion efficiency α ≃ 10^{−3}–10^{−1}, nonideal MHD models are also a promising avenue to reconcile the low turbulent activity measured in discs with their relatively high accretion rates.
Key words: accretion, accretion disks / protoplanetary disks / magnetohydrodynamics (MHD) / turbulence / planets and satellites: formation
© ESO 2018
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
1. Introduction
One of the most challenging problems 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 magnetorotational 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, 2013; Bai & Stone 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, nonideal MHD effects, such as ohmic diffusion, Hall effect and ambipolar diffusion prevail and tend to suppress any form of MRIdriven turbulence (Fleming et al. 2000; Sano & Stone 2002; Wardle & Salmeron 2012; Bai 2013, 2015; Lesur et al. 2014).
Observationally, the prevalence of a sustained and vigorous MRIturbulence in the regions ≳1 AU is also largely questioned. The high resolution imaging at optical, nearinfrared (VLT/SPHERE) and submillimetric wavelenght (ALMA) have provided precious information about the disc structure and turbulence, essentially by probing the spatial distribution of micrometre and millimetre dust grains. Few observations have been able to map the dust heightscale, by measuring the rings contrast (Pinte et al. 2016) or directly by imaging edgeon discs (Wolff et al. 2017). In particular, the survey of Pinte et al. (2016) in HL Tau suggests that the millimetre 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 one or two orders of magnitude smaller 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.
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). Although 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 dustdriftdriven 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 nonideal 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 nonideal effects lead in the nonlinear regime to the formation of selforganized 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 ringlike structures remains an open question. Since current submillimetre or nearinfrared observations are sensitive to dust, direct comparison with them requires to go one step forward and determine the dust dynamics in such nonideal 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; Fromang & Papaloizou 2006; Balsara et al. 2009; Zhu et al. 2015), it remains relatively unexplored in nonideal MHD simulations characterizing the disc outer regions (r ≳ 0.1–1 AU). Ruge et al. (2016) studied the case of zero net flux MRI with ohmic resistivity but did not include other relevant nonideal effects. Several simulations by Zhu et al. (2015) and Xu et al. (2017) 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 nonideal MHD effect to concentrate the grains into rings. The goal is to make comparison with recent ALMA observations for millimetre 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 nonideal MHD flows characterized 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 pressureless fluid made of different particle sizes. We explore different grain sizes from few micrometres 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 Sect. 2, we describe the model and review the main characteristics of nonideal physics and dustgas interaction. We also present the numerical methods used to simulate the dustgas dynamics. In Sect. 3, we perform numerical simulations without dust in order to characterize the properties of the flow subject to ambipolar diffusion (transport, turbulent levels, winds). In Sect. 4, we add the dust in the simulations and determine its vertical and radial distribution. We show in particular that the submillimetre 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 Sect. 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 Sect. 6 we compare the vertical and radial dust distributions with those inferred from ALMA observations (Pinte et al. 2016). We conclude in Sect. 7 by discussing the applications of our work on protoplanetary discs evolution and planet formation.
2. Model and numerical framework
2.1. 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 & LyndenBell 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 multifluid 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 , with c_{s} the uniform sound speed. The evolution of its density ρ, and total velocity field υ follows:(1) (2)
where the total velocity field can be decomposed into a mean shear and a perturbation 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 nonideal induction equation,(4)
The term ∇ × [η_{A}(J × e_{b}) × e_{b}] corresponds to ambipolar diffusion, with J = ∇ × B the current density, the unit vector parallel to the field line and η_{A} the ambipolar diffusivity. This diffusion is due to ionsneutral collisions and is assumed to be the dominant nonideal effect in the regions considered in this paper, that is 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 Sect. 2.4.
Finally, the dust is composed of a mixture of different species, characterizing different grain sizes. Each species is described by a pressureless fluid, with a given density ρ_{d} and velocity υ_{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 species are:(5) (6)
with γ_{g − > d} the drag acceleration induced by the the gas on a dust particle. If we label each dust species by a subscript k, we obtain by conservation of total angular momentum:(7)
In the next section, we specify the form of the drag term γ_{g − > d} and its dependence on the particle size.
2.2. 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:(8)
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(9)
This corresponds to the timescale on which frictional drag will cause an orderofunity 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(10)
If not precised, St denotes the Stokes number in the midplane. We note that the effective Stokes number in the disc atmosphere is larger than St, since it is inversely proportional to the density.
2.3. 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 where ρ_{0} is the midplane density and H = c_{s} / Ω is the disc scaleheight. Thus, combining these different relations, we obtain:(11)
To evaluate Σ, a first possibility is to assume that the disc can be described through the standard MMSN model (Hayashi 1981), with surface density(12)
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(13)
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:(14)(Wardle 2007), where γ_{i} = ⟨συ⟩_{i}/(m_{n} + m_{i}) and ⟨συ⟩_{i} = 1.3 × 10^{−9} cm^{3} s^{−1} is the ionneutral 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 , we also define the dimensionless ambipolar Elsasser number(15)
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, that is 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 horizontallyaveraged 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 sulphur 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} (PerezBecker & 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(16)
The profile of Am is then given by:(17)
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}.
2.5. Numerical methods
We use a modified version of the Godunovbased PLUTO code (Mignone et al. 2007) to simulate the gas and dust motions. This code provides a conservative finitevolume scheme that solves the approximate Riemann problem at each intercell boundary. Physical quantities are evolved in time with a second order RungeKutta scheme, and the orbitaladvection 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. 2016, 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 shearperiodic 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.
2.6. 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(18)
the betaplasma 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 nonideal MHD states. Their detailed setup and the range of Stokes numbers considered are given in Sect. 4.1.
2.7. Diagnostics
To analyse the statistical properties of the MHD flow, we define the box average:(19)
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:(20)
and its mean radial profile, averaged in y and z (over a size l_{z} ≤ L_{z}):(21)
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.(22)
3. AmbipolarMHD simulations without dust
A first and necessary step before computing the dustgas interaction is to characterize the properties of the flow in outer regions of PP discs. For that purpose, we start with nonideal 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 Sect. 4 to initialize simulations with dust.
3.1. Transport and turbulence properties
To characterize the properties of nonideal 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).
Fig. 1. Transport efficiency α as a function of time for different vertical magnetization and Am_{0} = 1. 
Figure 2 (top panels) shows the vertical profiles of horizontallyaveraged 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.
Fig. 2. Top: vertical profiles of horizontallyaveraged Reynolds and Maxwell stresses. Bottom: vertical profiles of horizontallyaveraged rms velocity fluctuations, computed in the three directions. The orange lines denote the absolute value of the mean wind ū(z) component. To help with reading, we add vertical dashed lines in the bottom panels which indicate the disc heightscale H. From left to right: β = 10^{5}, 10^{4} and 10^{3}. 
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 Sect. 6 that such result has consequences on the socalled 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 rms velocity u_{rms}(z) = in all directions. These quantities are important to quantify the diffusion of dust particles (see Sect. 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 − ū(z). The mean component ū(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
Selforganized axisymmetric structures, or commonly called “zonal flows” are found in a variety of nonideal MHD flows threaded by a net vertical field: Halldominated turbulence (Kunz & Lesur 2013), ambipolardominated flows (with or without ohmic diffusion, see Simon & Armitage 2014; Bai 2015) and flows combining all nonideal 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 spacetimediagram (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 pmodes (sound waves). These waves might be generated via the turbulent stress, through a process called “aerodynamic noise” (Lighthill 1952; Heinemann & Papaloizou 2009).
Fig. 3. Time evolution of the radial profiles of gas density (top panels) and vertical magnetic field (lower panels) averaged over the y − z plane within z = ±1.5H. From left to right: β = 10^{5}, 10^{4} and 10^{3}. 
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.
4. Simulations with dust and ambipolar diffusion
We now study the dust dynamics in the MHD flows described in Sect. 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 towards a quasisteady 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 structures and vertical dynamics. Finally, the presence of a mean radial pressure gradient in the disc is investigated.
4.1. 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(23)
with H_{d0} = 0.35H 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}ρ_{dk} / ρ ≪ 1), and whether the dust components are simulated altogether or separately.
Figure 4 shows the spacetime diagrams of (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 timeevolution of 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 towards 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 towards 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 towards 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 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).
Fig. 4. Timeevolution of the vertical profiles of dust density (in log space), averaged in x and y, for three different Stokes number in the case β = 10^{4} (left) and β = 10^{3} (right). The white dashed lines indicate the typical dust heightscale H_{d} measured using the relation (24). 
Fig. 5. Timeevolution of the horizontallyaveraged dust density for different Stokes numbers in the case β = 10^{4}. The top panel is for z = 0 (midplane) while the bottom panel is for z = 2H. 
4.2. Vertical profiles and dust scale height with respect to St
We now characterize in more detail the shape of the vertical equilibrium and estimate the typical dust scaleheight as a function 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.
Fig. 6. Vertical density profiles of dust grains for different Stokes number (in log scale) in the case β = 10^{4}. For comparison, the dashed blue line corresponds to the gas profile. 
To be more quantitative, we define the dust scaleheight H_{d} (St, β) as the altitude z such that(24)
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.
Fig. 7. Dust heightscale H_{d} as a function of the Stokes number for three different plasma β. The diamond markers correspond to direct measurement from numerical simulations, with the definition given by Eq. (25). 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; Fromang & Papaloizou 2006, and Sect. 5). The dashed/red horizontal line denotes the size of the numerical cells. 
By examining the scaling relations of Fig. 7, we found that(25)
for intermediate particle sizes, but saturates to H_{d} ≈ H in the limit of small St. The dependence on St^{−1/2} has been obtained in many other disc turbulence simulations coupling dust and gas dynamics (Fromang & Papaloizou 2006; Zhu et al. 2015) and can be understood, in a first and crude approach, within the framework of the classical diffusion theory (Morfill et al. 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 advectiondiffusion 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 Sect. 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 Sect. 5.
4.3. Effect of zonal flows
As we showed in Sect. 3, MHD flows dominated by ambipolar diffusion develop quasisteady zonal density structures (see Fig. 3). What are their impact on the radial distribution of grains and, more indirectly, on their vertical settling?
4.3.1. 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 bifluid assumption is not valid anymore when ρ_{d} / ρ ≃ 1.
Fig. 8. Timeevolution of the radial profiles of gas density (top panels) and dust densities (other panels). The second, third and fourth rows correspond to Stokes numbers of 0.1, 0.01 and 0.001. The densities are averaged over the y–z plane within z = ±1.5H. From left to right: β = 10^{5}, 10^{4} and 10^{3}. 
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.
4.3.2. 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; Gressel et al. 2015) 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 towards 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).
Fig. 9. Top panel: gas density (colourmap) and poloidal streamlines (cyan lines) showing the outflow topology. The line thickness accounts for the norm of the velocity. Center panel: dust density (colourmap) and poloidal streamlines for St = 0.01. The thickness of the arrows accounts for the norm of the mass flux ρ_{d}u_{d}. Lower panel: same for St = 0.001. 
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 towards 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 Sect. 5 for a more quantitative evidence of this result).
4.4. 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 characterizes the acceleration induced by the pressure gradient:(26)
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 TTauri discs suggest a disc aspect ratio H/R = 0.08 at R = 30 AU (see Fig. 8 of Pinte et al. 2016). We recall that in this configuration, the mean drift velocity Δυ = υ_{d} − υ is:(27)
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 maximize 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 with ĝ_{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.
Fig. 10. Right: dust vertical density profiles (averaged in x and y) at t = 220 Ω^{−1} for the case without mean radial pressure gradient (blue/plain line) and the case with ĝ_{e} = 0.05 (red/dashed line). Left: surface density profile in x for both cases (ρ_{d} averaged over y and z). 
In the limit St ≪ 1 and for perturbations wavelength , Jacquet et al. (2011) showed that the pressure gradient induces a streaming instability with growth rate σ:(28)
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.
5. Settling models
The aim of this section is to elaborate a settling model that accounts for the vertical density profiles obtained in Sect. 4.2. We start with the historical 1D diffusion model of Dubrulle et al. (1995) and attempt to apply it to our nonideal 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.
5.1. 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:(29)
By using such decomposition, the mass conservation equation (5), averaged in x and y, becomes:(30)
where Δυ = υ_{d} − υ is the drift velocity between dust an gas. The first term in the zderivative corresponds to the advectionstretching 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 Sect. 5.1.2 , it is possible to show that Eq. (31) takes the form of an advectiondiffusion equation:(31)(Morfill et al. 1985; Dubrulle et al. 1995; Schräpler & Henning 2004), with the diffusion coefficient and τ_{corr} the correlation time of the turbulent eddies. This equation and the relation liking D_{z} to the vertical rms fluctuations are derived in appendix B. To find equilibrium solutions of this equation, we assume the ergodicity of the system and . If the gas is in hydrostatic equilibrium , an analytical solution is(32)
For uniform diffusion coefficient D_{z} and z ≪ H, this gives:(33)
The distribution is Gaussian and the dust heightscale has a dependence on St^{−1/2} for St ≫ D_{z}/(ΩH^{2}). It tends towards unity in the limit of small St.
5.1.2. Check of the assumptions
To obtain Eq. (32) and the relation (34), one needs to make several assumptions (see Appendix B):

Gas density fluctuations small compared to the background density.

Vertical drift Δu_{z} = υ_{zd} − υ_{z} proportional to . This is called the “terminal velocity approximation”.

Wind advection negligible compared to the turbulent transport .

Turbulent transport in the form of a diffusion operator with D_{z} uniform diffusion coefficient (homogeneous and isotropic turbulence)
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 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 . 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.
Fig. 11. Vertical profiles of the flux terms appearing in the righthand side of Eq. (31), calculated directly form simulations and averaged in time, for β = 10^{4} and three different Stokes numbers. From left to right: St = 0.0025, 0.001 and 0.00025. Green, blue and red lines account respectively for the turbulent correlation term, mean wind and mean drift. The brown dashed line is the mean drift due to gravitational settling estimated in the terminal velocity approximation and assuming that (small gas density fluctuations). 
Finally, to check assumption (4), we plot in Fig. 12 the ratio(34)
Fig. 12. Turbulent diffusion coefficient computed numerically for different β and Stokes numbers (estimated from relation (35)). 
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 welldefined. In the midplane, within z ± 0.5H, we found that D_{z0} ≃ 1.9 × 10^{−4} for β = 10^{5} and D_{z0} ≃ 1.6 × 10^{−4} for β = 10^{4}. However D_{z} is not uniform in z and is multiplied almost by a factor of ten at z = 1.5H. The diffusion coefficient increases rapidly with altitude, fitting quitewell with an exponential function D_{z}(z)≃D_{z0}exp(1.6 z/H). These results are consistent with the theoretical prediction that with τ_{corr} ≃ 0.5 Ω^{−1}. In particular, we showed in Sect. 3.1 that the vertical rms fluctuations increase quasiexponentially 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 Eq. (34). However for St ≲ 0.001, this assumption does not hold; winds and nonuniform 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 downwards. This result suggests that the simple 1D diffusion theory is inappropriate to describe dust settling in this regime.
5.1.3. 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 Sect. 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_{z0} = 1.6 × 10^{−4} (blue lines) and a varying (orange lines) with h_{i} = 0.6 and . The theoretical profile in case of a nonuniform D_{z} is obtained analytically by integrating (33):(35)
Fig. 13. Comparison between numerical and theoretical density profiles for St = 0.01 and St = 0.00025 with Dubrulle’s diffusion prescription (β = 10^{4}). The blue line is derived from the simple model with uniform D_{z}, while the orange line is calculated with , h_{i} = 0.6. 
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_{z0}/(Ω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_{z0}/(Ω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 nonideal 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.
5.1.4. 1D models with a uniform and mean wind
We showed that for small particles (typically St < 0.001), the advection of the mean wind 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:(37)
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:(38)
Assuming a uniform D_{z} = D_{z0}, it is straightforward to show that for(39)
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 zdependence 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.
Fig. 14. Density profiles calculated in the 1D model of Dubrulle et al. (1995), extended to the case with a mean wind , for different St. They correspond to analytical solutions of Eq. (39) with D_{z}(z)=const.=1.6 × 10^{−4} and . 
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 xdependent selforganized and coherent axisymmetric flow induced by ambipolar diffusion. In particular, we showed in Sect. 4.3.2 that the wind is distributed nonuniformly 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 advectiondiffusion equation in 2D. This is the object of the next section.
5.2. 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 Sect. 5.1.2). Note that the 2D model described below applies also for the case β = 10^{4} in the regime of small dust particles.
5.2.1. 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 δυ:(40)
where denotes the average in time and y of a quantity X. Like in the 1D theory, we model the turbulent transport (crosscorrelation of density and velocity fluctuations) through a diffusion operator:(41) (42)
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:(43)
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:(44)
with A_{x} the amplitude of the radial density pattern. Instead of a uniform wind prescribed in 5.1.4, we consider here a xdependent steady axisymmetric outflow modelled by:(45)
where x_{i} = xcosi + zsini. 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 υ_{xd} = υ_{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:(46)
Fig. 15. Gas density (colourmap) and velocity streamlines (cyan arrows) of the 2D toy model. The lines thickness is proportional to the velocity ‖υ‖. 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 towards the centre of the 2D domain). 
5.2.2. 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 . The latter is calculated such that the horizontal average 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:(47) (48)
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.
5.2.3. 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 nonuniformly in x with δ = 0.7H. Both have nonzero inclination i = 40^{∘}. In the first case, the dust forms two layers off the midplane, located around z = 2 H, exactly as in the 1D case (see red line in Fig. 14 for comparison). We 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 Sect. 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 Ω_{w}^{−1} ≈ 200 Ω^{−1}. We 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.
Fig. 16. Equilibrium solutions of the 2D advectiondiffusion Eq. (44) with prescribed υ_{x} and υ_{z} (Eq. (46)) for two different grain sizes. The diffusion coefficients are uniform in z. The colourmaps indicate the dust density distribution and the blue lines are the velocity streamlines in the poloidal plane. The top panels correspond to a uniform wind with δ = ∞ and Ω_{w} = 6 × 10^{−3} while the bottom panels correspond to wind plume of width δ = 0.7H. In both cases, we fixed the inclination of the wind i = 40^{∘} and the amplitude of the zonal flow A_{x} = 0.5, as well as the two diffusion coefficients D_{x} = 2.5 × 10^{−4} and D_{z} = 3 × 10^{−4}. 
To check that our 2D model is compatible with simulations, we compare in Fig. 17 the horizontally averaged density profiles 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.
Fig. 17. Comparison between different settling models for St = 0.001 (top) and St = 0.00025 (bottom) in the case β = 10^{3}. The green line represents the vertical density profile obtained from the simulation. The red/dashed and orange lines are computed from the 2D semianalytical models with parameters detailed in Sect. 5.2.2. The red line is with nonuniform D_{z} in z while the orange corresponds to D_{z} = const. The blue and purple/dotted lines are analytical solutions derived from a 1D diffusion model (Dubrulle et al. 1995) respectively without wind and with uniform wind. All profiles are normalized to unity in the midplane to allow comparison. 
5.3. 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.
6. Discussion and link with observations
6.1. 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 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 nonzero net vertical flux) for which typical Schmidt number is never larger than 10 (Johansen et al. 2006; Fromang et al. 2006). We can interpret this result as a consequence of enhanced Maxwell to Reynolds stress ratio at large B_{z} (see Sect. 3.1). In fact, in nonideal 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 nonideal 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 (Fromang et al. 2013). We checked that for a box of size L_{z} = 6H, the vertical rms 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 overestimated 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.
6.2. 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 submillimetre ALMA observations of TTauri discs like HLTau (Pinte et al. 2016). These constraints, based on the measure of the rings contrast, suggest that submillimetre 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 = 100 AU (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 submillimetre dust particles measured by ALMA lies into the range St= 0.005–0.01 at 30 AU (see Sect. 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 (Fromang & Papaloizou 2006). 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 HLTau, in regions beyond 1 AU, is strongly affected by nonideal MHD processes such as ambipolar diffusion.
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 Fromang et al. (2006) 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 xaxis unit is in AU, the conversion is about 1H ≃ 2.5 AU at R = 30 AU. 
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 edgeon discs is probably a better and promising avenue to measure directly the dust scaleheight. Finally, although disc emission at λ = 1 mm 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.
We note that ongoing ALMA observations of edgeon 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 micrometre 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.
6.3. Rings and comparison with ALMA observations
The Table 1 of Pinte et al. (2016) summarizes the properties of the ringsgaps 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 submillimter 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 shearingboxes might influence the distance between the rings. Note also that other discs than HLTau, 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 HLTau) 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 HLTau. 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.
7. Conclusions
In summary, we characterized the dust dynamics in the outer regions of magnetized accretion discs 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 HLTau (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 (submillimetre 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 towards 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 nonuniformity 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 rms 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 υ_{turb} < 0.05c_{s} in the midplane. Nonideal 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 nonideal 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 (“radialdrift” 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. Nonideal 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 selfinduced 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 socalled pebble accretion (Ormel & Klahr 2010; Lambrechts & Johansen 2012), and therefore reduce the growth timescale of massive cores.
Acknowledgments
This work acknowledges funding from The French ANR under contracts ANR17ERC20007 (MHDiscs) and ANR16CE310013 (Planetformingdiscs). Part of this work was performed using HPC ressources from GENCIIDRIS under the allocation A0020402231 and using the Froggy platform of the CIMENT infrastructure (https://ciment.ujfgrenoble.fr), which is supported by the RhôneAlpes region (GRANT CPER0713 CIRA), the OSUG@2020 labex (reference ANR10 LABX56) and the Equip@Meso project (reference ANR10EQPX2901) of the programme Investissements d’Avenir supervised by the Agence Nationale pour la Recherche.
References
 Partnership, A. L. M. A., Brogan, C. L., Pérez, L. M., et al. 2015, ApJ, 808, L3 [NASA ADS] [CrossRef] [Google Scholar]
 Andrews, S. M., Wilner, D. J., Zhu, Z., et al. 2016, ApJ, 820, L40 [NASA ADS] [CrossRef] [Google Scholar]
 Bai, X.N. 2013, ApJ, 772, 96 [NASA ADS] [CrossRef] [Google Scholar]
 Bai, X.N. 2014, ApJ, 791, 137 [Google Scholar]
 Bai, X.N. 2015, ApJ, 798, 84 [Google Scholar]
 Bai, X.N., & Stone, J. M. 2013, ApJ, 767, 30 [NASA ADS] [CrossRef] [Google Scholar]
 Balbus, S. A., & Hawley, J. F. 1991, ApJ, 376, 214 [NASA ADS] [CrossRef] [Google Scholar]
 Balsara, D. S., Tilley, D. A., Rettig, T., & Brittain, S. D. 2009, MNRAS, 397, 24 [NASA ADS] [CrossRef] [Google Scholar]
 Béthune, W., Lesur, G., & Ferreira, J. 2016, A&A, 589, A87 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Béthune, W., Lesur, G., & Ferreira, J. 2017, A&A, 600, A75 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Birnstiel, T., Fang, M., & Johansen, A. 2016, Space Sci. Rev., 205, 41 [NASA ADS] [CrossRef] [Google Scholar]
 Carballido, A., Fromang, S., & Papaloizou, J. 2006, MNRAS, 373, 1633 [NASA ADS] [CrossRef] [Google Scholar]
 Chiang, E. & Youdin, A. N. 2010, Annual Review of Earth and Planetary Sciences, 38, 493 [Google Scholar]
 Dipierro, G., Ricci, L., Pérez, L., et al. 2018, MNRAS, 475, 5296 [NASA ADS] [CrossRef] [Google Scholar]
 Dubrulle, B., Morfill, G., & Sterzik, M. 1995, ICARUS, 114, 237 [CrossRef] [Google Scholar]
 Dullemond, C. P., & Dominik, C. 2004, A&A, 421, 1075 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Dullemond, C. P., & Penzlin, A. B. T. 2018, A&A, 609, A50 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Flaherty, K. M., Hughes, A. M., Rose, S. C., et al. 2017, ApJ, 843, 150 [NASA ADS] [CrossRef] [Google Scholar]
 Fleming, T. P., Stone, J. M., & Hawley, J. F. 2000, ApJ, 530, 464 [NASA ADS] [CrossRef] [Google Scholar]
 Flock, M., Dzyurkevich, N., Klahr, H., Turner, N. J., & Henning, T. 2011, ApJ, 735, 122 [NASA ADS] [CrossRef] [Google Scholar]
 Flock, M., Fromang, S., González, M., & Commerçon, B. 2013, A&A, 560, A43 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Flock, M., Ruge, J. P., Dzyurkevich, N., et al. 2015, A&A, 574, A68 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Fromang, S., & Papaloizou, J. 2006, A&A, 452, 751 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Fromang, S., Hennebelle, P., & Teyssier, R. 2006, A&A, 457, 371 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Fromang, S., Latter, H., Lesur, G., & Ogilvie, G. I. 2013, A&A, 552, A71 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gammie, C. F. 1996, ApJ, 457, 355 [NASA ADS] [CrossRef] [Google Scholar]
 Garufi, A., Benisty, M., Stolker, T., et al. 2017, The Messenger, 169, 32 [NASA ADS] [Google Scholar]
 Goldreich, P., & LyndenBell, D. 1965, MNRAS, 130, 125 [NASA ADS] [CrossRef] [Google Scholar]
 Gonzalez, J.F., Laibe, G., & Maddison, S. T. 2017, MNRAS, 467, 1984 [NASA ADS] [Google Scholar]
 Gressel, O., & Pessah, M. E. 2015, ApJ, 810, 59 [NASA ADS] [CrossRef] [Google Scholar]
 Gressel, O., Turner, N. J., Nelson, R. P., & McNally, C. P. 2015, ApJ, 801, 84 [NASA ADS] [CrossRef] [Google Scholar]
 Hawley, J. F., Gammie, C. F., & Balbus, S. A. 1995, ApJ, 440, 742 [NASA ADS] [CrossRef] [Google Scholar]
 Hayashi, C. 1981, Progress of Theoretical Physics Supplement, 70, 35 [NASA ADS] [CrossRef] [Google Scholar]
 Heinemann, T., & Papaloizou, J. C. B. 2009, MNRAS, 397, 64 [NASA ADS] [CrossRef] [Google Scholar]
 Isella, A., Guidi, G., Testi, L., et al. 2016, Phys. Rev. Lett., 117, 251101 [NASA ADS] [CrossRef] [Google Scholar]
 Jacquet, E., Balbus, S., & Latter, H. 2011, MNRAS, 415, 3591 [NASA ADS] [CrossRef] [Google Scholar]
 Johansen, A., & Klahr, H. 2005, ApJ, 634, 1353 [NASA ADS] [CrossRef] [Google Scholar]
 Johansen, A., Klahr, H., & Mee, A. J. 2006, MNRAS, 370, L71 [NASA ADS] [Google Scholar]
 Johansen, A., Youdin, A., & Klahr, H. 2009, ApJ, 697, 1269 [NASA ADS] [CrossRef] [Google Scholar]
 Kataoka, A., Tsukagoshi, T., Pohl, A., et al. 2017, ApJ, 844, L5 [NASA ADS] [CrossRef] [Google Scholar]
 Kunz, M. W., & Lesur, G. 2013, MNRAS, 434, 2295 [NASA ADS] [CrossRef] [Google Scholar]
 Laibe, G., & Price, D. J. 2011, MNRAS, 418, 1491 [NASA ADS] [CrossRef] [Google Scholar]
 Laibe, G., & Price, D. J. 2012, MNRAS, 420, 2345 [NASA ADS] [CrossRef] [Google Scholar]
 Lambrechts, M., & Johansen, A. 2012, MNRAS, 544, A32 [Google Scholar]
 Lesur, G., & Papaloizou, J. C. B. 2009, A&A, 498, 1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lesur, G., Kunz, M. W., & Fromang, S. 2014, A&A, 566, A56 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lighthill, M. J. 1952, Proc. Roy. Soc. London Philos. Trans. Ser. A, 211, 564 [Google Scholar]
 Mignone, A., Bodo, G., Massaglia, S., et al. 2007, ApJS, 170, 228 [Google Scholar]
 Miyake, T., Suzuki, T. K., & Inutsuka, S.I. 2016, ApJ, 821, 3 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Morfill, G. E. 1985, in in Birth and the Infancy of Stars, eds. R. Lucas, A. Omont, & R. Stora [Google Scholar]
 Okuzumi, S., Momose, M., Sirono, S.I., Kobayashi, H., & Tanaka, H. 2016, ApJ, 821, 82 [NASA ADS] [CrossRef] [Google Scholar]
 Ormel, C. W., & Klahr, H. H. 2010, A&A, 520, A43 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 PerezBecker, D., & Chiang, E. 2011, ApJ, 735, 8 [NASA ADS] [CrossRef] [Google Scholar]
 Pinte, C., Dent, W. R. F., Ménard, F., et al. 2016, ApJ, 816, 25 [NASA ADS] [CrossRef] [Google Scholar]
 Ruge, J. P., Flock, M., Wolf, S., et al. 2016, A&A, 590, A17 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Sano, T., & Stone, J. M. 2002, ApJ, 577, 534 [NASA ADS] [CrossRef] [Google Scholar]
 Schräpler, R., & Henning, T. 2004, ApJ, 614, 960 [NASA ADS] [CrossRef] [Google Scholar]
 Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337 [NASA ADS] [Google Scholar]
 Simon, J. B., & Armitage, P. J. 2014, ApJ, 784, 15 [NASA ADS] [CrossRef] [Google Scholar]
 Simon, J. B., Bai, X.N., Armitage, P. J., Stone, J. M., & Beckwith, K. 2013, ApJ, 775, 73 [NASA ADS] [CrossRef] [Google Scholar]
 Simon, J. B., Lesur, G., Kunz, M. W., & Armitage, P. J. 2015, MNRAS, 454, 1117 [NASA ADS] [CrossRef] [Google Scholar]
 Takahashi, S. Z., & Inutsuka, S.I. 2014, ApJ, 794, 55 [NASA ADS] [CrossRef] [Google Scholar]
 Venuti, L., Bouvier, J., Flaccomio, E., et al. 2014, A&A, 570, A82 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Wardle, M. 2007, Ap&SS, 311, 35 [NASA ADS] [CrossRef] [Google Scholar]
 Wardle, M., & Salmeron, R. 2012, MNRAS, 422, 2737 [NASA ADS] [CrossRef] [Google Scholar]
 Weidenschilling, S. J. 1977, MNRAS, 180, 57 [Google Scholar]
 Williams, J. P., & McPartland, C. 2016, ApJ, 830, 32 [NASA ADS] [CrossRef] [Google Scholar]
 Wolff, S. G., Perrin, M. D., Stapelfeldt, K., et al. 2017, ApJ, 851, 56 [NASA ADS] [CrossRef] [Google Scholar]
 Wünsch, R., Klahr, H., & Różyczka, M. 2005, MNRAS, 362, 361 [NASA ADS] [CrossRef] [Google Scholar]
 Xu, Z., Bai, X.N., & MurrayClay, R. A. 2017, ApJ, 847, 52 [NASA ADS] [CrossRef] [Google Scholar]
 Youdin, A., & Johansen, A. 2007, ApJ, 662, 613 [NASA ADS] [CrossRef] [Google Scholar]
 Youdin, A. N., & Goodman, J. 2005, ApJ, 620, 459 [NASA ADS] [CrossRef] [Google Scholar]
 Zhu, Z., Stone, J. M., & Bai, X.N. 2015, ApJ, 801, 81 [NASA ADS] [CrossRef] [Google Scholar]
Appendix A:
Tests of dust implementation in PLUTO
A.1 Numerical methods
The numerical integration of a given dust species 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; their Sect. 3.2):(A.1)
where the subscript k labels the different particle species. This prescription, derived from a Von Neumann analysis, guarantees the stability of the simple explicit Eulerscheme. It is then robust but perhaps not optimal for RungeKutta 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.
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 dustgas mixture along the radial direction. As a first step, we neglect the Coriolis force (no rotation) and the shear. Assuming a single dust species of size a and perturbations δρ, δρ_{d}, δυ_{x} and δυ_{xd} ∝ exp(ik_{x} − ωt), the equations governing the linear sound waves in the x direction are:(A.2) (A.3) (A.4) (A.5)
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:(A.6)
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.
Fig. A.1. Decay rates of 1D sound waves in a gasdust 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}. 
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; Laibe & Price 2012). Unlike SPH implementation, we did not experience any convergence problem when considering 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 dephasing 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.
A.3 Shearing waves with rotation
We performed similar tests for 2D nonaxisymmetric 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_{x0} + Sk_{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 timeevolution 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 wellcoupled 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.
Fig. A.2. Evolution of shearing waves kinetic energy in a gasdust mixture with χ = 0.4 and for three different stopping times. From top to bottom: τ_{s}Ω = 0.2, 2, and 20. The left panels correspond to the gas component while the right panels correspond to the dust component. The blue lines are the semianalytical evolution while other lines are the evolution computed numerically with PLUTO. 
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 Sect. 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 semianalytically by solving the linearised system given in Youdin & Goodman (2005). Figure A.3 shows the theoretical growth rate as a function of the stopping time τ_{s} (blue curve).
Fig. A.3. Growth rate of the streaming instability as function of the stopping time τ_{s}. The dust to gas ratio is χ = 0.2 and the radial normalized acceleration associated with the mean pressure gradient is . The blue line is the theoretical growth rate while the red diamond markers are the growth rates computed numerically with PLUTO. 
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 smallscale clumps with local dust to gas ratio greater than 1.
Appendix B:
Derivation of the 1D advectiondiffusion equation of Dubrulle et al. (1995)
In this appendix, we derive the 1D advectiondiffusion equation of Dubrulle et al. (1995; Eq. (32)), commonly used to model vertical dust settling in turbulent discs. In particular we aim to clarify the hypothesis behind this equation. Given the 1D decomposition of a quantity X into a mean component and a fluctuation δX(x, y, z), we showed already in Sect. 5.1.1 that the dust mass conservation equation expands to:(B.1)
We start by developing the mean drift , which appears in the third flux term. We use the “terminal velocity approximation” (Youdin & Goodman 2005) which assumes that the grain instantaneously responds to the gas drag (their inertia is negligible). This approximation is valid only if St ≪1 and if the typical eddies turnover time is larger than τ_{s}. In that case, the drift velocity at leading order is:(B.2)
where ρ^{⋆} = ρ + ρ_{d} is the total density. Yet, by considering the equation for the centremass vertical velocity , and neglecting the quadratic terms in Δυ, we have:(B.3)
The two inertia terms in the right hand side are generally negligible compared to the gravitational term and the gas wind advection. By taking g_{z} = −Ω^{2}z, we deduce from Eqs. (B.2) and (B.3) the mean vertical drift velocity:(B.4)
If the gas density fluctuations are small, we have .
We need then to find a closure relation to express the turbulent correlation term in Eq. (B.1). For that purpose, let us note χ = ρ_{d} / ρ ≪ 1 the dust concentration. By combining the mass conservation equation of both gas and dust species, we obtain(B.6)
By using the drift definition Δυ = υ_{d} − υ and simplifying the expression above, the dust to gas ratio obeys the following continuity equation:(B.7)
In the limit of small drifts which corresponds to small Stokes, the dust is a passive scalar and follows the pure advection equation:(B.8)
The horizontal average of the above equation writes:(B.9)
By expanding Eq. (B.8), we have also:(B.10)
Substracting Eq. (B.9) with Eq. (B.10) gives:(B.11)
To close the system of equations, one needs to neglect the term , which is justified in the limit of small turbulent correlation time. We obtain then:(B.12)
By considering the case without mean wind , and defining the correlation time τ_{corr} of the turbulent eddies such that ∂δχ/∂t ≃ δχ/τ_{corr}, we have then(B.13)
We finally use the fact that δρ/ρ ≪ 1 (anelastic approximation, so that ) to obtain a closure relation that links the turbulent term with the vertical gradient of :(B.14)
with τ_{corr} > 0. The turbulent correlation can be then reduced to a diffusive operator with coefficient D_{z}. Using Eq. (B.1), (B.4) and (B.14), and neglecting the mean wind term in Eq. (B.1), it is straightforward to show that:(B.15)
Appendix C:
The 2D advectiondiffusion solver
In Sect. 5.2, we study the 2D solutions ρ_{d}(x, z) that satisfy the following equation:(C.1)
where υ_{xd}, υ_{zd} and ρ are given 2D fields, and D_{x}, D_{z} are given functions of z. In this appendix, we describe the methods employed to calculate numerically such solutions. First, the equation is rewritten in the form:(C.2)
with A, B, C, D and E given functions of x and z. Then fields are discretized on a 2D grid of size (n_{x}, n_{z}) where first and second order partial derivative are approximated using secondorder centred differences. By defining a vector state X containing the values of ρ at each grid points, the discretized advectiondifffusion equation takes the form of a linear system LX = 0. A possible method to solve such system is to use a singular values decomposition (SVD) and seek for the kernel of L. However this method does not guarantee the positivity of the solution ρ_{d}(x, z) and generates a lot of noise, since the matrix is generally illconditioned. Another method we preferred is to add a temporal derivative ∂ρ_{d}/∂t in the left hand side of Eq. (C.1) and solve for the dynamical system (initial value problem). The solution is then obtained by convergence of the system toward a steady state. To integrate the PDE in time, we use a RungeKutta method of order 4 which guarantees the stability of the numerical scheme. The boundary conditions are periodic in x and with no gradient in z. The timestep is fixed and chosen to enforce the Courant–Friedrichs–Lewy condition. The solver is systematically initialized with the 1D Gaussian solution of Dubrulle et al. (1995). Convergence toward a steady state is obtained generally after few thousands of orbits, which corresponds in real time to few minutes with a single processor at resolution 64 × 64.
All Figures
Fig. 1. Transport efficiency α as a function of time for different vertical magnetization and Am_{0} = 1. 

In the text 
Fig. 2. Top: vertical profiles of horizontallyaveraged Reynolds and Maxwell stresses. Bottom: vertical profiles of horizontallyaveraged rms velocity fluctuations, computed in the three directions. The orange lines denote the absolute value of the mean wind ū(z) component. To help with reading, we add vertical dashed lines in the bottom panels which indicate the disc heightscale H. From left to right: β = 10^{5}, 10^{4} and 10^{3}. 

In the text 
Fig. 3. Time evolution of the radial profiles of gas density (top panels) and vertical magnetic field (lower panels) averaged over the y − z plane within z = ±1.5H. From left to right: β = 10^{5}, 10^{4} and 10^{3}. 

In the text 
Fig. 4. Timeevolution of the vertical profiles of dust density (in log space), averaged in x and y, for three different Stokes number in the case β = 10^{4} (left) and β = 10^{3} (right). The white dashed lines indicate the typical dust heightscale H_{d} measured using the relation (24). 

In the text 
Fig. 5. Timeevolution of the horizontallyaveraged dust density for different Stokes numbers in the case β = 10^{4}. The top panel is for z = 0 (midplane) while the bottom panel is for z = 2H. 

In the text 
Fig. 6. Vertical density profiles of dust grains for different Stokes number (in log scale) in the case β = 10^{4}. For comparison, the dashed blue line corresponds to the gas profile. 

In the text 
Fig. 7. Dust heightscale H_{d} as a function of the Stokes number for three different plasma β. The diamond markers correspond to direct measurement from numerical simulations, with the definition given by Eq. (25). 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; Fromang & Papaloizou 2006, and Sect. 5). The dashed/red horizontal line denotes the size of the numerical cells. 

In the text 
Fig. 8. Timeevolution of the radial profiles of gas density (top panels) and dust densities (other panels). The second, third and fourth rows correspond to Stokes numbers of 0.1, 0.01 and 0.001. The densities are averaged over the y–z plane within z = ±1.5H. From left to right: β = 10^{5}, 10^{4} and 10^{3}. 

In the text 
Fig. 9. Top panel: gas density (colourmap) and poloidal streamlines (cyan lines) showing the outflow topology. The line thickness accounts for the norm of the velocity. Center panel: dust density (colourmap) and poloidal streamlines for St = 0.01. The thickness of the arrows accounts for the norm of the mass flux ρ_{d}u_{d}. Lower panel: same for St = 0.001. 

In the text 
Fig. 10. Right: dust vertical density profiles (averaged in x and y) at t = 220 Ω^{−1} for the case without mean radial pressure gradient (blue/plain line) and the case with ĝ_{e} = 0.05 (red/dashed line). Left: surface density profile in x for both cases (ρ_{d} averaged over y and z). 

In the text 
Fig. 11. Vertical profiles of the flux terms appearing in the righthand side of Eq. (31), calculated directly form simulations and averaged in time, for β = 10^{4} and three different Stokes numbers. From left to right: St = 0.0025, 0.001 and 0.00025. Green, blue and red lines account respectively for the turbulent correlation term, mean wind and mean drift. The brown dashed line is the mean drift due to gravitational settling estimated in the terminal velocity approximation and assuming that (small gas density fluctuations). 

In the text 
Fig. 12. Turbulent diffusion coefficient computed numerically for different β and Stokes numbers (estimated from relation (35)). 

In the text 
Fig. 13. Comparison between numerical and theoretical density profiles for St = 0.01 and St = 0.00025 with Dubrulle’s diffusion prescription (β = 10^{4}). The blue line is derived from the simple model with uniform D_{z}, while the orange line is calculated with , h_{i} = 0.6. 

In the text 
Fig. 14. Density profiles calculated in the 1D model of Dubrulle et al. (1995), extended to the case with a mean wind , for different St. They correspond to analytical solutions of Eq. (39) with D_{z}(z)=const.=1.6 × 10^{−4} and . 

In the text 
Fig. 15. Gas density (colourmap) and velocity streamlines (cyan arrows) of the 2D toy model. The lines thickness is proportional to the velocity ‖υ‖. 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 towards the centre of the 2D domain). 

In the text 
Fig. 16. Equilibrium solutions of the 2D advectiondiffusion Eq. (44) with prescribed υ_{x} and υ_{z} (Eq. (46)) for two different grain sizes. The diffusion coefficients are uniform in z. The colourmaps indicate the dust density distribution and the blue lines are the velocity streamlines in the poloidal plane. The top panels correspond to a uniform wind with δ = ∞ and Ω_{w} = 6 × 10^{−3} while the bottom panels correspond to wind plume of width δ = 0.7H. In both cases, we fixed the inclination of the wind i = 40^{∘} and the amplitude of the zonal flow A_{x} = 0.5, as well as the two diffusion coefficients D_{x} = 2.5 × 10^{−4} and D_{z} = 3 × 10^{−4}. 

In the text 
Fig. 17. Comparison between different settling models for St = 0.001 (top) and St = 0.00025 (bottom) in the case β = 10^{3}. The green line represents the vertical density profile obtained from the simulation. The red/dashed and orange lines are computed from the 2D semianalytical models with parameters detailed in Sect. 5.2.2. The red line is with nonuniform D_{z} in z while the orange corresponds to D_{z} = const. The blue and purple/dotted lines are analytical solutions derived from a 1D diffusion model (Dubrulle et al. 1995) respectively without wind and with uniform wind. All profiles are normalized to unity in the midplane to allow comparison. 

In the text 
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 Fromang et al. (2006) 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 xaxis unit is in AU, the conversion is about 1H ≃ 2.5 AU at R = 30 AU. 

In the text 
Fig. A.1. Decay rates of 1D sound waves in a gasdust 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}. 

In the text 
Fig. A.2. Evolution of shearing waves kinetic energy in a gasdust mixture with χ = 0.4 and for three different stopping times. From top to bottom: τ_{s}Ω = 0.2, 2, and 20. The left panels correspond to the gas component while the right panels correspond to the dust component. The blue lines are the semianalytical evolution while other lines are the evolution computed numerically with PLUTO. 

In the text 
Fig. A.3. Growth rate of the streaming instability as function of the stopping time τ_{s}. The dust to gas ratio is χ = 0.2 and the radial normalized acceleration associated with the mean pressure gradient is . The blue line is the theoretical growth rate while the red diamond markers are the growth rates computed numerically with PLUTO. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.