The Validity of the Super-Particle Approximation during Planetesimal Formation

0.38 MB·English

by  Hanno Rein

Astronomy&Astrophysicsmanuscriptno.12870 c ESO2010
(cid:13)
January8,2010
The Validity of the Super-Particle Approximation during
Planetesimal Formation
HannoRein,GeoffroyLesur,andZoe¨ M.Leinhardt
UniversityofCambridge,DepartmentofAppliedMathematicsandTheoreticalPhysics,
CentreforMathematicalSciences,WilberforceRoad,CambridgeCB30WA,UK
e-mail:[email protected]
Submitted:11July2009-Revised:18December2009-Accepted:22December2009
0 ABSTRACT
1
0 The formation mechanism of planetesimals in protoplanetary discs is hotly debated. Currently, the favoured model involves the
2 accumulationofmeter-sizedobjectswithinaturbulentdisc,followedbyaphaseofgravitationalinstability.Atbest,onecansimulate
n a few million particles numerically as opposed to the several trillion meter-sized particles expected in a real protoplanetary disc.
a Therefore,singleparticlesareoftenusedassuper-particlestorepresent adistributionofmanysmallerparticles.Itisassumedthat
J small-scalephenomena donotplayaroleandparticlecollisionsarenotmodelled. Thesuper-particleapproximation isnotalways
validwhenappliedtoplanetesimalformationbecausethesystemcanbemarginallycollisional(oforderonecollisionperparticleper
8
orbit).Thesuper-particleapproximationcanonlybevalidinacollisionlessorstronglycollisionalsystem,although,inmanyrecent

