Particle-Trapping Eddies in Protoplanetary Accretion Disks

ICARUS

128, 213–229 (1997) IS975720

ARTICLE NO.

Particle-Trapping Eddies in Protoplanetary Accretion Disks H. Hubertus Klahr and Thomas Henning Research Unit ‘‘Dust in Star-Forming Regions,’’ Max Planck Society, Schillerga¨ßchen 3, D-07745 Jena, Germany E-mail: [email protected] Received August 16, 1996; revised January 9, 1997

In this paper, we present a new mechanism for number density enhancement and size segregation of particles in a nonlaminar accretion disk, e.g., the solar nebula. It can be shown that, if there is a vortical flow with a rotational frequency smaller than the local orbital frequency of the disk and a rotational axis parallel to the midplane of the nebula, particles will be trapped and concentrated in the gas flow. Due to the centrifugal force, a particle can be driven out of an eddy. It will be shown that this process is inhibited by the gravitational force induced by the protostar. Candidates for such slowly circulating flows are turbulent flow features in a size range close to the integral size of the turbulence as well as huge convection cells. The efficiency of the particle concentration depends on the coupling of the particle to the gas, i.e., the friction time. On account of the mass dependence of the friction time, a given eddy becomes a trap for particles of a characteristic size and causes a local change in the dust density. We calculate the maximal dust inhomogeneity due to this process. The strongest effect was observed for millimeter-sized particles, which can be concentrated by a factor of 100 within only 100 years. Our general estimates do not depend on special turbulence or convection models. In addition, we compare the analytical estimates with the results of numerical simulations of the dust motion. The effect will have an impact on the dust coagulation process and may be a link to the size distribution of chondrules.  1997 Academic Press

The first studies of this area have been performed by Adachi et al. (1976) and Weidenschilling (1977). Both assumed a laminar rotating accretion disk and discovered the general behavior of the particles to settle toward the midplane and to drift radially to the central object. The protoplanetary accretion disk is, however, assumed to be turbulent. Small particles will not be able to concentrate in the midplane as they are stirred up by the turbulentdriven diffusion (Vo¨lk et al. 1980, Cuzzi et al. 1993, Weidenschilling and Cuzzi 1993, Dubrulle et al. 1995). In these models a locally homogeneous distribution of particles was still assumed, but Wood (1963) had already suggested that there should be ‘‘dust gathered between eddies.’’ Squires and Eaton (1991) showed for the first time turbulent concentration of particles related to their aerodynamic stopping time (reviewed by Eaton and Fessler 1994). Dobrovolskis et al. (1993) and Cuzzi et al. (1996) calculated this particle concentration for the scenario of a protoplanetary accretion disk. They suggest that particles are concentrated in convergence zones between eddies in Kolmogorov-type homogeneous turbulence. In this work, only small and fast rotating eddies have been considered, and the influence of gravitation was neglected. In our work, we point out that for large and slowly rotating eddies, with a rotational axis parallel to the midplane of the nebula, e.g., convection cells, there is no tendency to collect grains between the eddies but, on the contrary, in their interior. One more model has appeared recently which produces particle concentration via vorticity, but operates by a completely different mechanism. Tanga et al. (1996) postulate a system of regular vortices driven by differential rotation of the disk, with rotational axes parallel to its angular momentum. The concentration mechanism then depends on Coriolis forces. We want to stress that by this mechanism particles of much larger size than in our work would be concentrated. The outline of this paper is as follows. In Section 2, we give a brief introduction to the topic of turbulence and particle motion in protoplanetary accretion disks. Section 3 presents the results of our analytical estimations, the derivation of which is in Appendix A. The comparison of

1. INTRODUCTION

In the framework of planetary formation, it is generally accepted that planetesimals form in a protoplanetary accretion disk by accumulation of presolar particles or nebula condensates (cf. review by Lissauer 1993). The dynamical behavior of micrometer-sized particles is dominated by the gas drag, while meter-sized bodies decouple from the gas and move on Keplerian orbits. However, there is a transition regime of particle sizes where the influence of the gas drag, the gravity of the central object, and the inertia forces of the particles are on the same order of magnitude. 213

0019-1035/97 $25.00 Copyright  1997 by Academic Press All rights of reproduction in any form reserved.

214

KLAHR AND HENNING

the analytical results to numerical simulations (see Appendix B), as well as an example for the expected particle concentration in convective zones, are given in Section 4. Section 5 summarizes and discusses our ideas. 2. FUNDAMENTALS

2.1. Nonlaminar Flow Pattern For a number of reasons, one expects a nonlaminar flow in protoplanetary accretion disks. Thermal convection has been suggested by several authors (Lin and Papaloizou 1980, Cabot and Pollack 1992) as a source of turbulent viscosity. However, in recent analytical and numerical investigations (Stone and Balbus 1996), it has been shown that this will not immediately lead to an effective angular momentum transport which is directed outward. Nevertheless, the occurrence of thermal convection is still possible and could be a source of disturbing a laminar flow pattern (Rozyczka et al. 1994, Kley et al. 1993). Another source perturbing the laminar flow is a sheardriven turbulence (Cuzzi et al. 1993, Tennekes and Lumley 1972). This possibility is only given in the late stages of grain evolution because a sufficiently high concentration of particles in the midplane of the accretion disk is needed. This concentration leads to a Keplerian rotational frequency of the midplane, while the upper and lower layer are still pressure supported and rotating with sub-Keplerian velocities. The most promising approach to turbulent viscosity in protoplanetary accretion disks is magnetic turbulence (Balbus and Hawley 1991, Stone et al. 1996). Stone et al. showed that the Balbus–Hawley instability will drive a magnetic turbulence fitting—at least on large scales—a Kolmogorov spectrum (Kolmogorov 1941). The average total stress can be scaled with the equatorial pressure by trf 5 aP(0). This equals a viscosity expressed in the Shakura–Sunyaev model (Shakura and Sunyaev 1973) with a # 1022. Our study of the grain motion is independent of any special model of turbulence in the disk. The only assumptions for the discussion on the order of magnitude effects are the following: The turbulence can be described by a velocity U and a length L at the largest scales, and we can derive the turbulent velocities on smaller scales by a power law as in the Kolmogorov case u(l ) 5 U

SD l L

1/3

.

(1)

The fundamental assumption is the existence of eddies, either as a macroscopic feature like a convective roll or as a special feature of a more general turbulent flow. The kinematical and dynamical properties of such flows can

be better understood by describing the flows in terms of individual events or streamline patterns (Wray and Hunt 1989, Hunt et al. 1988). These patterns can be classified in convergence zones, shear zones, streams, or eddies. Both eddies and shear zones have a high amount of vorticity, but the shear zones do not have a circulating flow pattern. In numerical simulations of incompressible isotropic turbulence about 20% of the volume was occupied by vortical flows, while about 30% was in a nonvortical state. The remaining volume could not be classified with the adopted criteria. Squires and Eaton (1991) and Eaton and Fessler (1994) investigated the effect of turbulence on the concentration of heavy particles in order to study their mixing behavior. Therefore, they studied the influence of different turbulent elements onto particles applying a flow pattern classification, as mentioned above, and ascertained a maximum particle enhancement in convergence zones, while eddies showed the strongest depletion of grains. Only statistical theories of particle dispersion have been adopted for diffusion and relative velocity calculations in protoplanetary accretion disks (Morfill 1983, Dubrulle et al. 1995, Cuzzi et al. 1993). Here the turbulence is treated in a purely statistical sense as a source of random forces applied to the particles. These theories always predict that initially uniformly distributed particles in a turbulent field will remain uniformly distributed (Squires and Eaton 1991). This prediction is only valid for scales larger than the scales of the turbulence itself. 2.2. Friction Time The friction time tf is the response time of a particle to the relative motion of the surrounding gas or liquid. In the absence of external forces, any relative velocity will be damped according to

