A&A 452, 751-762 (2006)
DOI: 10.1051/0004-6361:20054612

Dust settling in local simulations of turbulent protoplanetary disks

S. Fromang^1,2 - J. Papaloizou²

1 - Astronomy Unit, Queen Mary, University of London, Mile End Road, London E1 4NS, UK
2 - Department of Applied Mathematics and Theoretical Physics, University of Cambridge, Centre for Mathematical Sciences, Wilberforce Road, Cambridge, CB3 0WA, UK

Received 30 November 2005 / Accepted 25 February 2006

Abstract
Aims. In this paper, we study the effect of MHD turbulence on the dynamics of dust particles in protoplanetary disks. We vary the size of the particles and relate the dust evolution to the turbulent velocity fluctuations.
Methods. We performed numerical simulations using two Eulerian MHD codes, both based on finite difference techniques: ZEUS-3D and NIRVANA. These were local shearing box simulations incorporating vertical stratification. Both ideal and non ideal MHD simulations with midplane dead zones were carried out. The codes were extended to incorporate different models for the dust as an additional fluid component. Good agreement between results obtained using the different approaches was obtained.
Results. The simulations show that a thin layer of very small dust particles is diffusively spread over the full vertical extent of the disk. We show that a simple description obtained using the diffusion equation with a diffusion coefficient simply expressed in terms of the velocity correlations accurately matches the results. Dust settling starts to become apparent for particle sizes of the order of 1 to 10 centimeters for which the gas begins to decouple in a standard solar nebula model at 5.2 AU. However, for particles which are 10 centimeters in size, complete settling toward a very thin midplane layer is prevented by turbulent motions within the disk, even in the presence of a midplane dead zone of significant size.
Conclusions. These results indicate that, when present, MHD turbulence affects dust dynamics in protoplanetary disks. We find that the evolution and settling of the dust can be accurately modelled using an advection diffusion equation that incorporates vertical settling. The value of the diffusion coefficient can be calculated from the turbulent velocity field when that is known for a time of several local orbits.

Key words: accretion, accretion disks - magnetohydrodynamics (MHD) - methods: numerical - planets and satellites: formation

1 Introduction

Dust particles are very likely to be the basic building blocks that need to be assembled to make planets.

In the "core-accretion'' model for the formation of giant planets (Mizuno 1980; Pollack et al. 1996), the growth of their solid cores proceeds through the accretion of objects ranging from micron sized dust particles to planetesimals of radius $\sim$10 km to eventually reach about  $15~M_{\oplus}$. Dust is also one of the main observational tracers of the structure of protoplanetary accretion disks (D'Alessio et al. 2001; Adams et al. 1987). A detailed knowledge of its dynamics is therefore needed both in order to make theoretical models and to give the best interpretation of the observations.

The key ingredient that drives the evolution of dust particles is the drag force they feel from the gas (Weidenschilling 1977). Gas must be present if a giant planet is to form subsequently. Drag is responsible for example for the radial migration of solid bodies toward the inner parts of a disk for which the pressure decreases radially outwards. The timescale for this migration is so rapid for centimeter and meter size bodies that it has become a problem for standard planet formation theories (Weidenschilling 1977), a resolution of which might involve the particles accumulating close to pressure maxima possibly at the center of vortices.

Another important aspect of dust dynamics is the tendency to settle towards the midplane of the disk, which increases with particle size. Because of this, centimeter and metre sized bodies accumulate close to the disk midplane. As argued by Goldreich & Ward (1973), a gravitational instability in the dust sub-disk in triggered when the dust density becomes high enough (see also Safronov 1969; and, as remarked in a note in proof by Goldreich & Ward 1973; Gurevich & Lebedinsky 1950). At this point, larger bodies known as planetesimals form.

However, because the pressure force acts only on the gas, a velocity shear between the dust sub-disk and the upper layers of the gas disk develops. This can lead to the so called "shear-induced'' turbulence that may be able to mix the dust layer enough to prevent enough gravitational settling to satisfy the condition for gravitational instability (Gómez & Ostriker 2005; Garaud et al. 2004). However, the onset of this shear induced turbulence depends on the vertical profile of the dust sub-disk which is quite difficult to calculate. Indeed, this is determined as the result of the interplay between various physical processes, such as dust coagulation, Brownian motions and turbulent mixing (Dullemond & Dominik 2005). The purpose of this paper is to investigate the behaviour of the dust disk when just one of these operates, namely MHD turbulence, from first principles.

Turbulence in protoplanetary disks is believed to be the result of the magnetorotational instability (MRI). More than a decade of analytical and numerical work (see a review by Balbus & Hawley 1998) has shown that this is capable of leading to sustained turbulence, and magnetic field with associated stresses that transport angular momentum outwards. More and more of its properties have been analysed thanks to the increased computational power currently available. However, little work has focused on the effects of MHD turbulence on dust dynamics. Carballido et al. (2005) and Johansen & Klahr (2005) analysed the radial diffusion of small particles and calculated the effective Schmidt number, which is the ratio between the anomalous turbulent viscosity and the anomalous diffusion coefficient. They found different results which could be due to their different initial setup. Johansen et al. (2006) and Fromang & Nelson (2005) studied the radial migration of centimeter and meter size bodies in turbulent disk models and found that dust trapping in local pressure maxima could be an important confining effect. Very recently, Turner et al. (2006) investigated the dispersion of a passive scalar in stratified disk models by releasing particles at high altitude in the disk. They interpreted the effect of MHD turbulence in terms of a damped wave equation. However their analysis was restricted to infinitely small particles. They did not consider the vertical settling that begins to occur for larger particles. Thus in their case stratification is only significant for the gas, being felt by the dust particles only through their strong coupling to the gas.

A first attempt to study vertical settling numerically was made through global SPH calculations by Barrière-Fouchet et al. (2005). They confirmed earlier analytical results regarding settling in a quiescent disk (Garaud & Lin 2004). However, the question of the effect of MHD turbulence was left unanswered. In this paper, we tackle this problem by means of local MHD simulations that use Eulerian finite difference codes. Our goal is to relate the dust dynamics to the properties of the turbulence. and also to characterise the amount of settling as a function of the particle size. Because of the complexity of this problem, we neglect dust coagulation and only consider populations of particles with a single specified size.

The plan of the paper is as follows: in Sect. 2, we review the characteristics of dust-gas interaction and describe a simple model to account for the dispersion of a passive scalar in a turbulent medium with random velocity field. In Sect. 3, we present the numerical methods we used together with the basic disk model which exhibits sustained MHD turbulence. We go on to study of the diffusion of very small dust particles in Sect. 4. Their evolution is analysed in the context of the simple diffusion model discussed in Sect. 2. We then show in Sect. 5 that larger dust particles start to decouple from the gas and settle toward the midplane of the disk. In this case the evolution is found to be well described using an advection diffusion equation incorporating settling together with the same diffusion coefficient that applied to the very small particles. In Sect. 6, we investigate the effect a dead zone is likely to have on these processes by setting up non-ideal MHD calculations that result in a significant midplane centered dead zone. However, gas motions are effective at preventing complete settling for 10 cm sized particles in this case too. Finally, we discuss our results in Sect. 7.

   
2 Basic equations

We adopt a two fluid model of the dust and gas circulating in a protoplanetary accretion disk. The basic equations are those of MHD applied to the gas component together with those of hydrodynamics applied to the dust component. The model allows the two fluids to have different flow velocities and consequently exchange momentum through drag forces.

The basic equations for the gas component are the continuity, momentum conservation and induction equations. In a frame rotating with angular velocity ${\vec{\Omega}} = \Omega \vec{\hat{k}}$, with $\vec {\hat{k}}$ being the unit vector in the fixed direction of rotation, here called the vertical direction and $\Omega$ being the magnitude of the angular velocity, these take the form

    

$$\frac{\partial \rho}{\partial t} + \vec{\nabla}\cdot (\rho \vec{v}) = 0 ~,$$ (1)

