Table Of ContentMon.Not.R.Astron.Soc.000,000–000 (0000) Printed2February2008 (MNLATEXstylefilev2.2)
Kinematics of solid particles in a turbulent protoplanetary
disc
Augusto Carballido1,⋆ James M. Stone2,† Neal J. Turner1 ‡
8 1 Jet Propulsion Laboratory, MS 169-506, California Institute of Technology, Pasadena CA 91109, USA
0 2 Department of Astrophysical Sciences, Princeton University,Princeton NJ 08544, USA
0
2
n
a 2February2008
J
3
2 ABSTRACT
We perform numerical simulations of solid particle motion in a shearing box model
]
of a protoplanetary disc. The accretion flow is turbulent due to the action of the
h
p magnetorotational instability. Aerodynamic drag on the particles is modelled using
- the Epstein law with the gas velocity interpolated to the particle position. The effect
o of the magnetohydrodynamic turbulence on particle velocity dispersions is quantified
tr for solidsof differentstopping times ts,or equivalently,differentsizes.The anisotropy
s oftheturbulenceisreflecteduponthedispersionsoftheparticlevelocitycomponents,
a
withtheradialcomponentlargerthanboththeazimuthalandverticalcomponentsfor
[
particles larger than ∼ 10 cm (assuming minimum-mass solar nebula conditions at 5
1 AU).Thedispersionoftheparticlevelocitymagnitude,aswellasthatoftheradialand
v azimuthal components, as functions of stopping time, agree with previous analytical
6
results for isotropicturbulence. The relativespeed betweenpairs ofparticles with the
4 same value of ts decays faster with decreasing separation than in the case of solids
6
with different stopping time. Correlations in the particle number density introduce a
3
non-uniform spatial distribution of solids in the 10 to 100 cm size range. Any clump
.
1 of particles is disrupted by the turbulence in less than one tenth of an orbital period,
0
and the maximally concentrated clumps are stable against self-gravitationalcollapse.
8
0 Keywords: accretion,accretiondiscs–Solarsystem:formation–planetarysystems:
: protoplanetary discs – MHD – turbulence
v
i
X
r
a
1 INTRODUCTION caused by Brownian motion, other sources of relative ve-
locities in a protoplanetary nebula may come into play as
Inorderforplanetesimalstoforminacircumstellardisc,the
the interaction between solids and gas decreases. Differen-
parent dust grains must meet a variety of kinematic condi-
tial sedimentation, radial drift and turbulence can provide
tionsiftheyaretogrow inmass andsize. Theformation of
thenecessarykineticenergiestodrivethecollisionalmotion
kilometre-sized solid bodies by gravitational instability re-
of solids, in addition to changing their spatial distribution
quires, through the Toomre criterion, that the dust spatial
throughout the nebular disc. In this paper we focus on the
density be ρ >Ω2/πG, where Ω is the disc angular veloc-
d effectofaturbulentgas flowonsolid particle velocities and
ityandGthegravitationalconstant.Suchdensitiescouldbe
spatial arrangement.
achieved in a mid-plane dust layer after sedimentation has
occurred. On the other hand, growth by collisional sticking The formalism of V¨olk et al. (1980) makes it possible
maybeaccomplishedifthegrainrelativespeedsacquirecer- to calculate the velocity dispersions of particles in isotropic
tain values: experiments have shown that, for micron-sized turbulence.Theeffectofthefluctuatinggasvelocityonindi-
grains,collisionspeedslessthan∼0.2m/sleadtoaggregates vidualsolidsisapproximatedbyasumovertwotypesoftur-
withafractalgeometry(Blum&Wurm2000),andabove∼ bulenteddies:thosewithspatialfrequencies k andturnover
1 m/s disruption of thecolliding agglomerates takes place. times tk > ts, where ts is the particle stopping time, and
While binary collisions of the smallest dust grains are which merely advect the particle; and those with tk < ts,
whichdecaybeforetheparticlehashadtimetocrossthem.
Theresultisavelocitydispersionthatdecreaseswithparti-
⋆ E-mail:[email protected] cle stopping time as (Ωts)−1/2, with a sharp steepening for
† E-mail:[email protected] particleswithΩts ∼afewtimes10−1.Underminimummass
‡ E-mail:[email protected] solar nebula conditions at 5 AU, this corresponds to parti-
2
cle radii of tens of centimetres. The method is nicely sum- 2 METHOD
marisedandappliedbyCuzzi&Hogan(2003)tochondrule-
2.1 Numerical Method
sized objects.
The accumulation of a large number of solids in small Weuseathree-dimensionalversionoftheZEUScode(Stone
regionsoftheturbulentprotoplanetarydisccanleadtofeed- & Norman 1992a;b) in which the shearing box model of
backeffectsonthegasflowifthesolidmassdensityisgreater HGB is implemented. To simulate the solid particles, we
thanthelocalgasdensity.Suchparticleagglomerations can haveadded a module that solves theirequation of motion,
also lead to dust growth if the associated collision rates
are conducive to sticking. One recently suggested clumping 1
mechanism is a streaming instability that arises due to the f =− (vp−vg)−2Ω×vg+3Ω2xxˆ (1)
ts
relative motion between solids and gas in a Keplerian disk
(Youdin&Goodman2005). Thefeedbackofsolidsontothe where f is the force per unit mass on the particle, the first
gascangenerateexponentialgrowthofparticledensityper- termontheright-handsideistheEpsteindraglawforspher-
turbations as well as turbulence, and particle clumps can icalparticles,tsistheparticlestoppingtime,vpistheparti-
develop overdensities up to a factor of ∼ 50 relative to the clevelocity,vg isthevelocityofthebackgroundflowat the
mean gas density (Youdin & Johansen 2007, Johansen & particle’s position, and Ω=(0,0,Ω) is theangular velocity
Youdin2007). of the box. The last two terms of Eq. (1) are the Coriolis
and tidal forces, respectively, that arise as a result of the
The process of turbulent concentration (Maxey 1987)
non-inertialcharacteroftheshearingbox(xisthedistance
hasbeenexploredbyCuzzietal.(2001),whousetheratioof
in the radial direction xˆ). The particle module takes as in-
localparticledensitytoitsglobalaverageasthemeasureof
put thegas velocities v calculated by ZEUS on thefaces of
amultifractal distribution of concentrated chondrules.This
the computational grid cells. These velocities are obtained
distribution resembles that of the dissipation of the turbu-
by solving the ideal MHD equations for the shearing box,
lentkineticenergyontheKolmogorovscale(Chhabra1989).
together with an isothermal equation of state,
The connection between turbulent concentration and tur-
bulent dissipation suggests that the local structure of the
chondrule concentration field may be independent of flow ∂ρg +∇·(ρgv)=0 (2)
properties.Forlargerparticles,anaccuratecharacterisation ∂t
of their numberstatistics in theturbulent circumstellar en-
vironment is still lacking.
∂v 1 B2 (B·∇)B
Therecognition thatturbulenceduetothemagnetoro- +v·∇v=− ∇ P + + −2Ω×v+3Ω2xxˆ
∂t ρg „ 8π« 4πρg
tational instability (MRI; Balbus & Hawley 1991; Hawley, (3)
Gammie & Balbus 1995, henceforth HGB) can provide the
viscousstressesnecessarytoallowinwardaccretionindiscs,
∂B
has prompted investigations of different dynamical aspects =∇×(v×B) (4)
∂t
of dust evolution in the presence of magnetohydrodynamic
(MHD)turbulence.Planetesimalformationbygravitational
instability was addressed by Johansen et al. (2006, hence-
forthJKH),whosenumericalstudiesconcludedthatgroups P =ρgc2s (5)
of at least 1000 particles in their simulations could become where ρg is the gas density, B is the magnetic field, P is
gravitationally bound, assuming specific values of disc col-
umndensities(150and900gcm−2)andaratioofdiscscale thegas pressure, and cs is theisothermal sound speed. The
module interpolates each velocity component to the posi-
height toradiusof 0.04. Fromang& Nelson (2005) examine tion of each particle inside the respective computational
the accumulation of individual solids due to vortical struc-
grid cell, using a trilinear interpolation algorithm (Press et
turesinaglobalprotoplanetarydiscmodel.Theyfollowthe
al. 1992). Equation (1) is solved via a second-order Runge-
radial migration of up to 3000 bodies, whose drag interac- Kutta method. The back reaction of the particles on the
tionwiththegascorrespondstoobjectsofapproximately1
gas, inter-particle interactions, and particle coupling to the
m in size. Those particles able to migrate inwards do so at
magnetic field are not included.
different rates, with reductions in semi-major axes ranging
Despitethelikelypresenceofdeadzonesinprotoplane-
from a factor of 1.6 in approximately 50 orbits to a factor
tary discsas aresult of non-idealMHD processes (Gammie
of 2.2 in 200 orbits. The particles that become trapped in
1996), in this first study we neglect the effect of resistivity
vortices maintain an approximate constant radial position
and assume that our local disc model represents an active
untilthe end of thesimulation.
regioninthesurfacelayersoratalargeorbitalradius(Sano
In the following we present a study of the kinematics et al. 2000).
of solid particles in alocal MHDmodel of aprotoplanetary WeadoptthevaluesusedinHGBfortheinitialdensity,
disc. In Section 2 we describe the numerical method em- initialpressure,angularvelocity,plasmaparameterandbox
ployedtomodeltheaccretionflowandtheparticles.Section dimensions:ρ0=1,P0 =10−6,cs =Ω=10−3,β=400and
3 contains results obtained from measurements of particle H×2πH×H(whereHisthediscscaleheight),respectively.
velocities and spatial distribution. In Section 4 we discuss The initial magnetic field is vertical but, unlike HGB, it
theirsignificance in thecontext of particle growth in a pro- variessinusoidallyintheradialdirection,insuchawaythat
toplanetary nebula, and concluding remarks are presented the net magnetic flux is zero. The grid resolution used is
in Section 5. 84×180×84 zones.
Kinematics of solid particles in a turbulent protoplanetary disc 3
Simpletestsoftheparticleintegratorinvolvedtheevo- Also shown as plus signs are the ratios S ≡ σp/σg of the
v v v
lutionofasinglesolidparticleinazero-velocitybackground velocity magnitude dispersions σvp,g = h(vp,g−hvp,gi)2i
flow, with various values of the stopping time ts which, in (notice that S2 6= S2+S2+S2). The dqotted curve corre-
theEpstein regime, is given by(Weidenschilling 1977). v x y z
sponds to the analytical solution of V¨olk et al. (1980) and
Cuzziet al.(1993) for S .Thesolid anddashed curvescor-
ts= aρp (6) respond, respectively, tovthe radial and azimuthal solutions
csρg of Youdin and Lithwick (2007). With no gas drag, the par-
where a is the particle radius and ρp theparticle solid den- ticles would follow epicyclic paths with the ratio σvpx/σvpy
sity.Givenaninitialpositionandvelocity,theresponsetime of the radial and azimuthal velocity dispersions equal to 2.
of the particle to the drag force was found to adjust ac- Larger ratios σp /σp of around 2.3 measured in the MHD
vx vy
curately to the expected analytical behaviour for a range calculationareduetoweakdragforcesfromtheturbulence.
of values of the dimensionless stopping time Ωts between Separate numerical integration of the particle equation of
2×10−3 and 100. motion (1) with a sinusoidally time-varyinggas velocity re-
produces the anisotropy, given variation periods of one to
twoorbits.Fourieranalysisofthegasvelocitiesattheparti-
2.2 Simulations
clelocationsintheMHDcalculation indicatesthestrongest
Weperformed twosetsofsimulationsinwhich theshearing modes do have periods longer than one orbit. Our neglect
box is first evolved for 10 orbits (1 orbit=2π/Ω) to ensure of the vertical component of gravity appears to have little
thatturbulenceduetotheMRIhasdeveloped.Atthistime, effect on the ratios between particle velocity components,
N particles are introduced in the flow with random initial assimilaranisotropywasobservedinstratifiedshearingbox
positionsandwithinitialvelocitiessettozero.Insimulation calculations(Carballidoetal.2006),inwhichfortheΩts=10
set A, N = 80,000, and the total number of particles is case, σvpx ∼2.5σvpz.
divided in eight groups of 10,000 members. Each group is Itis important toremark thattheanalytical modelsof
given a specific value of Ωts. The values used were 100, 30, Fig. 2 were originally derived assuming an isotropic turbu-
10, 5, 2, 0.2, 0.02 and 0.002. The system was evolved for a lentvelocityfield,soinprincipleoneshouldnotexpectthat
further100 orbits. they be an accurate description of particle motion in MRI
Simulation set B, which was used to study the par- turbulence. Nevertheless, the particle data obtained from
ticle spatial distribution, consists of four groups of N = thesimulationsconformreasonablywelltothemodels.This
4,704,000 particles, with each group having one of Ωts=2, will bediscussed below.
0.2,0.125and0.05.Thecodewasrunforafurther30orbits.
3 RESULTS
3.3 Pair-wise relative velocities
3.1 Viscous evolution of disc
Fig. 3 shows the magnitude v of the relative velocity of
rel
Figure 1 shows the evolution of the magnetic (dotted line) pairsofparticlesthatarelocated insideacubicsub-volume
and hydrodynamic (dashed line) volume-averaged stresses inthebox,insimulationsetA,asafunctionofinter-particle
associated with our disc model, duringthelongest run per- separation r (in units of the box size H). The velocity has
formed. The stresses are normalized to the initial pressure. been normalised to the gas sound speed. The sub-volume
Their sum (solid line) gives an estimate of the volume- has dimensions 1×1×1, and extends along the y direc-
averaged α parameter of classical disc theory. The average tion between the values y = 2.64 and y = 3.64. At each
valueof αis 8×10−3, which falls in therange 10−4−10−2 snapshot, a sample of 700 particles, all with thesame stop-
inferred for T Tauri stars (e.g. Hartmann et al. 1998). ping time, is randomly selected in the sub-volume, and the
relative speed for each of the 244,650 possible pairs in the
sample is calculated. The relative separation between pair
3.2 Velocity dispersions
members is binned, with the bin size equal to the x-length
Simulation set A has been used to calculate the time- of a grid cell. The values of vrel are then averaged in each
averaged distributions of the particle velocity components. bin overa duration of 7 orbits. Thedata shown correspond
Tcohrerersepsounltdininggdgisapsevrseilooncsit,yσvpdiis(pie=rsxio,nys,zσ)vg,in,ohramvaelibseedenbymtehae- (tdoitahmeosntdosp)p,inagndtim2×es10Ω−ts3=(1fi0l0led(trciiarnclgelse)s.),T1o0c(asqlcuualraetse),v0re.l2,
sured as functions of particle stopping time. To calculate the contribution of the shear has been subtracted from the
thisratio,S ≡σp/σg ,thegasvelocitieswerecomputedat y component of theparticle velocities.
i vi vi
the position of the particles, at each snapshot and for each It can be seen that at small separations the relative
valueofΩts,andtheeffectofshearwassubtractedfromthe speed decreases appreciably. This is to be expected, since
particle and gas y-velocities. The ratios S thus obtained grains that are equally coupled to the gas move coherently
i
werethenaveragedoverthelast10orbitsofthesimulation, with small turbulent eddies. In the limit of particle stop-
aperiodinwhichthelargestparticles(thosewithΩts=100) ping times less than the lifetime of the smallest eddies, the
havereached dynamical equilibrium with the gas. laminar flow on scales below the dissipation scale does not
The results are summarised in Fig. 2, where the solid- contribute to the random relative grain velocity. Formally,
to-gas ratios of the x (squares), y (triangles) and z (filled thiscanbeexpressedforgrainsofstoppingtimests1 andts2
circles) components of the velocity dispersions are plotted. as (Weidenschilling 1984)
4
Figure1.Timeevolutionofthestressesinthelocaldiscmodel.Thestressesarenormalizedtotheinitialgaspressure.Thesumofthe
magneticandhydrodynamicstresses(solidline)providesanestimateoftheviscousαparameter,whoseaveragevalueis8×10−3.
3.4 Spatial distribution
Inordertoimprovethestatisticsofparticlecounts,wehave
taken as the unit of volume a group of 3×3×3 grid cells,
u
vrel ∼ τ|ts1 −ts2| (7) referred to as a 3-macrocell. In this way, the average oc-
cupation number is 100 particles per 3-macrocell, with an
associated standard deviation of 10%.
Theplacementofthedustparticlesintheshearingbox,
where u and τ are the velocity and lifetime of the smallest attimet=10orbits,israndom,andthereforefollowsaPois-
eddy,respectively. son distribution. This may no longer be true at subsequent
AsimilarcalculationasthatofFig.3ispresentedinFig. times, when the turbulence has mixed the particles in such
4, but for pairs of particles with different stopping times. a way that the initial stochastic distribution has been lost.
Each symbol represents a different pair of stopping times This situation is apparent in Fig. 5, where the left panel
(the error bars have been ommited for clarity, but give an shows the cumulative distribution function, as a solid line,
uncertainty of ∼ 10% for the data represented by the plus of the particle number n per 3-macrocell at the moment
signs, the diamonds and the triangles, and of ∼5% for the when the particles are introduced in simulation set B, for
filledcircles).Therelativespeedsdonotdecreaseassteeply theΩts=2case.Thedottedlineisthecorrespondingcumu-
with decreasing separation as in the case of like particles. lative Poisson distribution with the same mean (hni=100)
Again,fromtheorder-of-magnitudeestimate(7),vrelshould asthendistribution. Thetwoarealike, indicatingthat the
acquireafinitevalueinthelimitofsmallstoppingtimes.For initial particle positions are indeed drawn from a random
the particles considered, this value is a significant fraction sample.However,therightpanel,whichcorrespondstot=30
of the sound speed (∼12−15%). orbits, makes it evident that the spatial distribution of the
It is worthwhile to point out that as the inter-particle particles deviates from a Poisson distribution. The vertical
separationdecreases,theeffectsofnumericalresolutionmay dashed line in each panel indicates the point at which the
become important for the purpose of calculating the rela- distance d between the two distributions is greatest. This
tive velocities. The size ℓ of the smallest turbulent eddies, distance constitutes a measure of the deviation of one data
to which the smallest particles are coupled, varies with vis- set from the other. In the particular case shown in Fig. 5,
cosity ν as ℓ ∼ ν3/4. The viscosity in turn can be affected at t=10 orbits the maximum distance is d=0.026, whereas
by thegrid resolution (HGB). at t=30 orbits, d=0.55. The situation is similar at all other
Kinematics of solid particles in a turbulent protoplanetary disc 5
Figure 2. Ratio of particle to gas velocity dispersion, as a function of particle stopping time. Data shown correspond to the radial
(squares), azimuthal (triangles) and vertical (circles) components of velocity, as well as to the velocity magnitude (plus signs). The
curvesrepresentanalytical modelsderivedforisotropicturbulence.
timest>0,andforthethreeotherhydrodynamiccouplings non-Poissonian characteroftheparticlespatialdistribution
of set B. is also evident in Fig. 8, as the distributions possess ap-
Tomeasuretheamountofsolidparticleclumpinginthe preciable “tails” extending in thedirection of increasing C.
shearing box, we calculate the number of particles in each This is normally associated with spatial correlations in a
macrocell and divideby its mean: system of particles in which their positions are not statis-
tically independent (Landau & Lifshitz 1980). In the case
ofparticles immersed in turbulence,thesecorrelations arise
n
C ≡ (8) fromtheintermittencyofturbulentmixing.Thecorrelations
hni
aremanifestedbyanincreaseinthevarianceoftheparticle
counts, with respect to a purely random distribution.
Fig. 6 shows distributions of C for each particle-gas
Positive correlations can induce clustering at some
coupling, at t=20 orbits. These distributions are typical
range of spatial scales. One way to quantify this clustering
throughout the duration of the simulation. Also shown for
is by means of thepair correlation function
comparisonisthedistributionofgasclumping,ρg/hρgi.No-
tice that thesolids exhibit more over-densities with respect
tothemean(hCi=1) thanunder-densities,showingacon-
siderably high number of macrocells that have values of C hn(r)n(r+∆r)i
η(∆r)= −1 (9)
above the average occupation number hni. The highest C hn(r)i2
attained is Cmax=26.5, by the particles with Ωts=0.125, at
t=27.1 orbits.AhistoryofCmax isplottedinFig. 7,where
the time axis has been shifted to the instant at which the where the particle number density n is evaluated at two
particles are introduced. points separated by the distance ∆r. The measurement of
Thetime-averaged C distributionsareshown inFig. 8. this quantity is common in the study of cloud droplet clus-
Note the difference in the shapes of the distributions com- tering in the presence of turbulence in the Earth’s atmo-
pared to those in Fig. 6, particularly for Ωts = 0.2, 0.125 sphere (Kostinski & Jameson 2000). In order to calculate η
and 0.05 at C = 0. This can be understood by consider- forthepopulation of solids inourshearingbox,weproceed
ing that the time-averaging procedure effectively assigns a as follows. Given the anisotropic character of the turbulent
non-zero number of particles to each 3-macrocell, which at eddies present in the box (HGB), η is measured indepen-
a specific snapshot (but not all) is likely to be empty. The dentlyalong each spatial direction i (=1,2,3, or x,y,z):
6
Figure 3. Relative speed (in units of the sound speed) as a function of inter-particle separation (in units of scaleheight), for pairs of
particleswiththestoppingtimesshownontheupperleft-handcorner.Thenumberofparticlesperstoppingtimeis700.Datahavebeen
binned with respect to r, with a bin size equal to the size of a grid cell in the x direction (≈0.01H). The values of vrel decrease with
smallvaluesoftheseparation.
ally andz givesthevolume-averaged paircorrelation func-
tionhη (∆x)i.Thefunctionshη (∆y)iandhη (∆z)iareob-
x y z
hn(x )n(x +∆x )i
η (∆x )= i i i −1 (10) tained in an analogous manner.
i i hn(x )i2
i The volume-averaged pair correlation functions for the
Furthermore, to obtain a sufficient number of points Ωts=0.05 particles, at the time of their introduction in the
flow(t=10.0orbits),areidenticallyzero,reflectingthesta-
along each direction, at which to evaluate n, we consider
tistical independence of the initial particle numbers. This
2-macrocells instead of the 3-macrocells used to calculate
the clumping factor C. To avoid the effect of the periodic also holds for the cases Ωts=0.125, 0.2 and 2. The mea-
sured hη i at the end of the simulation are shown in Fig. 9.
boundary conditions in the y and z directions, as well as i
The columns correspond to the x, y and z directions, and
the quasi-periodic conditions along the radial direction x,
therows represent, from top to bottom, the stopping times
the separation ∆x between 2-macrocells is constrained to
i
be,atmost,half ofthecorrespondingboxlength.Thus,for Ωts =2, 0.2, 0.125 and 0.05. The error bars correspond to
thedata dispersion at each separation ∆x .
example, to evaluate η at separations of one macrocell, a i
x
horizontalrowofgridcellsischosen,sayonewithendcellsat One feature is immediately noticeable from Fig. 9: the
x=−0.5 and x=0.5, and with fixed y and z coordinates. vertical correlation function (third column) is very nearly
The particle number density n is measured at the first 2- flat for the four values of Ωts. This would seem to indicate
macrocell in this row, and also at the next macrocell. The that particles remain almost evenly mixed in the vertical
product of these two values of n is stored, and the process direction, without deviating much from the initial random
is repeated starting at the second macrocell. If m is the distribution. The particles show concentrations over radial
number of 2-macrocells in the row, the process ends with distancesofapproximately0.1H (forΩts=2)and0.05H (for
theproductofnevaluatedat macrocells m/2−1andm/2. the other stopping times). In the azimuthal direction the
Theaverageofalltheproductsiscalculatedandinsertedin respective correlation lengths are ≈ H and ≈ 0.3H. This
theright-hand side of Eq. (10). range of distances represents between 5% and 16% of the
The row of grid cells used in the calculation just de- box size.
scribed is located at specific y and z, and it is only one The largest correlation lengths belong to the Ωts=2
of 84 ×180 = 15,120 similar horizontal rows. With this particles, which also show the lowest average concentration
number of values for ηx(∆x), an average over all rows at (Figs. 7 and 8). Even taking the highest value of Cmax for
Kinematics of solid particles in a turbulent protoplanetary disc 7
Figure 4.SameasFig.3,butforpairsofparticleswithdifferentstoppingtimes.
Figure 5. Cumulative distributions of the number of particles per 3-macrocell, for Ωts = 2. The left panel represents the initial
distributionatt=10 orbits.Thesolidlinecorresponds totheparticledistributionn,whilethe dotted lineisthecorrespondingPoisson
distribution with the same mean. The vertical dashed line indicates the value of n at which the vertical separation d between the two
distributionsisgreatest.TherightpanelshowsthattheparticledistributionisnotPoissonianatt=30orbits.Thisistrueingeneralfor
t>0andthethreeothervaluesofΩts.
this stopping time (≈21), the surface density reached in a stoodtobeagroupofparticlesthathasasizesmallerthan
clump of dimensions hη i×hη i is ≈ 2.1×104H−2, com- the length δ of the diagonal of a 3-macrocell. The size of
x y
pared to 2.65×105H−2 for the most highly concentrated thegroupisdeterminedbycalculatingaclump“diameter”,
particles, those with Ωts=0.125. basedon themaximumandminimumvaluesoftheparticle
coordinates.
In order to constrain thelifetime of typical particle as-
sociations from thesimulation data,itisnecessary torefine ForeachvalueΩts=0.2,0.125and0.05,theclumpwith
the working definition of a clump. Here a clump is under- thehighest particle density was followed duringsubsequent
8
Figure 6. Distribution of the clumping number C ≡ n/hni at t =20 orbits, for all particle stopping times. All solids exhibit an
overabundance of 3-macrocells for which n is larger than the mean. All stopping times are more clumped than the gas (dash-tripple
dotted line).
snapshots. In all cases, the particles become detached from turbulence(Ormel&Cuzzi2007).Thoseresultsvalidatethe
the original clump after ∆t = 0.1 orbits, in such a way useof Eq.(7) in aregime wherethestoppingtimes of both
that a calculation of the new diameter yields a value larger particles are smaller than the lifetime τ of the smallest ed-
than δ. The fact that ∆t is the interval between snapshots dies.Fromameasurementofthecorrelation timeofthegas
indicates that the lifetimes are less than what is possible velocities,wefindτ ∼0.1orbits.Therefore,Eq.(7)willnot
to probe using the available data. The situation is similar apply to all the stopping times shown in Fig. 4, since the
if by clumps it is meant groups with horizontal dimensions stopping times for the Ωts=10 and Ωts=100 particles are,
corresponding tothe particle correlation length. respectively, ∼1.6 and ∼16 orbits. Nevertheless, a separate
calculation shows that relative velocities between pairs of
thethreesmallest stoppingtimes (Ωts=0.002, 0.02 and 0.2,
whichcorrespondto∼3×10−4,3×10−3and3×10−2orbits)
4 DISCUSSION
are well described by Eq. (7) at short separations, with the
The measurement of solid particle velocity dispersions in lowest velocity (v ≈ 0.09c ) corresponding to pairs with
rel s
MRI turbulence allows to test the analytic model of V¨olk Ωts=0.002 and 0.02.
et al. (1980) and Cuzzi et al. (1993), which adopted a Kol- For ts > tL, with tL being the overturn time of the
mogorov energy spectrum. The turbulence that develops in largest eddies (typically taken as one orbital period), the
the shearing box setup that we use follows roughly a Kol- model of Ormel & Cuzzi (2007) gives
mogorov scaling (HGB), and the dispersion of the particle
velocitymagnitude isgenerally consistentwiththeprevious
1 1
resulTtsh.e radial and azimuthal velocity dispersions were vr2el =vg2as„1+Ωts1 + 1+Ωts2« (11)
found to conform to the analytical results of Youdin and where vgas is the velocity of the turbulent flow, and it is
Lithwick (2007), who use a turbulence model in which gas assumed that ts1 > ts2 whitout loss of generality. Direct
velocitycorrelationshvxvyiareweakcomparedtohvx2+vy2i. measurementfrom oursimulation datayields vgas ≈0.12cs,
ThisisindeedthecaseintheMRI-generatedturbulentflow and so for the two largest stopping times studied (Ωts=10
present in our shearing box. and Ωts=100) we obtain vrel ≈ 0.04cs, in close agreement
Recent analytical calculations have provided closed- with thenumerical results of Fig. 4 (filled circles).
form expressions for the relative velocity between solids of The two-point correlation function calculated in Sec-
different sizes, as a function of stopping time, in isotropic tion 3.4 gives an estimate of the length scale over which
Kinematics of solid particles in a turbulent protoplanetary disc 9
Figure7.TimehistoryofthemaximumclumpingfactorCmax fortheparticlestoppingtimesusedinrunB.Thetimemarkedast=0
orbitscorrespondstotheinstantatwhichparticlesareintroducedintheflow.
particle clustering occurs. In a MMSN scenario at 5 AU, Eventhoughthestoppingtimesusedinthiscalculation
in which the scale height H of the disc is ∼ 4×1012 cm, andthosestudiedbyJKHaredifferent,wecanstillcompare
theminimumcorrelationlengthoftheparticlenumberden- our results for Ωts=2 with their Ωts=1 case. The main dif-
sity would be ξ ∼ 0.05H ∼ 2×106 km. A recent model of ference resides in the particle velocity dispersion, which we
chondrule formation, in which molten chondrules come to measureas0.036cs,approximately50%largerthantheyob-
equilibriumwith thegas evaporatedfrom otherchondrules, tain.IfourmaximallyconcentratedΩts=2groupweretobe
indicatesthattheseparticlesmusthaveoriginatedinregions gravitationally unstable,thevelocity dispersion would need
larger than 6,000 km in radius, or ξ∼6×10−4H (Cuzzi & to be a factor ∼ 2.5 smaller. We also note that the radial
Alexander2006). The smallest particles in our study of the drift incorporated by JKH leads to particle concentrations
spatial distribution correspond roughly to solids with a ra- thatarestrongerthaninthecasespresentedhere(Johansen
dius of 20 cm, two orders of magnitude larger than typical et al. 2007), making gravitational collapse more likely.
chondrulesizes.Itwouldbenecessarytoperformadditional
calculationswithparticlestoppingtimesthatcorrespondto
this size range, in order to further constrain the value of ξ.
5 CONCLUSIONS
However,itisunlikelythatcurrentcomputationalresources
wouldallowforthenecessaryresolutiontoprobesuchsmall We have performed numerical MHD calculations of solid
correlation lengths. particlevelocitiesandspatialdistributioninaturbulentpro-
toplanetary disc, in the context of the shearing box model.
We have estimated the Jeans radius R of the maxi- The results suggest that, even though the turbulent flow
J
mally concentrated clump for each stopping time, following has an anisotropic spectrum, the dispersion of the particle
calculationsbyJKH.Usingsimilardiscparameters(column velocity magnitude (normalised by the corresponding dis-
density Σ = 900 g cm−2, average dust-to-gas density ratio persion for the gas) as a function of particle stopping time
ζ = 0.01, and scale height-to-radius ratio H/r = 0.04) and is consistent with analytical models that assume isotropy.
a turbulent diffusion coefficient δt =α = 8×10−3 inferred Nevertheless,thedispersionsoftheindividualvelocitycom-
from our shearing box model, the values of R are smaller ponents for each stopping time are different. In the case of
J
thanthetypicalclumpsizebyfactorsbetween∼9and∼250. particles with Ωts & 1, the radial component of the solid-
The clump radius that is closest to R corresponds to the to-gas velocity dispersion ratio is larger than the azimuthal
J
Ωts=2 particles. Under these assumptions, our dense par- and vertical components typically by factors of ∼ 2.3− ∼
ticle associations would not be subject to a gravitational 2.5. The relative speed between pairs of particles with the
collapse. samestopping timedecreases towards small valuesat small
10
Figure8.Time-averageddistributionoftheclumpingnumberC,forallparticlestoppingtimes.Thetailsofthedistributions,although
notparticularlylong,areindicativeofthenon-Poissoniannatureofthenumberdensitycounts.
separations,whereasforpairsofdifferentstoppingtimesthe REFERENCES
relativespeedreachesafinitevaluethatdependsonthesizes
Balbus S.A., Hawley J. F. 1991, ApJ, 376,214
of the pair members.
Blum J., Wurm G. 2000, Icarus, 143, 138
The clustering induced on theparticle spatial distribu-
Carballido A., Stone J. M., Pringle J. E. 2005, MNRAS,
tion by the turbulent flow gives rise to short-lived clumps.
358, 1055
Under conditions of a minimum-mass solar nebula, these
Carballido A., Fromang S., Papaloizou, J. 2006, MNRAS,
clumps are at least nine times larger than their respective
373, 1633
Jeans radii, as a result of high particle velocity dispersions.
Chhabra A. B., Meneveau C., Jensen V. R., Sreenivasan
Gravitational collapse would therefore seem to be an un-
K. R.1989, Phys.Rev.A, 40, 5284
likely mechanism in the formation of more massive objects,
Cuzzi J. N., Dobrovolskis A. R., Champney, J. M. 1993,
undertheconditions that we assumed.
Icarus, 106, 102
Future calculations should increase the numerical res-
CuzziJ.N.,DobrovolskisA.R.&HoganR.C.1996,Chon-
olution used in the measurement of the relative velocities
drulesandtheprotoplanetary disk,R.Hewins,R.H.Jones,
between particles, in order to improve fittings to existing
and E. R.D. Scott, eds., Cambridge University Press.
analytical models.
Cuzzi J. N., Hogan R. C., Paque J. M., Dobrovolskis R.
2001, ApJ,546, 496
Cuzzi J. N.& Hogan R.C. 2003, Icarus, 164, 127
Cuzzi J. N.& AlexanderC. M. O. 2006, Nature,441, 483
Gammie C. F. 1996, ApJ,457, 355
Hartmann L., Calvet N., Gullbring E., D’Alessio P. 1998,
ApJ, 495, 385
ACKNOWLEDGMENTS
HawleyJ.F.,GammieC.F.,BalbusS.A.1995,ApJ,440,
This work benefited from discussions with Jim Pringle, 742 (HGB)
Richard Nelson, and Gordon Ogilvie. Suggestions from an Johansen A., Klahr H. 2005, ApJ, 634, 1353
anonymous referee were most valuable. AC acknowledges JohansenA.,KlahrH.&HenningT.2006,ApJ,636,1121
supportfrom CONACYTscholarship167912 andaCaltech (JKH)
Postdoctoral Scholarship. Part of this work was carried out Johansen A., Youdin A.N. 2007, ApJ, 662, 627
at the Jet Propulsion Laboratory, which is operated by the Johansen A., Oishi J. S., Mac Low M., Klahr H., Henning
California InstituteofTechnologyundercontracttoNASA. T. & Youdin,A.2007, Nature,448, 1022
