Subscriber Authentication Point
Free Access
 Issue A&A Volume 496, Number 3, March IV 2009 597 - 608 Astrophysical processes https://doi.org/10.1051/0004-6361/200811220 09 February 2009

## II. Dust settling

S. Fromang1,2 - R. P. Nelson3

1 - CEA, Irfu, SAp, Centre de Saclay, 91191 Gif-sur-Yvette, France
2 - UMR AIM, CEA-CNRS-Univ. Paris VII, Centre de Saclay, 91191 Gif-sur-Yvette, France
3 - Astronomy Unit, Queen Mary, University of London, Mile End Road, London E1 4NS, UK

Received 24 October 2008 / Accepted 15 January 2009

Abstract
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 1 micron up to a few millimeters.
Methods. We use a grid-based MHD code to perform global two-fluid 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 (Z>3H, where H is the scale-height). 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.

Key words: accretion, accretion disks - Magnetohydrodynamics (MHD) - methods: numerical - turbulence

## 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 mid-infrared 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 (Furlan et al. 2006; Apai et al. 2005; Kessler-Silacci et al. 2006). Because protoplanetary discs are optically thick, this mid-infrared emission mostly arises from the upper layers of the disc (van Boekel et al. 2003; Dullemond & Dominik 2008,2004), 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 solar-type, for example, the dust emission zone lies within a few tenths of an astronomical unit (AU) of the central star (Kessler-Silacci 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 scale-height 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 (Sicilia-Aguilar et al. 2007; Furlan et al. 2005), by combining Spitzer observations with observations at other wavelengths (Pinte et al. 2008) or even using completely different observational strategies (Rettig et al. 2006; Pinte et al. 2007). 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 scale-height 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 set-up as possible. This is not without consequences for their realism. Important issues related to the MRI, such as small scale dissipation (Lesur & Longaretti 2007; Fromang et al. 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 & Lynden-Bell 1965). A large number of them neglected density stratification in the vertical direction (Johansen et al. 2006b,a; Carballido et al. 2005) and thus only studied dust diffusion in an homogeneous environment. Others (Carballido et al. 2006; Fromang & Papaloizou 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 scale-heights). 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 set-up 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

In this paper, we will use both cylindrical and spherical coordinate systems. The former will be denoted by , and will be used mostly in the present Sect. 2 and in Sect. 3. Spherical coordinates will be used when we describe the numerical simulations and will use the notation  .

### 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 straight-forward analytic functions. The equation of state is locally isothermal: the sound speed, , only depends on R and is constant in time. Both and the disc midplane gas density, , obey power laws

 (1) (2)

where c0 and are the sound speed and the midplane gas density at a radius R0, 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 scale-height, H, is given by

 (5)

where is the disc scale-height at R0.