$$\frac{\partial \vec{v}}{\partial t} + ( \vec{v} \cdot \vec{\nabla... ...\vec{\nabla}{\times} \mbox{\boldmath {$B$ }} ) {\times} \mbox{\boldmath {$B$ }}$$

$$\hspace{3.6cm} -\nabla \Phi - \frac{\rho_{\rm d}}{\rho }{{\vec{F}_{\rm drag}}\over{m_{\rm p}}} ~,$$ (2)

$$\frac{\partial \mbox{\boldmath {$B$ }} }{\partial t} = \vec{\nabl... ...mbox{\boldmath {$B$ }} - \eta \vec{\nabla}{\times} \mbox{\boldmath {$B$ }} ) ~.$$ (3)

Here, P is the gas pressure, $\rho$ is the gas density, $\Phi$ is sum of the gravitational potential, here due to a central mass and the centrifugal potential $-\Omega^2 \vert\vec{r}{\times}\vec{\hat{k}}\vert^2/2$ with $\vec{r}$ being the position vector, ${ \mbox{\boldmath {$v$ }} }$ is the gas velocity and $\mbox{\boldmath {$B$ }}$denotes the magnetic field. The resistivity of the gas, $\eta,$ is taken to be nonzero only in Sect. 6 where we investigate the effect of non-ideal MHD. The last term in Eq. (2) is the acceleration of the gas due to drag on the dust component. The drag force acting on a single dust particle of mass $m_{\rm p}$ is ${\vec{F}_{\rm drag}}$ (see below) and  $\rho_{\rm d}$ is the dust density.

The equations for the dust are

   

$$\frac{\partial \rho_{\rm d}}{\partial t} + \vec{\nabla}\cdot (\rho_{\rm d} \mbox{\boldmath {$v$ }} _{\rm d}) = 0 ~,$$ (4)

$$\frac{\partial \mbox{\boldmath {$v$ }} _{\rm d}}{\partial t} + ( ... ...ath {$v$ }} _{\rm d} = -\nabla \Phi + {{\vec{F}_{\rm drag}}\over{m_{\rm p}}} ~,$$ (5)

with $\mbox{\boldmath {$v$ }} _{\rm d}$ being the velocity of the dust component.

2.1 The drag force

The dust interacts with the gas through drag. In this paper, we consider only dust particles that are small enough in size that they are in the Epstein regime (Weidenschilling 1977). The drag force  ${\vec{F}_{\rm drag}}$ acting on a single particle then takes the form

$${\vec{F}_{\rm drag}}=\frac{m_{\rm p}}{\tau_{\rm s}}( \mbox{\boldmath {$v$ }} - \mbox{\boldmath {$v$ }} _{\rm d}) .$$ (6)

This drag force is proportional to $\mbox{\boldmath {$v$ }} - \mbox{\boldmath {$v$ }} _{\rm d}$ which is the relative velocity between the gas and dust components, and  $\tau_{\rm s}$, the dust stopping time, defined by

$$\tau_{\rm s}=\frac{\rho_{\rm s} a}{\rho c_{\rm s}} ,$$ (7)

where a is the dust particle radius, $\rho_{\rm s}$ is the dust particle density and $c_{\rm s}$ is the local speed of sound which, as we shall take an isothermal equation of state, is independent of height.

When $\Omega$ is taken to be the Keplerian angular velocity at some disk radius, the size of the particles there can be measured in terms of the dimensionless parameter  $\Omega \tau _{\rm s}$ through the relation

$$a=\left( \frac{\rho}{\rho_{\rm s}} \right) ( \Omega \tau_{\rm s} ) H ~ ,$$ (8)

where $H = c_{\rm s}/\Omega$ is the disk scale height. Following Johansen et al. (2004), when we take $\rho_{\rm s}/\rho=10^{10}$ and H=10¹² cm, Eq. (8) gives

$$a=100 \left( \Omega \tau_{\rm s} \right)~{\rm cm} .$$ (9)

The above relation is also appropriate for the standard minimum mass solar nebula at 5.2 AU.

Using an approach based on Boltzmann averaging, Garaud et al. (2004) argued that the evolution of the dust particle distribution can be accurately determined by modelling it as a pressureless fluid as long as the dimensionless parameter  $\Omega \tau _{\rm s}$ is less than unity. This is the case for the dust particles over most of the disk domains considered in this paper. Accordingly we adopt the two fluid description in preference to, for example a much more computationally demanding N-body approach. We also assume throughout this paper that the dust particles remain electrically neutral so that they do not react to the magnetic field.

   
2.2 Single fluid with advection diffusion equation and settling

We here note that in the limit $\tau_{\rm s} \rightarrow 0,$ the two fluid description can be reduced to one of a single fluid with MHD together with an advection diffusion equation describing the evolution of the dust. To see this we introduce a density  $\overline{\rho}$ and a velocity  $\overline{ \mbox{\boldmath {$v$ }} }$ associated with the combined fluid. These are defined through

$$\overline{\rho} = \rho + \rho_{\rm d}$$ (10)

and

$$\overline{\rho} ~ {\overline { \mbox{\boldmath {$v$ }} }} = \... ...th {$v$ }} }+ \rho_{\rm d}{ \mbox{\boldmath {$v$ }} _{\rm d}}.$$ (11)

Adding the continuity Eqs. (1) and (4) and eliminating the drag force by combining Eqs. (2) and (5) and subsequently letting  $\tau_{\rm s} \rightarrow 0$ while assuming $\mbox{\boldmath {$v$ }} - \mbox{\boldmath {$v$ }} _{\rm d} = O(\tau_{\rm s})$gives

   

$${{\partial {\overline \rho}}\over {\partial t}} + \vec{\nabla}\cdot ({\overline \rho} ~ {\overline { \mbox{\boldmath {$v$ }} }}) = 0 ~,$$ (12)

