© ESO 2018

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 magneto-rotational instability (MRI, Balbus & Hawley 1991; Hawley et al. 1995) has long been considered as the main source of turbulence, among various other instabilities (of gravitational or thermodynamics nature) whose role is still debated. Local and global MRI simulations in the ideal limit suggest in particular that the rate of angular momentum transport is compatible with that inferred from observations (Flock et al. 2011, 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, non-ideal MHD effects, such as ohmic diffusion, Hall effect and ambipolar diffusion prevail and tend to suppress any form of MRI-driven turbulence (Fleming et al. 2000; Sano & Stone 2002; Wardle & Salmeron 2012; Bai 2013, 2015; Lesur et al. 2014).

Observationally, the prevalence of a sustained and vigorous MRI-turbulence in the regions ≳1 AU is also largely questioned. The high resolution imaging at optical, near-infrared (VLT/SPHERE) and sub-millimetric wavelenght (ALMA) have provided precious information about the disc structure and turbulence, essentially by probing the spatial distribution of 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 edge-on 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 dust-drift-driven viscous ring instability (Wünsch et al. 2005; Dullemond & Penzlin 2018), snow lines (Okuzumi et al. 2016), MHD turbulence with dead zones (Flock et al. 2015), zonal flows (Johansen et al. 2009) or secular gravitational instabilities in the dust (Takahashi & Inutsuka 2014).

A possible solution to account for all these observations is to consider the atypical dynamics induced by non-ideal MHD effects in protoplanetary discs, in particular the Hall effect and ambipolar diffusion. Indeed, if the latter tend to weaken considerably the disc turbulence, both are able to produce large scale horizontal magnetic fields within the disc, producing strong currents. These currents allow the development of vigorous outflows at the disc surface and a laminar Maxwell stress in the midplane (Lesur et al. 2014), both transporting significant angular momentum (Bai 2013; Simon et al. 2015; Gressel & Pessah 2015). Moreover, these non-ideal effects lead in the nonlinear regime to the formation of self-organized axisymmetric structures, such as zonal flows (Kunz & Lesur 2013; Bai 2015), corresponding to concentric rings in a global disc view (Béthune et al. 2016). These rings, which are very stable and coincide with gas density maxima, are likely to be ideal locations for dust accumulation and planet formation.

Whether the Hall effect and ambipolar diffusion are key physical processes accounting for the observed disc quiescence (in terms of turbulent motions) and ring-like structures remains an open question. Since current sub-millimetre or near-infrared observations are sensitive to dust, direct comparison with them requires to go one step forward and determine the dust dynamics in such non-ideal MHD flows. Ultimately it requires to compute its scattered and direct emission, but this remains out of the scope of the paper. Although the dust behaviour has been widely studied in local and global MRI disc simulations (Johansen & Klahr 2005; Carballido et al. 2006; Fromang & Papaloizou 2006; Balsara et al. 2009; Zhu et al. 2015), it remains relatively unexplored in non-ideal MHD simulations characterizing the disc outer regions (r ≳ 0.1–1 AU). Ruge et al. (2016) studied the case of zero net flux MRI with ohmic resistivity but did not include other relevant non-ideal effects. Several simulations by Zhu et al. (2015) 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 non-ideal 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 non-ideal 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 pressure-less 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³ to 10⁵. 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 non-ideal physics and dust-gas interaction. We also present the numerical methods used to simulate the dust-gas 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 sub-millimetre dust is highly sedimented and form ring structures with high density contrast. We also explore the effect of a mean radial pressure gradient in the disc. In 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 & Lynden-Bell 1965). This corresponds to a Cartesian patch of the disc, centred at r = R ₀, where the Keplerian rotation is approximated locally by a linear shear flow plus a uniform rotation rate Ω = Ωe_z. We note (x, y, z) respectively the radial, azimuthal and vertical directions. We adopt a multi-fluid approximation in which the ionized gas and the dust interact and exchange momentum through drag forces.

The gas is coupled to a magnetic field B and is assumed to be inviscid and isothermal, its pressure P and density ρ related by $P=\rho c_s^2$, with c_s the uniform sound speed. The evolution of its density ρ, and total velocity field υ follows:$$\begin{array}{c}\hfill {\displaystyle \frac{\partial \rho }{\partial t}}+\mathrm{\nabla }·\left(\rho \mathit{\upsilon }\right)=0,\end{array}$$(1) $$\begin{array}{c}\hfill \rho (\frac{\partial \mathit{\upsilon }}{\partial t}+\mathit{\upsilon }·\mathbf{\nabla }\mathit{\upsilon }+2\mathbf{\Omega }\times \mathit{\upsilon })=\rho \mathit{g}-\mathbf{\nabla }P+(\mathrm{\nabla }\times \mathit{B})\times \mathit{B}+\rho {\mathit{\gamma }}_{d->g},\end{array}$$(2)

where the total velocity field can be decomposed into a mean shear and a perturbation u:$$\begin{array}{c}\hfill \mathit{\upsilon }=-q\mathrm{\Omega }x{\mathit{e}}_y+\mathit{u},\end{array}$$(3)

with q = 3/2. The term g = −Ω²z e_z + 2q Ω² x e_x corresponds to the tidal gravity field in the local frame. The last term in the momentum equation (2) contains the acceleration γ_{d − > g} exerted by the dust drag force on a gas parcel (detailed below). The magnetic field B, which appears in the Lorentz force (∇×B)×B, is governed by the non-ideal induction equation,$$\begin{array}{c}\hfill \frac{\partial \mathit{B}}{\partial t}=\mathrm{\nabla }\times (\mathit{\upsilon }\times \mathit{B})+\mathbf{\nabla }\times [{\eta }_A(\mathit{J}\times {\mathit{e}}_{\mathrm{b}})\times {\mathit{e}}_{\mathrm{b}}].\end{array}$$(4)

