Contents

• Dust distribution in protoplanetary disks
• 1 Introduction
• 2 Dust dynamics
• 3 Description of the code
• 4 Simulations
• 5 Discussion
• 6 Conclusion
• References
• 7 Online Material
• Appendix A: Implicit schemes for the drag force

A&A 443, 185-194 (2005)
DOI: 10.1051/0004-6361:20042249

Dust distribution in protoplanetary disks

Vertical settling and radial migration

L. Barrière-Fouchet¹ - J.-F. Gonzalez¹ - J. R. Murray² - R. J. Humble³ - S. T. Maddison²

1 - Centre de Recherche Astronomique de Lyon (CNRS-UMR 5574), École Normale Supérieure de Lyon, 46 allée d'Italie, 69364 Lyon Cedex 07, France
2 - Centre for Astrophysics and Supercomputing, Swinburne University of Technology, PO Box 218, Hawthorn, VIC 3122, Australia
3 - Canadian Institute for Theoretical Astrophysics, University of Toronto, 60, St. George Street, Toronto, ON M5S 3H8, Canada

Received 25 October 2004 / Accepted 3 August 2005

Abstract
We present the results of a three dimensional, locally isothermal, non-self-gravitating SPH code which models protoplanetary disks with two fluids: gas and dust. We ran simulations of a 1 $M_\odot$ star surrounded by a 0.01 $M_\odot$ disk comprising 99% gas and 1% dust in mass and extending from 0.5 to ${\sim}300$ AU. The grain size ranges from 10^-6 m to 10 m for the low resolution ( ${\sim}25~000$ SPH particles) simulations and from 10^-4 m to 10 cm for the high resolution ( ${\sim}160~000$ SPH particles) simulations. Dust grains are slowed down by the sub-Keplerian gas and lose angular momentum, forcing them to migrate towards the central star and settle to the midplane. The gas drag efficiency varies according to the grain size, with the larger bodies being weakly influenced and following marginally perturbed Keplerian orbits, while smaller grains are strongly coupled to the gas. For intermediate sized grains, the drag force decouples the dust and gas, allowing the dust to preferentially migrate radially and efficiently settle to the midplane. The resulting dust distributions for each grain size will indicate, when grain growth is added, the regions when planets are likely to form.

Key words: stars: planetary systems: protoplanetary disks - hydrodynamics - methods: numerical

1 Introduction

Protoplanetary disks are by definition where we expect planets to form, with micron to millimetre size grains characteristic of the interstellar medium collecting and aggregating to form bodies thousands of kilometres in diameter. Grain growth initially occurs via individual grains sticking via collisions. However, once grains reach millimetre size, their collision velocities are sufficiently large to shatter the grains upon impact (Jones et al. 1996). One way of reducing the relative velocity of colliding grains is to increase their number density. Until recently, the dust dynamics of protoplanetary disks has mostly been studied from the viewpoint of calculating rates of radial migration into the central star. Classical work in this field was done by Weidenschilling (1977), who investigated the nature of the drag force between the dust and gas components of a non-turbulent protostellar disk and then derived equations describing the radial motions of individual dust particles. Weidenschilling concluded that the maximum radial velocity was independent of the drag law and was also insensitive to the nebula mass. Radial migration was fastest for meter sized objects with a velocity of ${\sim}10^4$ cm s^-1.

Stepinski & Valageas (1997,1996) added turbulence and grain coagulation and stressed the idea that the dust motion is decoupled from the gas for a given range of particle size. They found that grains with sizes 0.1 to 10⁴ cm have inward radial velocities larger than that of the gas, while larger particles have smaller velocities. Stepinski & Valageas (1997) modelled a Minimum Mass Solar Nebula (0.24 $M_\odot$) with radial extent 15 AU in which all the solids migrated onto the star, whereas a nebula of 0.002 $M_\odot$ and extending over 250 AU (which is close to our disk parameters) resulted in a distribution of solids close to that of our solar system. They did not investigate vertical settling to the midplane.

Rice et al. (2004) studied the concentration of dust in the spiral arms of marginally stable, self-gravitating protoplanetary disks. They followed the evolution of dust test particles added into their existing SPH code which models gas accretion disks. The test particles feel the gravitational effects of the star and gas disk, as well as gas drag, but do not influence the gas disk. They found that the dust density enhancement was significant only for particle sizes between 10 and 100 cm with their nebula parameters, suggesting that grain growth was possible due to the increased density in the dust layer and the decrease in the relative velocities of the dust grains to each other. Vertical settling was found to be insignificant inside the spiral arms.

The case of vertical settling was investigated by Garaud et al. (2004), who analytically derive fluid equations for the average motion of particles in a non-turbulent nebula. Starting from the momentum equation of a single particle in either Stokes or Epstein drag regime (for respectively large and small particles), a Boltzmann distribution is used to obtain the collective behaviour of dust particles which is consistent as dust particles do not interact with each other and the action of dust on the gas is neglected. Garaud et al. (2004) found that small particles move smoothly to the midplane, while large bodies oscillate about the midplane with decreasing amplitude.

It should be noted that the grain sizes discussed here are larger than those found in the interstellar medium. This is consistent with the spectral signature recently observed in the disk of CQ Tau (Testi et al. 2003) suggesting the presence of large grains (a few centimetres), as well as evidence of dust processing and grain growth found in other disks (e.g. Apai et al. 2004; Meeus et al. 2003).

In this article, we present our new three dimensional gas + dust code and compare the results of disk simulations briefly to the work of Weidenschilling (1977) and more extensively with Garaud et al. (2004) and Garaud & Lin (2004). The code allows us to follow both the radial migration and vertical settling of dust in high resolution, and the possibility of including additional physical processes such as turbulence, grain evolution, and radiative transfer in a simple approximation or detailed equation of state. It has already been applied in a preliminary study which uses scattered light as a diagnosis of dust settling (Barrière-Fouchet et al. 2004).