S D

t v 5 v0 exp 2 . tf

(2)

The drag force acting on a particle is Fd 5 2m

v 2 Vg , tf

(3)

where v is the particle velocity and Vg the gas velocity. If there is an acceleration a due to an external force, the particle will gain an equilibrium velocity ve 5 tf a,

(4)

because the friction and the external force cancel each other. We followed the formulation of drag forces by Weidenschilling (1977). In the Epstein regime, the mean

215

PARTICLE TRAPPING EDDIES

free path of the gas molecules l is much larger than the size of the dust particles (Epstein 1923). We consider this regime as long as the Knudsen number Kn 5 l /ap . 4/9. In the case of spherical objects, the friction time depends on the radius ap, the density of the particle rd, the thermal velocity of the gas vt, and the density of the gas rg :

tf 5

ap rd . vt rg

(5)

For a better physical understanding one can reformulate this relation as a dependence of the friction time on the mass mp and the cross-section area A of the particle. dP is a momentum transfer coefficient and characterizes different ways of reflection of the impinging gas molecules,

tf 5

3 mp . 4 dPAvt rg

(6)

In the case of spherical particles dP obtains values between dP 5 1 for specular reflection and dP 5 1 1 f/8 for diffuse reflection. Recent experimental and theoretical investigations suggest that this relation even holds for nonspherical particles and dust aggregates with dP P 1.2 ... 1.3 (Blum et al. 1996). For larger particles, i.e., a few centimeters at 1 AU, the Epstein regime is no longer valid for Kn , 4/9, and the friction time depends on the Reynolds number Re 5 2ap rgv/ h where h denotes the molecular viscosity of the gas. According to Whipple (1972), the friction time for a sphere can be written as 2rd a2p tf 5 for Re , 1. 9h

Here, V is the Keplerian frequency at radius r and z the height of the particle above the midplane using a cylindrical coordinate system with the origin at the central star. The gravitational acceleration is gz 5 2V2z.

Therefore, particles above or below the midplane settle within a typical time ts.

S D

t z 5 z0 exp 2 , ts

This corresponds to Stokes’ drag law. In what follows, we will use the fact that in this regime the friction time tf also does not depend on v. In our calculations Re , 1 is true at 1 AU as long as ap v , 3 105 cm2 sec21. This is valid up to meter-sized objects at 100 m sec21, and we are not interested in the treatment of larger bodies or higher velocities. In the analytical estimates of the particle motion, we assume that tf is constant, whereas in reality it is a function of vt and rg which primarily depend on the location in the disk and have special values for a selected disk region.

(10)

where ts is

ts 5

1 . tf V2

(11)

Furthermore, the rotation of the disk is supported by the gas pressure, leading to a sub-Keplerian rotation of the accretion disk. For this reason the particles undergo an orbital decay with velocity vr . For small particles with tf V ! 1, the radial velocity will be (Weidenschilling 1977) vr 5 22tf V DV.

(12)

There, DV is the difference between the Kepler velocity and the actual gas velocity. DV can be written as DV P

(7)

(9)

c 2s , vk

(13)

where cs is the speed of sound, and vk is the Keplerian velocity. The velocity difference DV is typically on the order of a few times 1023vk. Except for zones close to the midplane, the radial drift velocity is always smaller than the settling velocity. For objects fulfilling the criterion Vtf 5 1, the orbital decay rate has a maximum vr 5 2DV.

(14)

If the flow in the accretion disk is laminar all dust grains will settle toward the midplane of the disk with a velocity

In the inner part of the disk, these objects are meter-sized bodies with a radial speed of 104 cm sec21 (see Weidenschilling 1977). These bodies will be accreted by the central object within a few hundred years. Larger bodies, with Vtf @ 1, will have a slower orbital decay. They decouple from the surrounding gas only enduring some head wind. Their decay is given by

vz 5 2tf V2z.

vr 5 22 DV/(tf V).

2.3. Velocities without Turbulence

(8)

(15)

216

KLAHR AND HENNING

We want to stress again that all these calculations have been done assuming ... neither turbulent motions nor large scale circulations...’’ (Adachi et al. 1976). 3. GRAINS IN ACCRETION DISK EDDIES

In the standard framework of particle transport via a turbulent flow, eddies are always zones with a particle density smaller than any other flow feature (Squires and Eaton 1991, Eaton and Fessler 1994). Particles cannot follow a curved streamline but they spiral away from the center of the vortex (see Fig. 1a) due to centrifugal forces. Only in convergence zones (see Fig. 1b) are particles concentrated by the same mechanism. In the following, we will show that in protoplanetary accretion disks one might have an even higher particle concentration in the center of an eddy. Before presenting a mathematical treatment of this effect, we want to give a more vivid explanation of the concentration mechanism. 3.1. The General Mechanism The effect depends strongly on the gradient of the vertical component of the gravitational force. Due to this gradient, the particles are falling faster in the upper half of a turbulent eddy than in the lower half. If this concentration is stronger than the dispersion of grains due to centrifugal forces, then the particles are concentrated in an effective way. Therefore, the concentration effect will most likely occur in eddies of large extent and slow rotational velocity. In other words, for the right combination of time scales of settling and eddy turnover, there is a net motion toward an equilibrium point that is greater than the motion caused by the centrifugal force which tends to drive the grains outward. In Fig. 2 we show this effect by means of a force diagram. 3.2. Analytical Description A cylindrical coordinate system, such as in the previous section, will be used for describing the flow pattern of the eddy as well as the particle motion where r denotes the distance from the central object (the protostar) and z the height above the midplane of the accretion disk. In addition rotational symmetry is assumed; thus no quantity is dependent on u and ­u 5 0. We consider the most ‘‘general’’ eddy, being circular with radius l, bound to the r–z plane and rotating with frequency g. Its center is fixed at the radius r0 and at the height z0. We want to stress that the trapping capability does not depend on the circular shape of the eddies, but only on the existence of a closed vortical flow. Then the flow pattern of the eddy can be written as Vr 5 2(z 2 z0)g, Vz 5 (r 2 r0)g.

FIG. 1. Under laboratory-observed conditions particles (a) spiral out of the vortices (eddy), but (b) collect in convergence zones (strain fields).

As we derive in Appendix A.1.1, a particle of a certain size released in this eddy circulates in a similar way as the gas, and its motion can be described by r 5 a(t) sin(w (t)), z 5 a(t) cos(w (t)).

(16)

(17)

The radius of the circulation a varies with time as does

217

PARTICLE TRAPPING EDDIES

rs 5 tf V2z0 /g.

FIG. 2. The mechanism of particle concentration can be understood by considering the acting forces. Four particle positions are displayed. The dashed arrows indicate the local gas velocity. The solid arrows stand for the gas–particle relative velocities which are induced by the sum of gravitational and centrifugal forces (i.e., the acceleration (Fg 1 Fc)/mdust times the friction time). Thus, the solid arrows indicate the actual particle movement. On location 1, where the particle has the largest height above the midplane, the gravitation dominates over the centrifugal force and the particle lifting gas drag. Here, the particle moves closer toward the equilibrium point (indicated by a dashed cross). At point 2 as well as at point 4, the centrifugal force causes the particle to drift away from this equilibrium point, while gravitation and gas drag point in a tangential direction, not contributing to a radial motion. Finally location number 3 is dominated by the upward pushing gas drag, as the particle is the closest to the midplane and the gravitational force the smallest. At points 1 and 3 the particle is spiraling inward, while at points 2 and 4 it is spiraling outward. As long as the centrifugal effects are smaller than the effects caused by the gradient in the gravitational acceleration, the particle undergoes a net motion toward its equilibrium point.