$${\overline \rho}\left(\frac{\partial {\overline { \mbox{\boldmath... ...ath {$B$ }} ) {\times} \mbox{\boldmath {$B$ }} -{\overline \rho} \nabla \Phi ~.$$ (13)

Correct to first order in $\tau_{\rm s},$ we may set ${ \mbox{\boldmath {$v$ }} }_{\rm d} = {\overline { \mbox{\boldmath {$v$ }} }}$in the left hand side of Eq. (5) and then use that equation to find $\mbox{\boldmath {$v$ }} _{\rm d}$located on the right hand side. Substituting this into Eq. (4) and after making use of Eq. (13) gives rise to the equation

$$\frac{\partial \rho_{\rm d}}{\partial t} + \vec{\nabla}\cdot (\rho_{\rm d} { \mbox{\boldmath {$u$ }} }) = 0 ~$$ (14)

with

$${ \mbox{\boldmath {$u$ }} } = {\overline{ \mbox{\boldmath {$v$ }} }} + {\overline{ \mbox{\boldmath {$v$ }} }_{\rm s}},$$ (15)

where

$${\overline{ \mbox{\boldmath {$v$ }} }_{\rm s}}= -{\rho\tau_{... ...{\boldmath {$B$ }} ) {\times} \mbox{\boldmath {$B$ }} \right).$$ (16)

Using the above in Eq. (14) and setting $\rho_{\rm d} = X{\overline \rho},$ $\rho =(1- X){\overline \rho},$ results in the following equation for the dust mass fraction X

 

$${\overline \rho}\left( {{\partial X}\over {\partial t}} + {\overl... ...e \rho}} \right) \right]= \nabla\cdot (\rho c_{\rm s}^2\tau_{\rm s} X\nabla X).$$ (17)

This can be interpreted as an advection diffusion equation where the advection velocity  ${ \mbox{\boldmath {$u$ }} }_{\rm a}$ is the sum of the mean fluid velocity and a settling velocity such that

$${ \mbox{\boldmath {$u$ }} }_{\rm a} = {\overline{ \mbox{\bol... ...} \mbox{\boldmath {$B$ }} }{4\pi{\overline \rho}} \right)\cdot$$ (18)

We remark that for small X and no magnetic field, the vertical component of the settling velocity is $-g\tau_{\rm s}$with g being the vertical acceleration due to gravity. In addition there is a small diffusion coefficient arising through the second derivatives of X with respect to the coordinates such that  ${\cal D} = (1-X) Xc_{\rm s}^2\tau_{\rm s}.$ In practice, when turbulence is present, this is overwhelmed by the effect of the mean fluid velocity, that being particularly the case when X is small.

   
2.3 A simple theory for dust diffusion

The above analysis suggests that the dust mass fraction is advected as a passive scalar. When the gas is turbulent part of the velocity field will be turbulent. The characteristic correlation time associated with the turbulence is expected to be a fraction of an orbital period (see below), whereas the underlying settling takes place on a significantly longer timescale $\sim$ $1/(\Omega^2 \tau_{\rm s}).$Thus we separate the effects of turbulent diffusion and settling and consider the simplest case of the dispersion of a passive scalar in steady state, homogeneous and isotropic turbulence. This means that effects related to the imperfect coupling between dust and gas such as turbulent enhancement of grain-grain collisions are neglected (Morfill 1985; Voelk et al. 1980). It is well known (Batchelor 1949; Taylor 1921) that the mean square displacement of the scalar from its location at some initial time can be expressed in terms of the fluid velocity fluctuations. We adopt a Lagrangian approach and for simplicity consider only the vertical z-direction. Extension to consider the other coordinate directions is straightforward but uninformative in our case.

Let z(z₀,t) be the position of a particle which was located at z=z₀ at t=0. Its displacement Z(z₀,t) is given by

$$Z(z_0,t) = z(z_0,t)-z_0 = \int_0^t U(z_0,t') {\rm d}t',$$ (19)

where U(z₀,t) is the vertical component of the velocity of the particle as defined in a Lagrangian sense as seen along its path. It relates to its Eulerian equivalent v_z(z,t) through U(z₀,t)=v_z(z,t). We are interested in the mean square deviation of particles from their initial positions. Using the fact that U(z₀,t)and Z(z₀,t) are related by a time integral, the derivative of  Z²(z₀,t) with respect to t can be expressed in terms of the velocities as:

$$\frac{\partial Z^2(z_0,t)}{\partial t} 2\int_0^t U(z_0,t')U(z_0,t){\rm d}t'$$ =

$$2\int_0^t v_z(z(z_0,t'),t')v_z(z(z_0,t),t){\rm d}t' .$$ = (20)

To obtain a mean square particle displacement we take an ensemble average over many particle realisations, indicated by angular brackets, to obtain

 

$$\frac{{\rm d} \langle Z^2(z_0,t)\rangle}{{\rm d}t} 2\int_0^t S_{zz}(t,\tau) {\rm d}\tau$$ =

$$2\int_0^t\langle v_z(z(z_0,t-\tau),t-\tau)v_z(z(z_0,t),t)\rangle {\rm d}\tau ,$$ = (21)

where $S_{zz}(t,\tau)$ is the velocity correlation function which depends only on the properties of the turbulence. For turbulence in a statistically steady state, as there is no preferred time, it should depend only on the time difference $t-t' = \tau.$ Thus any value of t may be adopted in $S_{zz}(t,\tau)$ so that

$$S_{zz}(t,\tau)=S_{zz}(t=0,\tau)= S_{zz}(t=\tau,\tau) ,$$ (22)

which finally gives

$$S_{zz}(t,\tau)=\langle v_z(z(z_0,\tau),\tau)v_z(z_0,0)\rangle.$$ (23)

When $\tau \rightarrow 0$, S_zz becomes $\langle v_z^2\rangle$, i.e. it is a measure of the mean square turbulent velocity fluctuations. In the opposite limit, when $\tau \rightarrow \infty,$ S_zz tends to zero as the velocities become uncorrelated.

The calculation of the velocity correlation described above is greatly simplified if the integral in (21) can be performed at a fixed location using Eulerian velocity components, such that  $v_z(z(z_0,\tau),\tau)$ is replaced by $v_z(z_0,\tau)$ with z₀ being fixed. This is possible when the distance moved by a particle during the turbulence correlation time is small compared to the characteristic scale of the turbulence (Biferale et al. 1995) as would always be the case if the turbulent velocities were reduced in magnitude while retaining their spatial and temporal characteristics. This situation seems to be a reasonable approximation in the case of the MRI where the turbulence is subsonic and we note that it would be even better in dead zones where the turbulent velocity fluctuations are reduced in magnitude. In this case S_zz simply scales down as the square of the velocity amplitude. Furthermore for the turbulence we consider, there is no preferred location or particles in the above ensemble averages and all fixed locations can be assumed to get complete information about all particle realisations. Then it is reasonable to replace the ensemble average over all possible particle realisations as a function of time, by one over all possible fixed locations using the velocities there at that time. Therefore, from now on, we adopt

$$S_{zz}(\tau)=\langle v_z(z,\tau)v_z(z,0)\rangle ,$$ (24)

where, the ensemble average is over fixed locations z and over different realisations obtained from a set of models. We also define the quantity  $D_{\rm T}(\tau)$ to be

$$D_{\rm T}(\tau)=\int_{0}^{\tau} S_{zz}(\tau') {\rm d}\tau' .$$ (25)

2.4 Diffusion coefficient

As we shall see, for large $\tau,$ the quantity $D(\tau)$ defines an effective diffusion coefficient. For small $\tau$, S_zz takes on a finite value. Thus we expect $D_{\rm T}(\tau) \propto \tau$. This is known as the ballistic regime (Taylor 1921)

$$D_{\rm T}(\tau) \sim \langle v_z^2\rangle \tau \hspace{0.5cm} \textrm{(for small } \tau \textrm{).}$$ (26)

When $\tau$ becomes large, S_zz goes to zero, so we can postulate that $D(\tau)$ goes to a finite limit, meaning that  $\langle Z^2\rangle$ would scale like $\tau$. This is the diffusive regime. If $\tau_{\rm corr}$ is a typical correlation timescale for the turbulence, then we expect, from simple dimensional analysis

$$D_{\rm T}(\tau) \sim \langle v_z^2\rangle \tau_{\rm corr} \hspace{0.5cm} \textrm{(for large } \tau \textrm{).}$$ (27)

However, at this point we insert a cautionary note. Strict convergence as $\tau \rightarrow \infty$ would require ensemble averaging over an infinite number of realisations. In practical cases these will be finite in number and evaluating  $D_{\rm T}(\tau)$ for increasing $\tau$ corresponds to working out the average of a finite number of random walk displacements for increasing numbers of steps. Although this might appear reasonable for modest $\tau,$ it would not be expected to show strict convergence for very large $\tau$in practical cases. For this reason a time span of 8 orbits equivalent to about 50 correlation times is adopted below.

In this context we comment that the modelling of dust spreading using a diffusion coefficient obtained by the above procedure, on account of the averaging involved, of necessity only describes evolution occurring on time scales significantly longer than the correlation time associated with the turbulence. Behaviour occurring on the correlation time scale or faster cannot be meaningfully considered in this approach.

   
3 Numerical methods

3.1 The algorithms

In this paper, we use two well known finite difference codes: ZEUS-3D (Stone & Norman 1992a,b) and NIRVANA (Ziegler & Yorke 1997). Both are used here to solve the MHD equations in the shearing sheet approximation (Goldreich & Lynden-Bell 1965), including vertical stratification. In this approximation a small region of the disk is considered in the neighbourhood of some point rotating in circular orbit with the local Keplerian angular velocity. Local Cartesian coordinates are used with the x axis along the line connecting the origin to the central mass but pointing away from it, and the y axis directed in the direction of orbital motion. When dust is absent the equations solved are (1)-(3) with ${{\vec{F}_{\rm drag}}}=0$ adapted to the shearing sheet by taking

$$\nabla \Phi = -3\Omega^2 x\vec{\hat{i}} +\Omega^2z\vec{\hat{k}},$$ (28)

with $\vec{\hat{i}}$ being the unit vector in the x direction. The evolution of the magnetic field is calculated using the MOC-CT algorithm (Hawley & Stone 1995), such that the solenoidality constraint on the magnetic field is satisfied at all times.

In the two codes, we employ two different algorithms to describe the evolution of the dust particles that are the main focus of this paper. ZEUS-3D was extended to describe the dust particle evolution by means of a second, pressureless fluid (Fromang & Nelson 2005) using Eqs. (4), (5) together with appropriate drag forces in the low X limit. Because of the short stopping time of the dust particles, the effect of these is computed implicitly.

On the other hand, NIRVANA solved the MHD Eqs. (12), (13) together with (3). Consistent with the low X limit the mean flow velocity  $\overline{ \mbox{\boldmath {$v$ }} }$ was used in (3). The evolution of the dust mass fraction was calculated by solving Eq. (17) but with the use of an advection velocity limiter. This was applied such as to restrict the settling velocity ${{ \mbox{\boldmath {$u$ }} }}_{\rm a} - {{ \mbox{\boldmath {$v$ }} }}$ to be less than $0.1c_{\rm s}$ in magnitude. This is required in the low density upper layers where dust gas coupling becomes weak and where the theoretical description breaks down. However, because it applies only in the upper layers, calculations of settling dust are insensitive to this cut off that is required for numerical reasons.

   
3.2 The model properties

In this section, we describe the underlying disk model we used and highlight some of the properties of the MHD turbulence that are of importance for the dynamics of the dust that we go on to study.

The initial disk setup in the absence of dust is very similar to the model of Stone et al. (1996). The parameters are those of the standard shearing box. The equation of state for the gas is isothermal, $P=\rho c_{\rm s}^2$. Because of the vertical stratification, the initial density distribution is given by:

$$\rho \equiv {\overline \rho} =\rho_0 {\rm e}^{-z^2/2H^2},$$ (29)

where $\rho_0$ is the midplane density and $H = c_{\rm s}/\Omega$ is the disk scale height. The Cartesian box extends over the domains [-H/2,H/2] in x, $[0,2\pi H]$ in y and [-3H,3H] in z. We use the standard periodic boundary conditions in shearing coordinates in x(Hawley et al. 1995; Balbus & Hawley 1998) and periodic boundary conditions in y. One needs to be careful with the vertical boundary conditions. Because of the vertical density stratification, in particular the non vanishing vertical component of the gravitational force on the boundary, using simple periodic boundary conditions in that direction leads to unphysical fluctuations generated in the box.

In ZEUS-3D, we therefore made the vertical gravitational force continuous at the boundary by writing the gravitational potential as

$$V_{\rm grav}=\textrm{max}\left(-\frac{1}{2} \Omega^2 (z^2-H_0^2), 0\right) .$$ (30)

In practice, this means taking vertical gravity into account only in the region $z \in [-H_0,H_0]$. We found that taking H₀=2.4 H produces satisfactory results.

In NIRVANA, we further checked for the existence of numerically generated fluctuations by applying the additional procedure of making the vertical gravitational acceleration continuous throughout. To do this the usual acceleration was applied for $z \in [-2.25H,2.25H].$This was then reduced to zero in $\pm$ [2.25H ,2.4 H] by means of linear interpolation and set to zero for |z|> 2.4H. As indicated below, similar results were obtained with both codes.

In both cases the computational box is initially threaded by a magnetic field with zero net flux, being of the form

$$B_z=B_0 \sin (2\pi x/H) .$$ (31)

B₀ is calculated such that $\beta$, the ratio between the thermal pressure and the magnetic pressure is initially 100. Finally, the numerical resolution is  (N_x,N_y,N_z)=(32,100,192) for ZEUS and  ( 36, 100, 196) for NIRVANA. Models with no dust are initiated by imposing a small random velocity component and subsequently run for $\sim$100 orbits.

We describe here simulation results, focusing on those obtained with ZEUS-3D. The results obtained using NIRVANA were very similar and accordingly only a few need to be illustrated here. As expected, the initial set up is unstable to the MRI. Its early growth is illustrated in Fig. 1. Corresponding results obtained with NIRVANA are plotted in Fig. 2. Both figures show the time history of the volume averaged Maxwell stress (upper solid line) and Reynolds stress (lower solid line), normalised by the average midplane pressure P₀. They are respectively defined as

$$T_{r\phi}^{\rm Max} \int\!\!\!\int\!\!\!\int -{B_x B_y\over 4\pi} {\rm d}\tau ,$$ = (32)

$$T_{r\phi}^{\rm Rey} \int\!\!\!\int\!\!\!\int \rho (v_x-\bar{v_x})(v_y-\bar{v_y}) {\rm d}\tau,$$ = (33)

where $\bar{v_x}$ and $\bar{v_y}$ are respectively the vertically and azimuthally averaged radial and azimuthal velocities and d$\tau$ denotes the element of volume. The dashed line in Fig. 1 is the sum of the Maxwell and Reynolds stresses:

$$\alpha=\frac{T_{r\phi}^{\rm Max}+T_{r\phi}^{\rm Rey}}{P_0}\cdot$$ (34)

After a few orbits, it peaks at a maximum as the flow breaks down into turbulence and attains a saturated state. For the remaining part of the simulation, $\alpha$ varies noisily between 0.006 and 0.03, consistent with previous studies (Stone et al. 1996). We remark that there are quite large fluctuations in the stresses. At the present time, their origin is not clear. The origin of these large variations in our stratified zero net flux models should be addressed in the future using higher resolution calculations. In any case, the results show clear evidence that MHD turbulence is sustained for a very long time. The state of the disk is illustrated in Fig. 3, where we plot a typical snapshot of the density in the (r,z) plane. Turbulent density fluctuations superposed on the overall vertical stratification are obvious from this figure.

In Sect. 2.3, we argued that the diffusion of small dust particles depends on velocity fluctuations in the disk. In Fig. 4, we plot their vertical profile normalised by the sound speed. The dashed curve corresponds to the radial velocity fluctuations, while the azimuthal and vertical fluctuations are respectively given by the solid and dotted curves. The midplane values given by this plot are

$$(\delta v_x^2)^{1/2}=0.11 c_{\rm s} , ~ ~ ~ (\delta v_y^2)^{1/2}=0.08 c_{\rm s} \nonumber , ~ ~ ~ (\delta v_z^2)^{1/2}=0.07 c_{\rm s} \nonumber .$$

These values are consistent with numbers previously reported in the literature (Stone et al. 1996).

	Figure 1: Time history of the total stress ( dashed line), the sum of the Maxwell ( upper solid line) and the Reynolds stresses ( lower solid line) obtained with ZEUS-3D. All stresses are normalised by the midplane pressure. After the linear growth of the MRI, the total stress reaches a nonzero quasi-steady state, showing that the turbulence is sustained.
Open with DEXTER

	Figure 2: As in Fig. 1 but obtained with NIRVANA.
Open with DEXTER

	Figure 3: Gas distribution in the (r,z) plane. Density fluctuations across the disk are clearly visible, showing that the entire disk is turbulent at this stage.
Open with DEXTER

	Figure 4: Velocity dispersion as a function of height, normalised by the sound speed. The three curves corresponds to the radial ( dashed line), azimuthal ( solid line) and vertical ( dotted line) velocity dispersions.
Open with DEXTER

   
4 Vertical diffusion of small particles

In this section, we present the results we obtained for very small particles using ZEUS-3D. For a given disk model, simulations of dust evolution are defined by the value of  $\Omega \tau _{\rm s}$. This is actually a function of gas density. To obtain a fixed parameter defining a given simulation, $\Omega \tau _{\rm s}$ is evaluated using the initial uniform midplane disk gas density. Where this is not stated explicitly it should be taken as read. We focus here on particles for which  $\Omega \tau_{\rm s}=10^{-5}$ in the midplane of the initial disk. Equation (9) indicates that this corresponds to micron size particles. This small value of  $\Omega \tau _{\rm s}$ means that the dust particles are very well coupled to the gas and behave almost like a passive scalar.

Figure 5: The first four snapshots show the dust distribution in the (r,z) plane obtained with ZEUS-3D at times t=3, 5, 10 and 15 orbits after the dust is introduced ( from left to right). The last snapshot on the right represents the gas distribution at t=15 orbits. At that time, the micron-sized dust particles are well mixed in the entire disk. The dust density distribution is also seen to be well correlated with the gas density.

Open with DEXTER

The dust is initially distributed in a thin layer around the disk midplane. Initially, the vertical profile of $\rho_{\rm d}$ is taken to be a Gaussian with a thickness $H_{\rm d}=0.2H$. Under the effects of MHD turbulence, this small layer broadens with time. We found that its precise evolution depends largely on the particular time at which the dust particles are released in the disk. In order to improve the statistics, we computed nine simulations, by restarting the disk model described in Sect. 3.2 at times t=40, 45, 50, 55, 60, 65, 70, 75 and 80 orbits.

Figure 5 illustrates the typical evolution of such a model. From left to right, the first four snapshots on the left show the dust density distribution in the (r,z) plane at time t=3, 5, 10 and 15 orbits (measured after the dust has been introduced). The last snapshot shows the corresponding gas distribution at t=15 orbits. The dust is seen to spread quite rapidly. At 15 orbits, it is almost filling the entire vertical extend of the disk. At later times, this distribution does not change qualitatively. By comparing the dust and the gas distribution at t=15 orbits, one can also see that they are very well correlated. We also found that the initial distribution of the dust has very little effect on this final state: starting with an initial dust distribution such that the gas-to-dust ratio is uniform gives an end product almost indistinguishable from the last two snapshots of Fig. 5.

	Figure 6: Time history of the function $S_{zz}(\tau )$, averaged between the nine models that were run with ZEUS-3D. The solid line show the function as it is extracted from the models while the dotted line represents the function $S_0 \exp \left(-t/\tau_{\rm corr}\right)$, with $\tau _{\rm corr}=0.15$ orbits. This dotted curve is seen to match nicely the numerical result.
Open with DEXTER

In order to make a more quantitative comparison between the results of these simulations and the simple theory presented in Sect. 2.3, we computed the functions S_zz and $D_{\rm T}$according to Eqs. (24) and (25). To calculate the former, we first volume averaged the velocity product in radius and azimuth. The vertical average is taken only within Haround the disk equatorial plane. Indeed, the behaviour of the turbulence away from the midplane is likely to depart very much from being homogeneous and isotropic and may be influenced by the vertical boundary conditions. To reduce the statistical noise, the result is then averaged over the nine different models we computed. The resulting function $S_{zz}(\tau )$ is plotted in Fig. 6 using the solid line. The dotted line represents the time variation of the function

$$S_{zz}(\tau)=(0.07c_{\rm s})^2 {\rm e}^{-t/\tau_{\rm corr}} ,$$ (35)

where $\tau _{\rm corr}=0.15$ orbits. Given the good agreement between the solid and dotted lines in Fig. 6, we conclude that  $\tau_{\rm corr}$ is a measure on the typical timescale over which the fluid velocities become uncorrelated.

Integrating S_zz over time gives $D_{\rm T}(\tau)$. The left panel of Fig. 7 shows the time evolution of $D_{\rm T}/(c_{\rm s} H)$. Figure 8 gives results of similar calculations done from runs with NIRVANA (see below). $D_{\rm T}$ is first observed to rise at early times. This is in good agreement with Eq. (26). In fact, the prediction given by Eq. (26) is represented by the dashed line on the right panel of Fig. 7, which is an enlargement of the left panel at small times. It uses $(\delta v_z^2)^{1/2}=0.07 c_{\rm s}$, as obtained in Sect. 3.2. This early linear rise of the diffusion coefficient was also observed by Carballido et al. (2005). Here, we see that it is naturally understood in terms of the fluid velocity correlations. After this initial rise, $D_{\rm T}(\tau)$ is observed to reach a roughly constant value of $5.5 \times 10^{-3} c_{\rm s}H$. This value nicely compares with the naïve estimate of Eq. (27). Indeed, taking the value of $\delta v_z^2$and  $\tau_{\rm corr}$ derived above, one obtain $D_{\rm T} \sim 4.6 \times 10^{-3} c_{\rm s}H$, a value which is surprisingly close to the result of the numerical simulations. This analytical estimate is represented by the dashed line in Fig. 7.

If the dust distribution undergoes diffusive evolution as indicated by the simple theory, $\rho_{\rm d}$ should satisfy a diffusion equation of the form

$$\frac{\partial \rho_{\rm d}}{\partial t}=D \frac{\partial}{\p... ...{\partial z} \left( \frac{\rho_{\rm d}}{\rho} \right) \right],$$ (36)

where D is a constant diffusion coefficient. We tested this hypothesis by solving Eq. (36). To do so, we took an initial dust distribution identical to that of the numerical model and investigate three different values for the diffusion coefficient D. First, we used $D/(c_{\rm s}H)=D_{\rm th}=5.5 \times 10^{-3}$, i.e. the value given by our simple theoretical model in terms of the velocity fluctuations of the fluid. If this theory is correct, we expect this solution to be close to the numerical solution. In order to check its accuracy, we also computed the solution of Eq. (36) using $D/(c_{\rm s} H)=10^{-3}$ and $D/(c_{\rm s}H)=10^{-2}$.

The comparison between the analytical solution of Eq. (36) and the results of the numerical simulations is shown in Fig. 9. It is illustrated by nine panels. They correspond, from top left to bottom right, to times t=0.48, 0.80, 1.12, 1.44, 1.76, 2.08, 2.40, 2.72 and 3.04 orbits (measured after the dust was introduced into the disk). On each panel, there are four curves. The solid line shows the vertical profile of the dust-to-gas ratio, averaged between the nine models we ran and normalised by the value in the equatorial plane at t=0. The three dashed curves are the solutions of Eq. (36) using the three diffusion coefficients mentioned above. Of course, the smaller the value of D, the less the dust is spread across the disk at a given time.

The agreement between the numerical results and the simple model is fairly good. In all panels, the solid line is seen to have the best agreement with the middle dashed curve, calculated using  $D=D_{\rm T}$. The agreement is especially good at low altitude (typically z<H). This was to be expected: first, this is the region in space where we performed the volume average used to calculate the function S_zz. Second, the theory supposes that the turbulence is homogeneous and isotropic. This results in only possibly reasonable modelling near the midplane. The anisotropy of the turbulence is expected to increase away from the midplane, as the density stratification becomes stronger. From Fig. 9, there are some indications that the estimated diffusion coefficient might increase with height. Indeed, at late times, for z>H, the solid line shows better agreement with the dashed curve computed using $D/(c_{\rm s}H)=10^{-2}$ than at lower altitudes, indicating that dust particles are spread more efficiently than the simple theory suggests. This is also supported by results of runs performed with NIRVANA to test the simple diffusion theory. As with ZEUS-3D we considered runs restarted from the basic disk model after 33, 40, and 60 orbits. We evaluated the function  $D_{\rm T}(\tau)$ by performing ensemble averages both over grid points such that |z| < H, and for grid points such that H <|z| < 2Hand then averaging over the models. Results are shown in Fig. 8. The indication is, in agreement with ZEUS, that for |z| < H, $D_{\rm T}$ approaches $\sim$ $5\times 10^{-3} c_{\rm s} H,$ but for H < |z| < 2H, $D_{\rm T}$ is $\sim$4 times larger.

Figure 7: Diffusion coefficient as a function of time, normalised by $c_{\rm s}H$, obtained from the simulations done with ZEUS-3D. It is computed using Eq. (25) and combining an average over a part of the computational box (|z|<H) and over a set of models (see the description in the text). As shown on the left panel, the diffusion coefficient is observed to saturate to a well defined value $D\sim 5.5\times 10^{-3}c_{\rm s}H$ after an initial rise. The dashed line plots the analytical estimate computed from Eq. (27). The right panel shows an enlargement of the early evolution of $D_{\rm T}$, compared with its expected early behaviour as given by Eq. (26) for small time ( dashed line).

Open with DEXTER

Figure 8: As in Fig. 7 but for runs performed with NIRVANA. The lower curve applies for |z| < H, and the upper curve for H < |z| < 2H. In the former case the diffusion coefficient is observed to attain $D\sim 5\times 10^{-3}c_{\rm s}H$ after an initial rise in good agreement with the ZEUS results. In the latter case the estimated diffusion coefficient is less stable but is about 4 times larger. In that case the significance of the diffusion coefficient is less clear for the reasons given in the text.

Open with DEXTER

Figure 9: Vertical profile of the dust to gas ratio as a function of time. The different panels correspond, from top to bottom and from left to right, to t=0.48, 0.80, 1.12, 1.44, 1.76, 2.08, 2.40, 2.72 and 3.04 orbits. On each panels, the solid line is the result of the simulations performed with ZEUS-3d, averaged over the nine models described in the text. Three dashed curves are plotted. They show the solution of Eq. (36), computed using $D/(c_{\rm s} H)=10^{-3}$, $D/(c_{\rm s}H)=5.5\times 10^{-3}$ and $D/(c_{\rm s}H)=10^{-2}$. The smaller the value of D, the smaller the diffusion.

Open with DEXTER

   
5 Vertical settling of larger particles

In the previous section, we have studied the effect of MHD turbulence on very small dust particles, which behave essentially like passive scalars. In this section, we will focus on larger particles for which the settling processes described in the introduction are fast in the absence of turbulence.

We focus on three sizes, for which the parameter $\Omega \tau _{\rm s}$ is respectively 0.001, 0.01 and 0.1 in the midplane of the initial disk. Using Eq. (9), these values correspond respectively to 1 mm, 1 cm and 10 cm. In the absence of turbulence, the typical settling timescale can be written (Dullemond & Dominik 2004)

$$\frac{\tau_{\rm sett}}{T_{\rm orb}} \sim \frac{1}{2\pi} \frac{1}{\Omega \tau_{\rm s}} ,$$ (37)

where $T_{\rm orb}$ is the orbital time. In our three cases, $\Omega \tau _{\rm s}=0.001$, 0.01 and 0.1 correspond to $\tau_{\rm sett}=160$, 16 and 1.6 orbits respectively. Thus in the absence of turbulence, settling would be important within the duration of the simulations in the last two cases.

At t=40 orbits, we introduced the particles into the turbulent disk model with the same vertical distribution as for the small particles discussed in Sect. 4. Their evolution was then followed until the vertical profile of the dust-to-gas ratio reaches a steady state. Once again, by running models in which the initial dust-to-gas ratio is uniform, we checked that this final distribution does not depend on the initial conditions.

	Figure 10: Dust distribution in the (r,z) plane, obtained with ZEUS-3D, after a quasi-steady state has been reached. From left to right, the different panels correspond to $\Omega \tau _{\rm s}=0.001$, 0.01 and 0.1.
Open with DEXTER

The results obtained with ZEUS-3D are illustrated in Fig. 10. This shows the typical distribution of the dust particles in the (r,z) plane after they have reached a steady state. From left to right, the snapshots corresponds to  $\Omega \tau _{\rm s}=0.001$, 0.01 and 0.1. While the smallest dust particles are still efficiently spread across the disk (compare with the fourth panel of Fig. 5), there are some signs of vertical settling for the larger particles. This is obvious in particular for the largest particles. But note however that turbulence is quite efficient at preventing the complete settling of these particles. As discussed above, a very thin dust subdisk would form in just a few orbits in a quiescent disk. We have confirmed that very similar results for the degree of settling as a function of  $\Omega \tau _{\rm s}$ (evaluated in the midplane) are obtained with NIRVANA which employed an advection diffusion treatment of the dust mass fraction. Steady state distributions corresponding to those shown in Fig. 10 are plotted in Fig. 11. Again these are found to be initial condition independent for distributions initiated in the midplane regions.

	Figure 11: As in Fig. 10 but for results obtained with NIRVANA.
Open with DEXTER

Once again, we can make use of the simple theory developed in Sect. 2.3 to interpret these results. However, in this case the full advection diffusion Eq. (17) incorporating settling must be used with an anomalous diffusion coefficient. Neglecting Lorentz forces, in the low X or $\rho_{\rm d}/ \rho$ limit this gives (Schräpler & Henning 2004; Dubrulle et al. 1995).

$$\frac{\partial \rho_{\rm d}}{\partial t}-\frac{\partial}{\par... ...ial z} \left( \frac{\rho_{\rm d}}{\rho} \right) \right] \cdot$$ (38)

Assuming a steady state, this equation gives a simple first order ordinary differential equation for the vertical dust mass fraction profile. We calculate solutions for the same three values of the diffusion coefficient that were used in Sect. 4 and compare them with the numerical simulations. The results are shown in Fig. 12 for the case where  $\Omega \tau _{\rm s}=0.01$ and in Fig. 13 for the case $\Omega \tau _{\rm s}=0.1$. We do not consider the case $\Omega \tau _{\rm s}=0.001$ because the dust-to-gas ratio is almost uniform at the end of the simulation. In these figures the solid line shows the steady-state dust-to-gas ratio obtained at the end of the simulations while the three dashed curves are the corresponding solutions of Eq. (38) calculated with $D/(c_{\rm s} H)=10^{-3}$, $5.5\times 10^{-3}$ and 10^-2. As in Sect. 4, there is a good agreement between the solid line and the middle dashed curve, for which the magnitude of the diffusion coefficient corresponds to that estimated using the velocity fluctuations of the underlying model. The two other values of the diffusion coefficient can clearly be ruled out. There is also a marked difference between the solid curve and the middle dashed curve in the upper layers of the disk. This was also observed in Fig. 9. As pointed out in Sect. 4, this is due to the increase of the diffusion coefficient at high disk altitudes and also to the fact that our simple theory breaks down because the turbulence ceases to be something that can be profitably modelled as homogeneous and isotropic at these locations.

   
6 The effect of a dead zone

All the results presented so far in this paper suppose that gas and magnetic field are perfectly coupled throughout the entire vertical extent of the disk. However, protoplanetary disks are probably cold and dense enough that this perfect coupling is unlikely, at least in some regions of the disk. This situation has generated the "layered accretion'' paradigm (Gammie 1996) in which the gas is turbulent only in the upper layers of the disk, while a quiescent (or "dead'') zone exists around the midplane. The extent of this dead zone is very uncertain, since it depends on the ionising source and on the chemistry (Ilgner & Nelson 2006; Sano et al. 2000; Fromang et al. 2002), but its existence is likely. In this section, we investigate how the presence of a dead zone would influence the results described in the previous sections.

The problem of layered turbulence in local numerical simulations of stratified disks has already been studied by Fleming & Stone (2003). In this section, we reproduce one of their disk models by allowing the resistivity $\eta$ to be a function of position. We choose the vertical profile of $\eta$such that our model is identical to the "larger dead zone'' model of Fleming & Stone (2003):

$$\eta=\eta_0 \exp \left(-\frac{z^2}{2}\right) \exp \left( \fra... ...\sqrt{\pi}} \int_z^{\infty} {\rm e}^{-z'^2} {\rm d}z'\right) ,$$ (39)

where $\Sigma_0/\Sigma_{\rm CR}=30$ and $\eta_0$ is chosen such that the Reynolds magnetic number $\rm Re_m$, defined by

$${\rm Re_m}=\frac{c_{\rm s}H}{\eta} ,$$ (40)

equals 100 in the midplane of the disk. This disk model was run for 100 orbits. The evolution is the same as that found by Fleming & Stone (2003). In particular, we found that density waves are excited in the dead zone by the turbulent motions in the active layers. The velocity dispersions of the three velocity components are shown in Fig. 14 as a function of z. This figure should be compared with Fig. 4. Taken as a whole, there is a decrease of the velocity fluctuations compared to the fully turbulent model. The presence of the dead zone is apparent through a decrease in the azimuthal and vertical velocity fluctuations. The latter have a root mean square dispersion of roughly  $0.03 c_{\rm s}$, thereby showing a decrease by about a factor of two compared to the fully turbulent model. However, since it does not vanish, a nonzero diffusion coefficient can be expected.

	Figure 12: Steady state vertical profile of the dust-to-gas ratio when $\Omega \tau _{\rm s}=0.01$ ( solid line). It has to be compared with the vertical profile calculated using Eq. (36), shown with the dotted line for three different values of the dimensionless diffusion coefficient  $D/(c_{\rm s}H)$: 10-3, $5.5\times 10^{-3}$ and 10-2.
Open with DEXTER

	Figure 13: Same as Fig. 12, but for the case $\Omega \tau _{\rm s}=0.1$.
Open with DEXTER

Using this new underlying disk model and introducing dust particles, we used ZEUS-3D to recalculate the models described in Sect. 5 for which  $\Omega \tau _{\rm s}=0.01$ and  $\Omega \tau _{\rm s}=0.1$. After a few orbits, the dust distribution reaches a new equilibrium state. It is represented in Fig. 15 (for $\Omega \tau _{\rm s}=0.01$) and Fig. 16 (for $\Omega \tau _{\rm s}=0.1$). In both cases, the solid line shows the vertical profile of the dust-to-gas ratio in the nonideal case with a dead zone, while the dashed line corresponds to the fully turbulent case. As expected, the thickness of the dust sub-disk is smaller in the former case.

It is possible to compare these results to those obtained using the simple model presented in Sect. 2.3. To do so, we first note that close to the midplane, the stopping time  $\tau_{\rm s}$ is nearly constant. A steady-state solution to Eq. (38) can be written in that case as

$$\frac{\rho_{\rm d}}{\rho}=\left( \frac{\rho_{\rm d}}{\rho} \right)_0 {\rm e}^{-z^2/2\tilde{H_{\rm d}}^2} ,$$ (41)

where the dust scale height $\tilde{H_{\rm d}}$ is given by

$$\tilde{H_{\rm d}}=\sqrt{\frac{D}{\Omega^2 \tau_{\rm s}}} \cdot$$ (42)

We first focus on the case $\Omega \tau _{\rm s}=0.1$. When the disk is completely turbulent (or "active''), we found in Sect. 4 that the value of the diffusion coefficient was $D/(c_{\rm s}H)=5.5\times 10^{-3}$. Using Eq. (42), this would give a scale height

$$\tilde{H_{\rm d}}^{\rm active}=0.23$$ (43)

for the dust sub-disk. Next, we seek an estimate of the dust sub-disk scale height in the presence of a dead zone. By combining Eqs. (27) and (42), $\tilde{H_{\rm d}}$ can be related to the velocity fluctuations by

$$\tilde{H_{\rm d}} \propto (\delta v_z^2)^{1/2} .$$ (44)

Given the smaller value we obtained for $(\delta v_z^2)^{1/2}$ in the "larger dead zone'' model (see Fig. 14), we therefore expect

$$\tilde{H_{\rm d}}^{\rm dead} \sim \tilde{H_{\rm d}}^{\rm active}/2 \sim 0.11.$$ (45)

The two dust distributions computed with Eq. (41) using $\tilde{H_{\rm d}}^{\rm active}$ and  $\tilde{H_{\rm d}}^{\rm dead}$ (corresponding respectively to the fully turbulent case and to the "larger dead zone'' model) are plotted in Fig. 16 using dotted lines. Both are seen to match very accurately the solid and dashed curves that are their numerical analogues.

	Figure 14: Same as Fig. 4, but for the "larger dead zone'' model of Fleming & Stone (2003).
Open with DEXTER

	Figure 15: Steady state vertical profile of the dust to gas ratio in the ideal MHD case ( dashed line) and when a dead zone is present around the equatorial plane of the disk ( solid line). The parameter $\Omega \tau _{\rm s}$ equals 0.01 in that case.
Open with DEXTER

The same procedure was followed in the model for which $\Omega \tau _{\rm s}=0.01$. In that case, we found $\tilde{H_{\rm d}}^{\rm active}=0.77$ and $\tilde{H_{\rm d}}^{\rm dead}=0.37$. The dust-to-gas ratio vertical profiles derived using these two values are plotted in Fig. 15 with dotted lines. The agreement with the solid line is quite good. However, there is a poor agreement with the dashed line. This is because the dust is spread to higher altitudes in that case. The hypothesis that  $\tau_{\rm s}$ is a constant which we used to derive Eqs. (41) and (42) starts to break down, which explain the discrepancy with the numerical result.

   
7 Discussion

In this paper, we studied the effects of MHD turbulence on dust settling by means on local numerical simulations performed using ZEUS-3D and NIRVANA, being Eulerian MHD codes using finite differences.

	Figure 16: Same as Fig. 15, but for the case $\Omega \tau _{\rm s}=0.1$.
Open with DEXTER

We first investigated the case of very small particles which are strongly coupled to the gas. Turbulent velocity fluctuations were found to cause an initially thin dust sub-disk to spread. The time evolution of the vertical profile for the dust-to-gas ratio can be well modelled by a diffusion equation, with a diffusion coefficient Dthat can be expressed in terms of turbulent velocity correlations. We found that a simple analytical estimate of D can be obtained in terms of the mean square amplitude of the velocity fluctuations  $\delta v_z^2$ and their correlation time $\tau_{\rm corr},$ both of which are properties of the turbulence alone:

$$D=(\delta v_z^2)^{1/2}\tau_{\rm corr} .$$ (46)

Similarly to Turner et al. (2006), we found an increase of the diffusion coefficient with disk height. We also noticed that the ensemble averages used to calculate it show weaker convergence when the upper layers of the disk are included. While the diffusive description of dust spreading seems to work well in the neighbourhood of the midplane, it is less accurate at disk heights exceeding a few scale heights.

A standard procedure in this type of analysis is to determine the value of the Schmidt number $S_{\rm c}$, defined as the ratio between the anomalous viscosity and the diffusion coefficient. The standard approach in dust diffusion modelling is to take  $S_{\rm c}=1$(Ilgner et al. 2004; Schräpler & Henning 2004; Dullemond & Dominik 2004). In non zero net flux local simulations of radial dust diffusion, Carballido et al. (2005) found $S_{\rm c}=11$, while Johansen & Klahr (2005) reported $S_{\rm c}=1.5$ for vertical diffusion in zero net flux simulations. Turner et al. (2006) also reported a near unity Schmidt number in their calculations. It is worth comparing these values to the Schmidt number we can derive from our simulations. By averaging the total stress represented in Fig. 1 between 20 and 100 orbits, one obtains $\alpha=1.54 \times 10^{-2}$. This, together with the value of the diffusion coefficient obtained from the velocity fluctuations gives

$$S_{\rm c}=\frac{\alpha c_{\rm s} H}{D}=2.8 ,$$ (47)

an intermediate value between the measures of Johansen et al. (2006) and Carballido et al. (2005). However, it is important to stress here that the origin of a non zero diffusion coefficient is on account of the velocity fluctuations and not in the transport properties of angular momentum.

When dust particles grow to centimeter sizes, we found that they start to decouple from the turbulence and settle towards the midplane. The steady state profile of the dust-to-gas ratio is well approximated by the solution of an advection-diffusion equation. Even for particles as large as 10 cm, we found that the dust sub-disk is significantly spread since its semi-thickness $H_{\rm d}$ equals 0.23H, while the settling timescale in a quiescent disk is very short in that case (1.6 orbits). We note however that radial migration is important for particles of this size and the interplay between that migration and MHD turbulence in a stratified disk could lead to complex phenomena, such as local enhancement of the dust density (Fromang & Nelson 2005) that might affect this picture.

Because they are cold and dense, protoplanetary disks are unlikely to have adequate ionisation to be turbulent everywhere. Therefore we also investigated the effect of the presence of a dead zone around the midplane. As expected, we found thinner dust sub-disks in that case. Similarly to previous studies (Fleming & Stone 2003), we found the dead zone is able to maintain significant activity (excited by the turbulent velocity fluctuations of the active zone). This activity is able to prevent the complete settling of 10 cm size particles. However, we want to emphasise that for computational reasons, our analysis was limited to a case in which the mass of the dead zone roughly equals the mass of the active zone. We expect our result to be modified in cases where the mass of the dead zone is much larger than that of the active zone and therefore only apply them to dead zones that involve a modest fraction of the local surface density.

Nonetheless for conditions appropriate to a minimum mass solar nebula, the work presented here that considered MHD turbulence, taken together with that of eg. Gómez & Ostriker (2005) indicates that gravitational instability of the dust layer is unlikely and that the formation of objects of planetesimal size may depend on phenomena such as densification in vortices operating together with vertical settling.

For practical reasons, we neglected grain growth in this work. This is an important simplification, as dust particles are likely to grow at the same time as they settle toward the equatorial plane of the disk (Cuzzi et al. 1996; Dullemond & Dominik 2005). Simulations of the evolution of an entire dust population through grain growth and turbulent stirring, that take account of both radial and vertical disk structure, are very challenging with present day computational capabilities, but will have to be performed in the future.

Acknowledgements

Some of the simulations presented in this paper were performed on the QMUL High Performance Computing Facility purchased under the SRIF initiative.

References

• Adams, F. C., Lada, C. J., & Shu, F. H. 1987, ApJ, 312, 788 [NASA ADS] [CrossRef]
• Balbus, S., & Hawley, J. 1998, Rev. Mod. Phys., 70, 1 [NASA ADS] [CrossRef] (In the text)
• Barrière-Fouchet, L., Gonzalez, J. F., Murray, J. F., Humble, R. J., & Maddison, S. T. 2005, A&A, 443, 185 [EDP Sciences] [NASA ADS] [CrossRef] (In the text)
• Batchelor, G. K. 1949, Austr. J. Phys., 2, 437 [NASA ADS]
• Biferale, L., Crisanti, A., Vergassola, M., & Vulpiani, A. 1995, Phys. Fluids, 7, 2725 [NASA ADS] [CrossRef] (In the text)
• Carballido, A., Stone, J. M., & Pringle, J. E. 2005, MNRAS, 358, 1055 [NASA ADS] [CrossRef] (In the text)
• Cuzzi, J. N., Dobrovolskis, A. R., & Hogan, R. C. 1996, 35
• D'Alessio, P., Calvet, N., & Hartmann, L. 2001, ApJ, 553, 321 [NASA ADS] [CrossRef]
• Dubrulle, B., Morfill, G., & Sterzik, M. 1995, Icarus, 114, 237 [NASA ADS] [CrossRef]
• Dullemond, C. P., & Dominik, C. 2004, A&A, 421, 1075 [EDP Sciences] [NASA ADS] [CrossRef] (In the text)
• Dullemond, C. P., & Dominik, C. 2005, A&A, 434, 971 [EDP Sciences] [NASA ADS] [CrossRef] (In the text)
• Fleming, T., & Stone, J. M. 2003, ApJ, 585, 908 [NASA ADS] [CrossRef] (In the text)
• Fromang, S., & Nelson, R. 2005, MNRAS, 364, L81 [NASA ADS] (In the text)
• Fromang, S., Terquem, C., & Balbus, S. A. 2002, MNRAS, 329, 18 [NASA ADS] [CrossRef]
• Gammie, C. F. 1996, ApJ, 457, 355 [NASA ADS] [CrossRef] (In the text)
• Garaud, P., Barrière-Fouchet, L., & Lin, D. N. C. 2004, ApJ, 603, 292 [NASA ADS] [CrossRef]
• Garaud, P., & Lin, D. N. C. 2004, ApJ, 608, 1050 [NASA ADS] [CrossRef] (In the text)
• Goldreich, P., & Lynden-Bell, D. 1965, MNRAS, 130, 125 [NASA ADS] (In the text)
• Goldreich, P., & Ward, W. R. 1973, ApJ, 183, 1051 [NASA ADS] [CrossRef] (In the text)
• Gómez, G. C., & Ostriker, E. C. 2005, ApJ, 630, 1093 [NASA ADS] [CrossRef]
• Gurevich, L., & Lebedinsky, A. 1950, Seriya fizich, 14, 765 (In the text)
• Hawley, J., & Stone, J. 1995, Comput. Phys. Commun., 89, 127 [NASA ADS] (In the text)
• Hawley, J. F., Gammie, C. F., & Balbus, S. A. 1995, ApJ, 440, 742 [NASA ADS] [CrossRef]
• Ilgner, M., Henning, T., Markwick, A. J., & Millar, T. J. 2004, A&A, 415, 643 [EDP Sciences] [NASA ADS] [CrossRef]
• Ilgner, M., & Nelson, R. 2006, A&A, 445, 205 [EDP Sciences] [NASA ADS] [CrossRef]
• Johansen, A., Andersen, A. C., & Brandenburg, A. 2004, A&A, 417, 361 [EDP Sciences] [NASA ADS] [CrossRef] (In the text)
• Johansen, A., & Klahr, H. 2005, ApJ, 634, 1353 [NASA ADS] [CrossRef] (In the text)
• Johansen, A., Klahr, H., & Henning, T. 2006, ApJ, 636, 1121 [NASA ADS] [CrossRef] (In the text)
• Mizuno, H. 1980, Progr. Theor. Phy., 64, 544 [NASA ADS]
• Morfill, G. E. 1985, in Birth and the Infancy of Stars, 693
• Pollack, J., Hubickyj, O., Bodenheimer, P., et al. 1996, Icar, 124, 62 [NASA ADS]
• Safronov, V. S. 1969, Evoliutsiia doplanetnogo oblaka (In the text)
• Sano, T., Miyama, S. M., Umebayashi, T., & Nakano, T. 2000, ApJ, 543, 486 [NASA ADS] [CrossRef]
• Schräpler, R., & Henning, T. 2004, ApJ, 614, 960 [NASA ADS] [CrossRef]
• Stone, J. M., Hawley, J. F., Gammie, C. F., & Balbus, S. A. 1996, ApJ, 463, 656 [NASA ADS] [CrossRef] (In the text)
• Stone, J. M., & Norman, M. L. 1992a, ApJS, 80, 753 [NASA ADS] [CrossRef]
• Stone, J. M., & Norman, M. L. 1992b, ApJS, 80, 791 [NASA ADS] [CrossRef]
• Taylor, G. I. 1921, Proc. Lond. Math. Soc., 20, 196
• Turner, N. J., Willacy, G., Bryden, G., & Yorke, H. W. 2006, ApJ, 639, 1218 [NASA ADS] [CrossRef] (In the text)
• Voelk, H. J., Morfill, G. E., Roeser, S., & Jones, F. C. 1980, A&A, 85, 316 [NASA ADS]
• Weidenschilling, S. J. 1977, MNRAS, 180, 57 [NASA ADS] (In the text)
• Ziegler, U., & Yorke, H. W. 1997, Comput. Phys. Commun., 101, 54 [NASA ADS] (In the text)