Global MHD simulations of stratified and turbulent protoplanetary discs. II. Dust settling
Key Words.:
Accretion, accretion discs  MHD  Methods: numericalAbstract
Context:
Aims:The aim of this paper is to study the vertical profile of small dust particles in protoplanetary discs in which angular momentum transport is due to MHD turbulence driven by the magnetorotational instability. We consider particle sizes that range from approximately micron up to a few millimeters.
Methods:We use a grid–based MHD code to perform global twofluid simulations of turbulent protoplanetary discs which contain dust grains of various sizes.
Results:In quasi–steady state, the gravitational settling of dust particles is balanced by turbulent diffusion. Simple and standard models of this process fail to describe accurately the vertical profile of the dust density. The disagreement is larger for small dust particles (of a few microns in size), especially in the disc upper layers (, where is the scaleheight). Here there can be orders of magnitude in the disagreement between the simple model predictions and the simulation results. This is because MHD turbulence is not homogeneous in accretion discs, since velocity fluctuations increase significantly in the disc upper layer where a strongly magnetized corona develops. We provide an alternative model that gives a better fit to the simulations. In this model, dust particles are diffused away from the midplane by MHD turbulence, but the diffusion coefficient varies vertically and is everywhere proportional to the square of the local turbulent vertical velocity fluctuations.
Conclusions:The spatial distribution of dust particles can be used to trace the properties of MHD turbulence in protoplanetary discs, such as the amplitude of the velocity fluctuations. In the future, detailed and direct comparison between numerical simulations and observations should prove a useful tool for constraining the properties of turbulence in protoplanetary discs.
1 Introduction
Thanks to the Spitzer Space Telescope, the last few years have seen a dramatic improvement in our knowledge of dust emission features arising at midinfrared wavelengths from protoplanetary discs surrounding brown dwarfs, T Tauri stars and Herbig Ae/Be stars. The properties (size, composition) of these grains can be studied in great detail by analyzing the shape and strength of these emission features. The main result is that dust particles have been significantly processed compared to their interstellar medium cousins: grains are bigger (up to a few microns) and, for reasons not yet fully understood, crystallinity appears to be common among all observed spectral types (Apai et al. 2005; KesslerSilacci et al. 2006; Furlan et al. 2006). Because protoplanetary discs are optically thick, this midinfrared emission mostly arises from the upper layers of the disc (van Boekel et al. 2003; Dullemond & Dominik 2004, 2008), in the so–called “superheated” layer of Chiang & Goldreich (1997), and is attributed to the disc inner radii. In the case of T Tauri stars of solartype, for example, the dust emission zone lies within a few tenths of an astronomical unit (AU) of the central star (KesslerSilacci et al. 2007).
The presence of solid particles high above the disc midplane is the signature of the turbulent nature of the flow in protoplanetary discs. Turbulent velocity fluctuations lift up dust particles that would otherwise settle toward the disc midplane (Dubrulle et al. 1995; Dullemond & Dominik 2004) because of the vertical component of the central star’s gravitational potential. The vertical profile of the dust density is thus determined by the balance between gravitational settling toward the equatorial plane and upward lift due to turbulence. Its typical scaleheight is a function of grain size. Small particles are well coupled to the gas, and can be transported to higher altitudes (or smaller gas densities) above the equatorial plane than larger particles. Although the issue is still under discussion (Dullemond & Dominik 2008), this differential settling process can be inferred from a detailed analysis of the observations. This is the case for observations carried out by Spitzer alone (Furlan et al. 2005; SiciliaAguilar et al. 2007), by combining Spitzer observations with observations at other wavelengths (Pinte et al. 2008) or even using completely different observational strategies (Pinte et al. 2007; Rettig et al. 2006). Dust gravitational settling is also known to affect the spectral energy distribution of protoplanetary discs, as shown for example by Dullemond & Dominik (2004) or D’Alessio et al. (2006).
In general, such observational diagnostics of gravitational settling are inferred using simple parametric prescriptions of the effect of turbulence or of its consequence. Pinte et al. (2007, 2008) assume that the vertical profile of the dust density is a Gaussian (as is the case for the gas in the isothermal limit), leaving open the possibility that the dust disc scaleheight depends on the size of the particles. The point of their analysis is to demonstrate that this is indeed the case in GG Tau and IM Lupi. Furlan et al. (2005) use the models developed by D’Alessio et al. (2006) and assume a constant dust–to–gas ratio for two populations of solids, one consisting of small and well mixed particles, the other composed of large particles. Using Spitzer spectra, they demonstrate that the latter is depleted in the disc upper layers for the vast majority of a sample of 25 stars in Taurus. Dullemond & Dominik (2004) use a more physical approach in which turbulence is modeled as a diffusive process with a spatially constant diffusion coefficient (see also Dubrulle et al. 1995; Schräpler & Henning 2004). Using this formulation, one can derive analytic expressions for the dust density vertical profile. This is also the approach followed by Rettig et al. (2006) who used the strong settling (thin dust disc) limit of these formulae.
All of these approaches are very useful since they provide analytical formulae or a small number of model parameters to fit. They can thus be incorporated into radiative transfer tools, which can then generate a large grid of models to be used for interpreting the observations. All of them establish, at least qualitatively if not quantitatively, the basic result that differential gravitational settling occurs in protoplanetary discs. However, none of these derivations follows directly from first principles, but instead rely at best on a set of unchecked assumptions concerning the properties of the underlying turbulence. The origin of the latter, however, is thought to be magnetohydrodynamic in nature and driven by the magnetorotational instability (MRI, Balbus & Hawley 1991, 1998). It is now possible to perform global numerical simulations of turbulent protoplanetary discs, including solid particles, and to use these simulations as a test of the models described above. This is the purpose of the present paper. Despite the recent increase in computational resources, it is important to keep in mind that such simulations are still extremely challenging and need to use as simplified a setup as possible. This is not without consequences for their realism. Important issues related to the MRI, such as small scale dissipation (Fromang et al. 2007; Lesur & Longaretti 2007) and the possible presence of a dead zone (Gammie 1996) still have to be ignored. We will return to these issues in the last section of the paper when discussing the limits of our work.
Of course, the effect of MRI–induced MHD turbulence on dust dynamics and dust vertical settling itself has already been studied in numerical simulations, but this is the first time that both effects are incorporated in a single global simulation that takes vertical stratification into account. Barrière–Fouchet et al. (2005) performed global simulations of dust settling but neglected the effect of MHD turbulence. This prevents the system from reaching a quasi–steady state. Other global simulations, taking MRI–driven turbulence into account and including dust particles, neglected gaseous vertical stratification (Fromang & Nelson 2005) and instead focussed on the radial migration of larger bodies. Lyra et al. (2008) also neglected the vertical component of gravity acting on the gas, but included its effects on the dust particles. While this approach produces settling of the dust particles toward the midplane, its neglects the spatial inhomogeneity of the turbulence induced by vertical stratification of the gas. We shall see in the course of this paper that the latter is important when considering dust settling of small particles. The other published numerical simulations studying the effect of turbulence on dust dynamics were local simulations that use the shearing box model (Goldreich & LyndenBell 1965). A large number of them neglected density stratification in the vertical direction (Johansen et al. 2006a; Carballido et al. 2005; Johansen et al. 2006b) and thus only studied dust diffusion in an homogeneous environment. Others (Fromang & Papaloizou 2006; Carballido et al. 2006) included vertical density stratification but considered particles larger than one millimeter. Here, we want to concentrate on smaller particles that produce an observational signature at mid–infrared wavelengths. We note that such a simulation could in principle be done in the framework of the shearing box model. Indeed, this was done recently by Balsara et al. (2008), although the vertical extent of their shearing box is smaller than the simulations we present in this paper and thus less appropriate to study the dust distribution in the disc corona (i.e. above three scaleheights). However, since this work is intended to mark the beginning of an effort to compare observations and numerical simulations directly, it makes more sense to compute global models as these will be more readily comparable with the observations as they become more realistic. Our strategy in this paper will be simple. We will use exactly the setup presented by Fromang & Nelson (2006) for global simulations of turbulent and stratified protoplanetary discs. Dust particles of various sizes will be added to the disc and their subsequent evolution will be analyzed and compared with simple models of dust stratification in protoplanetary discs.
The plan of the paper is as follows. In Sect. 2, we introduce the model we use for the disc as well as convenient dimensionless parameters that appear in the problem. The relationship between these quantities and physical parameters in real systems will be outlined. In Sect. 3, we describe more quantitatively the different ways to model the quasi–steady state dust distribution resulting from the balance between dust settling and turbulent diffusion. These models are then compared with the results of our numerical simulations in Sect. 4. Finally, in Sect. 5 we discuss the implications and limitations of our work, and point the way toward future improvements.
2 Definitions
In this section, we describe the general properties of the disc model and the dust parameters we used. We also introduce the mathematical notation that will be required in the following sections.
2.1 Coordinate systems
2.2 Disc model
We consider a disc extending in radius between an inner radius and an outer radius . For simplicity and computational reasons we define the initial disc structure using straightforward analytic functions. The equation of state is locally isothermal: the sound speed, , only depends on and is constant in time. Both and the disc midplane gas density, , obey power laws
(1)  
(2) 
where and are the sound speed and the midplane gas density at a radius , respectively. The disc is initially axisymmetric and in radial and vertical hydrostatic equilibrium. The spatial distribution of density and angular velocity can thus be approximated by
(3)  
(4) 
in which . The disc scaleheight, , is given by
(5) 
where is the disc scaleheight at .
2.3 Integrated quantities
The density distribution can be used to work out the surface density of the disc:
(6) 
where the second relation serves as a definition for . The total disc mass follows by radially integrating the surface density between and :
(7) 
where we have assumed . This relation can be used to express the disc surface density at as a function of the disc parameters and :
(8) 
2.4 Dust size
In this paper, we shall study the effect of MHD turbulence on the dust. Gas affects solid body dynamics through the drag force it exerts on the dust particles. In the Epstein regime we are interested in, this force takes the simple form
(9) 
where and are the gas and dust velocities, respectively. is the dust stopping time. This is the typical time it takes for dust particles initially at rest to reach the local gas velocity. It depends on the dust particle mass density and its size through
(10) 
A relevant dimensionless parameter in the problem is the quantity . It compares the stopping time to the dynamical time. When , the stopping time is much smaller than the orbital period and the dust essentially follows the gas. When or larger, becomes comparable to and dust and gas start to decouple. As pointed out by Dullemond & Dominik (2004), this occurs at all radii for a given particle size provided is large enough. Using the disc parameters introduced above, can be expressed as a function of position according to
(11) 
where the parameter is the value of at in the disc midplane. It can be expressed in terms of the disc surface density at using Eq. (6):
(12) 
Using this expression along with Eq. (8), it is possible to express the dust size in term of , the disc parameters (mass and radius) and :
(13) 
2.5 Converting to physical units
Eq. (8) and (13) can be used to convert the dimensionless quantities describing the problem into physical quantities. When doing so, the numerical simulations we will describe in Sect. 4 should be thought of as covering a small fraction of the total disc, going from to an outer radius .
The dimensionless parameters describing the disc and dust particles, and the results of the numerical simulations presented below, can be rescaled to any physical system upon specifying the disc mass, the outer radius of the disc, and the value of (in astronomical units). For example, the disc surface density can be written as
(14) 
Likewise, taking gcm and using Eq. (13), we obtain an expression for the dust size:
(15) 
In Sect. 4, we will describe the results of three simulations. They are characterized by , and . In a disc of mass and AU in size, they would correspond to dust particles of size , and m, respectively, if we take AU. (This would mean that we consider the simulation to cover the radial extent to AU.) For the same disc mass and size, if we now take AU (i.e. we consider the simulation to cover radii ranging from to AU), the three sizes are , and m. Note, however, that these numbers are only illustrative. In general, for a given value of , the size of the dust particles decreases when the disc mass is decreased or its outer radius is increased.
3 Dust settling in turbulent discs
As pointed out in the introduction, a steady state is reached in a turbulent protoplanetary disc in which turbulent fluctuations oppose and balance against dust settling. In this section, we describe three approaches that can be used to model the vertical profile of the dust density.
3.1 A Gaussian profile
The simplest approach, and one that is commonly used when attempting to interpret observations (Pinte et al. 2008), is to assume that the dust density follows a Gaussian distribution, as the gas does, but with a different vertical scaleheight :
(16) 
where is the dust midplane density. In this approach, is different for each dust particle size . For example, while trying to model the dust properties in the protoplanetary disc orbiting the M dwarf IM Lupi, Pinte et al. (2008) found . One purpose of this paper is to establish whether such a description is supported by numerical simulations of turbulent protoplanetary discs.
3.2 Turbulence as a diffusive process
A more physical approach is to describe the transport of the dust particles by the turbulent fluctuations as a diffusion process. This has been used commonly in the literature (Dullemond & Dominik 2004). In this approach, the vertical evolution of the dust density can be described by the following partial differential equation (Schräpler & Henning 2004; Dubrulle et al. 1995):
(17) 
where is the dust particle density and is a diffusion coefficient that quantifies the turbulent diffusivity. This equation models the balance between vertical settling and turbulent diffusion. When looking for a steady state vertical profile for the density, the time derivative vanishes. Upon integrating once and rearranging terms, Eq. (17) gives
(18) 
The vertical integration of the last equation requires the knowledge of the diffusion coefficient as a function of . The simplest solution is to assume that it is constant. This is the approach described in the following subsection, while in Sect. 3.2.2 we outline a possible alternative.
A constant diffusion coefficient
When the dust diffusion coefficient is constant, Eq. (18) can be integrated to give
(19) 
where is a dimensionless diffusion coefficient defined as . The quantities and only depend on and are to be evaluated in the disc midplane. Note that, while deriving the last equation, we have assumed that the vertical distribution of gas remains Gaussian at all times, in agreement with local shearing box numerical simulations of the MRI (Miller & Stone 2000).
A common practice in this context is to express in units of the standard parameter introduced by Shakura & Sunyaev (1973). is then written as follows (Dullemond & Dominik 2004; Schräpler & Henning 2004):
(20) 
where Sc is the Schmidt number. In zero net flux MHD turbulence, the Schmidt number has been measured to be of order unity in local simulations of unstratified (Johansen & Klahr 2005; Johansen et al. 2006b) or stratified discs (Fromang & Papaloizou 2006; Ilgner & Nelson 2008). In the present paper, we will tune the Schmidt number in order to obtain the best agreement with the numerical simulations.
A vertically varying diffusion coefficient
Dust particles are diffused away from the disc midplane by the turbulent velocity fluctuations. Thus, the dust diffusion coefficient is intimately linked to the turbulence properties and particularly to the gas velocity fluctuations. It would then seem natural for to be constant in space if the turbulence was homogeneous. However, because of the vertical stratification, MHD turbulence is not homogeneous in protoplanetary discs and it is very likely that varies with at a given radius (even in the absence of a dead zone). Using local vertically stratified simulations of the MRI, Fromang & Papaloizou (2006) showed that the following simple expression gives a fairly good estimate to the numerically derived diffusion coefficient:
(21) 
In this equation, stands for the turbulent velocity fluctuations and is the correlation time of these fluctuations. Numerical estimates of both terms thus provide a path to calculating the value of . Their vertical variations will be investigated in Sect. 4.
Of course, the drawback of this approach is that it becomes impossible to explicitly integrated Eq. (18). We cannot provide an analytical expression for the vertical profile of the dust density and shall rely on a numerical integration once the vertical profile for is extracted from the numerical simulations.
4 Numerical simulations
4.1 Setup
The simulations presented in this paper are run using the code GLOBAL (Hawley & Stone 1995), which solves the ideal MHD equations using a spherical coordinate system as defined in Sect. 2.1. The setup we used is exactly that of model S5 in Fromang & Nelson (2006). Here, we describe it only briefly.
At the start of the simulation, the disc model presented in Sect. 2.2 is initialized on the grid. The units are such that , and (i.e. throughout the disc). The computational domain covers the range to in radius and the interval in . In the direction, the grid extends to scaleheights on both sides of the disc equatorial plane. The resolution for each simulation is . Following Fromang & Nelson (2006), a weak toroidal magnetic field is added to the disc, and small random velocity perturbations are also imposed. Time is measured in units of the orbital period at the inner edge of the computational domain in the following sections.
Because of the MRI, the presence of a weak magnetic field, together with the velocity perturbations, begins to drive MHD turbulence and angular momentum transport through the disc within a few orbits of the simulations starting. To reach a meaningful quasi steady state, however, the model is first evolved for orbits without dust particles. At that stage, the gas density is reset to its initial distribution and dust particles are introduced such that the dusttogas ratio is uniform through the computational domain (note that we neglect the back reaction of the solids onto the gas, so that the value of this ratio has no physical consequences). We ran three simulations for three different particle sizes. The three values of associated with these sizes are , and . All the simulations are further integrated for about orbits until the dust distribution itself reaches a steady state in which gravitational settling is balanced by turbulent diffusion. Examples of the disc structure at the end of such a run are illustrated in Fig. 1. Two snapshots of the gas (left panel) and dust (right panel) density in the plane are shown in the case at time . The dust disc appears thinner than the gas disc, indicating that significant settling has occurred.
In the following subsections, we will compare the dust distributions we obtained for the different sizes to the models described in Sect. 3. To do so, all relevant physical quantities will be averaged in time between and using snapshots so that they become statistically significant. We first start by describing the relevant properties of the turbulence that result from this procedure before concentrating on the degree of settling as a function of size.
4.2 Turbulence properties
The first relevant property of the turbulence is the vertically averaged angular momentum transport it generates. It is commonly measured using the parameter . Following Fromang & Nelson (2006), we calculated as a function of radius according to
(22) 
where and correspond respectively to the Reynolds and Maxwell stress contributions to . The overbar symbols denote densityweighted azimuthal and vertical averages (see Eq. (6) of Fromang & Nelson 2006). To reduce statistical noise, was further averaged in time between and . Its radial profile is shown in Fig. 2 for the model having . The dashed and dotted lines, respectively, show the variations of and , while the solid line represents , the sum of the two. As described in Fromang & Nelson (2006), presents large oscillations in space and time, but its time averaged value is well behaved, varying between and , in agreement with the results of Fromang & Nelson (2006). As is usually obtained in numerical simulations of this type, the Maxwell stress dominates over the Reynolds stress by a factor of about two to four. The volume averaged value of is also in agreement with the results of Fromang & Nelson (2006). Indeed, we obtained , , , , and respectively at times , , , , and . For the duration of the simulation, the turbulence is clearly in a quasi steady state.
In this paper, we present three simulations for different particle sizes. Because these simulations were obtained under different computing setups (i.e. using different numbers of CPUs), the details of the turbulent flows differ from one run to another. This is simply due to the chaotic nature of turbulence (Winters et al. 2003). However, the statistical properties of the turbulence are similar. This is shown on Fig. 3 where we compare the radial profile of for the models (solid line), (dotted line) and (dashed line). For each case, the curves are averaged in time between and . Fig. 3 demonstrates that we obtain very similar values of in the different simulations. Thus it is meaningful to compare the dust distributions in the three different models.
As explained in Sect. 3.2.2, the vertical velocity fluctuations and the correlation timescale of the turbulence are also of importance in affecting the diffusion of solid particles. In Fig. 4, we show the vertical profile of the vertical velocity fluctuations at , normalized by the speed of sound and averaged in time between and for the models (solid line), (dotted line) and (dashed line). Again the results are very similar for each model and in agreement with those of Fromang & Nelson (2006): the average is of order of the speed of sound in the disc midplane before rising up to values of the order of 20 – 30% in the disc corona, where weak shocks develop in locations where convergent flow speeds exceed the sound speed (Fromang & Nelson 2006). These supersonic turbulent motions are driven by magnetic stresses in regions where the Alfvén speed exceeds the sound speed. Fig. 4 shows that varies by a factor of about between the disc midplane and its corona.
As expressed by Eq. (21), the dust diffusion coefficient also depends on the turbulence correlation timescale . Although depends on the value of to the first power only, while it depends on the velocity fluctuations to the second power, vertical variations in the correlation timescale could still affect the numerical estimate of the diffusion coefficient. We thus calculated the value of at two locations, one in the disc midplane and the other in the disc corona. Following Fromang & Papaloizou (2006), can be evaluated by monitoring the time variation of the function
(23) 
where denotes an ensemble average. is expected to decrease toward zero from initially positive values in a time . To estimate the later, the model in which was restarted at time . We then calculated and averaged the function at seven different radii , , , , , and over a kernel of five cells in the radial direction. At all radii, the function was further averaged over two ranges in the meridional direction: to produce the function , and to produce the function . thus represents the velocity correlation function near the disc midplane while represents the correlation function in the upper layers including the corona. The two functions are plotted as a function of in Fig. 5, respectively, with a solid and a dashed line (both curves are normalized by their value at ). As expected, both display an initial decrease toward zero around which they stabilize after a few tenths of an orbit. This qualitative trend is in agreement with the results of Fromang & Papaloizou (2006). The dotted curves over plotted on the same figure represent the functions for , , and orbits. Although it is difficult to measure the correlation timescale precisely given the large fluctuations we obtained, these curves can still be used as a way to estimate the correlation time governing the functions and . They indicate that is in the corona and in the disc midplane. Both values are comparable to the value of orbit reported by Fromang & Papaloizou (2006). Although these two timescales are different, their ratio is at most two and certainly accounts for less variation in the diffusion coefficient than the variations in the vertical velocity fluctuations described above. Therefore, when plugging numerical estimates for the correlation time into Eq. (21), we will use orbits in the remaining of this paper. This value is such that and is in agreement with previously reported results in the literature (Johansen et al. 2006b) and with values used to model the effect of turbulence in theoretical studies of planet formation (see for example Weidenschilling 1984), thus giving additional support to such work.
4.3 Dust spatial distribution
For the three simulations we performed, we averaged the dust density in azimuth and over time between and orbits. The spatial distributions we obtained using this procedure are shown in Fig. 6 for the models having (left panel), (middle panel) and (right panel). Note that the radial extent of the snapshots is limited to . At larger radii, the outer buffer region we use (see Fromang & Nelson 2006) starts to affect the dust distribution. As expected, the smaller the dust particle sizes, the thicker the dust disc. This is simply because smaller particles are better coupled to the gas and can thus be lifted further away from the disc midplane by the turbulence before they decouple from the gas. The left panel also shows some sign of the dust disc flattening at large radii. This is simply because the strength of the turbulence, as measured for example by the parameter , decreases at large radii, as seen in Figs. 2 and 3. It should not be confused with the apparent dust disc flattening identified by Dullemond & Dominik (2004), which occurs because of self–shadowing in the presence of weak turbulence. Finally, Fig. 6 also highlights one of the limitations of our work: with the setup of the simulations presented in this paper, we cannot easily extend our parameter survey to smaller particles. Indeed, when , the dust disc already covers almost the entire computational domain in the vertical direction. Reducing further the size of the dust particles would require an increase in the size of the computational domain in the meridional region for the decoupling between gas and dust to occur within the computational grid. This would require an increase in the computing time required for a simulation performed with the same resolution. This will soon become possible as computational resources improve, but is currently very challenging.
To compare these results with the models presented in Sects. 3.1, 3.2.1 and 3.2.2, we computed the expected 2D dust distribution that each of them would predict. In doing so, we used the same disc and dust parameters as in the simulations.