### 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 R0 as a function of the disc parameters and R0:

 (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 a 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 a provided Z 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 R=R0 in the disc midplane. It can be expressed in terms of the disc surface density at R0 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 R0:

 (13)

### 2.5 Converting to physical units

Equations (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 R0 (in astronomical units). For example, the disc surface density can be written as

 (14)

Likewise, taking  g cm-3 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 , 10-3 and 10-4. In a disc of mass and 300 AU in size, they would correspond to dust particles of size 163, 16 and m, respectively, if we take R0=1 AU. (This would mean that we consider the simulation to cover the radial extent 1 to 8 AU.) For the same disc mass and size, if we now take R0=0.1 AU (i.e. we consider the simulation to cover radii ranging from 0.1 to 0.8 AU), the three sizes are 500, 50 and m. Note, however, that these numbers are only illustrative. In general, for a given value of R0, 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 scale-height :

 (16)

where is the dust midplane density. In this approach, is different for each dust particle size a. 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 (Dubrulle et al. 1995; Schräpler & Henning 2004):

 (17)

where is the dust particle density and D 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 z. 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.

#### 3.2.1 A constant diffusion coefficient

When the dust diffusion coefficient D is constant, Eq. (18) can be integrated to give

 (19)

where is a dimensionless diffusion coefficient defined as . The quantities and only depend on R 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 (Schräpler & Henning 2004; Dullemond & Dominik 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 et al. 2006b; Johansen & Klahr 2005) or stratified discs (Ilgner & Nelson 2008; Fromang & Papaloizou 2006). In the present paper, we will tune the Schmidt number in order to obtain the best agreement with the numerical simulations.

#### 3.2.2 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 D 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 D varies with Z 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 D. 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 D is extracted from the numerical simulations.

## 4 Numerical simulations

### 4.1 Set-up

 Figure 1: Snapshots of the gas density ( left panel) and the dust density ( right panel) at t=600 for the case . Open with DEXTER

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 set-up 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 GM=1, c0=0.1 and (i.e. H/R=0.1 throughout the disc). The computational domain covers the range to in radius and the interval in . In the -direction, the grid extends to 4.3 scale-heights 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 430 orbits without dust particles. At that stage, the gas density is reset to its initial distribution and dust particles are introduced such that the dust-to-gas 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 10-2, 10-3 and 10-4. All the simulations are further integrated for about 200 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 (R,Z) plane are shown in the case at time t=600. 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 t=550 and t=600 using 50 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

 Figure 2: Radial profile of ( solid line), ( dashed line) and ( dotted line) for the case . The data have been averaged in time between t=550 and t=600. Open with DEXTER

 Figure 3: Radial profile of (time averaged between t=550 and t=600) for the models ( solid line), 0.001 ( dotted line) and 0.0001 ( dashed line). Angular momentum transport is similar at all radii in the three models, despite differences arising because of the stochastic nature of MHD turbulence (see text for details). Open with DEXTER

 Figure 4: Vertical profile of the vertical velocity fluctuations (in units of the speed of sound) for the model having ( solid line), 0.001 ( dotted line) and 0.0001 ( dashed line) at R=2.93. The simulations data have been averaged between t=550 and t=600. Open with DEXTER

 Figure 5: Time variation of the correlation function of the vertical velocity fluctuations in the disc midplane ( solid line) and in the disc corona ( dashed line). The dotted lines represent functions decreasing exponentially toward zero with typical times ,0.1, 0.15 and 0.2 orbit. They can be used to estimate the typical correlation timescale of the turbulence. At both locations, it is not far from the value orbit we used when modeling the results of the numerical simulations. Open with DEXTER

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 density-weighted azimuthal and vertical averages (see Eq. (6) of Fromang & Nelson 2006). To reduce statistical noise, was further averaged in time between t=550 and t=600. 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 t=550, 560, 570, 580, 590 and 600. 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 set-ups (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 in Fig. 3 where we compare the radial profile of for the models (solid line), 0.001 (dotted line) and 0.0001 (dashed line). For each case, the curves are averaged in time between t=550 and t=600. Figure 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 R=2.93, normalized by the speed of sound and averaged in time between t=550 and t=600 for the models (solid line), 0.001 (dotted line) and 0.0001 (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 5% 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. Figure 4 shows that varies by a factor of about 5 between the disc midplane and its corona.

As expressed by Eq. (21), the dust diffusion coefficient also depends on the turbulence correlation timescale . Although D 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 t0=573. We then calculated and averaged the function at seven different radii R=2, 2.5, 3, 3.5, 4, 4.5 and 5 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.
 Figure 6: Dust density distribution in the (R,Z) plane in the numerical simulation for the model having ( left panel), ( middle panel) and ( right panel). For all cases, the raw data have been first azimuthally averaged and then time averaged between t=550 and t=600. Open with DEXTER
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 , 0.1, 0.15 and 0.2 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 0.05-0.1 in the corona and 0.1-0.2 in the disc midplane. Both values are comparable to the value of 0.15 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 D 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.

 Figure 7: Dust density distribution in the (R,Z) plane computed by fitting a Gaussian vertical profile to the simulation data (azimuthally and time averaged) at each radii. As is the case for Fig. 6, the left, middle and right panels respectively correspond to , 0.001 and 0.0001. Open with DEXTER

### 4.3 Dust spatial distribution

 Figure 8: Same as Fig. 7 using the model that assumes D=constant. Open with DEXTER

 Figure 9: Same as Fig. 7 using the model that assumes . Open with DEXTER

For the three simulations we performed, we averaged the dust density in azimuth and over time between t=550 and t=600 orbits. The spatial distributions we obtained using this procedure are shown in Fig. 6 for the models having (left panel), 0.001 (middle panel) and 0.0001 (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 set-up 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.

 Figure 10: Vertical profiles of the dust density, radially averaged between R=3 and R=5. As shown by the legend displayed on the upper panel, the different curves (normalized by their midplane values) were obtained using the numerical simulations ( solid line), by fitting a Gaussian profile to the simulation data ( dotted-dashed line), using a constant turbulent diffusion coefficient ( dashed line) and a vertically varying diffusion coefficient ( dotted line). The upper, middle and lower panels respectively correspond to , 0.001 and 0.0001. In the first case, all models successfully reproduced the simulation, while for the smaller dust sizes, only the later and more elaborate model produces satisfactory results. Open with DEXTER

To compare these results with the models presented in Sects. 3.13.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 et al. 2006b; Fromang & Papaloizou 2006; Johansen & Klahr 2005). 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 D 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 in Fig. 4 for the special case R=2.93) 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 scale-height is about 0.3H. 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 1.5H). 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 0.0001, 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 upper panel gathers the results corresponding to the case , the middle panel shows results obtained when . Finally, the lower 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 upper 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 lower 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 103 at Z=3H) and the model having a constant diffusion coefficient underestimates that density (the ratio at 3H is in excess of 106). 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 10, despite the fact that the dust density itself varies by more than 6 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 Sect. 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 in 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 Z=3H, 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 Z=3H. 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 Z=4H and increases with increasing Z.

## 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

 Figure 11: This figure compares the dust density vertical profile (shown with a solid line), in the case , with a set of Gaussian profiles with different scale-heights (shown with dotted lines). The different values of we used are 0.4H, 0.5H, 0.6H, 0.7H, 0.8H and 0.9H. None of the dotted lines provides an acceptable fit to the numerical dust density. This demonstrates that the dust vertical profile is not Gaussian. Open with DEXTER

 Figure 12: Dust disc scale-heights obtained through a Gaussian fit to the data are represented by the diamonds as a function of . The solid line displays the best power law fit to the points and indicates that . Open with DEXTER

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 , 0.5, 0.6, 0.7, 0.8 and 0.9. 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 scale-height is 0.9H while one would estimate it to be of order 0.6Hat . Both values are different from the dust disc scale-height 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 gas-dust 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 dust-gas ratio decreases as one moves further away from the midplane. These results emphasize the point that estimates of the dust disc scale-height 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 scale-height in our simulations using such a fit. We found , 0.84 and 0.40, respectively, for , 0.001 and 0.01. 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 -0.05, although it is clear that our results show a stronger relationship between the dust disc scale-height and particle size. Both values are, however, largely different from the value -0.5 which is often reported in the literature (Carballido et al. 2006; Dubrulle et al. 1995). 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 -0.5 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 a is decreased. In other words, the vertically averaged dust diffusion coefficient decreases with a. This is not unexpected, since the Schmidt number introduced in Sect. 3.2.1 is known to be an increasing function of (Youdin & Lithwick 2007; Cuzzi et al. 1993; Schräpler & Henning 2004). 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, Sc=1.5 already provides a good fit of the dust density in the case . For the cases and 10-4, we found that the vertically averaged Schmidt numbers that best fit the data are respectively Sc=0.4 and Sc=0.03. Although they indicate a scaling with close to linear, these results are not necessary in disagreement with the result of Youdin & Lithwick (2007), who found a quadratic scaling, because of the vertical average we made when doing such measurements. Note also that these fairly low values of the Schmidt number indicate that turbulent diffusion of the smallest dust particles can be much more efficient than angular momentum transport. This difference originates from their different physical origin: angular momentum is transported radially by the Maxwell and Reynolds stresses while dust diffuses away from the disc midplane because of the vertical velocity fluctuations of the gas. The former quantities decrease as one moves away from the disc midplane (see Fromang & Nelson 2006; Miller & Stone 2000) while the latter quantity increases away from the disc midplane. This means that angular momentum is inefficiently transported in the disc corona while small dust particles are efficiently diffused at the same location. As a consequence, the Schmidt number (i.e. the ratio of both diffusion coefficients) is much smaller than unity for the smallest of our particles.

### 5.3 A toy model

 Figure 13: Comparison between the dust density vertical profile (shown in both panels with the solid line) in the case with a simple toy model in which the turbulent velocity distribution is a step function (see text for details). In the upper panel, the velocity fluctuations are taken to be for |z|>2H, while ( dotted line), 0.05 ( dashed line) and 0.075 ( dotted-dashed line) for |z|<2H. In the lower panel, the velocity fluctuations are for |z|<2H and ( dotted line), 0.15 ( dashed line) and 0.30 ( dotted-dashed line) for |z|>2H. The results show that the dust density vertical profile is fairly insensitive to the midplane velocity fluctuations but more strongly depends on their amplitude in the disc upper layers. Open with DEXTER

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 in 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 , 0.05 and 0.075 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 , 0.05 and 0.1, 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 , 0.15 and 0.30. 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 CPU-intensive 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 (Simon et al. 2008; Fromang & Papaloizou 2007). 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 (Lesur & Longaretti 2007; Fromang et al. 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 a. 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 togenerate 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.

Acknowledgements
The simulations presented in this paper were performed using the QMUL High Performance Computing Facility purchased under the SRIF initiative. It is 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

... Sc=0.03
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).

## All Figures

 Figure 1: Snapshots of the gas density ( left panel) and the dust density ( right panel) at t=600 for the case . Open with DEXTER In the text

 Figure 2: Radial profile of ( solid line), ( dashed line) and ( dotted line) for the case . The data have been averaged in time between t=550 and t=600. Open with DEXTER In the text

 Figure 3: Radial profile of (time averaged between t=550 and t=600) for the models ( solid line), 0.001 ( dotted line) and 0.0001 ( dashed line). Angular momentum transport is similar at all radii in the three models, despite differences arising because of the stochastic nature of MHD turbulence (see text for details). Open with DEXTER In the text

 Figure 4: Vertical profile of the vertical velocity fluctuations (in units of the speed of sound) for the model having ( solid line), 0.001 ( dotted line) and 0.0001 ( dashed line) at R=2.93. The simulations data have been averaged between t=550 and t=600. Open with DEXTER In the text

 Figure 5: Time variation of the correlation function of the vertical velocity fluctuations in the disc midplane ( solid line) and in the disc corona ( dashed line). The dotted lines represent functions decreasing exponentially toward zero with typical times ,0.1, 0.15 and 0.2 orbit. They can be used to estimate the typical correlation timescale of the turbulence. At both locations, it is not far from the value orbit we used when modeling the results of the numerical simulations. Open with DEXTER In the text

 Figure 6: Dust density distribution in the (R,Z) plane in the numerical simulation for the model having ( left panel), ( middle panel) and ( right panel). For all cases, the raw data have been first azimuthally averaged and then time averaged between t=550 and t=600. Open with DEXTER In the text

 Figure 7: Dust density distribution in the (R,Z) plane computed by fitting a Gaussian vertical profile to the simulation data (azimuthally and time averaged) at each radii. As is the case for Fig. 6, the left, middle and right panels respectively correspond to , 0.001 and 0.0001. Open with DEXTER In the text

 Figure 8: Same as Fig. 7 using the model that assumes D=constant. Open with DEXTER In the text

 Figure 9: Same as Fig. 7 using the model that assumes . Open with DEXTER In the text

 Figure 10: Vertical profiles of the dust density, radially averaged between R=3 and R=5. As shown by the legend displayed on the upper panel, the different curves (normalized by their midplane values) were obtained using the numerical simulations ( solid line), by fitting a Gaussian profile to the simulation data ( dotted-dashed line), using a constant turbulent diffusion coefficient ( dashed line) and a vertically varying diffusion coefficient ( dotted line). The upper, middle and lower panels respectively correspond to , 0.001 and 0.0001. In the first case, all models successfully reproduced the simulation, while for the smaller dust sizes, only the later and more elaborate model produces satisfactory results. Open with DEXTER In the text

 Figure 11: This figure compares the dust density vertical profile (shown with a solid line), in the case , with a set of Gaussian profiles with different scale-heights (shown with dotted lines). The different values of we used are 0.4H, 0.5H, 0.6H, 0.7H, 0.8H and 0.9H. None of the dotted lines provides an acceptable fit to the numerical dust density. This demonstrates that the dust vertical profile is not Gaussian. Open with DEXTER In the text

 Figure 12: Dust disc scale-heights obtained through a Gaussian fit to the data are represented by the diamonds as a function of . The solid line displays the best power law fit to the points and indicates that . Open with DEXTER In the text

 Figure 13: Comparison between the dust density vertical profile (shown in both panels with the solid line) in the case with a simple toy model in which the turbulent velocity distribution is a step function (see text for details). In the upper panel, the velocity fluctuations are taken to be for |z|>2H, while ( dotted line), 0.05 ( dashed line) and 0.075 ( dotted-dashed line) for |z|<2H. In the lower panel, the velocity fluctuations are for |z|<2H and ( dotted line), 0.15 ( dashed line) and 0.30 ( dotted-dashed line) for |z|>2H. The results show that the dust density vertical profile is fairly insensitive to the midplane velocity fluctuations but more strongly depends on their amplitude in the disc upper layers. Open with DEXTER In the text