In a separate project, we incorporate dust physics into a parallel N-body +SPH code which calculates self-gravitational forces using a tree based data structure (Humble et al. 2005; Maddison et al. 2003). The data tree adds considerably to the algorithmic complexity of the code, but means we can consider gravitational stability problems.

In Sect. 2, we present the basic equations for dust hydrodynamics and in Sect. 3 we describe our code. In Sect. 4 we describe our simulations: the physical hypotheses, units and initial state, as well as the results. Finally in Sect. 5 we discuss our results in the light of the afore-mentioned references.

2 Dust dynamics

A single particle in circular orbit of radius r about a central body of mass $M_\star$ will move with the Keplerian velocity $v_{\rm k} = \sqrt{{\rm G}M_\star/r}$, where G is the gravitational constant. Gas moving in a disk around the same central body will typically orbit at a slower velocity due to the radial pressure gradient, which solid bodies moving in the gas disk do not experience and so orbit at a velocity different than that of the gas. The two phases interact via a drag which slows down the dust and forces it to migrate inwards to conserve angular momentum.

According to Whipple (1972,1973) and Weidenschilling (1977), small bodies are strongly coupled to the gas and have about the same velocity field, whereas large bodies follow marginally perturbed Keplerian orbits. Medium size particles experience the strongest perturbation to their movement, with increased accretion to the central object and vertical settling.

Before describing the equation of motion, we must explain the notations we have adopted for densities. We use the subscripts d and g to denote dust and gas respectively. We then define the intrinsic density ( $\rho_{\rm g}$, $\rho_{\rm d}$) of a fluid and the void fraction ( $\theta_{\rm g}$, $\theta_{\rm d}$) as the fraction of the volume occupied by a given phase, such that $\theta_{\rm d}+\theta_{\rm g}=1$. The density ( $\hat{\rho}_{\rm g}$, $\hat{\rho}_{\rm d}$) of a phase in a given volume of fluid is then given by

$$\left\{\begin{array}{l} \hat{\rho}_{\rm g} = \theta_{\rm g}~\... ...ho}_{\rm d} = \theta_{\rm d}~\rho_{\rm d} . \end{array}\right.$$ (1)

Dust particles have a high intrinsic density (in our case $\rho_{\rm d}=1$ g cm^-3), so a given mass of dust occupies a very small volume compared to the same mass of gas. Therefore, gas fills almost all the volume whereas dust occupies only a small portion of it, and:

$$\left\{\begin{array}{l} \hat{\rho}_{\rm {g}} \approx \rho_{\r... ...\ \hat{\rho}_{\rm {d}} \ll \rho_{\rm {d}}. \end{array}\right.$$ (2)

When the only force acting on the dust is drag, the equation of motion reads

$$\frac{{\rm d}\vec{v}_{\rm d}}{{\rm d}t}= -\frac{\vec{v}_{\rm ... ...}}{t_{\rm s}},\ {\rm with}\ t_{\rm s} = \frac{m v}{F_{\rm D}},$$ (3)

where $\vec{v}_{\rm d}$ is the dust velocity, $\vec{v}_{\rm g}$ the gas velocity, t the time, $t_{\rm s}$ the stopping time and $F_{\rm D}$ the drag force. The mass of the dust grain is

$$m=\frac{4 \pi}{3} \rho_{\rm d}~s^3$$ (4)

where we have assumed the grains are spherical with radius s. We use $v=\vert\vec{v}_{\rm d}-\vec{v}_{\rm g}\vert$ to denote the velocity difference between dust and gas phases at a particular point in space. The stopping time is the time it takes for a dust grain that starts with a velocity $\vec{v}_{\rm d}$ to reach the gas velocity $\vec{v}_{\rm g}$.

The functional form of the drag force is determined by the Reynolds number, Re, and by the ratio of the mean free path of gas molecules, $\lambda$, to the radius of the grains, s(see Stepinski & Valageas 1996; Weidenschilling 1977). The Reynolds number is given by ${\rm Re}=2s \rho_{\rm g} v/\eta$ where $\eta=\rho_{\rm g} \lambda C_{\rm s} /2$ is the gas molecular viscosity and $C_{\rm s}$ the sound speed. Thus ${\rm Re} = 4 s v / \lambda C_{\rm s}$.

If $\lambda<4s/9$, we are in the Stokes regime and the drag is due to the wake created by dust particles moving through the gas. The expression for the drag force varies according to the Reynolds number:

$$F_{\rm D}=\left\{\begin{array}{lll} 24~{\rm Re}^{-1}~F & {\rm... ...\\ 0.44~F & {\rm for} & {\rm Re} > 800 \\ \end{array}\right.$$ (5)

with

$$F=\pi s^2 \rho_{\rm g} v^2/2.$$ (6)

If $\lambda>4s/9$ and ${\rm Re}<1$, we are in the Epstein regime and the drag is due to thermal agitation and is given by:

$$F_{\rm D}=\frac{4 \pi}{3} \rho_{\rm g} s^2 v C_{\rm s}.$$ (7)

3 Description of the code

We use the Smoothed Particles Hydrodynamics (SPH) algorithm, a Lagrangian technique described by Monaghan (1992). The SPH equations and approximations have been rigorously established by Bicknell (1991).

Our code was based on that of Murray (1996), originally developed to study tidally unstable accretion disks in cataclysmic variables. We have made several major changes: the configuration is for that of a protoplanetary disk rather than a binary disk and we have modified the algorithm for finding near neighbours so that the code can better handle variable spatial resolution. Following the work of Monaghan & Kocharyan (1995) and Maddison (1998), we have added a second "dust'' phase in a full, self-consistent, fluid approach, contrary to other studies using only test particles to describe the dust phase (e.g. Rice et al. 2004). We take into account the pressure exerted by gas on dust and by dust on gas, as well as the drag force of gas on dust. In the results presented in this paper, we did not calculate the drag force of dust on gas because it is negligible in magnitude and computationally time consuming, but it is implemented into the code and can be turned on if required. We do not review the SPH method as it has been extensively described (see, for example, Monaghan 1992), but stress the changes made to the code written by Murray (1996) in the following subsections.

3.1 Density

As we are following the evolution of two fluids, we therefore have two independent density equations. The gas density is obtained by summation over only the gas neighbours and the dust density is obtained through summation over the dust neighbours. Using the subscripts a and b to refer to gas particles and i and j for dust particles, we find:

$$\left\{\begin{array}{l@{=}l} \hat{\rho}_a & \sum_b m_b W_{ab}\\ \hat{\rho}_i & \sum_j m_j W_{ij}, \end{array}\right.$$ (8)

where $W_{ij}=W(\vert\vec{r}_i-\vec{r}_j\vert)$ and W is the cubic spline kernel commonly used in SPH.

3.2 Smoothing length and link list

The SPH smoothing length is usually given by

$$h=h_0\left(\frac{\rho_0}{\rho}\right)^{1/3},$$ (9)

to ensure a roughly constant number of neighbours (in our case gas plus dust particles). Since our two types of particles have very different masses, to avoid numerical error we therefore chose to use the number density $n = \rho / m$ rather than the mass density $\rho$ in the calculation of h. This approach was extensively studied by Ott & Schnetter (2003) and found to be more accurate.

In order to calculate the number density and subsequently the smoothing length, we need to find all the close neighbours of a given particle. As our code does not include self gravity, rather that using a tree code to find neighbours (see Hernquist & Katz 1989), we use a less time consuming link list. In the original implementation of the code, the cells in the linked list were $2h_{\rm max}$ in size, where $h_{\rm max}$ was the maximum value of the smoothing length over the entire simulation (Murray 1996). With such a cell size, a given particle's neighbours laid either in the same link cell or in one immediately neighbouring it. However, with h varying by as much as a factor 1000, the number of particles in some cells became very large and the searches became inefficient. We implemented an alternative link-list scheme due to Speith (private communication) in which the size of the link cell is chosen to be sufficiently small so that each cell only contains a few particles. This means that, in order to find all the neighbours, the search has to be run over a large number of cells. We nevertheless find this to be efficient compared to the traditionally used link cell algorithm without introducing the programming complexity of a tree structure.

3.3 Drag Force

For the gas, the equation of motion is given by:

$$\frac{{\rm d}\vec{v}_a}{{\rm d}t}= \vec{P}_{ab}+\vec{M}_{aj}+\vec{D}_{aj}+\vec{G}_{a}.$$ (10)

$\vec{P}_{ab}$ is the usual SPH internal pressure term except that each fluid is scaled by the volume it occupies, hence

$$\vec{P}_{ab}=-\sum_b m_b \left(\frac{P_a \theta_a}{\hat{\rho}... ...a_b}{\hat{\rho}_b^2} + \Pi_{ab} \right) \vec{\nabla}_a W_{ab},$$ (11)

where $\Pi_{ab}$ is the SPH artificial viscosity as described in Murray (1996). The volume scaling is ensured by the multiplication by the void fraction $\theta$.

The mixed pressure term,

$$\vec{M}_{aj}=-\sum_j m_j \frac{P_a \theta_j}{\hat{\rho}_a \hat{\rho}_j }\vec{\nabla}_a W_{aj},$$ (12)

is the pressure exerted on one fluid by the other.

The drag force term is given by:

$$\vec{D}_{aj}=\sigma\sum_j\frac{m_j K_{aj}}{\hat{\rho}_j\hat{\rho}... ...rac{\vec{v}_{ja}\cdot\vec{r}_{ja}}{r_{ja}^2+\eta^2}\right) \vec{r}_{ja} W_{ja},$$

$${\rm with}\ K_{aj}=\frac{\rho_a \theta_j c_a}{s},$$ (13)

with $\sigma=1/D$, D being the number of dimensions, and c_a the sound speed associated with particle a. The drag exerted by dust on gas is negligible and time consuming so we have not included it in these computations.

Finally, $\vec{G}_{a}$ is the gravity exerted by the central star only, as we do not take self gravity into account. The derivation of the $\vec{M}_{aj}$ and $\vec{D}_{aj}$ terms can be found in Maddison (1998).

The dust equations are slightly different. The internal pressure term disappears because the dust fluid is made of large bodies (compared to gas molecules) that hardly ever encounter each other. These encounters are better described by an SPH shock treatment than by an internal pressure. Thus the dust equation of motion is given by:

$$\frac{{\rm d}\vec{v}_i}{{\rm d}t}=\vec{M}_{ib}+\vec{D}_{ib}+\vec{G}_{i},$$ (14)

where

$$\vec{M}_{ib}=-\theta_i \sum_b m_b \frac{P_b }{\hat{\rho}_i\hat{\rho}_b } \vec{\nabla}_i W_{ib}$$ (15)

is the mixed pressure term exerted on dust,

$$\vec{D}_{ib} = - \sigma \sum_b \frac{m_b K_{ib}}{\hat{\rho}_b... ...\cdot\vec{r}_{ib}}{r_{ib}^2+\eta^2}\right) \vec{r}_{ib} W_{ib}$$ (16)

is the drag force term exerted on dust, and G_i is the gravity exerted by the central star.

We use an implicit scheme which is iterative and hence a matrix inversion is not necessary to derive the drag force. The equations are given by Monaghan (1997) and can be found in Appendix A. Note that we use the simpler implicit scheme of Monaghan (1997), since our tests of the Tischer algorithm were very slow and did not give much better results.

3.4 Timestep

A timestep is defined for each physical force: the pressure ( $\delta t_{\rm p}$), gravity ( $\delta t_{\rm g}$), and drag ( $\delta t_{\rm d}$) timesteps. The pressure timestep is defined so as to verify the Courant condition and includes the viscosity timestep (see Maddison 1998):

$$\left\{\begin{array}{l} \delta t_{\rm p}=\min_a\displaystyle\... ...isplaystyle\sqrt{\frac{h}{\vert f_a\vert}}, \end{array}\right.$$ (17)

where subscripts a and b refer to two given SPH particles, c_a is the sound speed and f_a the acceleration (i.e. the net force per unit mass) of particle a, and $\zeta$ is the name we give to the SPH artificial viscosity parameter usually called $\alpha$ to avoid any confusion with the Shakura-Sunyaev $\alpha$ (Shakura & Sunyaev 1973). While $\delta t_{\rm p}$ is much larger than $\delta t_{\rm g}$, the gravity subroutine is much less time consuming than the pressure subroutine and so we therefore use operator splitting. This allows us to enter the gravity routine several times while only entering the pressure routine once. The pressure defines the global timestep $\delta t$, such that $\delta t = 0.4 \delta t_{\rm p}$, while $\delta t_{\rm g}$ gives the number of times the gravity routine is run after a single run of the pressure routine, following Murray (1996).

We also use operator splitting for the drag force because $\delta t_{\rm d} \gg \delta t_{\rm g}$ but $\delta t_{\rm d}$ can get much smaller than $\delta t_{\rm p}$ in case of strong drag (which happens close to the central object where the gas density is highest and for the smaller grain sizes). We choose "by hand'' the number of times the drag routine is run after a single run of the pressure routine: $\delta t_{\rm d}$ decreases with the grain size and we find that, down to 1 mm, we just need to run it once and then 5 times for 100 $\mu$m, 25 times for 10 $\mu$m and 100 times for 1 $\mu$m grains.

3.5 Limitations

The code treats turbulence in a very limited manner. The dissipation term in the equation of motion was chosen to reproduce the level of dissipation characteristic of the Shakura-Sunyaev turbulence model. However, the length scale on which dissipation occurs is of the order of the smoothing length and so a full turbulent cascade does not develop.

The resolution of the inner edge of the disk is also limited. In the calculations, the disk scale height, which is defined by $H=C_{\rm s}/\Omega$, where $\Omega=\sqrt{GM_\star/r^3}$ is the angular velocity, has to be larger than the smoothing length h otherwise the disk is not resolved. However, H<h at small radii and we therefore must set $H_{\rm code}=\max(H,h)$ to be able to carry on the calculations. It should be noted therefore that the disk is unresolved for r < 6 AU in our high resolution simulations and r< 20-40 AU in the low resolution simulations. Indeed, our disk extends over more than 2 orders of magnitude in radius making it difficult to resolve both the small and large radii. The solution is not, however, to simply truncate the disk at 6 AU (in the high resolution case) because this would result in boundary layer problems. Instead we know that we can trust the results from 6 AU. Hereafter, what we call the inner parts of the disk refer to the regions immediately outside the unresolved parts and not the actual disk inner edge.

4 Simulations

We model a protoplanetary disk of mass 0.01 $M_\odot$ with 1% of dust by mass around a 1 $M_\odot$ star. We consider dust particles ranging from 1 $\mu$m to 10 m and run a series of simulations that include only one grain size at a time. Our disk is vertically isothermal which means that the temperature is constant in the vertical direction but follows a radial power law ( $T \propto r^{-3/4}$). We implicitly assume that whatever heating process is acting within the disk, the subsequent heat is immediately radiated away.

We do not include the effects of photoevaporation, radiation pressure, Poynting Robertson flux, magnetic fields or grain coagulation, sublimation or condensation.

We scale our model so as to have numbers close to unity. The mass unit is one solar mass, the length unit 100 AU and the time unit is the orbital period of a particle at 100 AU around a 1 $M_\odot$object divided by $2\pi$ (1000 yr/$2\pi$). This choice of units leads to a gravitational constant set to 1.

We choose an initial state in near equilibrium conditions. In the following simulations, we let a gaseous non dusty disk relax from an initial distribution given by $\Sigma=\Sigma_0 r^{-3/2}$ and $H=C_{\rm s}/\Omega$. The initial velocity is Keplerian, given by $v_k=\sqrt{GM_\star/r}$. For the power laws of the temperature and initial surface density, we take the parameters of the Minimum Mass Solar Nebula (Hayashi et al. 1985).

Starting from that initial distribution, we let the disk relax for $\sim$8000 years, i.e. 8 orbits at 100 AU, which allows the pressure and artificial viscosity time to smooth out the velocity field. The gas velocity becomes slightly sub-Keplerian because of the pressure gradient. Once the gas disk has relaxed, we then add the dust particles on top of the gas particles with the same velocity and let the disk evolve for a further 8000 years.

We checked that the Reynolds number is always less than unity since the velocity difference between dust and gas particles stays subsonic. The mean free path is given by

$$\lambda=\frac{1}{n \sigma_{\rm0}} \approx \frac{m_{\rm H_2}} {\pi \rho_{\rm g} r^2_{\rm H_2}}\cdot$$ (18)

The right hand side is obtained using the approximations made by Stepinski & Valageas (1996), that hydrogen is considered as molecular and the cross section is assumed to be identical to the geometrical section. $\lambda$ decreases with increasing density, so if the criterion for the Epstein regime, $\lambda>4s/9$, is met in the densest regions (i.e. the centre of the disk) then it will be met everywhere. Because our disk is not massive ( $M_{\rm disk}=0.01~ {M_\odot}$) and is radially extended (300 AU), the highest gas volumetric density we get is of about $5\times10^{-11}$ g cm^-3, giving $9 \lambda / 4 > 11.5$ m. Therefore, for $10^{-6}~{\rm m} < s < 10$ m, we are in the Epstein regime, and

$$t_{\rm s}=\frac{s \rho_{\rm d}}{C_{\rm s} \rho_{\rm g}}$$ (19)

and the equation of motion for the dust is given by

$$\frac{{\rm d}\vec{v}_{\rm d}}{{\rm d}t}=-\frac{C_{\rm s} \rho_{\rm g}}{s \rho_{\rm d}}(\vec{v}_{\rm d}-\vec{v}_{\rm g}).$$ (20)

We studied a large range of dust grains sizes at low resolution ( ${\sim}25~000$ SPH particles) to see the effect of grain size on the dust settling morphology. The simulations proved to have evolved for long enough to see the most striking features of the dust settling. We then focused on intermediate size grains (10 cm to 100 $\mu$m) where the drag is most efficient and ran high resolution simulations ( ${\sim}160~000$ SPH particles). The results are shown in Figs. 1 and 2 for the low resolution, and Table 1 lists the simulation parameters.

 

 

Table 1: Simulation parameters.

$M_\star$	$1.0\ M_\odot$	s	10, 1, 10-1, 10-2, 10-3, 10-4, 10-5, 10-6 m
$M_{\rm gas}$	$0.01\ M_\star$	$C_{\rm s}$	C0 r-3/8, C0=0.1 code units
$M_{\rm dust}$	$0.01\ M_{\rm gas}$	$\Sigma$	$\Sigma_0~r^{-3/2}$
$R_{\rm disk}$	300 AU	$\alpha_{\rm SPH}$	$\zeta=0.1$
$\rho_{\rm d}$	1.0 g cm-3	$\beta_{\rm SPH}$	0.0

A high resolution single grain size computation requires an average of 15 000 CPU hours shared over 64 processors with OpenMP parallelisation on the SGI Origin 3800 of the French national computing centre CINES based in Montpellier. As the grain size decreases, the drag timestep gets smaller and the computation time increases. The 100 $\mu$m simulation required 22 400 CPU hours.

   
5 Discussion

Our simulations show that the dust in the solar nebula behaves qualitatively in the way of the analysis of Whipple (1972,1973) and Weidenschilling (1977). For example, the larger bodies ($s \ge 1$ m) are weakly coupled to the gas and follow marginally perturbed Keplerian orbits and the dust disk is flared except near the centre where the gas drag is the most efficient because of the high gas density. On the other hand, the smaller bodies ( $s \le 10\ \mu$m) are so strongly coupled to the gas that they follow its motion and again the dust disk is flared.

But for intermediate size particles (100 $\mu$m $\le s \le$ 10 cm), the drag strongly influences the dust dynamics and induces a motion that is very different from that of the gas. We first see a rapid stopping phase, where the dust settles into a region where the drag force dominates. This region is shown in Fig. 3 as the interior of the $\mu =1$ contours. This parameter $\mu$ was characterised by Garaud et al. (2004) as the ratio of the orbital time to the stopping time:

$$\mu=\frac{1}{\Omega}\frac{\rho_{\rm g}}{\rho_{\rm d}}\frac{C_{\rm s}}{s}\cdot$$ (21)

When $\mu \gg 1$, the Epstein drag is very strong. Once the larger grains reach the $\mu\ge1$ region, they experience rapid settling to the midplane as well as strong radial accretion and fall in the central star. The dust layer therefore gets very thin at small radii. For the smaller grains, the settling to the midplane is so efficient that radial migration cannot respond fast enough and there is a particle pileup and the dust layer gets thicker again. Unfortunately, we cannot quantitatively predict the thickness of the dust layer with these simulations, as dispersion in the velocities of individual SPH particles limits our resolution.

Figure 1: Density contours for dust particles for the disk seen face on at the end of the low resolution simulations, each quadrant representing a different size of particles. The vertical bar gives $\log\rho$ in g cm-3.

Figure 2: Density contours for dust particles of different sizes for the disk seen edge on at the end of the low resolution simulations. The outermost contour is the $\rho =1.9\times 10^{-22}$ g cm-3 contour for the gas. The vertical bar gives $\log\rho$ in g cm-3.

Figure 1 also shows the transition between regimes of weak and strong vertical settling, as well as weak and strong radial migration. The 10 m particles do not show any particular increase in density because the settling is weak, and their distribution stays very close to their initial configuration, which is that of the gas (as one can see in the edge-on plots in Fig. 2). Then for 1 m particles, the settling becomes efficient in the inner regions of the disk, but the radial migration is still weak. The particles therefore pile up towards the inner edge of the disk and the density increases. For 10 cm particles, the settling is even stronger and extends to larger radii, but then the radial migration becomes efficient and the particles are lost to the star and hence the density enhancement is weaker than for larger grains. For 1 cm particles, the settling is now so efficient that it extends over even larger radii and despite efficient accretion, too many particles arrive at the inner part of the disk at the same time and thus they pile up and we see a strong density enhancement. For 1 mm particles, the combined effects of settling and migration have proven so efficient that the outer parts of the disk are depleted of dust. Then the efficiency of the vertical settling decreases together with the radial migration, up to the point where particles are completely coupled to the gas as for 1 $\mu$m particles, for which the distribution is undistinguishable from their inital configuration.

In Fig. 4, one can see a large increase in the volume density ratio between dust and gas, rising from 10^-2 to 10^-1 for some radii. This increase is not as striking when we look at the surface density because the effect of the settling is hidden by the summation over the vertical scale height of the disk. Therefore 2D simulations will underestimate the density enhancement. The meaningful quantity to study is in fact the volumetric density, which has the potential to sufficiently damp the dust velocities to allow the grains to stick together via collisions instead of being shattered. It should be noted, however, that while we do not consider the drag of the dust on the gas because it is negligible when the dust to gas ratio is small, as this ratio approach one the dust could drag the gas with it into high density regions, thus resulting in a self regulation of the dust to gas density ratio.

Figure 3: $\mu =1$ contours (see text) at the end of the high-resolution simulations for 100 $\mu$m, 1 mm and 1 cm grains. For the 10 cm particles, this value was not reached and the contour is drawn for $\mu =0.5$ instead.

Figure 4: Surface density $\Sigma$ ( top panel) and volume density $\rho$( bottom panel) as a function of r at the end of the high-resolution simulations. In each case, the dust density is multiplied by 100 to better view the dust density increase.

Humble et al. (2005) show that, for values typical of the solar nebula, the distribution of solid particles evolves to an essentially stationary state and that the final radius of the dust disk can be characterised in terms of the particle size.

As in Weidenschilling (1977), we find that the radial migration rate is highest for a given grain size and drops quickly for larger and smaller particle sizes. We find this maximum is reached for 1 cm size particles, whereas Weidenschilling (1977) found it to be 1 m. This discrepancy comes from the difference in the parameters we use for the nebula and Weidenschilling (1977) showed that the gas density and dust intrinsic density have an impact on this maximum grain size.

To check the vertical evolution, we first compare our results with the self-similar behaviour of the dust phase during settling as described by Garaud & Lin (2004) and shown in their Figs. 1b and 2. Despite the noise due to the SPH technique, our results indeed show the same kind of behaviour, shown in Fig. 5. Note that our density increase is orders of magnitudes smaller than theirs because the computation is too time consuming to be run for as long an evolutionary time as given in Garaud & Lin (2004).

Figure 5: Volume density versus z at r=150 AU for the 100 $\mu$m dust at different times from dust injection (lightest curve) until the end of the high-resolution simulations (darkest curve).

We next validate the formula derived by Garaud et al. (2004) in the Epstein regime and reproduce their Fig. 5 with our conditions (see Fig. 6). Starting from the trajectory of a single particle and then using a Boltzmann averaging, Garaud et al. (2004) derived an analytical formula for the vertical velocity

$$v_z=-\frac{z}{\mu},$$ (22)

with $\mu$ defined in Eq. (21). To calculate the value of $\mu$and then of v_z, we use the gas density that comes out of the simulations and the sound speed prescription that we use in our code. We then compare this computed value of v_z with the one taken directly from the simulations in Fig. 6.

Figure 6: Contours of the vertical velocity for 1 mm ( top panels) and 100 $\mu$m ( bottom panels) particles for the high-resolution simulations 400 yr after dust injection. Left: data from our simulations, right: expected theoretical values from Eq. (22). The vertical bar gives vz in code units (see text).

The proportionality between vertical velocity and vertical position is most striking at the very beginning of the simulation when the vertical extension is maximal and the particle velocities towards the midplane are largest. The plots are shown 400 yr (0.4 orbits at 100 AU) after the addition of dust.

We can easily see the vertical stratification of the velocity for the 1 mm grains. It is less obvious for the 100 $\mu$m grains because the dust is more strongly coupled to the gas and settles more slowly towards the midplane. Furthermore, the noise from the simulation makes the structure more difficult to detect. At the beginning of the simulations for the 10 cm and 1 cm particles, the grains oscillate about the midplane following perturbed Keplerian orbits before the gas damps these oscillations efficiently as shown in Fig. 2 of Garaud et al. (2004).

Finally, looking at the sequence of edge-on images in Fig. 2, we see qualitatively the behaviour discussed by Garaud et al. (2004) that at a given height the region will be depleted of intermediate size particles.

6 Conclusion

6.1 Summary

In this article, we have shown the importance of the vertical settling of intermediate size grains (10 cm to 100 $\mu$m) in a protoplanetary disk with a central object of $M_\star$ = $1.0~M_\odot$, a gas disk mass of $M_{\rm gas}=0.01~M_\star$ and a dust disk mass of $M_{\rm dust}=0.01~M_{\rm gas}$.

We find that the dust distribution is completely different from the gas and the general approximation that the dust to gas ratio is constant throughout the disk is wrong for grains in this size range. This has an impact on the interpretation of observations made in the millimetre domain, for instance, as the opacity of millimetre size grains dominates. It is thus important to understand the specific behaviour of the dust to avoid misinterpretations such as an underestimate or overestimate of the gas density.

Our simulation results compare well with the radial migration seen in Weidenschilling (1977), the analytical work of Garaud et al. (2004) and the vertical settling of Garaud & Lin (2004).

6.2 Perspectives

The Shakura & Sunyaev (1973) turbulence modelling has been extensively used in the domain of protoplanetary disks and has proven to give consistent results. It should, however, be used cautiously and ideally a better prescription for turbulence should be used. Indeed, turbulence will likely reduce the dust settling velocities and lessen the density enhancements we observe in our simulations (see, e.g., Weidenschilling & Cuzzi 1993). The magneto-rotational instability has been consistently described by Balbus & Hawley (1991) and global numerical simulations for this instability are becoming more accurate (Fromang et al. 2004). It is now possible to have a consistent description of this kind of turbulence in a numerical simulation and a first step towards this goal was given by Monaghan (2002), who derived an SPH prescription of turbulence. We plan to incorpate SPH turbulence into our code in order to investigate how turbulence affects the dust settling.

While the locally isothermal approximation is roughly consistent with the temperature distribution in a protoplanetary disk, a lot of effort is currently going into the domain of radiative transfer in protoplanetary disks (e.g. Dullemond & Dominik 2004), suggesting that the energetics in such an object deserve a better treatment. Unfortunately, including full radiative transfer in our SPH code is not feasible because of the increased computation it would require. Nevertheless, we will be able to achieve a better description by implementing an adiabatic equation of state with cooling functions as one of the improvements we plan to add in our code.

This consistent description of the dust distribution in a protoplanetary disk is a first step towards planet formation. By implementing a consistent treatment of the grain growth, coagulation, and shattering, we will be able to better understand planetesimal formation and their distribution in the disk.

Acknowledgements

This work was funded by the French PNPS and L.B.-F. would like to thank the French CINES for the important amount of computing time that made these simulations possible. The authors are grateful to the French-Australian PICS for financial support for travel to Australia. L.B.-F. is also deeply grateful to Eduardo Delgado for comments and suggestions as well as to Roland Speith for the idea of the new link list routine.

References

 
• Apai, D., Pascucci, I., Sterzik, M. F., et al. 2004, A&A, 426, L53
   
• Balbus, S. A., & Hawley, J. F. 1991, ApJ, 376, 214
   
• Barrière-Fouchet, L., Pinte, C., Ménard, F., Gonzalez, J.-F., & Maddison, S. T. 2004, in SF2A-2004: Semaine de l'Astrophysique Française, ed. F. Combes, D. Barret, T. Contini, F. Meynadier, & L. Pagani, Paris, France, June 14-18, 2004 (Les Ulis: EDP Sciences), Conf. Ser., 63
   
• Bicknell, G. V. 1991, SIAM J. Sci. Atat. Comput., 12, 1198
   
• Dullemond, C. P., & Dominik, C. 2004, A&A, 421, 1075
   
• Fromang, S., Balbus, S. A., Terquem, C., & De Villiers, J.-P. 2004, ApJ, 616, 364
   
• Garaud, P., & Lin, D. N. C. 2004, ApJ, 608, 1050
   
• Garaud, P., Barrière-Fouchet, L., & Lin, D. N. C. 2004, ApJ, 603, 292
   
• Hayashi, C., Nakazawa, K., & Nakagawa, Y. 1985, in Protostars and Planets II, 1100
   
• Hernquist, L., & Katz, N. 1989, ApJS, 70, 419
   
• Humble, R. J., Maddison, S. T., Murray, J. R., Barrière-Fouchet, L., & Gonzalez, J.-F. 2005, MNRAS, to be submitted
   
• Jones, A. P., Tielens, A. G. G. M., & Hollenbach, D. J. 1996, ApJ, 469, 740
   
• Maddison, S. T. 1998, Ph.D. Thesis, Monash University
   
• Maddison, S. T., Humble, R. J., & Murray, J. R. 2003, in Scientific Frontiers in Research on Extrasolar Planets, ed. D. Deming, & S. Seager (ASP Conf. Ser.), 307
   
• Meeus, G., Sterzik, M., Bouwman, J., & Natta, A. 2003, A&A, 409, L25
   
• Monaghan, J. J. 1992, ARA&A, 30, 543
   
• Monaghan, J. J. 1997, J. Comp. Phys., 138, 801
   
• Monaghan, J. J. 2002, MNRAS, 335, 843
   
• Monaghan, J. J., & Kocharyan, A. 1995, Computer Physics Communications, 87, 225
   
• Murray, J. R. 1996, MNRAS, 279, 402
   
• Ott, F., & Schnetter, E. 2003 http://arxiv.org/abs/physics/0303112
   
• Rice, W. K. M., Lodato, G., Pringle, J. E., Armitage, P. J., & Bonnell, I. A. 2004, MNRAS, 355, 543
   
• Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337
   
• Stepinski, T. F., & Valageas, P. 1996, A&A, 309, 301
   
• Stepinski, T. F., & Valageas, P. 1997, A&A, 319, 1007
   
• Testi, L., Natta, A., Shepherd, D. S., & Wilner, D. J. 2003, A&A, 403, 323
   
• Weidenschilling, S. J. 1977, MNRAS, 180, 57
   
• Weidenschilling, S. J., & Cuzzi, J. N. 1993, in Protostars and Planets III, 1031
   
• Whipple, F. L. 1972, From plasma to planet, ed. A. Elvius (London: Wiley), 211
   
• Whipple, F. L. 1973, NASA SP, 319, 355
  

 

  
7 Online Material

   
Appendix A: Implicit schemes for the drag force

All the following equations are derived from Monaghan (1997). The momentum equations for gas and dust in Epstein regime read

$$\frac{{\rm d}\vec{v}_{\rm g}}{{\rm d}t}= -\frac{K}{\hat{\rho}_{\rm g}} (\vec{v}_{\rm g}-\vec{v}_{\rm d})$$ (A.1)

$$\frac{{\rm d}\vec{v}_{\rm d}}{{\rm d}t}= -\frac{K}{\hat{\rho}_{\rm d}} (\vec{v}_{\rm d}-\vec{v}_{\rm g}),$$ (A.2)

with

$$K=\frac{C_{\rm s} \hat{\rho}_{\rm g} \hat{\rho}_{\rm d}} {\rho_{\rm d} s},$$ (A.3)

and become in SPH notations (indices a and b refer to gas particles and i and j to dust particles and for a vector $\vec{q}$, $\vec{q}_{aj} = \vec{q}_a - \vec{q}_j$):

$$\frac{{\rm d}\vec{v}_a}{{\rm d}t}=\sigma \sum_j m_j \frac{K_{... ...t \vec{r}_{aj}}{r_{ja}^2 + \eta^2} \right) \vec{r}_{aj} W_{aj}$$ (A.4)

$$\frac{{\rm d}\vec{v}_j}{{\rm d}t}=\sigma \sum_a m_a \frac{K_{... ... \vec{r}_{aj}}{r_{aj}^2 + \eta^2} \right) \vec{r}_{ja} W_{ja}.$$ (A.5)

Then, it is possible to turn it to an implicit scheme without having to invert a matrix. We have

$$\frac{{\rm d}\vec{v}_a}{{\rm d}t}= \sum_j m_j s_{aj} (\vec{v}_{aj} \cdot \vec{r}_{aj}) \vec{r}_{aj}$$ (A.6)

$$\frac{{\rm d}\vec{v}_j}{{\rm d}t}= \sum_a m_a s_{aj} (\vec{v}_{aj} \cdot \vec{r}_{aj}) \vec{r}_{ja},$$ (A.7)

with

$$s_{aj}=\sigma \frac{K_{aj}}{\hat{\rho}_a \hat{\rho}_j} \left(\frac{W_{aj}}{r_{aj}^2 + \eta^2} \right)\cdot$$ (A.8)

With a timestep $\delta t$ and $\vec{v}_i^0$ and $\vec{v}_i^1$ respectively the old and new values of the velocity:

$$\vec{v}_a^1 = \vec{v}_a^0 - \delta t \sum_j m_j s_{aj} (\vec{v}_{aj}^1 \cdot \vec{r}_{aj})\vec{r}_{aj}$$ (A.9)

$$\vec{v}_j^1 = \vec{v}_j^0 - \delta t \sum_a m_a s_{aj} (\vec{v}_{aj}^1 \cdot \vec{r}_{aj})\vec{r}_{ja}.$$ (A.10)

Now, we will just consider one dust-gas particle pair at a time:

$$\vec{v}_a^1 = \vec{v}_a^0 - \delta t~m_j s_{aj} (\vec{v}_{aj}^1 \cdot \vec{r}_{aj})\vec{r}_{aj}$$ (A.11)

$$\vec{v}_j^1 = \vec{v}_j^0 - \delta t~m_a s_{aj} (\vec{v}_{aj}^1 \cdot \vec{r}_{aj})\vec{r}_{ja}.$$ (A.12)

Then taking the difference of the previous two equations and taking the scalar product with $\vec{r}_{aj}$, we get

$$\vec{v}_{aj}^1 \cdot \vec{r}_{aj} = \frac{\vec{v}_{aj}^0 \cdot \vec{r}_{aj}} {1 + \delta t (m_a + m_j) s_{aj} r_{ja}^2},$$ (A.13)

and

$$\vec{v}_a^1 = \vec{v}_a^0 - \frac{\delta t m_j s_{aj}\vec{r}_... ...^0 \cdot \vec{r}_{aj})}{1+\delta t (m_a + m_j) s_{aj}r_{ja}^2}$$ (A.14)

$$\vec{v}_j^1 = \vec{v}_j^0 - \frac{\delta t m_a s_{aj}\vec{r}_... ...dot \vec{r}_{aj})}{1+\delta t (m_a + m_j) s_{aj}r_{ja}^2}\cdot$$ (A.15)

The drag of dust on gas is very small, and to save computing time, we neglect it. So, at the end, only dust experiences drag with the following expression:

$$\vec{v}_j^1 = \vec{v}_j^0 - \sum_j \frac{\delta t m_a s_{aj}\... ...dot \vec{r}_{aj})}{1+\delta t (m_a + m_j) s_{aj}r_{ja}^2}\cdot$$ (A.16)

A.1 Tischer implicit scheme

Now, we subdivide each timestep into two, and we note $\vec{\tilde{v}}$the velocity at $\Delta = \delta t / 2$:

$$\vec{\tilde{v}}_a = \vec{v}_a^0 - \Delta\vec{r}_{aj} m_j s_{a... ...v}}_{aj}\cdot\vec{r}_{aj}+0.4\vec{v}_{aj}^0\cdot\vec{r}_{aj}),$$ (A.17)