For the first model (presented in Sect. 3.1), we fitted a Gaussian profile to the vertical profile of the dust density at each radius. This was done by performing a least squares fit to the dust density profile over the entire vertical extent of the disc. The resulting spatial distributions are plotted in Fig. 7 using the same color table and spatial domain as in Fig. 6.

For the second model (presented in Sect. 3.2.1), we used Eq. (19) to calculate the vertical profile of the dust density at each radius. The numerical value of the dimensionless diffusion coefficient at each radius was calculated using Eq. (20) in which we plugged the value of calculated according to Eq. (22). We used as this value turns out to provide the best fit to the simulations. Being of order unity, it is also in agreement with the results of previous local numerical simulations (Johansen & Klahr 2005; Johansen et al. 2006b; Fromang & Papaloizou 2006). The resulting spatial distributions of the dust density are plotted in Fig. 8.

For the third model (presented in Sect. 3.2.2), we integrated Eq. (18) numerically, using at each disc altitude the estimate for the diffusion coefficient provided by Eq. (21). In this equation, we used the azimuthally and time averaged vertical profile of the vertical velocity fluctuations at each radius (as shown on Fig. 4 for the special case ) and the correlation time orbit as explained in Sect. 4.2. The resulting spatial distributions of the dust density are plotted in Fig. 9.
The spatial distributions of the dust density obtained with the three different models are in rough agreement with the simulations and with naive expectations: all predict that the dust component of the discs thicken as dust particles decrease in size, in agreement with the MHD simulations. For the model having (i.e. for the largest particles), the agreement between the three models and the simulation is quite good. This can be understood easily: in this model, the dust scaleheight is about . In other words, the solid particles concentrate close to the equatorial plane, where the turbulence properties are fairly homogeneous (the vertical velocity fluctuations only start to rise significantly above ). Thus the diffusion coefficient calculated according to Eq. (21) is roughly constant and the second and third models (with constant and varying diffusion coefficient) give a similar result. Moreover, close to the equatorial plane (i.e. ), the leading order expansion to Eq. (19) turns out to be Gaussian, which explains the similarities between the three models for large particles.
On the other hand, models having smaller dust sizes, i.e. and , show differences between the different approaches. The Gaussian fit to the simulations always overestimates the dust density in the disc corona (), showing that in general the vertical profile of the dust density distribution is non Gaussian. On the other hand, the model having a constant diffusion coefficient always underestimates the dust density in the disc corona. Only the model taking the vertical variation of the diffusion coefficient into account gives a satisfactory fit to the simulations, especially in the disc upper layers.
These differences can be made more quantitative by comparing the vertical profile of the dust density between the simulations and the models. This is done in Fig. 10. The left panel gathers the results corresponding to the case , the middle panel shows results obtained when . Finally, the right panel shows results obtained when . In each panel, the solid curves plot the vertical profile of the dust density obtained in the MHD simulation by temporally, azimuthally and radially averaging the results. The radial averaging is performed over the range . The dot–dashed, dashed and dotted curves, respectively, correspond to the model using a Gaussian fit to the data, a constant diffusion coefficient and a vertically varying diffusion coefficient. Again, the left panel demonstrates the good agreement between all approaches in the case of large particles, as the curves used to represent the results of the different models are almost undistiguishable on that plot. It also provides, however, a first hint that the Gaussian fit increasingly overestimates the dust density as one moves away from the midplane, even for these larger particles. The middle and right panels confirm these results and show that the model having a vertically varying diffusion coefficient provides in general the best fit to the data in the disc corona. For example, in the case , the three models give a good fit to the data for . For , the Gaussian fit to the data overestimates the density obtained in the simulations (their ratio is about at ) and the model having a constant diffusion coefficient underestimates that density (the ratio at is in excess of ). The agreement between the simulations and the model having a vertically varying diffusion coefficient is better as the ratio between the predicted and observed dust densities is always bounded by , despite the fact that the dust density itself varies by more than orders of magnitude. It is interesting to point out, though, that the value of the density predicted by this model close to the equatorial plane seems to significantly underestimate the results of the simulations. This is most likely because the correlation timescale we used in calculating the diffusion coefficient underestimates the midplane correlation timescale of the turbulence (see section 4.2). Using a varying correlation timescale would certainly further improve the agreement with the numerical simulations, but such a level of refinement is probably meaningless given the other approximations involved in the simulations themselves. As shown on Fig. 10, the situation is similar in the case : the best fit to the numerical simulations is provided by the final and more elaborate model. Indeed, at , the dust density predicted by the model with a constant diffusion coefficient starts to underestimate significantly the results of the numerical simulations. Above , the difference becomes enormous, as in the case . The agreement, however, seems better with the model that uses a Gaussian fit to the data than in the case . The Gaussian fit indeed gives a good fit up to . This apparent agreement would most probably break down for , as the ratio with the simulated densities is already greater than two orders of magnitude at and increases with increasing .
5 Discussion and conclusions
In this paper, we have studied dust settling in turbulent protoplanetary discs using global MHD numerical simulations performed with the code GLOBAL using spherical coordinates. In this section, we summarize our main findings and discuss the limitations and future prospects of our work.
5.1 Dust density vertical profile
The first point that emerges is that Gaussian profiles fail badly to reproduce the data extracted from the simulations. This point is further illustrated by Fig. 11 where we compare the dust profile in the case (solid line) with a set of Gaussian profiles (dotted lines) computed according to the following equation,
(24) 
in which we used the values , , , , and . None of the dotted lines is an acceptable fit to the data. Close to the disc midplane, by eye inspection seems to suggest that the dust disc scaleheight is while one would estimate it to be of order at . Both values are different from the dust disc scaleheight returned by the least squares fit to the data, . The physical reason for this behaviour is clear: close to the midplane, where the gas density is high, the gasdust coupling is strong and the dust traces closely the Gaussian profile of the gas, whereas in the disc upper layers the coupling is weak and the dustgas ratio decreases as one moves further away from the midplane. These results emphasize the point that estimates of the dust disc scaleheight obtained using a Gaussian fit may lead to incorrect conclusions.
Nevertheless, in order to compare with results published previously in the literature, we measured the dust disc scaleheight in our simulations using such a fit. We found , and , respectively, for , and . The variation of as a function of is shown by the diamonds in Fig. 12. The dotted line displays the function
(25) 
which is the best fit to the data. Therefore, if we were to analyse our simulations assuming that the dust density is Gaussian, we would obtain since . Interestingly, this is not too different from the results of Pinte et al. (2008) who report an exponent equal to , although it is clear that our results show a stronger relationship between the dust disc scaleheight and particle size. Both values are, however, largely different from the value which is often reported in the literature (Dubrulle et al. 1995; Carballido et al. 2006). This is because the latter is obtained when solving Eq. (17) in the strong settling limit that is mostly relevant for large particles. In this case, the vertical variation of can be neglected, the vertical profile is Gaussian and the exponent is recovered. For the small particles we study here, the dust disc is thick and these vertical variations have to be taken into account, which leads to the more complicated expression given by Eq. (19) and departure from a Gaussian profile.
5.2 The Schmidt number
The second result that emerges from our work is that a model having a
constant diffusion coefficient increasingly underestimates the dust
density as the particle size is decreased. In other words, the
vertically averaged dust diffusion coefficient decreases with
. This is not unexpected, since the Schmidt number introduced in
Sect. 3.2.1 is known to be an increasing function of
(Cuzzi et al. 1993; Schräpler & Henning 2004; Youdin & Lithwick 2007). It is
unclear, however, whether such studies,
which assume homogeneous Kolmogorov–like turbulence,
are applicable to the highly magnetised flow of the corona. This is why we
did not attempt to make a direct comparison between our results and
these theories. For the sake of completeness, it is nevertheless
instructive to report here the vertically averaged Schmidt number we
measured in each case. As described above, already
provides a good fit of
the dust density in the case . For the cases
and , we found that the vertically
averaged Schmidt numbers that best fit the data are respectively
and
5.3 A toy model
The main result of this paper is the construction of a simple model that gives a reasonable fit to the simulations. In this model, the dust particles are diffused away by turbulence with a diffusion coefficient that scales with the square of the turbulent velocity fluctuations. Accordingly, it should be possible in principle to extract the vertical profile of the velocity fluctuations from the dust density vertical profile. This is, however, an inversion problem. As such, it is susceptible to being degenerate and we shall see in this section that this is indeed the case.
The relevant question to ask in this context is the following: provided we are able to measure the dust density vertical profile, what is it mostly sensitive to? Can we hope to constrain the velocity fluctuations in the disc midplane or is it mostly a consequence of the amplitude of the fluctuations in the upper layers. To answer that question, we designed the following toy model: guided by the vertical profile of the velocity fluctuation amplitudes shown in Fig. 4, we considered the following analytic vertical profile for the velocity fluctuations:
(26) 
where and stand for the turbulent velocity fluctuations in the disc midplane and in the disc corona. As shown on Fig. 4, typical numerical values are and . To investigate the sensitivity of the results to the midplane velocity fluctuations, we calculated the dust density vertical profile in the case by numerically integrating Eq. (18), using Eq. (21) and (26) with , and and (i.e. we used in the disc corona the value suggested by the simulation data) . The results are summarized on the upper panel of Fig. 13, where the dust density profile in the case is shown with the solid line, while the numerically integrated profiles are represented by dotted, dashed and dot–dashed lines for , and , respectively. The last three curves are very similar in this plot and all give a fairly good fit to the simulations, especially in the disc upper layers. This shows that the dust density vertical profile is fairly insensitive to the midplane velocity fluctuations. To estimate the sensitivity of the dust density profile to the velocity fluctuations in the upper layers, we repeated the same analysis using (i.e. we used in the disc midplane the value suggested by the simulation results) and , and . The results are shown on the bottom panel of Fig. 13 with the same conventions as for the upper panel. In this case, the dust density vertical profile is seen to be much more sensitive to the upper layer velocity fluctuations and only the dashed curve, for which (i.e. in rough agreement with the numerical data), is in good agreement with the data.
The implications of these results are twofold. First, the simple toy model described by Eq. (26) would be suitable to use when trying to fit the observations as it reproduces the numerical data fairly well if we choose the values and , compatible with the simulations. But these results also show that such a fit would only provide sensitive information about the turbulent velocity fluctuations in the disc upper layers. The physical reason for this last point is that such a fit is mostly sensitive to the properties of the region where the gas and dust decouple. For the small particles studied in this paper (and observed using the Spitzer telescope), this region turns out to lie in the disc upper layers. When observing larger particles at longer wavelengths (for example with ALMA), it will become possible to constrain the turbulent velocity fluctuations of the disc closer to the midplane as such particles will settle and decouple deeper in the disc.
5.4 Limitations and future prospects
Of course, there are strong limitations to our work due to the complex and CPUintensive nature of global simulations of protoplanetary discs. On the purely numerical side, the limited resolution we used is of course an issue, as was pointed recently in a number of studies (Fromang & Papaloizou 2007; Simon et al. 2008). Proper simulations should include explicitly microscopic diffusion coefficients (viscosity and resistivity), as the latter have been shown to be important in determining the saturation level of the turbulence (Fromang et al. 2007; Lesur & Longaretti 2007). However, the resolutions required to include these processes in global simulations are currently out of reach and one must instead rely on the subgrid model provided by numerical dissipation to carry out global numerical simulations. On a more physical side, there are also limitations due to the simple disc model we used. The locally isothermal equation of state we used is not appropriate for the disc inner parts that we are simulating as the gas there is optically thick. This could have numerous effects. For example, Dullemond (2002) and D’Alessio et al. (1998) report a temperature increase in the inner disc upper layers. Such an increase would strengthen the coupling between gas and grains in the disc corona (through a decrease in the local value of the parameter ). This would cause small particles to settle less than reported in this paper and could change the relationship between and . Obviously, the significance of the comparison we tentatively made between our simulations and the observations of Pinte et al. (2008) should be taken with care. Another topic of concern in our simulations is the assumption of ideal MHD. It is indeed well known that protoplanetary discs are so poorly ionized because of their large densities and low temperatures that parts of the flow, refered to as dead zones, remain laminar (Gammie 1996). We completely ignored the effects of dead zones in the present paper. Clearly, future work should improve the thermodynamic treatment of the gas, possibly including radiative transfer and dead zones.
Nevertheless, we have shown in this paper that dust observations can be used in principle to constrain the properties of MHD turbulence in discs. We have found that even the simplest simulations provide disagreements with previously used fits and diffusion models because of the nature of disc turbulence. This illustrates even further the need to compare directly observations and numerical simulations. It will be important in the future to generate a grid of more realistic discs models (varying the disc parameters, including dead zones, flaring discs, non isothermal discs) and produce synthetic observations that could be compared in the next few years with multiwavelengths observations. This comparison would provide diagnostics of disc turbulence, the existence (or not) of dead zones and thus constrain planet formation models. The recent work of Pinte et al. (2008) shows that such multiwavelengths observations are starting to become feasible. When combined with future instruments like Herschel and ALMA, large samples will become available and will provide a wealth of constraints on disc structure and properties when combined with appropriate global numerical simulations of protoplanetary discs.
Acknowledgments
The simulations presented in this paper were performed using the QMUL High Performance Computing Facility purchased under the SRIF initiative. It is also a pleasure to acknowledge fruitful discussions with J–C.Augereau and P–Y. Longaretti. We also thank an anonymous referee whose comments significantly improved the paper.
Footnotes
 offprints: S.Fromang
 The fact that we obtain Schmidt numbers lower than one, in contrast to Cuzzi et al. (1993), Schräpler & Henning (2004) and Youdin & Lithwick (2007) is due to our definition being different to that used in these theoretical studies, as discussed in detail by Youdin & Lithwick (2007)