The term ∇ × [η_A(J × e_b) × e_b] corresponds to ambipolar diffusion, with J = ∇ × B the current density, ${\mathit{e}}_{\text{b}}=\mathit{B}/\Vert \mathit{B}\Vert$ the unit vector parallel to the field line and η_A the ambipolar diffusivity. This diffusion is due to ions-neutral collisions and is assumed to be the dominant non-ideal effect in the regions considered in this paper, that is R₀ = 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 pressure-less 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:$$\begin{array}{c}\hfill {\displaystyle \frac{\partial {\rho }_d}{\partial t}}+\mathrm{\nabla }·\left({\rho }_d{\mathit{\upsilon }}_d\right)=0,\end{array}$$(5) $$\begin{array}{c}\hfill {\rho }_d({\displaystyle \frac{\partial {\mathit{\upsilon }}_d}{\partial t}}+{\mathit{\upsilon }}_d·{\mathbf{\nabla }\mathit{\upsilon }}_d+2\mathbf{\Omega }\times {\mathit{\upsilon }}_d)={\rho }_d(\mathit{g}+{\mathit{\gamma }}_{g->d}),\end{array}$$(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:$$\begin{array}{c}\hfill {\mathit{\gamma }}_{d->g}=-{\displaystyle \frac{1}{\rho }}\sum _k{\rho }_{d_k}{\mathit{\gamma }}_{g->d_k}.\end{array}$$(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:$$\begin{array}{c}\hfill {\mathit{\gamma }}_{g->d}={\displaystyle \frac{1}{{\tau }_s}}(\mathit{\upsilon }-{\mathit{\upsilon }}_d).\end{array}$$(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$$\begin{array}{c}\hfill {\tau }_s={\displaystyle \frac{{\rho }_{s}a}{\rho c_s}}·\end{array}$$(9)

This corresponds to the timescale on which frictional drag will cause an order-of-unity change in the momentum of the dust grain. This is a direct measure of the coupling between dust particles and gas. A useful dimensionless quantity to parametrize this coupling is the Stokes number$$\begin{array}{c}\hfill \text{St}=\mathrm{\Omega }{\tau }_s.\end{array}$$(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 $\mathrm{\Sigma }={\rho }_{0}H\sqrt{2\pi }$ where ρ₀ is the midplane density and H = c_s / Ω is the disc scaleheight. Thus, combining these different relations, we obtain:$$\begin{array}{c}\hfill \text{St}=a_d\left({\displaystyle \frac{{\rho }_s\sqrt{2\pi }}{\mathrm{\Sigma }}}\right).\end{array}$$(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$$\begin{array}{c}\hfill \mathrm{\Sigma }=1700{(R_0/1\text{AU})}^{-3/2}\text{g}{\text{cm}}^{-2}.\end{array}$$(12)

This gives Σ = 10.34 g cm⁻² at R₀ = 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⁻² at R₀ = 30 AU (see Pinte et al. 2016, in the HL Tau disc for instance). By assuming ρ_s = 2.5 g cm³, we obtain the following conversion$$\begin{array}{c}\hfill \text{St}\simeq 0.3-0.6\left({\displaystyle \frac{a}{\text{1 cm}}}\right){\left({\displaystyle \frac{R_0}{30\text{AU}}}\right)}^{3/2}·\end{array}$$(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:$$\begin{array}{c}\hfill {\eta }_A={\displaystyle \frac{{\mathit{B}}^2}{{\gamma }_i{\rho }_n{\rho }_i}}\end{array}$$(14)(Wardle 2007), where γ_i = ⟨συ⟩_i/(m_n + m_i) and ⟨συ⟩_i = 1.3 × 10⁻⁹ cm³ s⁻¹ is the ion-neutral collision rate. ρ_n and ρ_i are respectively the neutral and ion mass density; m_n and m_i their individual mass. By introducing the Alfvén speed ${\upsilon }_A=B/\sqrt{\rho }$, we also define the dimensionless ambipolar Elsasser number$$\begin{array}{c}\hfill \text{Am}={\displaystyle \frac{{\upsilon }_A^2}{\mathrm{\Omega }{\eta }_A}}={\displaystyle \frac{{\gamma }_i{\rho }_i}{\mathrm{\Omega }}}·\end{array}$$(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 horizontally-averaged density, integrated from the vertical boundary, we assume that FUV are blocked for Σ_i(z) > Σ_ic = 0.1 g cm⁻². Given that carbon and sulphur atoms are fully ionized by FUV and immune to recombination on grains, the ionization fraction x_i = ρ_im_n/ (ρm_i) above z_io is 10⁻⁵ (Perez-Becker & Chiang 2011). To avoid any discontinuity of this quantity at z = z_io, we use a smooth function to make the transition between the midplane regions and the ionized layer, so that the ionization profile is$$\begin{array}{c}\hfill x_i(z)={10}^{-5}exp[-{\left({\displaystyle \frac{{\mathrm{\Sigma }}_i(z)}{{\mathrm{\Sigma }}_{\mathit{ic}}}}\right)}^4]·\end{array}$$(16)

The profile of Am is then given by:$$\begin{array}{c}\hfill \text{Am}={\text{Am}}_0+3.3\times {10}^{12}{x}_i(z){\displaystyle \frac{\rho }{{\rho }_0}}{(R_0/1 \text{AU})}^{-5/4},\end{array}$$(17)

where Am₀ 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.5. Numerical methods

We use a modified version of the Godunov-based PLUTO code (Mignone et al. 2007) to simulate the gas and dust motions. This code provides a conservative finite-volume scheme that solves the approximate Riemann problem at each inter-cell boundary. Physical quantities are evolved in time with a second order Runge-Kutta scheme, and the orbital-advection algorithm FARGO is used to reduce numerical dissipation. The implementation of ambipolar diffusion is described in Lesur et al. (2014; Appendix B) and has been used in various studies (Lesur et al. 2014; Béthune et al. 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 defines our unit of time and H = c_s/Ω the disc scaleheight, our unit of length. Boundary conditions are periodic in y and shear-periodic in x. In the vertical direction, we use a standard outflow condition for both the gas and dust velocities but compute a hydrostatic balance in the ghost cells for the gas density. For the magnetic field, we adopt the so called “vertical field” boundary with B_x = B_y = 0 at z = ±L_z/2. The presence of magnetic fields produces outflows, which depletes the mass contained in the box. Then, at each time step, we artificially multiply ρ in each cell by a uniform factor such that the total mass in the box is maintained constant. We checked that the mass injected at each orbital period is negligible compared to the total mass (≈1%) and does not affect the results. Note that the mean B_z is conserved in the box but not the mean horizontal magnetic field, whose evolution depends on the electromotive forces (EMFs) at the boundary.

2.6. Initialization and main parameters

We fix R₀ = 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$$\begin{array}{c}\hfill \beta ={\displaystyle \frac{2c_s^2{\rho }_0}{B_z^2}},\end{array}$$(18)

the beta-plasma parameter with ρ₀ = 1 the midplane density. The midplane ambipolar Elsasser number is fixed to Am₀ = 1 in all cases. Simulations with dust are systematically initialized from developed non-ideal MHD states. Their detailed setup and the range of Stokes numbers considered are given in Sect. 4.1.

2.7. Diagnostics

To analyse the statistical properties of the MHD flow, we define the box average:$$\begin{array}{c}\hfill {\overline{X}}^B=\frac{1}{L_x{L}_y{L}_z}{\displaystyle {\int }_{V}X\mathrm{d}V,}\end{array}$$(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:$$\begin{array}{c}\hfill \overline{X}={\overline{X}}^Z={\displaystyle \frac{1}{L_x{L}_y}}{\displaystyle \int {\displaystyle \int X\mathrm{d}x\mathrm{d}y,}}\end{array}$$(20)

and its mean radial profile, averaged in y and z (over a size l_z ≤ L_z):$$\begin{array}{c}\hfill {\overline{X}}^R={\displaystyle \frac{1}{l_z{L}_y}}{\displaystyle {\int }_{-l_z/2}^{l_z/2}{\displaystyle \int X\mathrm{d}y\mathrm{d}z}}.\end{array}$$(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_xu_y and Maxwell stresses M_xy = −B_xB_y, divided by the box average pressure.$$\begin{array}{c}\hfill \alpha ={\displaystyle \frac{{\overline{H}}_{\mathit{xy}}^B+{\overline{M}}_{\mathit{xy}}^B}{{\overline{\rho }}^B{c}_s^2}}·\end{array}$$(22)

3. Ambipolar-MHD simulations without dust

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

3.1. Transport and turbulence properties

To characterize the properties of non-ideal MHD flows with ambipolar diffusion, we perform three different simulations (assuming Am₀ = 1) with β = 10⁵, 10⁴ and 10³ 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 quasi-steady 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 β₀ = 10⁵ and α ≃ 0.02 for β₀ = 10⁴, which is identical to what Simon et al. (2013) found for a similar setup with the ATHENA code (see their Table 1).

Figure 2 (top panels) shows the vertical profiles of horizontally-averaged Maxwell and Reynolds stress. Each of them contributes positively to the radial transport efficiency α (except H_xy in the midplane for β = 10³), although the Maxwell stress always dominates over the Reynolds stress. For all β, both quantities peak around z ≃ 2H where the ionization is sufficiently large to excite the MRI. Their profiles show a dip in the midplane where the ionization is weak and the MRI absent.

However, as the vertical field increases, a laminar stress due to the large scale B_x and B_y starts to be important in the midplane and induces a shallower dip in Maxwell stress. In particular, the ratio of Maxwell to Reynolds stress sharply increases with magnetization, from 2.5 at β = 10⁵ to ∼100 at β = 10³, and thus deviates significantly from that of ideal MRI. We will see in Sect. 6 that such result has consequences on the so-called Schmidt number, ratio between the gas effective viscosity and dust turbulent mixing. Note that a transition to a magnetic and quasi-laminar stress occuring around β ≃ 10³ 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) = ${(\overline{{\mathit{u}}^2})}^{1/2}$ 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⁵ and 10⁴, radial and vertical velocity dispersions in the midplane are comparable in both cases (≈2 ×10⁻²) 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³) , this wind component is stronger and becomes similar in amplitude to the fluctuating part for z ≳ 1.5H. Note that the azimuthal velocity dispersion increases significantly with B_z in the midplane. This is a consequence of the emergence of large scale radial zonal flows whose amplitude is particularly enhanced at large B_z (see next section).

3.2. Zonal flows, pressure maximas and sound waves

Self-organized axisymmetric structures, or commonly called “zonal flows” are found in a variety of non-ideal MHD flows threaded by a net vertical field: Hall-dominated turbulence (Kunz & Lesur 2013), ambipolar-dominated flows (with or without ohmic diffusion, see Simon & Armitage 2014; Bai 2015) and flows combining all non-ideal effects (Lesur et al. 2014; Bai 2015). It has been pointed out that ambipolar diffusion enhances these features (Bai 2015) and can lead to zonal flows with density amplitude of 10– 20% the background disc density (Simon & Armitage 2014). These zonal flows become large scale rings in the global configuration, as revealed by the simulations of Béthune et al. (2017).

To check whether these axisymmetric structures exist in our stratified simulations, we show in Fig. 3 the spacetime-diagram (t, x) of the mean gas density (top) and B_z (bottom), averaged in y and z, for our three different β. Note that the average is done between z = −1.5H₀ and z = 1.5H₀ 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⁵, the density structures are very faint and do not exceed 5% of the mean background density. For β = 10⁴ and 10³, 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⁴ and 10⁵, the density plots of Fig. 3 reveal the presence of large scale acoustic waves propagating radially in the disc. By computing the Fourier transform of ρ in time and space (along the x direction), we checked that the signal frequencies match the dispersion relation of p-modes (sound waves). These waves might be generated via the turbulent stress, through a process called “aerodynamic noise” (Lighthill 1952; Heinemann & Papaloizou 2009).

The bottom panels show that for β = 10⁴ and 10³, 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 quasi-steady state. We then characterize their vertical distributions and typical scaleheight with respect to St and β. We analyse the effect of zonal flows (pressure bumps) on their radial 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$$\begin{array}{c}\hfill {\rho }_d(t=0)={\rho }_{0}exp(-{\displaystyle \frac{z^2}{2H_{d_0}^2}}),\end{array}$$(23)

with H_d0 = 0.35H and ρ₀ 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⁵, 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⁴, 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³, 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 space-time diagrams of ${\overline{\rho }}_d(z,t)$ (horizontally averaged vertical density distribution) obtained from the simulations, for β = 10⁴ (left) and β = 10³ (right), and for different Stokes numbers. The density profiles seem to converge to a steady state depending on the Stokes number after a certain period of time. To quantify more precisely this timescale, we plot in Fig. 5 the time-evolution of ${\overline{\rho }}_d$ at z = 0 (midplane) and z = 2H, for β = 10⁴. 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 Ω⁻¹/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 Ω⁻¹) 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³, 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).

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⁴ 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 = ρ₀exp(−z/λ_d) with λ_d ≃ H for the minimum grain size.

To be more quantitative, we define the dust scaleheight H_d (St, β) as the altitude z such that$$\begin{array}{c}\hfill {\rho }_d(z=H_d)={\rho }_d(z=0)e^{-\frac{1}{2}}\simeq 0.6{\rho }_{d0}.\end{array}$$(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⁻⁴, the distribution is probably not fully converged and one might expect that H_d is slightly underestimated.

By examining the scaling relations of Fig. 7, we found that$$\begin{array}{cc}\hfill H_d/H& \simeq 0.013{\text{St}}^{-\frac{1}{2}}\text{for}\beta ={10}^{4-5},\hfill \\ \hfill & \simeq 0.019{\text{St}}^{-\frac{1}{2}}\text{for}\beta ={10}^3,\hfill \end{array}$$(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 advection-diffusion equation is solved in the vertical direction which gives the relation (25) in the limit of H_d ≪ H and a flatter dependency in the limit of small St (see Sect 5.1.1). For a comparison with the numerical data, we superimposed in Fig. 7 the predicted scaleheight from this model for the three different magnetizations (see Sect. 5 for details about their calculation).

Although the turbulent diffusion theory seems to reproduce quite well the numerical dependency for β = 10⁵ and for β = 10⁴, it does not explain the flared shape of the density profiles at small St (Fig. 6). In the case β = 10³, 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³, 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 quasi-steady zonal density structures (see Fig. 3). What are their impact on the radial distribution of grains and, more indirectly, on their vertical settling?

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⁵ (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⁴ and β = 10³), dust rings also form, with density contrast and concentration much higher than in the case β = 10⁵ 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³) 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 Ω⁻¹. The evolution of such structure with high dust concentration is hard to predict from our simulations, since the bi-fluid assumption is not valid anymore when ρ_d / ρ ≃ 1.

We check that in most cases, the location where the dust is trapped corresponds to a gas pressure maximum (for comparison, we plotted the radial density profiles of the gas in the top panels of Fig. 8). However, for small St, some rings develop outside gas pressure maxima, like the upper one in the right panels of Fig. 8, for St = 0.01 (β = 10³). 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³. 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 Ω⁻¹. 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).

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:$$\begin{array}{c}\hfill {\widehat{g}}_e=-{\displaystyle \frac{1}{H{\mathrm{\Omega }}^2}}{\displaystyle \frac{1}{\rho }}\left({\displaystyle \frac{\partial P}{\partial R}}\right)\approx \left({\displaystyle \frac{H}{R}}\right),\end{array}$$(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⁻²(R/AU)^2/7, which gives g_e ≃ 0.07 at 30 AU. Observations of T-Tauri discs suggest a disc aspect ratio H/R = 0.08 at R = 30 AU (see Fig. 8 of Pinte et al. 2016). We recall that in this configuration, the mean drift velocity Δυ = υ_d − υ is:$$\begin{array}{c}\hfill \mathbf{\Delta }\mathit{\upsilon }=-c_s{\widehat{g}}_e[{\displaystyle \frac{\text{St}}{1+{(f_g\text{St})}^2}}{\mathit{e}}_x+{\displaystyle \frac{f_g{\text{St}}^2}{2+2{(f_g\text{St})}^2}}{\mathit{e}}_y],\end{array}$$(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 Ω²ĝ_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 Ω⁻¹. When a radial drift is present, the dust layer is slightly wider in z, although the difference remains very small. Axisymmetric rings are still formed but the dust concentration inside is weaker. We checked that the gas radial density structures as well as the dust rings remain at the same location as in the the case ĝ_e = 0.

In the limit St ≪ 1 and for perturbations wavelength $\lambda \gg 2\pi {\widehat{g}}_e{\text{St}}^2$, Jacquet et al. (2011) showed that the pressure gradient induces a streaming instability with growth rate σ:$$\begin{array}{c}\hfill \sigma /\mathrm{\Omega }\simeq f_p{\text{St}}^3{\widehat{g}}_e^2{\left({\displaystyle \frac{k{k}_{x}H}{k_z}}\right)}^2,\end{array}$$(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 non-ideal MHD flows. We examine in particular the influence of a mean wind component in the model. In a second part, we propose a 2D model particularly adapted to the case β = 10³, 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:$$\begin{array}{c}\hfill {\rho }_d=\overline{{\rho }_d}+\delta {\rho }_d;\mathit{\upsilon }=\overline{\mathit{\upsilon }}+\delta \mathit{\upsilon };{\mathit{\upsilon }}_d=\overline{{\mathit{\upsilon }}_d}+\delta {\mathit{\upsilon }}_d.\end{array}$$(29)

By using such decomposition, the mass conservation equation (5), averaged in x and y, becomes:$$\begin{array}{c}\hfill {\displaystyle \frac{\partial \overline{{\rho }_d}}{\partial t}}+{\displaystyle \frac{\partial }{\partial z}}(\overline{{\rho }_d}\overline{{\upsilon }_z}+\overline{\delta {\rho }_d\delta {\upsilon }_{z_d}}+\overline{{\rho }_d}\overline{\mathrm{\Delta }{\upsilon }_z})=0,\end{array}$$(30)

where Δυ = υ_d − υ is the drift velocity between dust an gas. The first term in the z-derivative corresponds to the advection-stretching of dust by the mean vertical gas flow (wind) that we leave behind here. The second term is due to the correlation of turbulent fluctuations which is reduced to a diffusion operator in Dubrulle’s theory. The third term accounts for the mean vertical drift motion of dust due to gravitational settling. By making some assumptions, listed in Sect. 5.1.2 , it is possible to show that Eq. (31) takes the form of an advection-diffusion equation:$$\begin{array}{c}\hfill {\displaystyle \frac{\partial \overline{{\rho }_d}}{\partial t}}={\displaystyle \frac{\partial }{\partial z}}\left(z{\mathrm{\Omega }}^2\overline{{\tau }_s}\overline{{\rho }_d}\right)+{\displaystyle \frac{\partial }{\partial z}}\left[D_z\overline{\rho }{\displaystyle \frac{\partial }{\partial z}}\left({\displaystyle \frac{\overline{{\rho }_d}}{\overline{\rho }}}\right)\right]\end{array}$$(31)(Morfill et al. 1985; Dubrulle et al. 1995; Schräpler & Henning 2004), with $D_z\simeq ⟨\overline{{\upsilon }_z^2}⟩{\tau }_{\text{corr}}>0$ the diffusion coefficient and τ_corr the correlation time of the turbulent eddies. This equation and the relation liking D_z to the vertical rms fluctuations are derived in appendix B. To find equilibrium solutions of this equation, we assume the ergodicity of the system and $\langle \partial \overline{{\rho }_d}/\partial t\rangle \simeq 0$. If the gas is in hydrostatic equilibrium $\rho ={\rho }_0{e}^{-\frac{z^2}{2H}}$, an analytical solution is$$\begin{array}{c}\hfill {\rho }_d(z)={\rho }_{d_0}exp(-{\displaystyle \frac{z^2}{2H^2}})exp(-{\displaystyle \int {\displaystyle \frac{\text{St}\mathrm{\Omega }z{e}^{\frac{z^2}{2H^2}}}{D_z(z)}}\mathrm{d}z})·\end{array}$$(32)

For uniform diffusion coefficient D_z and z ≪ H, this gives:$$\begin{array}{c}\hfill {\rho }_d(z)\simeq {\rho }_{d_0}exp(-{\displaystyle \frac{z^2}{2H_d^2}})\text{with}{\displaystyle \frac{H_d}{H}}={\displaystyle \frac{1}{\sqrt{1+{\displaystyle \frac{\text{St}\mathrm{\Omega }{H}^2}{D_z}}}}}·\end{array}$$(33)

The distribution is Gaussian and the dust heightscale has a dependence on St^−1/2 for St ≫ D_z/(ΩH²). 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):

1. Gas density fluctuations small compared to the background density.

2. Vertical drift Δu_z = υ_zd − υ_z proportional to $z\text{St}/\overline{\rho }$. This is called the “terminal velocity approximation”.

3. Wind advection ${\overline{u}}_z\overline{{\rho }_d}$ negligible compared to the turbulent transport $\overline{\delta {\rho }_d\delta {\upsilon }_{z_d}}$.

4. Turbulent transport in the form of a diffusion operator $-D_z\overline{\rho }{\displaystyle \frac{\partial }{\partial z}}\left({\displaystyle \frac{\overline{{\rho }_d}}{\overline{\rho }}}\right)$ with D_z uniform diffusion coefficient (homogeneous and isotropic turbulence)

The aim of this paragraph is to discuss the validity of these assumptions in the different regimes probed by our simulations.

The first assumption (1) is verified for β = 10⁵ and β = 10⁴. Indeed, in both cases, the density fluctuations never reach more than 30% of the background density. However it is not valid for β = 10³ 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⁴. The turbulent flux $\overline{\delta {\rho }_d\delta {\upsilon }_{z_d}}$ is computed directly from the 3D numerical data, with perturbations δρ, δu obtained by subtracting to each quantity its mean value. We superimposed in the same figures the approximated settling term $z\text{St}\overline{{\rho }_d}/\overline{\rho }$. 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⁻⁴. Regarding hypothesis (3), Fig. 11 shows that the flux associated with wind advection (blue/plain lines) remains small compared to gravitational settling and turbulent transport, provided that St ≳ 0.0025. However for St < 0.001, it is comparable to other terms and carries a large fraction of the dust vertically. By neglecting the wind component, the 1D model of Dubrulle et al. (1995) is probably incomplete for small particles with St ≲ 0.001.

Finally, to check assumption (4), we plot in Fig. 12 the ratio$$\begin{array}{c}\hfill D_z=-{\displaystyle \frac{⟨\overline{\delta {\rho }_d\delta {\upsilon }_{z_d}}⟩}{\langle \overline{\rho }{\displaystyle \frac{\partial }{\partial z}}\left({\displaystyle \frac{\overline{{\rho }_d}}{\overline{\rho }}}\right)\rangle }},\end{array}$$(34)

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 Ω⁻¹. For β = 10⁵ and β = 10⁴ (top panel), D_z is positive and well-defined. In the midplane, within z ± 0.5H, we found that D_z0 ≃ 1.9 × 10⁻⁴ for β = 10⁵ and D_z0 ≃ 1.6 × 10⁻⁴ for β = 10⁴. 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 quite-well with an exponential function D_z(z)≃D_z0exp(1.6 |z|/H). These results are consistent with the theoretical prediction that $D_z\simeq ⟨\overline{{\upsilon }_z^2}⟩{\tau }_{\text{corr}}$ with τ_corr ≃ 0.5 Ω⁻¹. In particular, we showed in Sect. 3.1 that the vertical rms fluctuations increase quasi-exponentially with altitude and have similar values for both β = 10⁵ and β = 10⁴. 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 non-uniform D_z have to be included in the model to obtain acceptable solutions. For β = 10³ (bottom panel of Fig. 12), we found however that the 1D diffusion coefficient is not well-defined, 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⁴ (the case β = 10³ 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⁻⁴ (blue lines) and a varying $D_z=D_{z_0}exp(|\tilde{z}|/h_i)$ (orange lines) with h_i = 0.6 and $\tilde{z}=z/H$. The theoretical profile in case of a non-uniform D_z is obtained analytically by integrating (33):$$\begin{array}{c}\hfill \overline{{\rho }_d}={\overline{\rho }}_{d_0}exp(-{\displaystyle \frac{z^2}{2H^2}})exp(-{\displaystyle \frac{\text{St}\mathrm{\Omega }{H}^2}{D_{z_0}}}\xi (z)),\end{array}$$(35)

with$$\begin{array}{c}\hfill \xi (z)={\displaystyle \frac{1}{h_i}}\sqrt{{\displaystyle \frac{\pi }{2}}}{e}^{-{\displaystyle \frac{1}{2h_i^2}}}[\text{erfi}\left({\displaystyle \frac{h_i\tilde{z}-1}{\sqrt{2}{h}_i}}\right)-\text{erfi}\left({\displaystyle \frac{-1}{\sqrt{2}{h}_i}}\right)]+e^{{\displaystyle \frac{{\tilde{z}}^2}{2}}-{\displaystyle \frac{\tilde{z}}{h_i}}}-1.\end{array}$$(36)

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²), 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⁻⁴ ≃ D_z0/(ΩH²), both solutions fail to reproduce the numerical vertical distributions. If the heightscale calculated from the theoretical solution is not too far from the simulated one, the shape of the profile differs significantly. In particular, the 1D model does not account for the flaring at large z.

In conclusion, we showed that the 1D diffusion model of Dubrulle et al. (1995) applied to a non-ideal MHD turbulent flow works pretty well for β = 10⁵ and β = 10⁴, 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 $\overline{{\rho }_d}\overline{{\upsilon }_z}$ in Eq. (31) becomes comparable to the gravitational settling. Then, assumption (3) needs to be relaxed. A first and naive approach is to model the wind by a linear function of z but uniform in x:$$\begin{array}{c}\hfill \overline{{\upsilon }_z}={\mathrm{\Omega }}_{w}z,\end{array}$$(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⁻⁴, 7 × 10⁻⁴ and 6 × 10⁻³ Ω for β = 10⁵, 10⁴ and 10³. The modified equilibrium, including the effect of the wind, is:$$\begin{array}{c}\hfill {\rho }_d(z)={\rho }_{d_0}exp(-{\displaystyle \frac{z^2}{2H^2}})exp(-{\displaystyle \int {\displaystyle \frac{(\text{St}\mathrm{\Omega }{e}^{\frac{z^2}{2H^2}}-{\mathrm{\Omega }}_w)z}{D_z(z)}}\mathrm{d}z})·\end{array}$$(38)

Assuming a uniform D_z = D_z0, it is straightforward to show that for$$\begin{array}{c}\hfill \text{St}\le {\displaystyle \frac{{\mathrm{\Omega }}_w}{\mathrm{\Omega }}}-{\displaystyle \frac{D_{z_0}}{\mathrm{\Omega }{H}^2}},\end{array}$$(39)

the dust density increases with altitude in the vicinity of z = 0. In case β = 10⁴ with Ω_w = 7 × 10⁻⁴, this occurs for a critical St = 5.4 × 10⁻⁴. We plot in Fig. 14 the solutions of Eq. (39) for different Stokes numbers and β = 10⁴. For St = 2.5 × 10⁻⁴, the density profile exhibits two bumps at z ≃ 1.5H where the dust settles. This configuration corresponds to “floating” or “levitating” dust layers. Particles are trapped between z ≃ H where the wind dominates and regions of higher altitudes where gravitational settling dominates (the latter is an exponential function of z while the wind advection is only proportional to z). We checked that assigning a z-dependence of D_z leads to similar results. These floating layers are actually reminiscent of the simulations by Miyake et al. (2016), who found that magnetorotational driven winds can concentrate dust grains with size 20–40 μm around four disc scaleheights by a similar mechanism.

In conclusion, the simple 1D wind model does not reproduce at all the profiles computed numerically. Indeed, none of them exhibit such distribution with floating particles. The reason behind that is intimately related to the baseline assumptions of the 1D model and particularly to the decomposition (41). While the gas density fluctuations remain small compared to the background density, it is not necessarily true for the gas vertical velocity (see Fig. 2 for example). The term “fluctuations” here is actually misleading since it does include both random turbulent motions and the x-dependent self-organized and coherent axisymmetric flow induced by ambipolar diffusion. In particular, we showed in Sect. 4.3.2 that the wind is distributed non-uniformly in x and rather active in the gas density minima, where the dust concentration is low. By averaging quantities in x, we thus overestimate the wind effect. Moreover, the presence of zonal flows in x leads to a large scale meridional circulation (see Fig. 9) which is not taken into account in the 1D model. To better understand the dust distribution in magnetized discs with zonal flows (typically for β = 10³ and marginally for β = 10⁴), one needs to treat the advection-diffusion 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³ 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⁴ 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 δυ:$$\begin{array}{c}\hfill {\rho }_d=⟨{\overline{{\rho }_d}}^a⟩+\delta {\rho }_d;\mathit{\upsilon }=⟨{\overline{\mathit{\upsilon }}}^a⟩+\delta \mathit{\upsilon };{\mathit{\upsilon }}_d=⟨{\overline{{\mathit{\upsilon }}_d}}^a⟩+\delta {\mathit{\upsilon }}_d,\end{array}$$(40)

where $⟨{\overline{X}}^a⟩(x,z)=⟨\int X\mathrm{d}y⟩$ denotes the average in time and y of a quantity X. Like in the 1D theory, we model the turbulent transport (cross-correlation of density and velocity fluctuations) through a diffusion operator:$$\begin{array}{c}\hfill ⟨{\overline{\delta {\rho }_d\delta {\upsilon }_{z_d}}}^a⟩=-D_z⟨{\overline{\rho }}^a⟩{\displaystyle \frac{\partial }{\partial z}}\left({\displaystyle \frac{⟨{\overline{{\rho }_d}}^a⟩}{⟨{\overline{\rho }}^a⟩}}\right),\end{array}$$(41) $$\begin{array}{c}\hfill ⟨{\overline{\delta {\rho }_d\delta {\upsilon }_{x_d}}}^a⟩=-D_x{\displaystyle \frac{\partial ⟨{\overline{{\rho }_d}}^a⟩}{\partial x}},\end{array}$$(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:$$\begin{array}{c}\hfill {\displaystyle \frac{\partial ({\rho }_d{\upsilon }_{x_d})}{\partial x}}+{\displaystyle \frac{\partial ({\rho }_d{\upsilon }_{z_d})}{\partial z}}={\displaystyle \frac{\partial }{\mathrm{d}z}}\left[D_z\rho {\displaystyle \frac{\partial }{\partial z}}\left({\displaystyle \frac{{\rho }_d}{\rho }}\right)\right]+{\displaystyle \frac{\partial }{\mathrm{d}x}}\left[D_x{\displaystyle \frac{\partial {\rho }_d}{\partial x}}\right]·\end{array}$$(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:$$\begin{array}{c}\hfill \widehat{\rho }={\rho }_0[1-A_{x}cos\left({\displaystyle \frac{2\pi x}{L_x}}\right)]exp(-{\displaystyle \frac{z^2}{2H^2}}),\end{array}$$(44)

with A_x the amplitude of the radial density pattern. Instead of a uniform wind $\overline{{\upsilon }_z}={\mathrm{\Omega }}_{w}z$ prescribed in 5.1.4, we consider here a x-dependent steady axisymmetric outflow modelled by:$$\begin{array}{c}\hfill \{\begin{array}{c}{\upsilon }_z={\mathrm{\Omega }}_{w}z\mathcal{G}(x)\text{with}\mathcal{G}(x)={\sum }_{n=-\infty }^{\infty }exp(-{\displaystyle \frac{{(x_i+n{L}_{x}cosi)}^2}{2{\delta }^2}})\hfill \\ {\upsilon }_x=-{\upsilon }_{z}tani\hfill \end{array}\end{array}$$(45)

where x_i = xcosi + zsini. This simple prescription mimics the inclined outflow obtained for β = 10³ 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:$$\begin{array}{c}\hfill {\upsilon }_{d_z}={\upsilon }_z-{\displaystyle \frac{\text{St}\mathrm{\Omega }z}{\widehat{\rho }}}·\end{array}$$(46)

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 ${\mathrm{\Omega }}_w=6\times {10}^{-3}{L}_x/(\delta \sqrt{2\pi })$. The latter is calculated such that the horizontal average $\overline{{\upsilon }_z}(z)$ matches with the 1D prescription of 5.1.4. The amplitude of the zonal flow A_x is fixed to 0.5. The two diffusion coefficients remain to be estimated from simulations. We remind that their definition in the 2D formalism are:$$\begin{array}{c}\hfill D_x(z)=\int {\displaystyle \frac{⟨\int \delta {\rho }_d\delta {\upsilon }_{x_d}\mathrm{d}y⟩}{{\displaystyle \frac{\partial }{\mathrm{d}x}}⟨\int {\rho }_d\mathrm{d}y⟩}}\mathrm{d}x,\end{array}$$(47) $$\begin{array}{c}\hfill D_z(z)=\int {\displaystyle \frac{\langle \int \delta {\rho }_d\delta {\upsilon }_{z_d}\mathrm{d}y\rangle }{\langle \int \rho \mathrm{d}y\rangle {\displaystyle \frac{\partial }{\mathrm{d}z}}\left({\displaystyle \frac{\langle \int {\rho }_d\mathrm{d}y\rangle }{\langle \int \rho \mathrm{d}y\rangle }}\right)}}\mathrm{d}x.\end{array}$$(48)

The numerical calculation for β = 10³ gives D_x ≃ 2.5 × 10⁻⁴ and D_z ≃ 3 × 10⁻⁴ 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 non-uniformly in x with δ = 0.7H. Both have non-zero inclination i = 40^∘. In the first case, the dust forms two layers off the midplane, located around z = 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²/D_x ≈ 4000 Ω⁻¹ is much longer than the outflow timescale Ω_w⁻¹ ≈ 200 Ω⁻¹. 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.

To check that our 2D model is compatible with simulations, we compare in Fig. 17 the horizontally averaged density profiles $\overline{{\rho }_d}$ of the 2D solutions (red dashed lines and orange lines) with those obtained in the shearing box (green lines). For both Stokes numbers St = 0.001 and St = 0.00025 our model is able to accurately reproduce the shape of the vertical profiles. For St = 0.00025, the result is independent on whether the diffusion coefficients are uniform or not. However this is not necessarily true for St = 0.001 for which the 2D model with uniform D_z fails to reproduce the simulations data. We also superimpose in Fig. 17 the profiles obtained with the 1D Dubrulle’s model (blue lines) and the 1D model with wind (purple dotted line). Clearly the 2D models are better to explain the simulations. This indicates that the dust poloidal circulation induced by the zonal flows is crucial to model the dust vertical equilibrium.

5.3. Conclusions on settling models

In the regime of β ≳ 10⁴ 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⁻⁴. This corresponds to values of the Schmidt number S_​c = αΩH² / D_z between 10 for β = 10⁵ and 200 for β = 10³. Hence, the gas transport is always much larger than the vertical dust transport. This in contrast with ideal MRI turbulence simulations (either zero or non-zero net vertical flux) for which typical Schmidt number is never larger than 10 (Johansen et al. 2006; 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 non-ideal MHD flows, angular momentum is mainly extracted through coronal winds and laminar Mawell stress, while dust is more sensitive to turbulent gas motions. Note that the large S_​c found in our simulations, as well as the the small H_d / H and τ_corrΩ ≃ 0.5, seem incompatible with the recent non-ideal MHD simulations of Zhu et al. (2015) with ambipolar diffusion. Indeed in their case, they found S_​c ≲ 1 and diffusion coefficients 10 to 20 times larger than ours. The main difference between their simulations and ours is the vertical and radial box size (they used L_x = 4H, L_z = 6H. while we use a larger domain with L_x = 8H, L_z = 12H). We ran few simulations by varying both L_x and L_z and found that the discrepancy is mainly due to the vertical extent L_z. The reason is that the wind properties, in particular its strength and mass loss rate, depend strongly on the the vertical extent of the box (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 over-estimated by our simulations. To avoid such bias, it is necessary to model the wind in a global way and characterize the dust settling in global simulations, which will be the object of a future paper.

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 sub-millimetre ALMA observations of T-Tauri discs like HL-Tau (Pinte et al. 2016). These constraints, based on the measure of the rings contrast, suggest that sub-millimetre particles are contained in a geometrically thin disc with scale height between 0.7 and 2 AU at 100 AU. If we assume that the ratio H_d / H does not vary strongly with R and consider disc aspect ratio H/R ≃ 0.1 at R = 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 sub-millimetre 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³, 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⁻³) is 20 times larger than those obtained in our simulations with ambipolar diffusion. All these results are indications that the dynamics of discs like HL-Tau, in regions beyond 1 AU, is strongly affected by non-ideal MHD processes such as ambipolar diffusion.

There are however many caveats associated with such comparison. First the uncertainties on the gas surface density in the outer regions can be more important than those estimated. The current MMSN models and the observations can differ by one order of magnitude at most. We have also assumed that the dust to gas scaleheight ratio H_d / H does not vary with radius, which is a crude approximation. In reality, it may encounter large variations between 30 AU and 100 AU, either due to radial variation of turbulence strength or local change in the gas to dust ratio (Pinte et al. 2016). The estimation of H_d from the rings contrast is also largely arguable, since it depends on the dust opacities and radiative transfer used in the synthesized model of the disc. High resolution imaging of edge-on discs is probably a better and promising avenue to measure directly the dust scaleheight. Finally, although disc emission at λ = 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 edge-on discs (in particular in HH30, Menard et al., in prep), combined with previous HST surveys, are going to provide accurate measure of both millimetre and 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³. 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 rings-gaps in HL Tau estimated from the millimetre emission. In particular information about the gap width and depth is given. The width associated with the gap at 30 AU is about 10 AU, although it varies moderately from one gap to another. In our simulations, the typical width of the sub-millimter dust rings (St ≃ 0.01), averaged over the different structures, is respectively 0.4H, 1.7H and ≳3.5H for β = 10⁵, 10⁴ and 10³. Given that H ≃ 2.5 AU in our case, there is again a good agreement between the numerical result and the observation for β = 10³. We stress however that the boundary conditions and the lack of global radial structure in shearing-boxes might influence the distance between the rings. Note also that other discs than HL-Tau, like Elias 24, exhibit much wider rings or gaps with size of few tenth of AU. To explain such systems, our simulations are probably not relevant anymore; the presence of an embedded planet with a mass comparable to that of Jupiter might be required to form such gaps (Dipierro et al. 2018).

As for the gap depth, which measures the ratio between the average dust density in the ring and the density at the minimum, direct comparison is much less obvious and has to be taken with caution. Indeed, due to the instrumental resolution (≃3–5 AU for HL-Tau) the contrast between the rings cannot be measured accurately. At best, a lower bound for the gaps depth can be estimated, which is around 20 for HL-Tau. In the simulations, this quantity is also extremely difficult to estimate since it drastically changes between the different gaps. In the case β = 10³, 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⁶, 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³. 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³, the rings separation (or gaps width) is of order 3–4 H, which is in agreement with measures in HL-Tau (Pinte et al. 2016). We checked that the presence of a mean pressure gradient in the box does not change drastically these results, although our simulations are not enough resolved to trigger the streaming instability and thus determine its impact on the dust distribution.

For small Stokes number (sub-millimetre to micron size particles), we showed on contrary that the dust is less sedimented in the vertical direction, with H_d / H ≈ 1 and forms much fainter rings. In case of a strong magnetization (β = 10³), 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 non-uniformity of the wind in the radial direction. Indeed, models assuming a uniform outflow (inclined or not) lead to “floating” dust layers which are not seen in the present simulations.

We discuss finally the implications of this work on accretion and planet formation scenario in protoplanetary discs. Since they are cold and dense, such objects are inevitably subject to non-ideal MHD effects when threaded by a magnetic field. Our simulations suggest that both accretion efficiency (up to α ≃ 10⁻²) 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. Non-ideal MHD effects are therefore a possible avenue to reconcile the accretion theory with the observations. A next step is to perform global simulations, including a more realistic dust size distribution and the three non-ideal effects. In particular the role of the Hall effect (Lesur et al. 2014; Bai 2014) on the dust is still unknown and has to be investigated.

This work is also directly relevant for planet formation theory and may be important to understand the growth of grains to pebbles (“radial-drift” barrier) and from planetesimals to planets (“fragmentation” barrier) (Birnstiel et al. 2016). Indeed, the confinement of solids in pressure maxima (here the rings) is known to be a very efficient way to overcome these barriers. In addition to stop their radial migration and accelerate their growth, such trapping regions also reduce their relative velocities, which prevents grain fragmentation (Gonzalez et al. 2017). The high dust to gas ratio in these regions can ultimately lead to gravitational runaway of the solids. Non-ideal MHD effects appear to be an alternative to other trapping mechanisms such as vortices, whose stability and origin remain unclear (Lesur & Papaloizou 2009). It may be also more generic and efficient at collecting dust particles than the streaming instability (and related self-induced traps) whose excitation requires very specific conditions. The high dust to gas ratio in the midplane due to the strong sedimentation induced by MHD zonal flows could also enhance the so-called pebble accretion (Ormel & Klahr 2010; Lambrechts & Johansen 2012), and therefore reduce the growth timescale of massive cores.

Acknowledgments

This work acknowledges funding from The French ANR under contracts ANR-17-ERC2-0007 (MHDiscs) and ANR-16-CE31-0013 (Planet-forming-discs). Part of this work was performed using HPC ressources from GENCI-IDRIS under the allocation A0020402231 and using the Froggy platform of the CIMENT infrastructure (https://ciment.ujf-grenoble.fr), which is supported by the Rhône-Alpes region (GRANT CPER07-13 CIRA), the OSUG@2020 labex (reference ANR10 LABX56) and the Equip@Meso project (reference ANR-10-EQPX-29-01) of the programme Investissements d’Avenir supervised by the Agence Nationale pour la Recherche.

References

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):$$\begin{array}{c}\hfill \mathrm{\Delta }t<\underset{k}{min}\left({\displaystyle \frac{\rho {\tau }_{s_k}}{\rho +{\rho }_{d_k}}}\right),\end{array}$$(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 Euler-scheme. It is then robust but perhaps not optimal for Runge-Kutta solvers. Finally, in order to avoid strong numerical diffusion, we implement a version of the FARGO algorithm for the dust components, which treats the shear advection separetely from the other flux and source terms.

A.2 Sound waves

To check our implementation, we first perform a series of tests involving the propagation of simple axisymmetric sound waves in a dust-gas mixture along the radial direction. As a first step, we neglect the Coriolis force (no rotation) and the shear. Assuming a single dust 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:$$\begin{array}{c}\hfill -i\omega \delta \rho +i{k}_x\rho \delta {\upsilon }_x=0,\end{array}$$(A.2) $$\begin{array}{c}\hfill -i\omega \delta {\rho }_d+i{k}_x{\rho }_d\delta {\upsilon }_{x_d}=0,\end{array}$$(A.3) $$\begin{array}{c}\hfill (-i\omega +{\displaystyle \frac{\chi }{{\tau }_s}})\delta {\upsilon }_x+i{\displaystyle \frac{k_x{c}_s^2}{\rho }}\delta \rho -{\displaystyle \frac{\chi \delta {\upsilon }_{x_d}}{{\tau }_s}}=0,\end{array}$$(A.4) $$\begin{array}{c}\hfill (-i\omega +{\displaystyle \frac{1}{{\tau }_s}})\delta {\upsilon }_{x_d}-{\displaystyle \frac{\delta {\upsilon }_x}{{\tau }_s}}=0,\end{array}$$(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:$$\begin{array}{c}\hfill {\omega }^3+{\displaystyle \frac{i{\omega }^2}{{\tau }_s}}(1+\chi )-k_x^2{c}_s^2\omega -i{\displaystyle \frac{k_x^2{c}_s^2}{{\tau }_s}}=0.\end{array}$$(A.6)

The solution is composed of one standing wave and two oscillatory modes with frequencies close to k_xc_s (in the limit of small τ_s). Because the dust and gas exert a mutual drag, the oscillatory modes decay in time with a rate that depends on the dust to gas ratio and τ_s. In Fig. A.1, we plotted in red/dashed line the theoretical decay rate as function of τ_s, calculated analytically from the dispersion relation, for a fixed dust to gas ratio χ = 1.

To check that our PLUTO implementation is correct, we simulate a linear sound wave for different τ_s between 0.1 and 1000, and for χ = 1. The wave corresponds initially to a pure eigenmode of the linear system with an amplitude of 10⁻³. 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 de-phasing between gas and dust is small, of order c_sτ_s, and must be resolved numerically in order to capture properly the physics of the wave dissipation. Our tests in PLUTO suggest however that the dephasing length does not need to be resolved to obtain correct dissipation. Such length is actually purely geometrical and is not associated with any energy or momentum transfer.

A.3 Shearing waves with rotation

We performed similar tests for 2D non-axisymmetric waves of the form exp(ik_x + ik_y − Ωt) with k_y ≠ 0, in the presence of rotation and shear. Since these waves have a wavenumber k_x = k_x0 + Sk_yt 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⁻³. We tested three different stopping time by fixing the dust to gas ratio χ = 0.4. The equation of state is isothermal and the gas is unmagnetized. The results are shown in Fig. A.2, where the blue curves represent the time-evolution of kinetic energy integrated with our linear solver. The other curves correspond to PLUTO simulations with different resolutions. Our code reproduces quite well the desired solution during the first shearing times but the gas component dissipates at longer times (t > 10 Ω⁻¹) as the wave is strongly sheared out and damped by numerical diffusion. We checked in particular that increasing the numerical resolution makes the evolution of the waves closer to the theoretical prediction. Note that for τ_sΩ = 0.2, the dust component tends to follow the gas oscillations within the first orbits, since τ_s is smaller than the oscillation period. Dust and gas are well-coupled and the wave is rapidly damped due to the strong drag between the two components. However for τ_s = 20 Ω⁻¹, 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.

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 semi-analytically 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).

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⁻³. 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 Ω⁻¹, however we do not reproduce the correct growth rate. Actually we found that this discrepancy is not due to our implementation but to the presence of unstable harmonics with much higher growth rates, probably excited by the initial numerical noise. To make sure that our implementation is robust, we also simulate the nonlinear evolution of such mode in a bigger box (L_x = L_z = 6H) with higher resolution (N_X = 512, N_Z = 256) during 1000 orbits. We recover the classical result of Youdin & Johansen (2007) that dust concentrates into small-scale clumps with local dust to gas ratio greater than 1.

Appendix B:

Derivation of the 1D advection-diffusion equation of Dubrulle et al. (1995)

In this appendix, we derive the 1D advection-diffusion 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 $\overline{X}(z)$ and a fluctuation δX(x, y, z), we showed already in Sect. 5.1.1 that the dust mass conservation equation expands to:$$\begin{array}{c}\hfill {\displaystyle \frac{\partial \overline{{\rho }_d}}{\partial t}}+{\displaystyle \frac{\partial }{\partial z}}(\overline{{\rho }_d}\overline{{\upsilon }_z}+\overline{\delta {\rho }_d\delta {\upsilon }_{z_d}}+\overline{{\rho }_d}\overline{\mathrm{\Delta }{\upsilon }_z})=0,\end{array}$$(B.1)

We start by developing the mean drift $\overline{\mathrm{\Delta }{\upsilon }_z}$, 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:$$\begin{array}{c}\hfill \mathbf{\Delta }\mathit{\upsilon }={\tau }_s({\displaystyle \frac{\mathbf{\nabla }P}{{\rho }^{\star }}}-{\displaystyle \frac{(\mathrm{\nabla }\times \mathit{B})\times \mathit{B}}{{\rho }^{\star }}}),\end{array}$$(B.2)

where ρ^⋆ = ρ + ρ_d is the total density. Yet, by considering the equation for the centre-mass vertical velocity ${\upsilon }_z^{\star }={\displaystyle \frac{\rho {\upsilon }_z+{\rho }_d{{\upsilon }_d}_z}{{\rho }^{\star }}}$, and neglecting the quadratic terms in Δυ, we have:$$\begin{array}{c}\hfill {[{\displaystyle \frac{\mathbf{\nabla }P}{{\rho }^{\star }}}-{\displaystyle \frac{(\mathrm{\nabla }\times \mathit{B})\times \mathit{B}}{{\rho }^{\star }}}]}_z=g_z-\frac{\partial {\upsilon }_z^{\star }}{\partial t}-{\mathit{\upsilon }}^{\star }·\mathbf{\nabla }{\upsilon }_z^{\star }.\end{array}$$(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 = −Ω²z, we deduce from Eqs. (B.2) and (B.3) the mean vertical drift velocity:$$\begin{array}{c}\hfill \overline{\mathrm{\Delta }{\upsilon }_z}=-{\overline{\tau }}_s{\mathrm{\Omega }}^{2}z,\end{array}$$(B.4)

with$$\begin{array}{c}\hfill {\overline{\tau }}_s={\displaystyle \frac{a{\rho }_s}{c_s}}\overline{{\rho }^{-1}}.\end{array}$$(B.5)

If the gas density fluctuations are small, we have $\overline{{\rho }^{-1}}\simeq {\overline{\rho }}^1$.

We need then to find a closure relation to express the turbulent correlation term $\overline{\delta {\rho }_d\delta {\upsilon }_{z_d}}$ 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$$\begin{array}{c}\hfill {\displaystyle \frac{\partial \chi }{\partial t}}=-{\displaystyle \frac{1}{\rho }}(\mathrm{\nabla }·\left(\chi \rho {\mathit{\upsilon }}_d\right)-\mathrm{\nabla }·\left(\rho \mathit{\upsilon }\right)\chi ).\end{array}$$(B.6)

By using the drift definition Δυ = υ_d − υ and simplifying the expression above, the dust to gas ratio obeys the following continuity equation:$$\begin{array}{c}\hfill \rho ({\displaystyle \frac{\partial \chi }{\partial t}}+\mathit{\upsilon }·\mathrm{\nabla }\chi )=-\mathrm{\nabla }·\left(\chi \rho \mathbf{\Delta }\mathit{\upsilon }\right).\end{array}$$(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:$$\begin{array}{c}\hfill {\displaystyle \frac{\partial \chi }{\partial t}}+\mathit{\upsilon }·\mathrm{\nabla }\chi =0.\end{array}$$(B.8)

The horizontal average of the above equation writes:$$\begin{array}{c}\hfill {\displaystyle \frac{\partial \overline{\chi }}{\partial t}}+\overline{{\upsilon }_z}·{\displaystyle \frac{\partial \overline{\chi }}{\partial z}}+\overline{\mathit{\upsilon }·\mathrm{\nabla }\delta \chi }=0.\end{array}$$(B.9)

By expanding Eq. (B.8), we have also:$$\begin{array}{c}\hfill {\displaystyle \frac{\partial \chi }{\partial t}}+\overline{{\upsilon }_z}·{\displaystyle \frac{\partial \overline{\chi }}{\partial z}}+\delta {\upsilon }_z{\displaystyle \frac{\partial \overline{\chi }}{\partial z}}+\mathit{\upsilon }·\mathrm{\nabla }\delta \chi =0.\end{array}$$(B.10)

Substracting Eq. (B.9) with Eq. (B.10) gives:$$\begin{array}{c}\hfill {\displaystyle \frac{\partial \delta \chi }{\partial t}}+\delta {\upsilon }_z{\displaystyle \frac{\partial \overline{\chi }}{\partial z}}+\overline{\mathit{\upsilon }}·\mathrm{\nabla }\delta \chi +\delta \mathit{\upsilon }·\mathrm{\nabla }\delta \chi -\overline{\delta \mathit{\upsilon }·\mathrm{\nabla }\delta X}=0.\end{array}$$(B.11)

To close the system of equations, one needs to neglect the term $\delta \mathit{\upsilon }·\mathrm{\nabla }\delta \chi -\overline{\delta \mathit{\upsilon }·\mathrm{\nabla }\delta \chi }$, which is justified in the limit of small turbulent correlation time. We obtain then:$$\begin{array}{c}\hfill {\displaystyle \frac{\partial \delta \chi }{\partial t}}+\delta {\upsilon }_z{\displaystyle \frac{\partial \overline{\chi }}{\partial z}}+\overline{{\upsilon }_z}·{\displaystyle \frac{\partial \delta \chi }{\partial z}}=0.\end{array}$$(B.12)

By considering the case without mean wind $\overline{{\upsilon }_z}=0$, and defining the correlation time τ_corr of the turbulent eddies such that ∂δχ/∂t ≃ δχ/τ_corr, we have then$$\begin{array}{c}\hfill \delta \chi \simeq -{\tau }_{\text{corr}}\delta {\upsilon }_z{\displaystyle \frac{\partial \overline{\chi }}{\partial z}}.\end{array}$$(B.13)

We finally use the fact that δρ/ρ ≪ 1 (anelastic approximation, so that $\delta \chi \simeq \delta {\rho }_d/\overline{\rho }$) to obtain a closure relation that links the turbulent term with the vertical gradient of $\overline{\chi }$:$$\begin{array}{c}\hfill \overline{\delta {\rho }_d\delta {\upsilon }_{z_d}}=-D_z\overline{\rho }{\displaystyle \frac{\partial }{\partial z}}\left({\displaystyle \frac{\overline{{\rho }_d}}{\overline{\rho }}}\right),\end{array}$$(B.14)

with $D_z\simeq \overline{{\upsilon }_z^2}$τ_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 $\overline{{\rho }_d}\overline{{\upsilon }_z}$ in Eq. (B.1), it is straightforward to show that:$$\begin{array}{c}\hfill {\displaystyle \frac{\partial \overline{{\rho }_d}}{\partial t}}={\displaystyle \frac{\partial }{\partial z}}\left(z{\mathrm{\Omega }}^2\overline{{\tau }_s}\overline{{\rho }_d}\right)+{\displaystyle \frac{\partial }{\partial z}}\left[D_z\overline{\rho }{\displaystyle \frac{\partial }{\partial z}}\left({\displaystyle \frac{\overline{{\rho }_d}}{\overline{\rho }}}\right)\right].\end{array}$$(B.15)

Appendix C:

The 2D advection-diffusion solver

In Sect. 5.2, we study the 2D solutions ρ_d(x, z) that satisfy the following equation:$$\begin{array}{c}\hfill {\displaystyle \frac{\partial ({\rho }_d{\upsilon }_{x_d})}{\partial x}}+{\displaystyle \frac{\partial ({\rho }_d{\upsilon }_{z_d})}{\partial z}}={\displaystyle \frac{\partial }{\mathrm{d}z}}\left[D_z\rho {\displaystyle \frac{\partial }{\partial z}}\left({\displaystyle \frac{{\rho }_d}{\rho }}\right)\right]+{\displaystyle \frac{\partial }{\mathrm{d}x}}\left[D_x{\displaystyle \frac{\partial {\rho }_d}{\partial x}}\right],\end{array}$$(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 re-written in the form:$$\begin{array}{c}\hfill A{\displaystyle \frac{{\partial }^2{\rho }_d}{\partial x^2}}+B{\displaystyle \frac{\partial {\rho }_d}{\partial x}}+C{\displaystyle \frac{{\partial }^2{\rho }_d}{\partial z^2}}+D{\displaystyle \frac{\partial {\rho }_d}{\partial z}}+E{\rho }_d=0,\end{array}$$(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 second-order centred differences. By defining a vector state X containing the values of ρ at each grid points, the discretized advection-difffusion 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 ill-conditioned. 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 Runge-Kutta 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 time-step 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