numericalsimulationsthisisnotthecase.
]
P In this work, we present new results from numerical simulations of planetesimal formation via gravitational instability. A scaled
E systemisstudiedthatdoesnotrequiretheuseofsuper-particles.Thissystemissimplifiedforcomputationalpracticalityandproper
identificationofimportantprocesses:1)theevolutionofparticlesisstudiedinalocalshearingbox;2)theparticle-particleinteractions
.
h suchasgravity,physicalcollisions,andgasdragaresolveddirectlyassumingaconstantbackgroundshearflowwithoutanyfeedback
p fromtheparticles.Wefindthatthescaledparticlescanbeusedtomodeltheinitialphasesofclumpingifthepropertiesofthescaled
- particlesarechosensuchthatallimportanttimescalesinthesystemareequivalenttowhatisexpectedinarealprotoplanetarydisc.
o
Constraintsaregivenforthenumberofparticlesneededinordertoachievenumericalconvergence.
r
Wecomparethisnewmethodtothestandardsuper-particleapproach.Wefindthatthesuper-particleapproachproducesunreliable
t
s resultsthatdependonartifactssuchasthegravitationalsofteninginboththerequirementforgravitationalcollapseandtheresulting
a
clumpstatistics.Ourresultsshowthatshort-rangeinteractions(collisions)havetobemodelledproperly.
[

Keywords.planetformation–planetesimals–N-bodysimulations–convergence–accretiondisc
2
v
9
0 1. Overview (Weidenschilling 1977b). Even if km-sized planetesimals were
1 able to form, they would be in danger of being ground down
Extrasolar planets have been observedaround a variety of par-
0 againbymutualcollisions(Stewart&Leinhardt2009).
. entstars frompulsarsto solar-typestars to M-dwarfs(see e.g.,
1
Chauvinetal. 2004; Wolszczan&Frail 1992) indicating that Gavitational instability is often considered to be a solution
0
planetformationiscommonandsuccessfulinabroadrangeof to most of these problems because the intermediate sizes are
0
environments. However, the process of planet formation itself avoidedalltogether.Inthistheory,thedustlayerbecomesdense
1
is not directly observable, leaving theory and numerical simu- enough for the Keplerian shear and velocity dispersion of the
:
v lationstofillintheblanksbetweenobservationsofhotcircum- dust particles to be unable to support the dust against its own
Xi stellardiscsaroundyoungstarsandplanetsorbitingmaturestars. self gravity. The dust then collapses into clumps that eventu-
Oneofthemostimportantunansweredquestionsinthetheoryof ally cool via drag forces and mutual collisions into planetesi-
ar planetformationiswhatisthemechanismforplanetesimalfor- mals.Manydifferentauthorshaveworkedonthissubject.Until
mation,i.e.,theprocessbywhichthebuildingblocksofplanets recently, the focus has been on quiet, non-turbulent, and low
areformed. densityregionsofthe accretiondisc(seee.g.,Michikoshietal.
There are two main theories for planetesimal formation: 2009, 2007; Tangaetal. 2004). However, the turbulent gas in
mutual collisions (e.g., Hayashietal. 1977) and gravitational theprotoplanetarynebulastirsthedust,whichincreasestheve-
instability in the dust layer (Goldreich&Ward 1973). In the locity dispersion of the dust particles. Several ideas have been
first hypothesis, dust particles grow as the result of accretion- proposedtoovercometheturbulence-inducedmixingofthedust
dominated collisions. Although the formation of planetesimals particlesandcreatelocalizedclumps.Forexample,Cuzzietal.
by mutualcollisionsis consistentwith meteoritic evidence,the (2008) suggest that the same turbulence that stirs the dust on
collision speed between dust particles (or aggregates) must be largerscalesmayalsocollectthedustparticlesonsmallscales.
much slower than the typical velocity dispersion in a standard AsimilarideawasproposedbyJohansenetal.(2007),inwhich
protostellar disc to avoid destructive collisions (Blum&Wurm dust particles are localized into clumps using both turbulence
2008). In addition, the planetesimal formation process is so and the streaming instability (Youdin&Goodman 2005). The
slow that meter-sized particles are in danger of spiraling into clumpsthenbecomegravitationallyunstable.Athirdhypothesis
the star before growing large enoughto decouplefrom the gas suggests that large structures, such as vortices, may be able to
2 HannoReinetal.:Thevalidityofthesuper-particleapproximation
collectandprotectdustparticlesfromtheturbulentbackground 2. Ordersofmagnitude
(Barge&Sommeria1995).Inthispaper,wefocusonthegravi-
2.1.Definitionoftimescales
tationalcollapseinaverydenseandturbulentregionofthepro-
toplanetarydisc.Theoverdensityinthedustlayer,whichisap- The aim of this work is to simulate the dynamics of dust (or
proximately two orders of magnitude higher than the standard particles)interactinggravitationallyinsideaprotoplanetarydisc.
minimum-masssolarnebula,mighthavebeencreatedbyanyof Eachparticleissubjecttofivedifferentphysicalprocesses,each
theprocessesdescribedabove. operatingonvarioustypicaltimescales:
Numericalsimulationsmustbeusedtotestmodelsofplan- Stoppingtimeτ
s
etesimal formation. However, the separation of scales in the Each particle feels the effects of the surrounding gas
problemishuge.Ameter-sizedobjectismorethan12ordersof through a linear drag force. This force can be written as
magnitudesmallerthantheprotoplanetarydiscthicknessforcing F = mp(v v ), where τ is the stopping time of a
numericiststousesuper-particles.Thesesuper-particleshavea drag τs g − p s
particle of mass m , v is the velocity of the gas, and v is
massmuchhigherthanthemassofasingledustparticle,butthe p g p
thevelocityoftheparticle(Weidenschilling1977a).
forces(e.g.,dragforces)actingonthemareequivalenttothose
ofasingledustparticle.Tangaetal.(2004)wereabletousethis
Physicalcollisiontimescaleτ
approximationtosuccesfullymodelthegravitationalcollapsein Dust particles will suffer acphysical collision with another
aregimeinwhichcollisionsareunimportant.However,whenthe dust particle on a timescale of τ = (σ v¯ n) 1, where σ
c c p − c
surfacedensityisincreasedbytwoordersofmagnitude,asmen-
is the geometrical cross-section of the particles, v¯ is the
tionedabove,thesystemdoesbecomemarginallycollisionaland p
particle velocity dispersion, and n is the number density of
thereforethesuper-particleapproachbreaksdown.
particles.
Michikoshietal.(2009)performasimilarsetofsimulations,
alsoinalowdensity,non-turbulentregionofthedisc.Although Orbitaltimescaleτe
the authors describe their particles as super-particles, the way Thetimescaleassociatedwiththeparticlesorbitingthecen-
in which they accurately treat the collisions is the same as the tral object in a local shearing patch is the epicyclic period
method presented in this paper. Other approaches use Monte τe = 2πΩ−1, where Ω is the angular velocity at the semi-
Carlotypemethodstoresolvestatisticallytheoutcomeofover- majoraxisoftheshearingbox.
lappingparticles(seee.g.,Lithwick&Chiang2007).
Inthispaper,we considertwolimitsforthe gravitationalinter-
Inthefollowing,welookcarefullyatthenumericalrequire- actionsbetweendustparticles,long-rangeandshort-rangeinter-
ments of modelling gravitationalinstability accurately and test actions. The long-range interaction can be seen as a collective
the validity of using super-particles in a high density region. process involving interactions between clouds of particles. On
For simplicity, we begin with a system without turbulenceand the other hand, the short-range interaction is important for re-
assume that the surface density has already been increased by solvingclose approachesbetweenpairsofparticles(i.e.,gravi-
turbulence or some equivalent process. We find that the super- tationalscattering).Forthesakeofclarity,weseparatethesetwo
particle approach gives erratic results when collisions become processes.
important.
Gravitationalcollapsetimescaleτ
Gl
We thereforego on andpresenta methodthatdoes notuse Weintroducealengthscaleλ δr,whereδr istheaverage
≫
super-particles.Wesimulateascaledsysteminwhichthenum- distance between neighbouring dust particles. In this limit,
ber of particlesis dramaticallylower thanin a real protoplane- which we call the long-range interaction, the dust particle
tarydisc.Tokeepthesurfacedensitythesame,the massofin- distribution can be approximated by a continuous density
dividualparticleshasto be increased.The dragforceproduced ρ, and one can define a gravitational timescale, which is
by the surrounding gas is scaled accordingly, so that the stop- of the order of the free fall time of the system, defined by
pingtimeremainsconstant.Furthermore,weincreasethephys- τGl =1/√GρwhereGisthegravitationalconstant.
ical size of each particle. This allows us to keep the collision
timescales in simulations of different particle numbers exactly GravitationalscatteringtimescaleτGs
thesame.Theserequirementsleadtoresultsthatareindependent Short-rangegravitational interactions can be seen as an in-
ofthenumberofparticles,whichistheonlyfreeparameterinthe teractionbetweenapairofsingleparticlesthatcanresultin
simulation,andwecallthosesimulationsconverged.Theresults ascatteringevent.Aswithphysicalcollisions,thetimescale
agreewithothersimulationsperformedinthecontextofSaturn’s fora gravitationalinteractiondependson the cross-section,
rings(Goldreich&Tremaine1982;Daisaka&Ida1999).Inthis τGs =(σGv¯pn)−1,whereσGisthegravitationalcross-section
paper,wearguethatthismethodshouldalsobeusedincurrent of the particle. Using Kepler’s laws, one finds that σG
simulationsofplanetesimalformation. G2m2/v¯4.Thegravitationalcross-sectionisvelocitydepen∼-
p p
dent, which makes it qualitatively different from physical
The outline of this paper is as follows. In Sect. 2, we de-
(billardball)collisionsbetween two particles. We note that
fine the importanttimescales in the problemand show that the
gravitationalfocusing can also lower the physical collision
super-particleapproachbreaksdowninhighdensityregions.In
timescale.
Sect. 3, we discussthe numericalmethodandtheinitialcondi-
tionsofoursimulations.Resultsofnumericalsimulationsusing As mentioned earlier (Sect. 1), the particles might also be af-
boththesuper-particleapproachandourscalingmethodarepre- fected by excitation from a turbulent background state on a
sented in Sect. 4. We conclude the paper with resolution con- timescale defined by the turbelence itself. In this paper, we
straints for numerical simulations and discuss the implications model this excitation in a simplified way, described in Sect.
fortheplanetesimalformationprocessthroughgravitationalin- 3. In our case, the turbulence is scale independent and has no
stabilityinSect.5. timescaleassociatedwithit.
HannoReinetal.:Thevalidityofthesuper-particleapproximation 3
2.2.Therealphysicalsystem super-particles undergo gravitational scattering events, which
is unphysical since these events never occur in a real system.
Toidentifythedominantphysicalprocessesinaprotoplanetary
Individualparticle-particlecollisionsarenotmodelled.
disc, one has to quantify and compare the relevant timescales
There are various examples where this approach is used
asdefinedabove.Inthefollowingdiscussionandinthenumer-
successfully. For example, smoothed particle hydrodynamics
ical simulations, we assume R = 1 AU, a gas disc thickness
0 (SPH) uses the super-particle approximation to simulate many
of H/R = 0.01, and a minimum mass solar nebula (MMSN)
0 gasmolecules(Lucy1977).Thesesystemsareoftenassumedto
withsurfacegasdensityΣ = 890gcm 2.Oneusuallyassumes
− be strongly collisional to ensure thermodynamicalequilibrium
asolidtogasratioofρ /ρ =0.01.We,however,assumeasolid
s g inside each super-particle. One can therefore assign collective
to gas ratio of unity. As mentioned above, this value is justi-
propertiesto clouds of particles such as pressure and tempera-
fied by numericalsimulations(Johansenetal. 2007, 2009) that
ture. Another example is the evolution of galaxies. When two
demonstrate overdensities of the order of 101 102 can easily
− galaxiescollide,individualparticles(stars)willusuallynotun-
occur because of the interaction with a turbulent gas disc, the
dergo gravitational scattering events or physical collisions. In
streaming instability, or vertical settling. We note that signifi-
thatcase,thesuper-particleapproachmodelsacollisionlesssys-
cantly differentmodelsof the solar nebulahave been proposed
teminwhichcollectivedynamicsistheonlyimportantphysical
(Desch 2007). However, all our results can easily be scaled to
process.Inbothcases,thesuper-particleapproximationisavalid
differentscenarios.
approachfor simulating the system numerically,but will break
We simplify the system by assuming that all solid compo-
downassoonasthesystemismarginallycollisional.
nents of the disc are in meter-sized particles. We assume that
Asanexample,weconsidertwocloudsofparticlesundergo-
thesebouldershaveatypicalvelocitydispersioncausedbygas
inga“collision”.Inthestronglycollisionalcase,thecloudswill
turbulence v¯ 0.05c 30m/s, where c is the local sound
p ∼ s ∼ s slowly merge and the thermodynamicvariables (e.g., tempera-
speed, in accordance with numerical results (Johansenetal. ture) will diffuse between the clouds, the particles inside each
2007).Sincev¯ c ,theparticlestendtosedimenttowardthe
p ≪ s cloudfollowingarandomwalktrajectoryduetonumerouscol-
midplane,forming a finite thicknessdust layer due to the non-
lisions. On the other hand, in the collionless regime, the two
zerovelocitydispersion.
clouds simply do not see each other because there is no short-
Formeter-sizedboulders,thephysicalcross-sectionisσ
1m2,whereasthegravitationalscatteringcross-sectionisσc ∼ rangeinteractionpresentbetweentheparticles.Ina marginally
10 21m2 showing clearly that gravitational scattering is irGre∼l- collisional regime, some particles will collide with particles of
− theothercloud,leadingtoapartialthermalisationofthevelocity
evant to the dynamics of dust particles embedded in a disc.
distribution,butsomeotherparticleswillnothaveanycollision
However, all other timescales are roughly equivalent. Meter-
atallandwillfollowapproximatelyastraightline.Evidently,the
sized boulders are weakly coupled to the gas with a stopping
outcomeofthiseventcannotbedescribedusingasuper-particle
timeτ τ Ω 1 (Weidenschilling1977a).Thephysicalcol-
lisiontsim∼escea∼leis−τ =7.3Ω 1andthelong-rangegravitational approach.
c −
interactiontimescalesisτ Ω 1.
Gl −
∼
This physical system is expected to become gravitationally 2.4.Scalingmethod
unstable,accordingtotheToomrecriterion(Toomre1964).The
Theideaofourscalingmethodisthatoneshouldkeepallimpor-
instability occurswhen the gravitationalcollapse timescale τ
Gl tanttimescalesinanumericalsimulationasclosetothoseofthe
is shorter than the transit time due to random particle motions
realphysicalsystemaspossibleandmodelallparticlecollisions
λ/v¯ and the orbitaltimescale. Assuming thatthe particlescan
p explicitly.
bemodelledbyanisotropicgas1withasoundspeedv¯ ,thesys-
p The numerical system consists of N particles in a box
tembecomesunstablewhen num
with length H, simulating a patch of the disc at a fixed radius
v¯ Ω R . The particle mass is m = ΣH2/N (H is the box size
p 0 num num
Q < Q 1. (1)
≡ πGΣ crit ≃ andonescaleheight).Thus,thedensityintheboxandthelong-
rangegravitationalinteractionsare unchangedcomparedto the
Inthatcase,themostunstablewavelengthisgivenby realsystem.
The gravitational scattering cross-section of the numerical
2π2GΣ
λ = . (2) systemisthengivenby
T Ω2
G2Σ2H4
Inthefollowing,wecomparethephysicalparametersfromthis σ = . (3)
section to their counterpart in numerical simulations. We first G,num Nn2umv¯4p
summarize the super-particle approach before presenting the
Fortheinitialsurfacedensityandvelocitydispersionusedinthe
scalingmethod.
simulation(Sect.2.2),onefinds
H2 0.0001AU2
2.3.Super-particleapproximation σ . (4)
G,num ≃ N2 ≃ N2
Asuper-particlerepresentsmanysmallerparticles.Thedynam- num num
icsofthesmallparticlesarenotcalculatedexactly.Itisassumed The physical collision cross-section in the simulation is
thatallparticlesbehavesimilarlyandinacollectivemanner.To σ = πa2 where a , is the radius of the particles in
c,num num num
simulate gravitational interaction between super-particles, one thesimulation.Wederivetwoconstraintsfromthephysicaland
has to use a softened potential. Without this, the gravitational gravitationalcollisioncross-sections:
scatteringcross-sectionofsuper-particlesbecomestoolargeand
1. Samemeanfreepathinsimulationandrealphysicalsystems
1 Thiswouldbetrueforastronglycollisionalsystem,butisnotfor- That means Nnumσc,num = Nσc, where σc and N are the
mallyvalidinthestudiedregime. geometrical cross-section and the number of particles in a
4 HannoReinetal.:Thevalidityofthesuper-particleapproximation
regionofsizeH3inarealdisc,respectively.Inotherwords, The interactionterm F is divided into componentsrelated to
int
this condition ensures that the physical collision timescale self gravity, physical collisions between particles, or drag and
inthesimulationisexactlythesameasintherealsystem. excitationforces:
2. Negligiblegravitationalscattering Fint = Fgrav+Fcol+Fdrag+Fturb. (7)
Whentwoparticlesapproacheachother,theoutcomeshould
Wesolvetheresultingequationsofmotionwithaleapfrog(kick
beaphysicalcollision,whichmeansthatσ σ .
c,num ≫ G,num drift kick) time stepping scheme. In the following subsections,
The gravitational scattering timescale remains long com-
we describe the numerical methods used to compute each of
paredtothephysicalcollisiontimescale.
thesetermsandtheirphysicalrelevance.
Thefirstconditionplacesaconstraintontheparticlesizeinthe The box size of 0.01 AU was chosen such that there are
simulation.WithN = 4.7 1018andaparticleradiusofa =1m always several unstable modes in the box, when the system is
as found in a real disc w·ithin a box of volume H3, using the pushedintotheregimeofgravitationalcollapse,asestimatedby
parametersfromabovefindsthat Eq.2.
σ Nσ 0.014
a = c,num = c AU. (5)
num r π rNnumπ ≃ √Nnum ghostboxes
Wenotethatoncethisconditionissatisfiedintheinitialcondi-
y
tionsitwillbeautomaticallysatisfiedatalltimes.Insimulations,
wetypicallyfindthatthecollisiontimescaleisreducedbymore
thanoneorderofmagnitudeduringthegravitationalcollapse.
ghostboxes
Thesecondconditionisthensatisfiedbychangingthenum-
ber of particles. An interesting result is that if N > 10, the
num
second condition is easily satisfied if the first one is satisfied.
We note however that σ depends strongly on the velocity
G,num
dispersion. In particular, a velocity dispersion 5 times smaller
x
tostar
than the initial value (as found in some simulations) leads to
an increase by a factor of 600 in the gravitational scattering
cross-section.Therefore,we suggestusinga largesafety factor
(N >105)toensurethatthesecondconditionisalwayssatis-
num
fied,evenforsignificantlysmallerv¯ .Thisconditionalsoallows
p
ustoestimatewhenourapproachbreaksdown,namelywhenthe
numberofclumpsinthesystemissolow(N < 10 100)that
∼
gravitationalscatteringbecomesimportant.
Althoughitwouldbehelpfulasafurthersimplification,itis
not possible to perfom the above mentioned simulation in two
dimensions. In that case, the filling factor, which is defined to
be the ratio of the volume (or area) of all particles to the total Fig.1.Shearingbox,simulatingasmallpatchoftheprotoplan-
volume(orarea),isoftheorderofonefortheaboveparameters. etary disc. The grey box represents the shearing box of inter-
Wenotethatin2Danincreaseinparticlenumber(anddecrease est,thedashedboxessurroundingthecentralboxrepresentone
inparticlesizeasrequiredbyEq.5)doesnotdecreasethefilling ghostring.
factorifthecollisionallifetimeremainsconstant.
3. Methods
3.1.Selfgravity
Twodifferentkindsofsimulationareconsideredinthispaper:
In the N-body problem that we consider, we have to solve
1. Particlesareassumedtobepointmassesandhavenophysi-
Newton’sequationsof universalgravitationfor a large number
calsize,thegravitationalfieldoftheparticlesbeingapprox-
ofparticlesN.Thegravitationalforceonthei-thparticleisgiven
imated with a smoothing length to avoid numerical diver-
by
gences(seeSect.3.1).
2. Particleshaveaphysicalsizeand,therefore,nogravitational mm
F = G i j e , (8)
scmluodoedth.ing is required but physical collisions must be in- grav Xj,i rij+b 2 ij
(cid:16) (cid:17)
We referto theparticlesassuper-particlesandscaledparticles, where b is the smoothing length used to avoid divergences in
respectively. numericalsimulations(insimulationsthatincludephysicalcolli-
We performoursimulationsinacubicboxwithshearperi- sionsb=0)ande istheunitvectorinthedirectionofthegrav-
ij
odicboundaryconditionsintheradial(x)andperdiodicbound- itational force between the i th and j th particle. Calculating
ary conditions in the azimuthal (y) and vertical (z) directions, the gravity for each particle−from each−other particle results in
as illustrated in Fig.1 (Wisdom&Tremaine 1988). In the local O N2 operations.Toreducethenumberofoperations,onecan
approximation,the force per mass on each particle is a sum of use(cid:16) di(cid:17)fferent approximations to Eq. 8. In this paper, we use a
thecontributionsfromHill’sequationsandtheinteractionterms.
Barnes-Hut(BH)treecode(Barnes&Hut1986).
Hill’sequationcanbewrittenas
TheBHtreeinthreedimensionsdividestheoriginalboxinto
F = 2Ωe v +3Ω2 xe Ω2ze . (6) eightsmallercellswithhalfthelengthoftheoriginalbox.This
hill z p x z
− × −
HannoReinetal.:Thevalidityofthesuper-particleapproximation 5
processis continuedrecursivelyuntil there is only one particle casescenariobecausealltimescalesareequivalent,weincrease
percellleft.Thedepthofthetreeinahomogeneousmediumis theradiusbyafactorof4inallsimulations.
approximatelylog N.Onethencalculatesthetotalmassandthe
8
centreofmassforeverycellateverylevelofthetree.
3.3.Dragforce
To calculate the force acting on a particle, one starts at the
topof the tree anddescendsinto the tree asfar as necessaryto Eachparticlefeelsadragforce.Thebackgroundvelocityofthe
achieveagivenaccuracy.Ifthecurrentcellisfarawayfromthe gasv isassumedtobeasteadyKeplerianprofile
g
particleforwhichtheforceiscalculated,thenthedetailedden-
sitystructurewithinthiscellisnotimportant.Allthatmattersis 3
v = Ωxe . (9)
the box’s monopole moment (total mass and centre of mass). g −2 y
One therefore does not have to descend into this tree branch
Althoughaccretionflowsareturbulent,wechoosetouseasim-
anyfurther.TheBH tree reducesthe numberof calculationsto
plevelocityprofilesothatthebehaviourofself-gravitatingpar-
O NlogN .
ticles can be understood in conditions that can be easily con-
(cid:0)We use(cid:1)8 rings of ghost boxes in the radial and azimuthal
trolled.
direction(Fig.1showsonering).Aghostboxissimplyashifted
copyof the main box.The gravity on each particle is then cal-
culated by summing over contributionsfrom each (ghost) box. 3.4.Randomexcitation
Thissetupapproximatesa mediumofinfinite horizontalextent
and avoidslargeforce discontinuitiesat the boundaries.We do The system of particles described above is dynamically unsta-
not need any ghost boxes in the vertical direction because the bleevenwithoutself-gravity,asthecoefficientofrestitutionand
discisstratified. thestoppingtimedonotdependontheparticlevelocities2.The
particles can either settle down in the midplane and create a
Afinitenumberofghostringscanonlyactasanapproxima-
razor-thindiscif thestoppingtime istooshort,ortheycan ex-
tionofaninfinitemediumandthegravitationalforcewilltendto
pandverticallyforeverifthestoppingtimeisabovesomecritical
concentrateparticleshorizontallyin the centre of the box.Due
value and the excitation mechanism is provided only by colli-
tothe8ghostrings,theasymmetryofthegravitationalforcebe-
sions. Dustparticlesinaccretiondiscsare foundin the settling
tweenthecentreandthefacesoftheboxisreducedbyafactor
regime.However,a completesettling neveroccursin accretion
of20comparedtousingnoghostringsatall.Wenote,however,
disksbecausethebackgroundflowisalwaysturbulentduetothe
that a small asymmetry remains in our simulations and leads
Kelvin-Helmoltzinstability (Johansenetal. 2006), the MRI, or
to higherconcentrationsin theboxcentre.Inothersimulations
such as Tangaetal. (2004), this effect is not seen because the other hydrodynamic instabilities (see e.g., Lesur&Papaloizou
2009)whichdiffuseparticlesvertically(Fromang&Papaloizou
system is initially gravitationally unstable and integrated only
2006). This stirring process of turbulence is not present in the
for approximatelyone orbit. The system that we are interested
gasvelocityfield ofthesimulationspresentedhere(see Eq.9).
inismarginalgravitationallyunstableandthisslightasymmetry
To approximate the turbulent mixing, we added a random ex-
will become importantafter several orbits. In the beginningof
citation (white noise in space and time) to the particles in the
oursimulations,weintegrateformanyorbitsinordertoreacha
simulation. This allows us to have a well defined equilibrium
stable equilibrium.We thenpush the system intoa gravitation-
in whichthesystem isstable ratherthanstartingfromunstable
allyunstableregime(seeSect.3.5fordetails).Duringthestable
initialconditionsthatmightinfluencethefinalstate.
phase,horizontalover-concentrationsaretobeexpectedandare
indeedobservedbecauseofthisasymmetry.However,sinceina We perturb the velocity components of each particle after
realprotoplanetarydiscthegravitationalinstabilityisconsidered each timestep on a scale ∆vi = √δtξ, where δt is the current
toappearlocally(e.g.,insideavortexorasimilarstructure),this timestep and ξ is a random variable with a normaldistribution
effectofhavingapreferredlocationforgravitationalinstability around0andvariances.Thisexcitationmechanismheatsupthe
is not unphysicaland does not affect any conclusions made in particlesandallowsustohaveawelldefinedthreedimensional
thispaper.Wetestedthisbyapplyingalinearcut-offtothegrav- equilibrium as shown in Sect. 4. In our simulations, we use a
itationalforceatadistanceofoneboxlength(A.Toomre,private value of s = 1.3 1015 m2s−3 and s = 8.9 1014 m2s−3 for the
· ·
communication).Inthiscase,therewasnopreferredlocationin simulationswithoutandwithcollisions,respectively.Theseval-
theboxandthesimulationevolvedinexactlythesameway. ueswerechosensuchthattheequilibriumstateisapproximately
equivalentforbothtypesofsimulations.
3.2.Physicalcollisions
3.5.Initialconditions
Physical collisions between particles are treated in the follow-
ing way. After each timestep, we check whether any two par- All particles, which have equal mass, smoothing length, and
physical size, are placed randomly inside the box in the x
ticles are overlapping.Using the already existing tree structure
−
y plane. We note that particles have either a physical size or a
from the gravity calculation, one can again reduce the compu-
smoothing length associated with them, dependingon the type
tationalcostsfromO N2 to O NlogN . Inthese simulations,
ofsimulationweperform(Sect.3).Inthez-direction,theparti-
a collision between p(cid:16)artic(cid:17)les is(cid:0)defined(cid:1)as an overlap between clesareplacedinalayerwithaninitiallyGaussiandistribution
twoparticlesthatare approachingeachother.Whenacollision
aboutthemidplaneandastandarddeviationof0.05H.Weallow
isdetected,itisresolvedassumingenergyandmomentumcon- thesystemtoreachequilibriumbyintegratingituntilt=30Ω 1
servation(perfectlyelasticcollisions).Wenotethatthetimestep −
(4.8yrs).
hastobesmallenoughsuchthatnocollisionismissedandthe
overlapisalwayssmallcomparedtotheparticlesize.
2 Thisinstabilityis qualitatively similar to the instabilitydescribed
TheparticleradiusisgivenbyEq.5.Inordertokeepthecol- byGoldreich&Tremaine(1978)forSaturnrings,althoughinthelatter
lisiontimescaleτcclosetounity,whichisnumericallytheworst casethedragforcesareabsent.
6 HannoReinetal.:Thevalidityofthesuper-particleapproximation
lengthactsasaneffectivesmoothinglength.Ingeneral,tocheck
N=80000
25 m/s N=160000 thenumericalresolution,theparticlenumber Nnum isincreased
N=320000
N=160000 LS and the smoothing length b is reduced independently.A simu-
lation is resolved when the result is independent of both N
num
and b. This turns out to be impossible in the present situation.
The main reason is that the smallest scale in a gravitational
n 20 m/s collapse will ultimately depend on the smoothing length. By
o
ersi varying both parameters at the same time, an empirical scal-
disp ingof b 1/√Nnum worksfairlywellif Nnum is largeenough.
ocity However,∼this procedureis not justified and is unphysical.The
vel smoothing length was introduced to avoid divergentterms and
15 m/s not to model any small-scale physical process. Therefore, it
shouldnothaveanyimpactonthephysicaloutcomeofthesim-
ulation.Incidentally,theexistenceofthisdependencyindicates
thatshort-rangeinteractionsareimportantfortheresultofthese
simulationsandshouldbemodelledwithcare.
10 m/s
30 Ω-1 35 Ω-1 40 Ω-1 45 Ω-1 50 Ω-1
InFig.2,weplotthevelocitydispersionasafunctionoftime.
time
Thesimulationsbeginfromthe stable equilibriumdescribedin
Fig.2.Super-particles:Velocitydispersionasafunctionoftime Sect. 3. After t = 30 Ω 1, we switch off the excitation mecha-
infourdifferentruns.ThesimulationlabeledLShasatentimes −
nism.Becausethesystemcontinuestocool,itbecomesgravita-
largersmoothinglength.Thecurvesdonotoverlapbecausethe
tionally unstable within one orbit, as estimated by Eq. 1. Once
simulationsarenotconverged.
clumpingoccurs,thevelocitydispersionbeginstoriseagain.All
simulationsusetheparticleradiusgivenbyEq.5asasmoothing
35 m/s N=160000 lengthexcepttheoneslabeledLS,whichuse a tentimeslarger
N=320000
N=640000 value.The largersofteninglengthis approximatelya 128-thof
theboxlengthandillustratesthekindofevolutionexpectedfrom
30 m/s anFFTmethodusinga1283grid.
Snapshotsoftheparticledistributionoftwosimulationsare
n
persio 25 m/s pshaorawmneitnerFs,ige.x4c.eBptotthhesismmuoloatthioinngsluesneg1th6,0a0re00thpeasratimclee.sTahnedtaolpl
dis
velocity 20 m/s rwohweriesatshethseimbuoltatotimonrwowithuasessmaoostmhionogthliennggtlhenggivthenthbayt Eisq.te5n,
times larger. The simulations correspond to the blue (medium
dashed)andred(solid)curveinFig.2.Weshowthesnapshotsto
illustratetheimportanceofthesmoothinglengthinsimulations
15 m/s
withoutphysicalcollisions.Onecanseethatthesimulationsdif-
fer already before clumps form. Stripy structures appear on a
scale given by Eq. 2. As soon as we enter the unstable regime
10 m/s
30 Ω-1 35 Ω-1 40 Ω-1 45 Ω-1 50 Ω-1 at t 32Ω 1 (i.e., the velocity dispersion begins to rise, see
−
time also F∼ig.2) the simulationsevolveverydifferently.Thesimula-
Fig.3.Scaled-particles:Velocitydispersionasafunctionoftime tiononthetoprowofFig.4formsmanyclumpsatanearlytime,
inthreedifferentrunsthatincludephysicalcollisions.Allcurves whereasthebottomrowsimulationformsonlyafew,moremas-
overlapbecausethesimulationsareconverged. siveclumpsatlatertimes.
Thiscanbeconfirmedbylookingatspectraofthesametwo
Wetestedvariouswaysofpushingthesystemintotheunsta- simulationsasshowninFig.6.Thesespectraweregeneratedby
bleregime,eitherbyreducingtheturbulencestirringorbyshort- mappingthe particlesonto a 128 128 grid in the x y plane
eningthestoppingtime,butdidnotseeanyqualitativedifference andcomputingthefastFouriertra×nsformofthemapp−inginthe
aslongaswestartfromawelldefinedequilibriumstate.Inthe x direction. The resulting spectra were finally averaged in the
simulationspresentedinthispaper,weswitchofftheturbulence y direction (see Tangaetal. 2004, for a complete description
stirring.Followingthismodification,thesystembecomesgravi- of the procedure). As suggested by the snapshots, the nonlin-
tationallyunstableandboundclumpsformwithinafeworbits. ear dynamics of the gravitational instability strongly depends
on the smoothing length used. In particular, one observes that
smaller scales are amplified more slowly for larger smoothing
4. Results lengths. The resulting spectra at t = 40Ω 1 also differ signif-
−
icantly. With a small enough smoothing length, the spectrum
4.1.Superparticles
looks almost flat, whereas a large smoothing length introduces
We first presentsimulationswithoutphysicalcollisions(super- acutoffatk/2π 10.Wealsonotethatthesmoothinglengthin
≃
particles) which rely on the smoothing length b to avoid any thelattercaseisoftheorderofthegridsize(k/2π 100).The
divergencies. Although we use a tree code and therefore a cutoff observedclearly demonstratesthat the smoo∼thinglength
smoothed potential for the force calculation is assumed, the modifiesthe dynamicson scales up to 10 times larger than the
results are equivalent to an FFT-based method where the grid smoothinglengthitself.
HannoReinetal.:Thevalidityofthesuper-particleapproximation 7
Fig.4.Super-particles:Snapshotsoftheparticledistributioninthe xyplane.Bothsimulationsuse160000particleswithdifferent
smoothinglengths.Thesimulationonthebottomusesatentimeslargersmoothinglengththantheoneonthetop.Thesnapshots
weretaken(fromlefttoright)att = 0,30,37,40Ω 1.Withalargesmoothinglength,theoutcomelooksverydifferent,thesystem
−
ismorestable,morestripystructurecanbeseenandclumpsformlater,ifatall.
Fig.5. Scaled particles: Snapshots of the particle distribution in the xy plane. The simulations (from top to bottom) use 40 000,
320000,640000particleswiththeirphysicalsizegivenbyEq.5.Inallsimulations,thecollisiontimescalewaskeptconstant.The
snapshotsweretaken(fromlefttoright)att = 0,30,37,40Ω 1.Thesimulationwith40000particleshasalargefillingfactorby
−
thelastframe,whichpreventsclumpsfromforming.Theintermediateresolutionsimulation(middlerow)andthehighestresolution
simulation(bottomrow)havemoreandsmallerparticlesandthusasmallerfillingfactor.Theresults(i.e.,numberofclumpsinthe
lastframe)areverysimilarintheintermediateandhighresolutionsimulations.
8 HannoReinetal.:Thevalidityofthesuper-particleapproximation
4.2.Scaledparticles
The velocity dispersion evolution of simulations with N =
0.01
160000,320000,640000particlesincludingphysicalcollisions
is presentedin Fig.3. In these runs, no gravitationalsmoothing N=160000
t=38
lengthisneeded.Again,thesimulationsstartfromastableequi-
libriumandaftert=30Ω−1,wechangethesameparametersas 0.001
in the non-collisionalcase to push the system into the gravita-
tionallyunstableregime.Snapshotsoftheparticledistributions er Nt==31660000 LS Nt==31660000
w
forthreerunsarealso plottedin Fig.5. The middleandbottom Po
rowsofsnapshotscorrespondtothepurple(mediumdashed)and 0.0001
darkblue(longdashed)linesofFig.3.
N=160000 LS
One can see in Fig.3 as well as Fig.5 that the intermedi- t=38
ate and high resolution runs produce very similar results. We
N=160000 LS
call those simulations converged,as a change in particle num- 1e-05 N=160000
t=30
ber does not change the outcome. An additional, lower resolu-
tion run (N = 40000)is notconvergedbecausethe filling fac-
toris toolargein denseregions,preventingclumpsfrombeing 1 10 100
gravitationally bound.The problem does not occur in the non- Wavelength 1/L
collisionalcasesbecauseparticlesareallowedtooverlap.Since Fig.6. Super-particles: Density distribution spectrum with
the filling factor scales with particle number as 1/√N , this 160000particlesandtwodifferentsmoothinglengthsthatdiffer
num
issue is resolved for runs with more than a few hundred thou- bya factor10attimest = 30(red,solidcurve),t = 36(green,
sandparticles.Wenotethatthevelocitydispersionmightdiffer longdashedcurve),andt =38(blue,shortdashedcurved).The
slightly at later stages for the convergedruns because the final spectraatsimilartimesarenotontopofeachotherbecausethe
clump radius is still determined by the physical particle radius simulationsarenotconverged.
andthereforebytheparticlenumber.
WeshowthedensitydistributionspectruminFig.7asafunc-
tionofthenumberofparticles.Theresultingspectradonotde-
pend on the number of particles in the system, as expected. In 00..0011
comparisonwithFig.6,thesimulationswithcollisionsarecon-
vergedsincethedynamicisnotcontrolledbythenumericalpa-
rameterN,whichistheonlyfreeparameterinthesystem.
At very late times (t & 38Ω 1), the spectra begin to show 00..000011 tt==3388
−
asystematic trendtowardsmorestructureonsmallerscalesfor
erer
runswithlargerN.Thisisexpectedandduetothefillingfactor, PowPow tt==3366
already mentioned above. As soon as one clump in the simu-
00..00000011
lation is onlydeterminedbythe size ofthe individualparticles
it consists outof, our approachbreaksdown. One can also see
these clumps forming by looking at the velocity dispersion in tt==3300
NN==116600000000
Fig.3. The particles inside a clump begin to dominate the ve- 11ee--0055
locity dispersionoverthe backgroundaftert 38Ω 1. A clear NN==332200000000
−
indication of this is the spiky structure with∼a typical correla- NN==664400000000
tiontimeof Ω 1correspondingtodifferentclumpsinteracting
∼ − 11 1100 110000
andmergingwitheachother.Onewayaroundthisissue,which WWaavveelleennggtthh 11//LL
willbeconsideredinfuturework,istoallowparticlestomerge
Fig.7. Scaled particles: Density distribution spectrum with
(Michikoshietal. 2009). Using this accretion model, the mass
160000,320000,and640000particles(colorsareequivalentto
oftheclumpcanbeusedasanupperlimit.
Fig.6).Ateachtime,thespectraareconverged,i.e.,thespectra
are on top of each other and independent of the particle num-
ber, besides the noise level on very small scales. At late times
5. Discussionandimplications
(t & 38,bluecurve),onebeginsto see morestructureonsmall
In this paper,we have shown thatconvergencein N-bodysim- scalesinsimulationswithlargernumberofparticles(seetext).
ulations of planetesimal formation via gravitational instability
can be achieved when taking into accountall relevantphysical
processes.Itisabsolutelyvitaltosimulategravity,damping,ex-
citation, and physical collisions simultaneously, as the related collapse.Simulationswithandwithoutphysicalcollisionswere
timescales are of the same order and therefore all effects are studiedforarangeofparticlenumberstotestconvergence.
stronglycoupled. Incaseswithoutphysicalcollisions,convergencecannotbe
Asetofsimulationsisdefinedtobeconvergedwhenthere- achieved.However,itwaspossibletochangemultiplefreesim-
sultsdonotdependontheparticlenumberoranyotherartificial ulationparametersatthesametime,namelytheparticlenumber
numericalparametersuchasasmoothinglength.Asatestcase,a N andthesmoothinglengthb,suchthatinspecialcasesthere-
boxwithshearperiodicboundaryconditionswasused,contain- sults did not depend on the particle number. We note however
ing hundreds of thousands of self-gravitating (super-)particles thatthereisafreeparameterinthesimulation(i.e.,thesmooth-
ina stratifiedequilibriumstate.Theparticleswerethenpushed inglength)thateffectstheoutcomeandmakesit impossibleto
intoaselfgravitatingregimethateventuallyledtogravitational findtherealphysicalsolution.
HannoReinetal.:Thevalidityofthesuper-particleapproximation 9
Intheprotoplanetarydisc,physicalcollisionsdominateover lessthanunityatt=0isgivenby
gravitational scattering. In the case of planetesimal formation,
4
thesystemismarginallycollisional.Physicalcollisionswillpar- π(a )3 N C L2H, (11)
num num f
tially randomise the particle distribution, whereas a smoothing 3 · ≤
lengthdoesnotbecausethesofteninglengthisverylargecom- whereC isasafetyfactor.AlthoughC . 1wouldbeenough
f f
paredtothegravitationalcross-section.Evenifthegravitational toresolvecollisionsinitially,itisinsufficienttosimulatethecol-
particle-particle scattering becomes important, it can be seen lapsephase.Duringthecollapse,thefillingfactorrisesrapidly.
from Eq. 3 that the velocity dependence of the cross-section Assuming that one wishes to resolve collisions correctly when
σG differs fundamentally from the velocity independentcross- theparticleshavecontractedbyoneorderofmagnitude,onehas
sectionσC.We thereforearguethatthesuper-particleapproach to include a safety factor of Cf = 10−3. If the collapse occurs
shouldnotbeusedwhencollisionsbecomeimportant. mainlyintwodimensions,asinthesimulationspresentedhere,a
With collisions, the simulation outcome is both quantita- factorofC =10 2issufficient.Furthermore,onehastoensure
tivelyand qualitativelyverydifferentfromsimulationswith no thattheboxf conta−insatleastafewunstablemodes(seeEq.2).
physicalcollisions.Theinitialsizeoftheclumpsislargerandthe Allthistogetherplacestightnumericalconstraintsonnumerical
numberof clumpsis smaller. As there is no smoothinglength, simulationsofplanetesimalformationviagravitationalinstabil-
thereiseffectivelyonefewerfreeparameter.Changingthepar- ity.
ticlenumberwhilekeepingthecollisiontimeconstantdoesnot In future work, we will expand the discussion to include
change the outcome. We therefore call these simulations con- more physics, namely a proper treatment and feedback of the
verged. background gas turbulence. This will allow us to further clar-
Additionaltestsincludinginelasticcollisionswithanormal ify the behaviourof planetesimals in the early stages of planet
coefficientofrestitutionof0.25havebeenperformedtoconfirm formationandeventuallytestthedifferentformationscenarios.
thattheresultsdonotdependonthewayphysicalcollisionsare
modelled.Asexpected,thequalitativeoutcomeintermsofnum- Acknowledgements. We thank J. Papaloizou, G. Ogilvie, S.-I. Inutsuka, A.
ber of clumps and clump size is different because the physical Johansen, Y. Lithwick, and A. Youdin for helpful discussions and com-
ments. Hanno Rein is supported by an Isaac Newton Studentship and St
propertiesofthesystemhavebeenchanged.However,oncein- John’sCollegeCambridge.GeoffroyLesuracknowledges supportfromSTFC.
dividual collisions are resolved the results converge in exactly Zoe¨ M. Leinhardt is an STFC postdoctoral fellow. Simulations were per-
thesamemannerasinthosesimulationspresentedabove(which formed on the Astrophysical Fluids computational facilities and on Darwin,
haveacoefficientofrestitutionof1). the Cambridge University HPC facility. The authors are grateful to the Isaac
NewtonInstituteforMathematicalSciencesinCambridgewherethefinalstages
Wehavealsoexploredothergravitysolversbycomparinga
of this work were carried out during the Dynamics of Discs and Planets re-
fastFourier(FFT)gravitysolverwiththeBHtreecode(usinga
search programme. A simplified version of the simulation used in this pa-
smoothinglength)andconfirmedourpreviousresults.Thegrid per can be downloaded free of charge as an application for the iPhone at
in the FFT code acts as a smoothing length and a direct com- http://itunes.com/apps/gravitytree.
parisonbetweenthetwogravitysolversgaveanapproximateef-
fectivesmoothinglengthofaquarterofthegridlength.Because
References
onehastosolveindividualparticle-particleinteractions,thegrid
sizehastobesmallerthanthephysicalsizeoftheparticles.As Barge,P.&Sommeria,J.1995,A&A,295,L1
a result, the FFT method is always slower than a tree code. In Barnes,J.&Hut,P.1986,Nature,324,446
addition,thetreestructurecanbereusedtosearchefficientlyfor Blum,J.&Wurm,G.2008,ARA&A,46,21
collisions. Chauvin, G., Lagrange, A.-M., Dumas, C., et al. 2004, A&A,
This paper focused on the numericalrequirementsto study 425,L29
gravitational instability and the formation of planetesimals in Cuzzi,J.N.,Hogan,R.C.,&Shariff,K.2008,ApJ,687,1432
protoplanetary discs. The initial conditions were chosen such Daisaka,H.&Ida,S.1999,Earth,Planets,andSpace,51,1195
thattheequilibriumisawelldefinedstartingpointforthecon- Desch,S.J.2007,ApJ,671,878
vergence study. The collision, damping, collapse and orbital Fromang,S.&Papaloizou,J.2006,A&A,452,751
timescales are all of the same order, which is effectively the Goldreich,P.&Tremaine,S.1982,ARA&A,20,249
worstcasescenario.Toprovideanyconstraintsonplanetesimal Goldreich,P.&Tremaine,S.D.1978,Icarus,34,227
formationitself,onewouldhavetoproperlysimulatetheturbu- Goldreich,P.&Ward,W.R.1973,ApJ,183,1051
lentbackgroundgasdynamicswhichisbeyondthescopeofthis Hayashi,C.,Nakazawa,K.,&Adachi,I.1977,PASJ,29,163
paper. Johansen,A.,Henning,T.,&Klahr,H.2006,ApJ,643,1219
However, we were able to determine the numerical resolu- Johansen,A., Oishi, J. S.,Low, M.-M.M.,et al. 2007,Nature,
tion needed to resolve dust particles in a protoplanetary disc. 448,1022
Assuming that one wishes to simulate dust particles realisti- Johansen,A.,Youdin,A.,&MacLow,M.2009,ApJ,704,L75
cally, including collisions, and that the particles are uniformly Lesur,G.&Papaloizou,J.C.B.2009,ArXive-prints
distributed in a box of base L2 and height H, the size of each Lithwick,Y.&Chiang,E.2007,ApJ,656,524
particle is determinedby the collision rate in the real disc (see Lucy,L.B.1977,AJ,82,1013
Eq.5) Michikoshi, S., Inutsuka,S.-i., Kokubo,E., & Furuya, I. 2007,
ApJ,657,521
N Michikoshi, S., Kokubo, E., & Inutsuka, S.-i. 2009, ApJ, 703,
a = a, (10)
num rN 1363
num
Stewart,S.T.&Leinhardt,Z.M.2009,ApJ,691,L133
whereaandNarethesizeandnumberofparticlesinthevolume Tanga,P.,Weidenschilling,S.J.,Michel,P.,&Richardson,D.C.
L2H oftherealdisc. Thenumberofparticlesin the simulation 2004,A&A,427,1105
N has to be large enough to ensure that the filling factor is Toomre,A.1964,ApJ,139,1217
num
lowuntilclumpsform.Therequirementthatthefillingfactoris Weidenschilling,S.J.1977a,MNRAS,180,57
10 HannoReinetal.:Thevalidityofthesuper-particleapproximation
Weidenschilling,S.J.1977b,Ap&SS,51,153
Wisdom,J.&Tremaine,S.1988,AJ,95,925
Wolszczan,A.&Frail,D.A.1992,Nature,355,145
Youdin,A.N.&Goodman,J.2005,ApJ,620,459