References
 Apai, D., Pascucci, I., Bouwman, J., et al. 2005, Science, 310, 834
 Balbus, S. & Hawley, J. 1991, ApJ, 376, 214
 Balbus, S. & Hawley, J. 1998, Rev.Mod.Phys., 70, 1
 Balsara, D. S., Tilley, D. A., Rettig, T., & Brittain, S. A. 2008, ArXiv eprints
 Barrière–Fouchet, L., Gonzalez, J. F., Murray, J. F., Humble, R. J., & Maddison, S. T. 2005, A&A, 443, 185
 Carballido, A., Fromang, S., & Papaloizou, J. 2006, MNRAS, 373, 1633
 Carballido, A., Stone, J. M., & Pringle, J. E. 2005, MNRAS, 358, 1055
 Chiang, E. I. & Goldreich, P. 1997, ApJ, 490, 368
 Cuzzi, J. N., Dobrovolskis, A. R., & Champney, J. M. 1993, Icarus, 106, 102
 D’Alessio, P., Calvet, N., Hartmann, L., FrancoHernández, R., & Servín, H. 2006, ApJ, 638, 314
 D’Alessio, P., Canto, J., Calvet, N., & Lizano, S. 1998, ApJ, 500, 411
 Dubrulle, B., Morfill, G., & Sterzik, M. 1995, Icarus, 114, 237
 Dullemond, C. P. 2002, A&A, 395, 853
 Dullemond, C. P. & Dominik, C. 2004, A&A, 421, 1075
 Dullemond, C. P. & Dominik, C. 2008, A&A, 487, 205
 Fromang, S. & Nelson, R. 2005, MNRAS, 364, L81
 Fromang, S. & Nelson, R. P. 2006, A&A, 457, 343
 Fromang, S. & Papaloizou, J. 2006, A&A, 452, 751
 Fromang, S. & Papaloizou, J. 2007, A&A, 476, 1113
 Fromang, S., Papaloizou, J., Lesur, G., & Heinemann, T. 2007, A&A, 476, 1123
 Furlan, E., Calvet, N., D’Alessio, P., et al. 2005, ApJ, 628, L65
 Furlan, E., Hartmann, L., Calvet, N., et al. 2006, ApJS, 165, 568
 Gammie, C. F. 1996, ApJ, 457, 355
 Goldreich, P. & LyndenBell, D. 1965, MNRAS, 130, 125
 Hawley, J. & Stone, J. 1995, Comput. Phys. Commun., 89, 127
 Ilgner, M. & Nelson, R. P. 2008, A&A, 483, 815
 Johansen, A. & Klahr, H. 2005, ApJ, 634, 1353
 Johansen, A., Klahr, H., & Henning, T. 2006a, ApJ, 636, 1121
 Johansen, A., Klahr, H., & Mee, A. J. 2006b, MNRAS, 370, L71
 KesslerSilacci, J., Augereau, J.C., Dullemond, C. P., et al. 2006, ApJ, 639, 275
 KesslerSilacci, J. E., Dullemond, C. P., Augereau, J.C., et al. 2007, ApJ, 659, 680
 Lesur, G. & Longaretti, P.Y. 2007, MNRAS, 378, 1471
 Lyra, W., Johansen, A., Klahr, H., & Piskunov, N. 2008, A&A, 479, 883
 Miller, K. A. & Stone, J. M. 2000, ApJ, 534, 398
 Pinte, C., Fouchet, L., Ménard, F., Gonzalez, J.F., & Duchêne, G. 2007, A&A, 469, 963
 Pinte, C., Padgett, D. L., Menard, F., et al. 2008, ArXiv eprints, 808
 Rettig, T., Brittain, S., Simon, T., et al. 2006, ApJ, 646, 342
 Schräpler, R. & Henning, T. 2004, ApJ, 614, 960
 Shakura, N. I. & Sunyaev, R. A. 1973, A&A, 24, 337
 SiciliaAguilar, A., Hartmann, L. W., Watson, D., et al. 2007, ApJ, 659, 1637
 Simon, J. B., Hawley, J. F., & Beckwith, K. 2008, ArXiv eprints
 van Boekel, R., Waters, L. B. F. M., Dominik, C., et al. 2003, A&A, 400, L21
 Weidenschilling, S. J. 1984, Icarus, 60, 553
 Winters, W. F., Balbus, S. A., & Hawley, J. F. 2003, MNRAS, 340, 519
 Youdin, A. N. & Lithwick, Y. 2007, Icarus, 192, 588