(19)

Only in the midplane of the accretion disk z0 5 0 both rotational motions have the same center. This shift in particle motion reduces the eddie’s capabilities in trapping dust grains. As one can see from Eq. (19) and Fig. 3, the particle path can possibly cross the boundary of the eddy and thus the particle gets lost to the eddy trap. Thus the trapping depends on the initial location of the particle within the eddy. In a ‘‘normal’’ eddy (i.e., an eddy in an earth-bound system, without the special gravitational field of a Keplerian accretion disk) this circling dust would be spiraling outward, due to centrifugal forces. However in the accretion disk, the gradient of the gravitational acceleration leads to a stronger settling of the dust in the upper half of the circle than in the lower part. If the gradient is strong enough in comparison to the rotational frequency of the eddy, then we expect the particles to spiral inward. Hence this is a new possibility of concentrating dust in a protoplanetary accretion disk. All the particles of interest are in the regime of small Stokes numbers St 5 tfg ! 1. The capability of the eddy for collecting particles depends only on the relation between the Keplerian and the eddy frequency. Particles with Stokes numbers larger than unity are moving on Keplerian orbits perturbed by gas drag and can be neglected for our investigation. The general particle motion is described by a circulation with an increasing or decreasing radius a a 5 a0 exp(tf ct),

(20)

the frequency w˙ (t). This variation of a will be called radial motion throughout our work, while radial drift will always stand for the orbital decay. An eddy will prevent dust from settling toward the midplane as long as the dust has a settling velocity smaller than the rotational velocity of the eddy. This can be expressed as a relation between the friction time, the height of the eddy center above the midplane, the radial position in the accretion disk, or the orbital Keplerian frequency V, eddy size, and frequency (see derivation of Eq. (50) in the Appendix),

tf , lg/(V2z0).

(18)

Particles fulfilling condition (18) will rotate with the gas as long as the eddy lives. But (as it is explained by Eq. (53)) eddy and particle rotate around different centers. The rotation of the grain is shifted by an amount rs ,

FIG. 3. The possibility of getting trapped depends on the initial position of the particle within the eddy. If the eddy center is not in the midplane of the accretion disk, then the equilibrium point (or rotational center) of the particle is no longer at the center of the eddy (for details see text).

218

KLAHR AND HENNING

where c (see Appendix A for the derivation) is 1 c 5 g2 2 V2. 2

(21)

The smaller the rotational frequency of the eddy is, the stronger is the concentration capability. As long as V @ g the drift velocity toward the equilibrium point is proportional ptf . The eddy criterium for particles spiraling inward is 1 V. g, Ï2

(22)

It is obvious that the effect is independent of tf . Here, the dependence of tf on z has been neglected. However, as tf increases with its height above the midplane, the velocity relative to the gas does not increase linearly anymore, but as vz p zb with b . 1, thus supporting the general effect. As soon as an eddy is rotating significantly faster than the critical value for concentrating particles, one can neglect V. By this means it is possible to investigate the particle concentration between eddies, as it arises due to the decreasing particle density in vorticity zones. Until now, this has been investigated only numerically (Dobrovolskis et al. 1993, reviewed by Eaton and Fessler 1994). Equation (21) shows that for small Stokes numbers (tf c) is increasing with ptf while in the case of high Stokes numbers it is decreasing: (tf c) p tf21/3. Thus the maximum effect occurs for Stokes numbers around unity, as for gtf 5 1/ Ï2 both relations are giving the same result. The Stokes number St 5 1/ Ï2 for the maximal velocity out of the eddy, fits very well to the numerically obtained Stokes numbers for maximal concentration between the eddies. 3.3. Interpretation Now, we consider typical eddies in a protoplanetary accretion disk, in order to find scenarios for particle concentration. In the simplest case, there are large convective rolls rotating with a quite moderate frequency. Our 2D hydrodynamical calculations showed that these frequencies are 40 times smaller than the Keplerian frequency. That means that dust grains will be collected and concentrated in these huge eddies. The velocities in the macroscopic flow exceed 103 cm sec21. Hence, quite large ‘‘grains’’ of up to centimeters in size are able to be trapped and transported high above the midplane. Regarding the overall turbulent motion, one must consider the turbulence properties at the large scales L, i.e., the turbulent velocity U and the frequency gL 5 U/L. The turbulence follows a cascade down to small length scales l where the turbulent energy is dissipated by the molecular

friction (Richardson 1926). As a first approximation, a Kolmogorov spectrum has been adopted for the frequency distribution over the length scales, i.e., isotropic and homogeneous turbulence with a constant energy transfer rate « on all scales. This assumption is consistent with recent calculations of magneto-hydrodynamical turbulence in accretion disks (Stone et al. 1996). With « 5 U 3 /L 5 u3 /l

(23)

one obtains an eddy-size dependence of the rotational frequency of g(l ) 5 gL

SD l L

22/3

.

(24)

The large scale is assumed to be on the order of the pressure scale height, and the frequency of the large scale turbulence is some orders smaller than the Keplerian frequency. This depends on the a value which can be interpreted as the relation between turbulent speed and speed of sound. All large eddies fulfill the criterion for collecting particles. However, on smaller scales, as the frequency increases, this effect becomes less. At the critical value lcrit , eddies start to drive the dust out of their interior. The combination of Eq. (22) and Eq. (24) with the assumption gL 5 aV yields lcrit 5 L (aÏ2)2/3.

(25)