$$\vec{\tilde{v}}_j = \vec{v}_j^0 - \Delta\vec{r}_{aj} m_j s_{a... ...v}}_{aj}\cdot\vec{r}_{aj}+0.4\vec{v}_{aj}^0\cdot\vec{r}_{aj}),$$ (A.18)

and

$$\vec{v}_a^1 = 1.4\vec{\tilde{v}}_a - 0.4\vec{v}_a^0 - 0.6\Delta\vec{r}_{aj} m_j s_{aj}(\vec{v}_{aj}^1 \cdot \vec{r}_{aj}),$$ (A.19)

$$\vec{v}_j^1 = 1.4\vec{\tilde{v}}_j - 0.4\vec{v}_j^0 - 0.6\Delta\vec{r}_{aj} m_a s_{aj}(\vec{v}_{aj}^1 \cdot \vec{r}_{aj}).$$ (A.20)

For simpler notations, we define $A=\Delta r^2_{aj} (m_j+m_a) s_{aj}$ and $B=\Delta r_{aj} m_a s_{aj}$. Then

$$\vec{\tilde{v}}_{aj}\cdot\vec{r}_{aj}= \frac{1-0.4A}{1+0.6A}\vec{v}^0_{aj}\cdot\vec{r}_{aj},$$ (A.21)

$$\vec{v}^1_{aj}\cdot\vec{r}_{aj}= \frac{1-0.8A}{(1+0.6A)^2}\vec{v}^0_{aj}\cdot\vec{r}_{aj},$$ (A.22)

and

$$\vec{v}^1_{j}=\vec{v}^0_{j}+B\vec{r}_{aj}\frac{2+0.36A}{(1+0.6A)^2} \vec{v}^0_{aj}\cdot\vec{r}_{aj}.$$ (A.23)