The smaller the a value, the wider is the range of eddy scales on which grains are collected. For an a as low as ˙ 5 1029 M( year21 0.5 3 1024 and a mass accretion rate of M even the turbulent eddies on the Taylor micro scale (i.e., the scale that dominates vorticity) are not able to drive dust out of their interior. Although the above used parameters are rather low for steady accretion, in the quiescent phase of FU-Orionis these values are not unrealistic (Bell and Lin 1994). 3.4. Efficiency The concentration of particles is limited by three mechanisms. These are the limited lifespan of the eddies, turbulent diffusion, and radial drift of the particles. 3.4.1. Lifespan of the eddies. Global convective rolls are assumed to be very stable, even when they are changing their position. The closed flow pattern remained in our calculations for many rotational periods. Turbulent eddies have a lifetime smaller than that of the convective rolls, because they transport energy from large scales to smaller ones and have no steady energy support by the temperature gradient like the convection. The eddy dissipates,

PARTICLE TRAPPING EDDIES

219

while the energy is transported to the smaller scales. This time is called integral time scale of the turbulence, which is in our case the lifespan of an eddy. According to Safronov (1969), this time is Y 5 1/V. According to Cuzzi et al. (1993) it is Y Q 1/g. We assume an eddy to rotate for about one period before the rotational frequency has changed significantly. In Appendix A.2.1, we derive the density concentration factor F per revolution, i.e., the relation between initial particle number density rp(0) and the number density after one eddy rotational period rp(Y ), F 5 exp(ftf V2 /g),

(26)

neglecting the diffusion of particles due to turbulent motion on length scales smaller than the eddy. As long as g is several times smaller than V, the density increases rapidly. For convective cells with (according to our hydrodynamical calculations) V/g 5 40, this factor is huge. Thus large areas are empty, while the grains are drifting toward an accumulation point. In Fig. 4 we plotted F for spherical particles of 0.1 cm radius at 1 AU over different length scales l of turbulence, where l0 refers to the largest scales: l0 5 HD . We show in ascending order the concentration after 0.5, 1, 1.5, and 2 eddy revolutions. The concentration is the strongest on the largest scale. After only one rotation of the eddy, there is a particle enhancement by a factor of 10 and after one more revolution this factor is above 100. On smaller scales and faster rotational frequencies the concentration capability of the eddies quickly decreases, as the centrifugal force gains importance. We stopped the process after two revolutions for a couple of reasons.

FIG. 5. Relative particle number density in comparison to the initial density versus different particle radii under the assumption of spherical particles. In ascending order after 0.5, 1, 1.5, and 2 eddy periods. This is a plot of the analytical solution without diffusion. See Eq. (26). The eddy size is one scale height. Its location is at 1 AU.

• First, there is no diffusion included. But with a particle enhancement by a factor of hundred and more, a steep gradient in the particle number density arises and diffusion becomes important (see Section 3.2.2). • Second, if particles are concentrated by such a huge factor, their influence on the gas dynamics cannot be neglected anymore. It must be checked by multicomponent fluid models, whether the eddy pattern will be destroyed, amplified, or stays the same. • Third, after the time for two eddy revolutions the eddy flow pattern can have changed significantly. Future work will have to implement all these features. Figure 5 shows the concentration factor for different particle sizes in large scale eddies. The plots refer again to 0.5, 1, 1.5, and 2 eddy revolutions. While the effect is the strongest for millimeter-sized particles, it decreases quickly with smaller grains. We did not plot results for larger objects, though the concentration effect would have been even stronger, because for these particles other effects become important and they cannot be trapped and concentrated qualitatively (see Section 3.2.3). 3.4.2. Turbulent diffusivity. Turbulence on smaller length scales than the eddy will drive a diffusion process. However, turbulence on larger length scales can be neglected for our considerations. This is obvious for two reasons:

FIG. 4. Relative particle number density in comparison to the initial density over different eddy sizes after (in ascending order) 0.5, 1, 1.5, and 2 eddy periods. This is a plot of the analytical solution without diffusion. See Eq. (26). The radial position in the disk is 1 AU.

• The velocity fluctuations on larger scales disappear on smaller scales. • The time for a velocity variation of the larger scale turbulence is larger than the typical time for the concentration process, i.e., the rotational period of the eddy.

220

KLAHR AND HENNING

Thus any large scale turbulence can be regarded as a superimposed systematic and static flow not affecting small scale processes. In Appendix A.2.2, a maximal concentration which will be obtained after the mechanisms of diffusing and advecting have come to an equilibrium was derived. The solution is a distribution around the center of the particle circulation in the form of a Gaussian function

S

D

2tf (2c) 2 a . 2D

rp(a) 5 rmax exp

(27)

Here, D denotes the diffusivity due to turbulence on subeddy scales. The maximal concentration is

rmax 5 r ini

S

tf (2c)l 2

S

2D 1 2 exp 2

DD

tf (2c) 2 l 2D

,

(28)

where rini denotes the initially homogeneous distribution and rmax the maximal achievable concentration in the center of the distribution. The concentration factor increases with tf and with c. In a first approximation the diffusivity D is proportional to the kinematic viscosity D 5 n /Sc,

(29)

where Sc is the Schmidt number which gives the coupling efficiency of a particle with a certain friction time to the turbulence. The maximum coupling occurs for Sc 5 1 (small-particle case), i.e., particles with friction times smaller than the time scale of the turbulent motion. The kinematic viscosity can be assumed as the product of a typical length and a typical velocity of the turbulent motion,

n 5 lt vt 5 l 2tgt .

(30)

For the assumption of turbulent diffusion, gt is expected to be significantly larger then g, and lt smaller than the eddy size l. Parameterizing this subeddy turbulence by gt 5 «g, lt 5 «

(31)

23/2

l,

(32)

with « . 1 determines the diffusion constant to be D 5 «22l 2g/Sc.

(33)

While a determines the absolute efficiency of the large scale turbulence by scaling the large scale turbulent veloc-

FIG. 6. Relative particle number density in comparison to the initial density over the radius a of the concentration field for different particle radii ap . In ascending order 0.01, 0.1, and 1.0 cm. Including the diffusion process. The eddy size is one scale height. This is the steady-state solution after unlimited time, due to competition between concentration and diffusion.

ity with the speed of sound, « divides the fluctuating part of turbulence from the nonfluctuating part of turbulence, at least over one eddy period. Figure 6 shows the maximal possible number density enhancement for different particle sizes. The radial position in the disk is 1 AU and « is assumed to be 10; thus one order of magnitude lies between the typical time for advection and diffusion. The local concentration factors for 0.01-, 0.1-, and 1-cm grains are plotted over the radius of the eddy, thus showing the increase of concentration toward the equilibrium point, or in this case the center of the eddy. Millimeter-sized particles are concentrated by a factor of up to 100. At radii larger than 25% of the eddy size the particle density is strongly depleted. Thus about 90% of the eddy volume are free of millimeter-sized particles. Centimeter-sized particles would be concentrated even stronger, but their concentration is limited by effects other than diffusion, as we already mentioned. If the particles are only 0.1 mm in size, their concentration is limited to a factor of 10. Their concentration is already limited by the lifetime of eddies, as we have shown in the previous section. It would take at least 10 eddy rotations to produce a 10-fold concentration and we do not know, whether we can assume such a long lifespan for vortical flows. If one considers only the largest eddies, where l 5 HD and g 5 aV, Eq. (28) simplifies to F5

«2 tf V . 4 a

(34)

PARTICLE TRAPPING EDDIES

221

This result is similar to the results of Dubrulle et al. (1995) F 5 Ïtf V/a, which gives the maximum number density enhancement of particles, if turbulence is working against settling toward the midplane of the accretion disk. However, this effect is one dimensional. Thus the possible concentration is smaller than for collecting particles in eddies. 3.4.3. Orbital decay. The radial drift toward the protostar becomes important for compact particles in the size of centimeters at 1 AU. Only grains with a radial drift velocity smaller than the rotational speed of the eddy can be trapped inside. This is the same mechanism preventing particles from settling toward the midplane. Using the formulation for the radial drift velocity in Eq. (12), we can modify Eq. (18) and get

tf , lg/(2V DV ).

(35)

In Section 4.1.3 and Fig. 8 we will show tracks of particles of different size, which will elucidate this effect. 3.5. Brownian Motion Brownian diffusion has the same effect as turbulent diffusion. Both diffusion processes can be mathematically treated the same way by a constant D denoting the diffusivity. Nevertheless, in our case, the Brownian diffusivity is always by orders of magnitude smaller than the turbulent diffusivity. Only in the case of scales below the turbulent micro scale, Brownian motion gains its importance. 4. NUMERICAL MODELING

In order to test our analytical estimates, we solved the full set of differential equations (90)–(95) using numerical integration schemes (see Appendix B). In the numerical calculations, it is possible to include the dependency of tf on r and z, sub-Keplerian rotation, Coriolis forces, and ‘‘real’’ flow pattern.

FIG. 7. (a) Track of a 0.1-cm spherical particle in a large scale eddy at 1AU. (b) Normalized radial velocity a˙ /a; the mean value is given by the solid line and the theoretical prediction by the dotted line.

4.1. Ideal Eddy Flow Pattern The flow pattern used in this part is not the result of a hydrodynamical calculation, but an idealized eddy, as given by Eq. (16). This is necessary to check the criteria for collecting particles given by Eq. (22). All other disk quantities are the results of 2D hydrodynamical calculations of ˙ 5 1027 M( year21; a protoplanetary accretion disk (M 22 a 5 10 ) and stored in a table. The gas density in the midplane has typical values of r P 3 3 10210 g/cm3, the temperature is T P 600 K, and the surface density is o P 300 g/cm2 at 1 AU. The radial and vertical gas velocities are calculated according to Eq. (16). The following parameters are used: g0 5 0.01 V, l0 5 Hd , r0 5 1 AU, and z0 5 0.

4.1.1. Particle spiraling inward. The eddy parameters in this calculation are g 5 g0, l 5 l0. This represents a large scale eddy, where one expects to see the particletrapping in the most perfect way. Figure 7 shows the trace of a grain with a radius of 0.1 cm which was started at a0 5 l/2 circling around its stable point and spiraling inward. The simulation follows the particle motion during three eddy rotations. This equals to 300 years. 4.1.2. Radial motion. Figure 7a states the qualitative behavior of the grain only. In order to test the particle motion quantitatively, we calculated the normalized radial component of the particle motion in respect to the equilib-

222

KLAHR AND HENNING

FIG. 8. Track of 0.5- to 3.0-cm particles in a large scale eddy. At about 1 cm in radius, the radial drift starts to dominate over the collecting process. The smaller particles are the spiraling particles at the right. The larger particles are those which move to the left.

rium point. This can be interpreted as the relative radius change per eddy period T. Figure 7b compares the time evolution of this velocity (a˙ /a)T with the mean value ka˙ /alTT and the analytical prediction. The quantity a˙ /a fluctuates strongly due to the varying gravitational field, additional effects of nonconstant tf , the Coriolis force, and a sub-Keplerian orbital frequency. Nevertheless, the timeaveraged value is very close to the analytical prediction. 4.1.3. Orbital decay. The radial drift toward the protostar limits the possibility of trapping particles. Here (Fig. 8) particle tracks are calculated for grains from 0.5 to 3 cm in size. They were all started at the same point a0 5 l/2. While the smaller particles are trapped in the eddy, the larger particles escape in the direction toward the central object (left). The equilibrium point is shifting with increasing friction time of the particles, to a situation, where this point would be outside of the eddy. 4.1.4. Eddy variation. The capability of particle trapping depends on the eddy frequency. Here, particles of 0.1 cm in size are used and put in different eddies with an increasing frequency. Figure 9 shows the time-averaged ka˙ /alTT for a wide range of eddy sizes and the analytical values from Eq. (20). The eddy size ranges from Hd to 1023 Hd . The frequencies are calculated via a Kolmogorov spectrum aV , g , V, with a 5 1022. The radial velocity has its maximum at the large scales, and it drops very fast with decreasing scales. If the eddy has a Keplerian frequency, the velocity switches the sign, denoting a drift out of the eddy. Analytical estimates (dotted line) and the numerical results (solid line) are in good agreement.

FIG. 9. Normalized radial velocity a˙ /(ag) for a range of eddy sizes; the mean value is given by the solid line and the theoretical prediction by the dotted line.

4.1.5. Particle variation. The particle trapping also depends on the particles friction time (see Eq. (20)). Here, one large scale eddy is used as in (4.1.1), and the particle size is varied from micrometers to centimeters (see Fig. 10). Analytical prediction (dotted line) and numerical simulation (solid line) of the mean radial velocity are in good agreement for particle sizes up to a few millimeters in size. Larger particles approach their equilibrium point in a time shorter than one eddy period; thus the process of time

FIG. 10. Normalized radial velocity a˙ /(ag) for a range of particle sizes; the mean value is given by the solid line and the theoretical prediction by the dotted line.

223

PARTICLE TRAPPING EDDIES

averaging is no longer correct and the predicted value is an overestimation. 4.2. Combination with 2D Hydrodynamics Several models simulated thermal convection in a protoplanetary accretion disk. As example for 2D radiation hydrodynamical simulations we want to mention the recent studies by Rozyczka et al. (1994) and Kley et al. (1993). So far the only 3D calculation, but without proper dissipation and correct radiative transfer, was presented by Stone and Balbus (1996). This simulation suggests that the flow pattern calculated in 2D simulations is also an element of the more general 3D flow pattern, and the 2D restriction of the problem seems to be validated. For our calculations, we ourselves used the 2D code of Kley (1987, 1989). For a comprehensive discussion concerning the simulation of the thermal convection, we refer to Kley et al. (1993) where the details of the code as well as the physical background are extensively described. The interested reader will find a brief description of the computational method in Appendix B of our paper. The input parameters of our calculation were the accre˙ 5 1028 M( /year and the efficiency of the viscostion rate M ity a 5 1023. We simulated a segment of the accretion disk reaching from 0.85 out to 1.15 AU. We chose a resolution of 122 radial grid points and 60 grid cells in the vertical direction. The initial condition was obtained as a self-similar solution assuming local thermal equilibrium by using a 1D vertical accretion disk model. This led to a gas density in the midplane of r P 8 3 10210 g/cm3, a temperature of T P 350 K, and a surface density of o P 500 g/cm2. Figure 11a shows the resulting flow pattern after more than 600,000 computational steps, which equals to a physical time of almost 1400 years or more than 1000 orbital periods at the outer radius. This long time is necessary to damp the initial distortions with the small amount of viscosity assumed to be present in the cold disk and reach a stationary solution with an almost constant accretion rate over the entire disk radius. From the theory of thermal convection in protoplanetary accretion disks and our chosen parameters, we already expected to simulate a transition zone between a nonconvective inner part (left) and a convective outer part (right). Thus, the inner part is dominated by a viscously determined inflow–outflow feature, while with increasing radius the typical cells for thermal convection gain prominence. For this flow pattern we calculated the traces of 1000 randomly distributed spherical particles of 0.1 cm radius (see Fig. 11b). For our nebula parameters these particles will have a typical friction time of about tf P 104 s. The particle motion was calculated over a period of 160 years, which corresponds to about four revolutions of the convec-

tive cells. This clearly indicates the fulfillment of the trapping condition g , V (i.e., g 5 5 3 1029 and V 5 2 3 1027). The term ftf V2 /g (see Eq. (72)) characterizing the particle concentration per eddy revolution is in the order of unity which gives cause for the expectation of effective particle trapping. The trapping condition Eq. (35) is also fulfilled as the radial drift velocity on the order of 10 cm/ sec and the sedimentation velocity at z 5 HD on the order of 103 cm/sec while the convective gas velocity is on the order of 102 –105 cm/sec, depending on the height in the disk. In Fig. 11c, the traces of 500 particles are shown (for technical reasons). One can distinguish between sedimentation in the nonconvective region and circulation in the convective region. The particles collect in the convective cells, but not in the convergence zones at the midplane, which would be expected in earth-like turbulence. Figure 11d shows the position of the 1000 grains after 160 years. The particles tend to concentrate at their equilibrium points in the convective cells far above the midplane. This simulation demonstrates clearly that dust trapping eddies need not be perfectly circular in shape. To accomplish this simulation we varied the particle sizes. First we used particles of 0.1 mm in size. In this case the trapping condition is fulfilled; nevertheless the typical concentration parameter is 10 times smaller, as it is proportional to tf , thus indicating a much weaker concentration process. We used the same initial distribution as in the 1-mm case. Figure 12a shows the location of the 1000 particles after 160 years, which is the same time we stopped the simulation in the case of the larger particles. One clearly sees the faint but almost negligible particle concentration around almost the same knots as in Fig. 11d. As one treats the behavior of larger particles, i.e., of 1 cm radius, one expects the particles not to participate in the concentration process for most initial particle locations. The radial drift velocity is on the order of 102 cm/sec and the sedimentation velocity at z 5 HD on the order of 104 cm/sec while the convective gas velocity is still on the order of only 102 –105 cm/sec. Only a few particles get the chance to be trapped by a locally stronger flow, and as one can see in Fig. 12b most particles tend to sediment toward the midplane. This situation and the effects of turbulent diffusion working against sedimentation were already described in several studies (Vo¨lk et al. 1980, Cuzzi et al. 1993, Weidenschilling and Cuzzi 1993, Dubrulle et al. 1995), making the high concentration of particles close to the midplane an artifact as we neglected the assumed turbulence in the nonconvective part of the disk. 5. DISCUSSION

We have shown that there is a very effective mechanism for concentrating particles in a nonlaminar accretion disk:

224

KLAHR AND HENNING

FIG. 11. (a) Mass flux in the convective region of a protoplanetary accretion disk. This is the result of a 2D radiation hydrodynamical simulation. (b) The homogeneous initial distribution of 0.1-cm grains. (c) The traces of a part of the grains during 160 years. (d) The position of the grains after 160 years.

1. Particles cannot only be concentrated by settling, thus, forming a particle subdisk in the midplane of the accretion disk. They can also collect in eddies far above the midplane. By this mechanism, convective regions of the protoplanetary accretion disk become very interesting places for particle processing. 2. The concentration capability increases exponentially with tf up to a critical value, where the radial drift and/or the settling velocity exceed the velocity of the non-laminar flow. Thus the concentration is strongly size segregating, which leads to a narrow size distribution in convective cells. 3. While on small scales, the particles behave like in earth-like turbulence (i.e., particles tend to leave the eddy and concentrate in convergence zones); this behavior is diminished and reversed with increasing scales. Thus the

picture of particles interacting with a turbulent flow pattern must be checked for the case of the primordial solar nebula. We gave four equations which determine the qualitative and quantitative properties of this effect. The first two equations (see Eq. (18) and Eq. (35)) indicate whether a particle can be trapped at all. Additionally Eq. (18) describes the dependency of the effect on the eddy frequency and Eq. (26) the dependency of the strength of the effect on the Stokes number which can be rewritten in the form F 5 exp(fStV2 /g) 5 exp(fSta22).

(36)

This will have implications on the processes in the solar nebula in various ways:

225

PARTICLE TRAPPING EDDIES

A disadvantage of our work is the fact that the mechanism of turbulence for a protoplanetary accretion disk is still not understood. Neither numerical nor analytical models are able to deal with high Reynolds numbers on the order of 1010. Computer power must be increased by orders of magnitude before an accretion disk without ad hoc a turbulence can be calculated. Future work might adopt more sophisticated analytical models of turbulence. We will direct our interest toward the convective cells, because there is little known about their three-dimensional shape and dynamical evolution. Recent calculations performed by Stone et al. (1996) indicated that convective cells are of limited azimuthal extension and of a spiral shape. The consequences of these properties of convection will be investigated as soon as the hydrodynamical calculations deliver reliable results. APPENDIX A

Analytical Treatment A.1. Analytical Treatment of Particles in Eddies The most general eddy is bound to the r–z plane with a radius l and a rotation frequency g. The center of the vorticity is located at a radius r0 and a height z0 . The flow pattern of the eddy is described by Vr 5 2zg, Vz 5 rg.

(37)

The general equation of motion for a particle in a two-dimensional gas flow is given by mr¨ 5 2m(r˙ 2 Vr)/ tf 1 Fr , mz¨ 5 2m(z˙ 2 Vz)/ tf 1 Fz .

FIG. 12. Here a calculation as for Fig. 11 was repeated with the same flow pattern and the same initial position of particles, but with different particle sizes. (a) This is the position of 0.1-mm grains after 160 years. The particle concentration is very faint. (b) For larger particles, i.e., 1.0 cm radius, the sedimentation and radial drift dominates over the concentration process for most of the initial positions of the grains. Again we show the position of the grains after 160 years.

Here r˙ and z˙ denote the particle velocities, the first term of the right side of the equations represents the friction force, and F can denote some additional external forces. A.1.1. No gravitation. The equation of motion for the particle in the absence of external forces reads r¨ 5 2(r˙ 1 zg)/ tf , z¨ 5 2(z˙ 2 rg)/ tf .

1. Grain growth: The number density enhancement by a factor of 100 will strongly influence the time scales of coagulation. 2. Chondrule formation: Only grains of the size of chondrules (or with an equivalent friction time) concentrate in convective zones. 3. Charge separation: If grains of a different size have distinct charges, the size segregation can lead to strong electric fields. In this case, the charge separation is not driven by sedimentation but by convection; thus the energy provided for this process is much larger.

(38)

(39)

A small grain with the Stokes number St 5 tf g ! 1 is coupled to the gas motion; i.e., the particle rotates with the same frequency g as the eddy. The particle positions are r 5 a(t) sin(gt) z 5 a(t) cos(gt),

(40)

where a denotes the distance between the particle position and the eddy center. Due to the centrifugal force a¨ 5 2a˙ / tf 1 g2a,

(41)

the dust grain spirals outward. Neglecting the second derivative (adiabatic elimination of fast variables) one gets

226

KLAHR AND HENNING a˙ 5 tf g2a,

(42)

and hence,

one gets rid of the constant term: r¨9 5 2(r˙9 1 zg)/ tf ,

a 5 a0 exp(tf g2t).

(54)

z¨ 5 2(z˙ 2 r9g)/ tf 2 V z. 2

(43)

(55)

There is a unstable equilibrium point at the center of the circulation (a0 5 0 leads always to a(t) 5 0).

These equations are valid for any height z0 of the eddy center above the midplane. After introducing polar coordinates

A.1.2. Constant gravitation. Including a gravitational field in the equation of motion

r9 5 a(t) sin(w),

(56)

z 5 a(t) cos(w),

(57)

r¨ 5 2(r˙ 1 zg)/ tf ,

(44)

z¨ 5 2(z˙ 2 rg)/ tf 1 g.

(45)

According to the previous section, there is an equilibrium point at the origin of the coordinate system. Introducing the gravitational field, this equilibrium point shifts to rs 5 2gtf /g.

(46)

In this point, the drag of the streaming gas and the gravitation cancel each other. Shifting the coordinate system to this new point r R r9 2 rs, Eq. (45) results in z¨ 5 2[z˙ 2 (r9 2 gtf /g)g]/ tf 1 g 5 2(z˙ 2 r9g)/ tf .

(47)

Now, it is obvious that all considerations for the case without gravitation are still valid for the case with constant gravity. The dust grains are spiraling outward with a velocity dependent on their friction time and the rotational velocity of the eddy. In contrast to the case without external forces, the dust grains no longer co-rotate with the gas but rotate around their equilibrium point, depending on the particle mass. Due to the finite size of an eddy, there is a critical mass where the particles leave the eddy. If a dust grain has a friction time of

tf . lg/(2g),

and neglecting the second derivative of a, one obtains a˙ 5 tf w˙ 2a 2 tf V2z cos w.

The cosine term arises from the fact that only the radial component of the gravitational acceleration acts on a˙. Replacing z R a cos(w) gives a˙ 5 atf (w˙ 2 2 V2 cos2 w).

S

D

1 ka˙lw 5 atf w˙ 2 2 V2 . 2

(60)

For the small grains with gtf ! 1 and w˙ 5 g, the critical value for diffusion or collecting the particles is

(48)

A.1.3. Harmonic potential. In the framework of accretion disk eddies, the gravitational force cannot be considered to be constant. The gravitational acceleration depends on the distance z to the midplane of the disk.

(59)

The term within the brackets can be positive as well as negative, corresponding to outward or inward spiral motion. As we are only interested in the net change of the rotational radius a, during a whole cycle, the integration over w results in a mean value for

gcrit 5

its equilibrium point would be outside the eddy. This is a sufficient condition for a grain to leave an eddy, but not a necessary one. Depending on their initial positions, grains that meet these conditions can still escape.

(58)

1 Ï2

V.

(61)

The smaller the rotational frequency of the eddy is, the stronger is the concentration ‘‘capability.’’ The critical frequency does not depend on the friction time for small particles, In the following, the abbreviation a˙ 5 tf ca

(62)

will be used, where c stands for gz 5 2V2z

(49)

represents an harmonic potential. The condition (48) for a dust grain to be able to leave the eddy changes to

tf . lg/(V2z0),

(50)

where z0 denotes the height of the eddy center above the midplane. Here the same remarks as for Eq. (48) are valid. The equations of motion are now

1 c 5 g2 2 V2. 2

A.2. Estimate for the Maximal Concentration A.2.1. No diffusion. If diffusive processes working against the particle concentration are neglected, the concentration depends only on the lifetime of the eddy. Starting with the continuity equation

r˙ p 1 =(rpvp) 5 0,

r¨ 5 2(r˙ 1 zg)/ tf ,

(51)

z¨ 5 2(z˙ 2 rg)/ tf 2 V2z 2 V2z0 .

(52)

(63)

(64)

introducing a polar coordinate system and assuming rotational symmetry, this equation can be written as

Shifting the coordinate system r R r9 1 rs rs 5 tf V2z0 /g,

(53)

1 r˙ p(a) 1 ­a (arp (a)vp(a)) 5 0. a

(65)

227

PARTICLE TRAPPING EDDIES

rp 5 krpl 1 rp9 , vp 5 kvpl 1 v9p .

The velocity vp(a) 5 a˙ can be expressed as in Eq. (62) 1 r˙ p (a) 1 ­a (rp (a)tf ca2) 5 0 a

(66)

Expanding and Reynolds averaging of the continuity equation yields kr˙ pl 5 2=(krpl kvpl 2 krp9 v9pl).

which leads to

r˙ p (a) 5 22tf crp (a) 2 tf ca­a rp (a).

(67)

Assuming that the initial distribution is homogeneous (­arp (a, t 5 0) 5 0), one gets for t 5 0 1 dt

r˙ p (a, 0 1 dt) 5 2rp (a, 0) exp(22tf cdt),

(68)

and the distribution will remain homogeneous as every rp(a) is increased by the same amount. Setting rp(a, 0) 5 rp ini

rp (t) 5 rp ini exp(22tf ct).

(69)

For small particles in large eddies, g can be neglected:

rp (t) 5 rp ini exp(tf V2t).

(70)

During the lifetime T Q g21 there is a concentration by a factor F(t) 5 rp(t)/ rp ini F 5 exp(tf V /2g). 2

(71)

This factor must be regarded as an upper limit. The effect of particles to get lost by crossing the eddy boundary increases with their friction time. Also eddies that are far away from the central plane will be less effective at collecting grains because particles tend to fall out of them. Another useful expression is the concentration per revolution of the eddy which is given by F 5 exp(ftf V2 /g).

(72)

Thus, for V . g, all particles in an eddy of unlimited life time in the center of the disk could collect in the same point. The density, while homogeneous over a decreasing area, diverges. A.2.2. Diffusion. An effect acting against this high concentration is turbulent and Brownian diffusion. First one must estimate the size of the diffusion coefficient for both cases. The diffusion coefficient D can be assumed as the product of a typical length and a typical velocity. Brownian motion. For Brownian motion one can use the result derived by Einstein (1905) D 5 kT/(6fhap),

(76)

(77)

Here, the gradient diffusion hypothesis (GDH) is adopted, which assumes that the particle flux jp resulting from the fluctuating term is proportional to the mean gradient of the particle density (see Cuzzi et al. 1993): k jpl 5 krp9 v9pl 5 D=krpl.

(78)

The diffusion coefficient is assumed to be proportional to the kinematic viscosity nt , D 5 nt /Sc,

(79)

where Sc is a dimensionless coefficient called the Schmidt number, expressing the coupling of the particle to the turbulence. As only mean values remained, the bracket notation will no longer be used and (77) simplifies to

r˙ p 5 2=(rpvp 2 D=rp).

(80)

In a polar coordinate system and assuming rotational symmetry this can be written as 1 r˙ p 5 2 ­a (arpvp 2 aD­arp). a

(81)

Using the relation vp 5 a˙ according to Eq. (62) yields 1 r˙ p 5 2 ­a (tf crpa2 2 aD­arp). a

(82)

The term in brackets is easily identified as the local particle flux jo and we get

r˙ pa 5 o˙ p 5 2­a( jo).

(83)

Looking for a stationary solution, we set r˙ p 5 0, yielding ­a ( jo) 5 0.

(84)

Therefore, jo is a constant. If one now assumes a closed system, where there is no flux over the eddy boundary jo(l ) 5 0, jo is always zero, and we get

(73) 0 5 tf crpa2 2 aD­a rp .

(85)

or, in the special case of the Epstein regime, D 5 kT

@S

D

4 frg a2p vt . 3

For a . 0 this is (74)

But for the particles and the length scales we are interested in, this diffusivity is always smaller than the diffusivity due to turbulent motion, and Brownian motion can be neglected further on. Turbulent motion. Starting with the continuity equation

r˙ p 1 =(rpvp) 5 0,

(75)

the particle density and the velocity are split up in a mean part and a fluctuating part (designated by a prime):

05

tf (2c) 2arp 1 ­a rp . 2D

(86)

Note that, as long as c is negative, the solution of Eq. (86) is a Gaussian function:

r(a)p 5 C exp

S

D

2tf (2c) 2 a . 2D

(87)

Multiplying with 2fa and integrating over a from 0 to l yields the column density Np of the dust distribution

228

KLAHR AND HENNING Np 5 C

S

S

tf (2c) 2 2fD l 1 2 exp 2 tf (2c) 2D

DD

.

(88)

Setting Np 5 fl 2 rp ini , where rp ini denotes the initial particle density, and identifying C with rp max, i.e., the maximal possible particle density, the concentration factor F can be determined to be F5

rp max 5 rp ini

S

tf (2c)l

2

S

DD

tf (2c) 2 l 2D 1 2 exp 2 2D

.

(89)

Due to an adaptive step-size procedure, the same results were obtained from both integrators.

ACKNOWLEDGMENTS The authors thank P. Artymovicz, K. R. Bell, J. Blum, J. N. Cuzzi, A. Heines, S. Kempf, W. Kley, R. Mucha, S. Pfalzner, R. Sablotny, and S. Weidenschilling for supporting us with their specific knowledge and their merciless criticism. A part of the numerical calculations has been performed at the computing center of the Friedrich Schiller University Jena and at the super power computing center HLRZ Ju¨lich.

APPENDIX B REFERENCES

Numerical Solution of the Equation of Motion In order to treat the particle motion in a realistic way, one first must calculate an accretion disk model to get spatially distributed values for density, temperature, and flow pattern.

B.1. Hydrodynamical Part The 2D hydrodynamical code we use was originally developed by Kley (1987, 1989). The set of equations to be solved consists of six partial differential equations for density, momentum, energy of the gas, and radiation energy. The numerical solution is obtained by applying a finite difference method on a fixed Eulerian grid which is of second-order accuracy in space and semi-second order in time. The program is independent of a special coordinate system; i.e., two-dimensional problems can be solved in Cartesian as well as in spherical or cylindrical coordinates. The flux-limited diffusion part is treated with SOR. This code calculates all the physical information we need for the dust tracer. The data will be stored in tabulated form and can be called by the particle program. The relevant values are density and temperature for the determination of the friction time, as well as a three-dimensional velocity vector describing the local gas flow.

B.2. Particle Part In the numerical code, no cylindrical coordinate system as in the analytical part is used but a spherical coordinate system (r, u, w) is. The reasons for this are that we use the same discretization in the hydrodynamical code, and the calculation of the gravitational force needs the smallest numerical effort in the spherical coordinate system. The equations of motion including gravitation and friction are given by

Adachi, I., C. Hayashi, and K. Nakazawa 1976. The gas drag effect on the elliptical motion of a solid body in the primordial solar nebula. Prog. Theor. Phys. 56, 1756–1771. Balbus, S. A., and J. F. Hawley 1991. A powerful local shear instability in weakly magnetized disks. I. Linear analysis. II. Nonlinear evolution. Astrophys. J. 376, 214–233. Bell, K. R., and D. N. C. Lin 1994. Using FU Orionis outbursts to constrain self-regulated protostellar disk models. Astrophys. J. 427, 987–1004. Blum, J., G. Wurm, S. Kempf, and Th. Henning 1996. The Brownian motion of dust particles in the solar nebula. An experimental approach to the problem of pre-planetary dust aggregation. Icarus 124, 441–451. Cabot, W., and J. B. Pollack 1992. Direct simulations of turbulent convection. II. Variable gravity and differential rotation. Geophys. Astrophys. Fluid Dynam. 55, 97–133. Cuzzi, J. N., A. R. Dobrovolskis, and J. M. Champney 1993. Particle-gas dynamics in the midplane of a protoplanetary nebula. Icarus 106, 102–134. Cuzzi, J. N., A. R. Dobrovolskis, and R. C. Hogan 1996. Turbulence, chondrules, and plantesimals. In Chondrules and the Protoplanetary Disk (R. Hewins, R. Jones, and E. R. D. Scott, Eds.). Cambridge Univ. Press, Cambridge. Dobrovolskis, A. R., J. N. Cuzzi, and R. C. Hogan 1993. Particle sorting and segregation in a preplanetary nebula. Bull. Am. Astron. Soc. 25, 1122. Dubrulle, B., G. Morfill, and M. Sterzik 1995. The dust subdisk in the proto-planetary nebula. Icarus 114, 237–246. Eaton, J. K., and J. R. Fessler 1994. Preferential concentration of particles by turbulence. Int. J. Multiphase. Flow 209, 169–209. ¨ ber die von der molekularkinetischen Theorie der Einstein, A. 1905. U Wa¨rme geforderten Bewegung von in ruhenden Flu¨ssigkeiten suspendierten Teilchen. Ann. Phys. 17, 549–560.

dr pr 5 , dt m

(90)

p du 5 u , dt mr2

(91)

pf df 5 , dt mr 2 cos2 u

(92)

M m mVgas 2 pr dpr 5 31 3 2 c s2 1 , dt mr mr cos2 u r t

(93)

p2f sin u rmGgas 2 pu dpu 52 2 1 , dt mr cos3 u t

(94)

Kley, W. 1989. Radiation hydrodynamics of the boundary layer in accretion disks. Astron. Astrophys. 208, 98–110.

dpf rm Hgas 2 pf 5 . dt t

(95)

Kley, W., and G. Hensler 1987. Two-dimensional numerical models of boundary layer of accretion disks in cataclysmic variables. Astron. Astrophys. 172, 124–142.

Depending on the strength of the friction term, it is useful to switch between an explicit integrator like Runge–Kutta for larger particles and an implicit method like the backward differencing scheme for small grains.

Kley, W., J. C. B. Papaloizou, and D. N. C. Lin 1993. On the momentum transport associated with convective eddies in accretion disks. Astrophys. J. 416, 679–689.

p2u

p2f

Epstein, P. 1923. On the resistance experienced by spheres in their motion through gases. Phys. Rev. 22, 710–733. Hunt, J. C. R., A. A. Wray, and P. Moin 1988. Eddies, streams, and convergence zones in turbulent flows. In Proc. of the Summer Program 1988, pp. 193–208. Center for Turbulence Research, Stanford University, Stanford, CA.

PARTICLE TRAPPING EDDIES

229

Kolmogorov, A. N. 1941. The local structure of turbulence in incompressible viscous fluid for very large Reynolds number. Dokl. Akad. Nauk SSSR 30, 9–13 [Reprinted in Proc. R. Soc. London A 434, 9–13 (1991)]

Stone, J. M., J. F. Hawley, C. F. Gammie, and S. A. Balbus 1996. Three dimensional magnetohydrodynamical simulations of vertically stratified accretion disks. Astrophys. J. 463, 656.

Lin, D. N. C., and J. Papaloizou 1980, On the structure and evolution of the primordial solar nebula. Mon. Not. R. Astron. Soc. 191, 37–48.

Tanga, P., A. Babiano, B. Dubrulle, and A. Provenzale 1996. Forming planetesimals in vortices. Icarus 121, 158.

Lissauer, J. J. 1993. Planet formation, Ann. Rev. Astron. Ap. 31, 129–174. Morfill, G. 1983. Physics and chemistry in the primitive solar nebula. Presented at Les Houches Summer School, Aug. 1983. Richardson, L. F. 1926. Atmospheric diffusion shown on a distance neighbor graph. Proc. R. Soc. London A 110, 709. Rozyczka, M., P. Bodenheimer, and K. R. Bell 1994. A numerical study of viscous flows in axisymmetric a-accretion disks. Astrophys. J. 423, 736–747. Safronov, V. S. 1969. Evolution of the Protoplanetary Cloud and Formation of the Earth and Planets. [NASA TTF-677] Shakura, N. I., and R. A. Sunyaev 1973. Black holes in binary systems. Observational appearance. Astron. Astrophys. 24, 337–355. Squires, K. D., and J. K. Eaton 1991. Preferential concentration of particles by turbulence. Phys. Fluids A 3, 1169–1178. Stone, J. M., and S. A. Balbus 1996. Angular momentum transport in accretion disks via convection. Astrophys. J. 464, 364.

Tennekes, H., and J. L. Lumley 1972. A First Course in Turbulence. MIT Press, Cambridge, MA. Vo¨lk, H. J., F. C. Jones, G. E. Morfill, and S. Ro¨ser 1980. Collisions between grains in a turbulent gas. Astron. Astrophys. 85, 316–325. Weidenschilling, S. J. 1977. Aerodynamics of solid bodies in the solar nebula. Mon. Not. R. Astron. Soc. 180, 57–70. Weidenschilling, S. J., and J. N. Cuzzi 1993. Planetesimal formation in the protoplanetary nebula. In Protostars and Planets III (E. H. Levy and J. L. Lunine, Eds.). Univ. of Arizona Press, Tucson. Whipple, F. L. 1972. On certain aerodynamic processes for asteroids and comets. In From Plasma to Planet, Proceedings of the Nobel Symposium 21 (A. Elvius, Ed.). Wiley, New York. Wood, J. A. 1963. Chondrites and chondrules. Sci. Am. Oct. 1963. Wray, A. A., and J. C. R. Hunt 1989. Algorithms for classification of turbulent structures. In Proc. IUTAM Conf. of Topology of Fluid Mechanics, pp. 95–104. Cambridge Univ. Press, Cambridge, UK.