Eccentricity growth of planetesimals in a self-gravitating protoplanetary disc

0.26 MB·English

by  M. Britsch

Mon.Not.R.Astron.Soc.000,000–000(0000) Printed2February2008 (MNLATEXstylefilev2.2)
Eccentricity growth of planetesimals in a self-gravitating
protoplanetary disc
M. Britsch1,2,4, C. J. Clarke1 and G. Lodato3,1
8
0 1InstituteofAstronomy,MadingleyRoad,Cambridge,CB30HA
0 2ITA,Zentrumfu¨rAstronomiederUniversita¨tHeidelberg,Albert-Ueberle-Str.2,69120Heidelberg
2
3DepartmentofPhysicsandAstronomy,UniversityofLeicester,Leicester,LE17RH

n 4Presentaddress:Max-Planck-Institutfu¨rKernphysik,Saupfercheckweg1,69117Heidelberg,Germany
a
J

4
2February2008
1

]
h ABSTRACT
p We investigate the orbital evolution of planetesimals in a self-gravitating circumstellar disc
- in the size regime (∼ 1−5000 km) where the planetesimals behave approximately as test
o
particles in the disc’s non-axisymmetricpotential. We find that the particles respond to the
r
t stochastic,regenerativespiralfeaturesin the disc byexecutinglargerandomexcursions(up
s
a toafactoroftwoinradiusin∼1000years),althoughtypicalrandomorbitalvelocitiesareof
[ orderonetenthoftheKeplerianspeed.Thelimitedtimeframeandsmallnumberofplanetes-

imalsmodeleddoesnotpermitustodiscernanynetdirectionofplanetesimalmigration.Our
1
chiefconclusionisthatthehigheccentricities(∼ 0.1)inducedbyinteractionwithspiralfea-
v
turesinthediscislikelytobehighlyunfavourabletothecollisionalgrowthofplanetesimals
8
0 inthissizerangewhilethediscisintheself-gravitatingregime.Thusif,asrecentlyarguedby
1 Riceetal2004,2006,theproductionofplanetesimalsgetsunderwaywhenthediscisinthe
2 self-gravitatingregime (either at smaller planetesimal size scales, where gas drag is impor-
. tant,or via gravitationalfragmentationof the solid component),then theplanetesimalsthus
1
producedwouldnotbeabletogrowcollisionallyuntilthediscceasedtobeself-gravitating.
0
8 It is unclear, however, given the large amplitude excursions undergone by planetesimals in
0 theself-gravitatingdisc,whethertheywouldberetainedinthediscthroughoutthisperiod,or
: whethertheywouldinsteadbelosttothecentralstar.
v
i Key words: accretion, accretiondiscs: gravitation:instabilities: stars: formation:planetary
X
systems:formation:planetarysystems:protoplanetarydiscs
r
a
1 INTRODUCTION Once the disc is no longer self-gravitating in the dominant
(gas)component,itisevidentlyimpossibletoformplanetsthrough
Itiscurrentlyunclearwhethertheprocessofplanetformationbe-
either gas phase Jeans instability or dust accumulation in spi-
ginsduringtheearliestphasesofstarformation(i.e.incircumstel-
ral arms and thus the most successful models invoke the colli-
lardiscsthat arestrongly self-gravitating).Thehigh discmassat
sional growth of dust (Pollacketal. 1996) followed possibly by
this stage (with gas mass ∼ 10% of the central star’s mass (e.g.
the accretion of a gaseous envelope. Despite the lower surface
Eisner and Carpenter 2006) is favourable to this hypothesis and
densities during later phases (estimates of gas disc mass from
providestheonlyopportunityforgasgiantstoformvia(gasphase)
sub-millimetre dust measurement suggest typical masses that are
gravitationalinstability(Boss2000).Amoresubtleeffect-again
around an order of magnitude lower than in the self-gravitating
stronglyfavouredduringtheself-gravitatingcase-arisesthrough
case; Andrews&Williams 2005), the great advantage for planet
theinterplaybetweenpressuregradientsinspiralarmsandgasdrag
formation issimply the longer duration of this later phase. Discs
onsolidparticlesandresultsintheaccumulationofdustinspiral
with sufficient gas and dust to form planets typically survive
arms(Riceetal.2004).Thisdustmaythenbeabletoaccumulate
for 106 − 107 years (Haischetal. 2001; Armitageetal. 2003;
further through the action of collisions and/or self-gravity in the
Wyattetal. 2005), whereas the existence of dusty debris discs
dustphase(Riceetal.2006).1
aroundolderstarssuggeststhat(potentiallyplanetbuilding)colli-
1 Notethatsuchrapidaccumulation,whichmightresultintherapidforma-
tionoflargesizeplanetesimals,isrequiredinorderforthegrowingplan- jecttoveryfastmigrationdownintothehottestinnerdisc(Weidenschilling
etesimals to escape the dangerous meter-size range, where they are sub- 1977).
2 M. Britsch,C. J. ClarkeandG. Lodato
sionsbetweenrockyplanetesimalsproceedsfor108 yearsormore axisymmetric structure in a disc that is subject to the magneto-
(Kenyon&Bromley 2004). These numbers should be contrasted rotational(MRI)instability.Inthisstudyitwasfoundthattheplan-
with the brief window of opportunity (∼ 105 years) in the self- etesimalsunderwentstochasticmigration,butthatitwasnotpossi-
gravitatingphase. bletodiscernthemeanmagnitude(orevensign)ofthemigration
In summary, then, it is unclear without detailed calculation rateoverthelimitedtimeframeofthesimulation(afewhundredor-
whichoftheseregimesisfavouredinthetrade-offbetweenlifetime bits).Inourexperimentalso,conductedoverasimilartimeframe
anddiscmass. toNelson(2005),weareunabletodiscernanetdirectionofmigra-
Theviabilityofplanetformationintheself-gravitatingregime tion,butwefind(perhapsunsurprisinglygiventhefactthatourdisc
howeverrequiresthatanumberofconditionsaresatisfiedbeyond isabout anorderof magnitudelargerinmass)thattheamplitude
theinitialphasesofaccumulationoutlinedabove.Inparticular,we ofthestochasticvariationsintheorbitalradiusismuchlarger,with
needtounderstand thedynamical evolutionof theprotoplanetary severalplanetesimalsundergoingradialexcursionsofafactortwo
components(whetheratthescaleofdustaggregates,planetesimals ormoreinbothdirections.Evidently,thisbehaviourhasimportant
orproto-giantplanets)withaviewtoansweringthefollowingques- implicationsfortheabilityofplanetesimalstogrowbycollisions
tions:(i)doesorbitalmigrationresultinsuchbodiesspirallinginto duringtheself-gravitatingphase.
thecentralstar,orcantheyberetainedinthediscatleastoverthe Thestructureofthepaper isasfollows.InSection2wede-
initialself-gravitatingphase, and(ii)(inthecaseof dust or plan- scribethemodelingoftheself-gravitatingdisc,whichismaintained
etesimals) arethekinematics of suchbodies conducive tofurther ina’self-regulated’(i.e.nonfragmenting)statethroughtheimpo-
collisionalgrowth? sition of an ad hoc cooling law for which the cooling time is a
As afirst step inanswering thisquestion we here undertake suitably large multiple of the local dynamical time. In Section 3
pilot simulations of a small number of test particles in a self- weanalysetheorbitalresponseof12testparticlesplacedinadisc
gravitatingdisc,withaviewtobothconsideringtheissueoftheir ofmass0.1and0.5M⊙.Wealsoanalysethelocationofthedomi-
orbital migration and to establishing the velocity dispersion that nanttorquecontributionontheplanetesimalsandassesstherange
is set up through interaction with spiral structure in the disc. It ofplanetesimalmassandsizescalesforwhichweexpect’testpar-
should be stressed that the calculation differs in several impor- ticle’behaviourandforwhichgasdragcanbeignored.Insection
tant ways from most calculations of orbital migration and plan- 4wediscusstheimplicationsofournumericalfindingsforthemi-
etesimalkinematicsreportedintheliterature.Firstly,conventional grationandcollisionalgrowthof planetesimalsinself-gravitating
planetary migration, whether classified as Type I or Type II, is discs.Section5containsourchiefconclusions.
the consequence of the gravitational torque exerted on the planet
(or planetesimal) by non-axisymmetric structure induced by the
planet/planetesimals. Sincethis torque contribution can have sig-
nificant contributions at the size scale of the Hill radius (or be- 2 NUMERICALSIMULATIONS
low: Bateetal. 2003), it is evidently important that this scale is
2.1 Numericalcode
well resolved innumerical calculations, which makes the model-
ingofthemigrationoflowmassobjectsparticularlychallenging, In this paper we simulate the interaction of a gaseous self-
especially with codes such as Smoothed Particles Hydrodynam- gravitatingdiscwithasmallnumber ofplanetesimalsusingSPH,
ics(SPH)(deVal-Borroetal.2006).Inthepresentcase,however, a Lagrangian particle based method (Lucy 1977; Benz 1990;
weareinterestedinaregimewherethedischaspronouncednon- Monaghan 1992). This allows the simulation of gaseous and N-
axisymmetric structure in the absence of the planet/planetesimal bodyparticles,usingindividualtime-steps(Bateetal.1995),thus
(duetothepresenceofself-gravitatingspiralmodesinthedisc)and resultinginalargesavingincomputationaltimewhenalargedy-
wearethereforenotobliged toresolvetheHillradiusinorderto namicrangeoftimestepsisinvolved.
capturethedominanttorquecontribution.Evidently,asthemassof The general setup of our simulations is similar to
theplanetesimalisincreased,thecontributiontothetorquewhich Lodato&Rice(2004,2005).Weconsideramassivegaseousdisc
it experiences from self-induced spiral structure becomes signifi- orbitingarounda1M⊙centralstar.Thediscistakentoextendini-
cant and we will assess this a posteriori in order to set an upper tiallyfrom0.25to25AUandwehaveconsideredtwocases,where
mass scale on the applicability of our calculations. Secondly, in thediscmassisequalto0.1M⊙(standardruns)and0.5M⊙(massive
conventionalcalculationsofthekinematicsofplanetesimalswarms runs),respectively.
(Kokubo&Ida 2000; Thommesetal. 2003), the velocity disper- OursimulationsemployanumberN =250,000particlesand
sion of planetesimals is set by an equilibrium between what is arethusintermediate-highresolutionsimulations.Previoussimula-
sometimestermed’viscousstirring’(i.e.thetransferofkineticen- tions(Lodato&Rice2004,2005)atthesameresolutionhavebeen
ergy from larger bodies to the planetesimal swarm via two-body showntoreproducereasonablywelltheinternaldiscdynamics.In
scattering) and eccentricity damping by gas drag. In the present particular,oncethediscreachesaquasi-steadyself-regulatedstate
case,however,thevelocitydispersionofourparticlesissetbythe (seebelow),thediscthicknessH isresolvedwithanaverageof2
leveltowhicheccentricityispumpedthroughinteractionwiththe smoothinglengthsforthelowdiscmasscaseandwith4smoothing
fluctuatingpotentialoftheself-gravitatingdisc.Inthispilotcalcu- lengthsforthehotter,highdiscmasscase.
lationweomitgasdrag,whichplacesalowerlimittothesizescale Asmentionedabove, ourcodefollowsthedynamicsofboth
ofobjectstowhichthecalculationisapplicable;again,wesetthis SPHparticles and of point masses, which are able to accrete gas
lowerlimitaposteriori. particles, ifthey come closer thana given accretion radius tothe
Althoughwehavestressedthequalitativedifferencesbetween point mass. For the central object, such accretion radius is set to
thiscalculationandthatreportedinthelargebodyofliteratureon 0.25AU.Attheendofthesimulationstypicallyafractionnolarger
planetarymigrationandplanetesimalkinematics,wepointoutthat than10%oftheinitialdiscmasshasbeenaccretedontothecen-
it bears close comparison with the calculation of (Nelson 2005), tralstar. Theplanetesimals haveanominal mass equal to thede-
whichstudiestheresponseofplanetesimalstothefluctuatingnon- fault gas particle mass, that is 4×10−7M⊙ for the standard runs
Eccentricitygrowthof planetesimalsin a self-gravitatingprotoplanetarydisc 3
and2×10−6M⊙ forthemassiverun.Thischoiceensuresthatthe
gravitational force of individual planetesimals does not influence
theoveralldiscstructureandthedynamicsoftheplanetesimals.In
ordertopreventplanetesimalsfromaccetinggaseousparticles,we
havesettheiraccretionradiustoaverysmallvalue.
The gas particles are allowed to heat up via pdV work and
artificialviscosityandtocooldownwiththefollowingsimplepre-
scription
du u
=− (1)
dt(cid:12)(cid:12)cool tcool
(cid:12)
whe(cid:12)(cid:12)rethecoolingtime-scaletcool = βΩ−1 and βisfixedinspace
andtime.Whilethisisnotarealisticprescriptionfor thecooling
inactualprotostellardiscs,itprovidesaconvenientwayofparam-
eterisingthecoolingphysicsinorderthatonecanensurethatone
ismodelingdiscsthatremainintheself-gravitatingregimewithout
fragmenting.Inparticular,ifthevalueoftheparameterβislarge
enough(β & 5−6;Gammie2001;Riceetal.2005),itallowsthe
realizationofafeedbackloopwherethediscevolvesintoaquasi-
steady,self-regulatedcase,andismaintainedclosetomarginalsta-
bility. For smaller values of β the disc undergoes fragmentation.
Inoursimulationswehavesetβ = 7.5.Wenotethatsincethisis
closetothecriticalvalueofβ,weexpectsuchadisctoberelatively
Figure1. Initialpositionsoftheplanetesimalsintheevolvedquasi-steady
‘active’(inthesensethatthemaintenanceofthermalequilibrium
spiral structure of the 0.1 M⊙ gas disc at the beginning ofthe run. The
requiresspiralmodestobeofrelativelylargeamplitude).Werevisit
dimensions are60AUineachdirection andthelogarithmic colourscale
thispointinSection4. showsthesurfacedensityΣandgoesfrom10to105gcm−2.
2.2 Initialconditions
Initiallythegasinthediscisrelativelyhot,withatemperaturepro-
file T(R) ∝ R−1/2 and a minimum value of the stability parame-
ter (Toomre 1964) Q = 2, attained at the outer disc edge. Since
marginal stability corresponds to Q = 1, the whole disc is ini-
tiallygravitationallystable.Weinitiallyfollowtheevolutionofthe
gaseous disc only, as it cools down and becomes gravitationally
unstable.Theinitialevolutionisverysimilartotheonedescribed
in several other papers for example (Lodato&Rice 2004, 2005;
Mejiaetal. 2005): when Q ≈ 1, the discdevelops a spiral struc-
ture, inducing dissipation through compression and mild shocks
andthereforeheatingupthediscandestablishingafeedbackloop
whichkeepsthediscclosetomarginalstability,sothattheToomre
Qparameterisclosetounityovertheradialrange∼1−20.Thisini-
tialevolutiontakesroughly60outerdynamicaltimesforthestan-
dardrunsand45outerdynamicaltimesforthemassiveruns.Inthis
state,thesurfacedensitydistributionisdominatedbypronounced
(oforderunityinamplitude)spiralfeatures,which,althoughregen-
erative,areindividuallytransientstructures(i.e.withlifetimesthat
areoforderthelocaldynamicaltimescale).
Anensembleoftwelveplanetesimalswasthenintroducedbe-
tween10.8and18.8AU,alignedinfourequispacedradialspokes
(seeFigure1forthe0.1M∗discandFigure2forthe0.5M∗disc). Figure2. SameasFig1butforthe0.5 M⊙ gasdisc.Notethewarmer
The default velocities were chosen to be Keplerian with a small
conditions in the massive disc in the self-regulated state lead to broader
correction for the gravitational potential of the disc, although we spiralfeaturesthaninthestandardcase(Fig1)
alsoinvestigatedthecasewheretheplanetesimalswereintroduced
withanon-zeroinitialeccentricity.
planetesimals are initially set in nearly circular orbits. Although
someplanetesimalsundergolargeamplitudeexcursions(factortwo
3 RESULTS or more), theaverageorbital radius of theensemble variesby no
morethan10%duringthesimulation,indicatingthatplanetesimals
3.1 Orbitalevolutioninthe0.1solarmassdisccase
inagivenregionareapparentlyaslikelytogoinwardsasoutwards.
InFig.3weshowthetimeevolutionoftheorbitalradiiofallthe Thethickdashedlineinfig.4showsthemeaneccentricityofthe
planetesimals for the case where M = 0.1M and where the planetesimalsasafunctionoftime,showingthate∼ 0.05istypi-
disc ⋆
4 M. Britsch,C. J. ClarkeandG. Lodato
0.4 3e+39
0.35 2.8e+39
0.3 2.6e+39
s
citie 0.25 2.4e+39
eccentri 0.2 - E [J] 2.2e+39
an 0.15 2e+39
e
m
1.8e+39
0.1
1.6e+39
0.05
1.4e+39
0 6e+47 6.5e+47 7e+47 7.5e+47 8e+47 8.5e+47 9e+47
0 50 100 150 200 250 300 J [kg m2 / s]
t [yr]
Figure4. Timeevolutionofthemeaneccentricityoftheplanetesimalsfor Figure6. Energyandangularmomentum(intheinertialframe)areplotted
variousinitialeccentricities: e(t = 0) ≈ 0(bolddashedline),e(t = 0) = forthreeplanetesimalsinthemassivedisccaseoverthetimeinterval550−
0.06(boldsolidline), e(t = 0) = 0.03(thinsolidline)e(t = 0) = 0.37 700yr.Thenearlinearrelationship,withthesameslope,suggeststhatthese
(thindashedline).Duringthecourseofthesimulationsaquasi-steadymean threeplanetesimals areinteracting withamodewithpatternspeedΩP =
eccentricitye≈0.07isestablishedinallcaseswheretheinitialeccentricity 4−5×10−9s−1,i.e.withco-rotationat∼14A.U..
islessthan0.1.
beconserved, andthatsince E isrelatedtotheenergy, E, inthe
J
cal(althoughindividualplanetesimalsmaytemporarilyattainstill inertialframe,andangularmomentum,L,via
highereccentricities,∼0.3).Thusalthoughmuchofthefinestruc-
E =E−LΩ (2)
tureinFigure3issimplyaconsequenceofellipticalorbitalmotion, J p
itisalsoevidentthatinadditionsomeplanetesimalsundergoabrupt (whereΩ isthepatternspeed),wecantestforsuchbehaviourby
p
and/orlargescalechangesinorbitalradius(forexample,oneplan- plotting E versus L. This exercise indicates that there are indeed
etesimalmigratesover20AUin4000years,andanotherfollowsa periodswhencertainplanetesimalsappeartobeexhibitingthisbe-
periodofrelativeorbitalquiescencebymigrating∼5AUinabout haviour.ThiscanbeevidencedbyalinearrelationbetweenEand
200years).Suchmigrationratesimplyradialvelocitieswhichare L,asshowninFigure6)forthreeofthesimulatedplanetesimals.
attimesasignificantfractionoftheKeplerianspeed. Fromtheanalysisofthisrelation,wecanthendeterminethepattern
InFigure4,wealsoshowtheresultsofstartingtheplanetes- speedofthedrivingpattern.TheplanetesimalsillustratedinFigure
imal ensemble with non-zero initial eccentricity, and see that the 6appeartobeinteractingwithamodewithco-rotationataround14
results are unchanged for initial eccentricity less than ∼ 0.1. On AU:inspectionofananimatedseriesofsnapshotsofthesimulation
theotherhand,ifalargerinitialeccentricityisemployed,themean indicatesthatthesethreeplanetesimalsareindeedroughlycorotat-
eccentricityoftheensembleremainsaroundthisvalueforthedu- ingwithaspiralfeatureduringthisperiod,butthat,asthefeature
rationoftheexperiment. dissolvesshortlythereafter,theythenmigrateandstallclosetoan-
otherspiralarmatadifferentradius.Thisbehaviourisnotseenin
thestandard(0.1M⊙disc)simulations,whichisconsistentwiththe
3.2 Orbitalevolutioninthe0.5solarmassdisccase expectation(see,forexample,Lodato&Rice2005)thatthemore
massive discproduces rather morecoherent and long-lived spiral
Figure5illustratestherichvarietyoforbitalhistoriesofplanetes-
features.Nevertheless,wefindonlyoccasionalevidenceforplan-
imals in the massive disc case. The amplitude of orbital migra-
etesimalsthat are being driven by a dominant quasi-steady spiral
tionissomewhat larger than inthe standard case, although again
mode sothat theanalysis of orbital familiesin thissituation(see
the change in average orbital radius for the ensemble is small,
e.g.Vorobyov&Shchekinov2006)isoflimitedapplicabilityhere.
< 20%. We note that, as in the standard case, the planetesimals
undergostochasticmigrationeventsevenwhentheyarelocatedat
radii(> 20AU)whereQisgenerallysomewhatlargerthanunity.
3.3 Planetesimal-planetesimalencountersandgasdensities
Theseeventscoincidewithspiralstructurestemporarilyextending
closetotheplanetesimals
outintotheserelativelyquiescentregions.
The mean eccentricity of the ensemble is about double its Inspectionofanimationsofthesimulationssuggestthatplanetes-
value in the standard case ∼ 0.15, with individual planetesimals imals tend to spend time preferentially in spiral structures in the
temporarilyattainingeccentrcitiesofupto0.5.Althoughsomeof disc.Assuchstructuresdissolveonaroughlydynamicaltimescale,
thequasi-periodicstructureinFigure5issimplyattributabletoel- theplanetesimalsthenmigrateandfindfurthertemporarylodging
lipticalorbitalmotion(i.e.atroughlyconstantangularmomentum), in a new spiral structure. (The strongly fluctuating nature of the
we also see planetesimals in which, temporarily, the energy and spiral potential means that, asdiscussed above, theplanetesimals
angular momentum aresubject tocorrelatedquasi-periodicvaria- usually have insufficient time to attain steady orbits in the frame
tions.Suchbehavioursuggeststhattheseplanetesimalsareatthis co-rotatingwiththelocalspiralpattern).Thispredilectionforbe-
stagebeingdrivenbyaquasi-steadyrotatingnon-axisymmetricpo- inginthespiralarmssuggeststhattheplanetesimalsareexpected
tential.Wehaveinvestigatedthispossibilitybynotingthatinthis tocomeclosertoeachotherandtosamplehighergasdensitiesthan
casetheJacobiconstant(E )oftheplanetesimal’smotion(i.ethe inthecaseofparticlesorbitinginanaxisymmetricdisc.
J
totalenergymeasuredintheframeoftherotatingpattern)should We illustrate the increased tendency for planetesimal-
Eccentricitygrowthof planetesimalsin a self-gravitatingprotoplanetarydisc 5
35
30
25
20
U]
A
R [ 15
10
5
0
0 500 1000 1500 2000 2500 3000 3500 4000
t [yr]
Figure3.EvolutionofthedistancetothestarofthreerepresentativeplanetesimalsfortheMdisc=0.1M⊙case.Whileindividualplanetesimalsmightmigrate
significantlyoverthecourseofthesimulation,theaverageorbitalradiusonlyvariesbylessthan10%.
30
25
20
U]
A 15
R [
10
5
0
0 200 400 600 800 1000 1200 1400 1600 1800
t [yr]
Figure5.EvolutionofthedistancetothestarofthreerepresentativeplanetesimalsfortheMdisc=0.5M⊙case.
d<0.5 0.5<d<1 1<d<1.5 1.5<d<2 η= 2(ρp−ρR), (3)
ρ +ρ
p R
0.1 M⊙-disc 10 17 17 19
whereρ isthegasdensityinthevicinityoftheplanetesimaland
non-self-grav. 0 17 13 16 p
ρ istheazimuthallyaverageddensityattheinstantaneousorbital
R
Table1.Numberofencounterswithclosestapproachsmallerthand(mea- radiusoftheplanetesimal.Thisplanetesimalappearstospendmost
suredinAU)betweenplanetesimalsduringthefirst15000timeunitsinthe ofitstimeinthespiralarmsandthusatlargergasdensities.This
0.1M∗discrunandthenon-self-gravitatingdiscrun. istypical,althoughsomeplanetesimalsalsoundergophaseswhere
they sample lower gas densities as can be seen in Fig. ??. Such
periodsinunder-denseregionsarelessapparentinthecaseofthe
planetesimal encounters in Table 1, where we contrast the num- massivediscsimulation.
berofencountersbetweentheplanetesimalsinvariousseparation WeseefromFigure7thatalthoughthisplanetesimalsamples
rangeswiththecorrespondingstatisticsfromacontrolruninwhich gas densities larger than the local azimuthal average, η isalways
theplanetesimalsorbitinasmooth, non-self-gravitatingdisc. We lessthanoroforderunity.Thisistobeexpected:theamplitudeof
see that the number of encounters with closest approach smaller thespiralfeaturesinthegasisonlyoforderunityandthusevenif
than 0.5 AU is significantly enhanced in the self-gravitating disc planetesimalsspent all theirtimeinspiralarms, theywouldonly
case.Asexpected,theeffectofplanetesimalslingeringnearspiral experienceadensitythatexceededthelocalmeanbyamodestfac-
featuresistodecreassetheminimumdistance betweenplanetesi- tor.
malsinthesimulation.Westressthat,giventheverylow masses
oftheplanetesimals(i.e.equalinmasstoasinglegasparticle),the
3.4 Numericalissues
effectofmutualencountersontheplanetarydynamicsisnegligible
-asintended-evenfortheclosestencountersinthesimulation. Firstly,ourexperimentshaveaimedtostudytheresponseofeach
We illustrate the effect of non-axisymmetric structure in the planetesimal tothefluctuating disc potential and haveonly mod-
disconthegasdensityfieldsampledbytheorbitingplanetesimals eledtheevolutionofanensemble ofplanetesimalsforreasonsof
inFig.7.Weplot -foragivenplanetesimal -thedensityexcess, computationaleconomy.Therefore,weneedtobesatisfiedthatthe
whichwedefineas stochastic behaviour we see is a result of planetesimal-disc and
6 M. Britsch,C. J. ClarkeandG. Lodato
0.5 1
0.4
0.3 0.8
0.2
0.1 que 0.6
η -0 .01 ated tor
-0.2 egr 0.4
nt
i
-0.3
-0.4 0.2
-0.5
-0.6
0 500 1000 1500 2000 2500 3000 3500 4000 0
0 2 4 6 8 10 12 14 16 18 20
t [yr] R [AU]
Figure7.Comparisonbetweenthedensityinthevicinityofoneplanetes- Figure9.Integrated zcomponentofthetorqueasafunctionofdistance
imalsinthe0.1M∗discwiththemeangasdensityatthecurrentradiusof fromtheplanetesimal,normalisedtothetorquefromwithin20AUofthe
theplanetesimal. planetesimal.
M 1/3
0.4 Rh=Rp 3Mp ! , (4)
⋆
0.2
where R is the radial coordinate of the planetesimal and M the
0 p p
massoftheplanet,sincethisistheregionfromwhichthebulkof
-0.2 thetorquesontheplanetesimalarise.
-0.4 Inthepresentcase,bycontrast,whereweexaminemigration
η -0.6 drivenbyspiralstructureintheself-gravitatingdisc,itisnolonger
important toresolvetheHillradius andwecan thusexamine the
-0.8
migration of low mass planetesimals for which resolution of the
-1
Hillradiuswouldbeimpracticable.Nevertheless,weneedtocheck
-1.2 thatourassumptioniscorrect,i.e.thatthedominanttorquesindeed
-1.4 ariseatrelativelylargedistancesfromtheplanet.Fig.9illustratesa
-1.6 typicalsnapshotofthecumulativedistributionforthezcomponent
0 500 1000 1500 2000 2500 3000 3500 4000 ofthetorqueofthediscontheplanetesimal asafunctionofdis-
t [yr]
tancefromtheplanetesimal.Weseethatthedominantcontribution
Figure8. Evolutionofthedensityexcessη(seeequation(3)asinFigure tothetorqueisaround1AUfromtheplanetesimalandthatthecon-
7,butforadifferentplanetesimal. tributiontothetorquefromthepoorlyresolvedregionaroundthe
planetesimal(i.e.fromwithinatypicalparticlesmoothinglengthin
thatregion,around0.2AU)issmall(<20%).Wearethussatisfied
not planetesimal-planetesimal interactions. We find, in fact, that thattheorbitalevolutionisnotdrivenbythedisc’sresponsetothe
the torques experienced by the planetesimals due to other plan- planetesimalandthattheplanetesimalshouldbebehavingapprox-
etesimalsaretypicallytwotothreeordersofmagnitudelessthan imatelyasatest particle. Weexamine thisassumption further by
those arising from the disc. This is to be expected. The closest investigatinghow themigration depends on planetesimal massin
planetesimal-planetesimal interactions occur at distances of ∼ 1 Section3.5.
AU:giventhatatypicalsmoothinglengthinthediscis∼ 0.2AU Ontheotherhand,inthecaseofthecomponentofthetorque
and that a smoothing kernel contains 50 particles, it follows that perpendicular toz,i.e.actingontheinclination,thecontributions
evenatclosestapproachthereareabout6000gasparticlesthatare fromclosetotheplanetesimal actuallycanmakeupalargefrac-
closer to a planetesimal than its nearest planetesimal neighbour. tionofthetorque(i.e.around40%ofthetorqueoriginatingfrom
Giventhattheplanetesimalmassesinthesimulationsareequalto within a smoothing length). Thus the inclination development of
thegasparticlemasses,itfollowsthatthegravitationaltorquesex- thesimulationscannotbetrusted.Neverthelesstheinclinationsstay
ertedbetweenplanetesimalsarenegligible. relativelysmall(typicallylessthan0.02radiansinthestandardand
Secondly, we are interested here in examining the regime 0.03forthemassivedisc).Thisissmallcomparedwiththeaxisra-
where the dominant torques experienced by the planetesimal de- tioofthediscH/R ∼ 0.1andsotheplanetesimalsremainmainly
rive from the non-axisymmetric density distribution that the disc confinedwithinthedisc.
wouldhaveintheabsenceoftheplanetesimal(i.e.thatduetothe
spiralstructureinthedisk)ratherthanthenon-axisymmetricstruc-
3.5 Validregime
tureinducedinthediskbytheplanetesimal.Notethatinnon-self
gravitatingdiscs,andthusinalldiscussionsofplanetesimalmigra- Wenowdiscusstherangeofplanetesimalmassesforwhichthecur-
tionintheliterature(withtheexceptionofthatinNelson2005)itis rentcalculationisvalid.AsnotedinSection2.1,wearemodelling
thelattereffectthatisresponsibleforplanetmigration.Inthiscase, theregimeinwhichtheplanetesimal behavesasatestparticlein
itisessentialthatthecalculationproperlyresolvestheHillradius thegravitationalpotentialofthestar-discsystemandwherethere-
oftheplanetesimal sultsshouldbeindependentofplanetesimalmass.Inreality,atlow
Eccentricitygrowthof planetesimalsin a self-gravitatingprotoplanetarydisc 7
massestheeffect ofgasdrag(omittedinthepresent simulations) 4
wouldbreakthismassindependence, whereasathighmasses,the
discresponsetothepresenceoftheplanetesimalaffectsthetorques 3
ontheplanetesimal.Weconsidereachoftheseissuesinturn.
Alowerlimitforthemassrangeofoursimulationcanbeob- nits] 2
tainedbycalculatingthemassscaleatwhichthestoppingtimedue y u
todragiscomparablewiththetimescaleonwhichtheplanetesimal bitrar 1
undergoesstochasticmigrationduetointeractionwithspiralstruc- ar
tureinthedisc.FrominspectionofFigure3weseethatplanetesi- M [Pl 0
malstypicallyreversetheirdirectionofmigrationontimescalesof T/
about103years.Thusifwedemandthatthetimescaleforgasdrag -1
tosignificantlyperturbtheorbitshouldexceedsay104years(∼58
orbits),weshouldsafelybeintheregimewheregasdragplaysa
-2
minorroleinthedynamicsofthesystem. 0 20 40 60 80 100 120 140 160
t [yr]
ThedragforceonsolidbodieswithinaprotostellardiscF is
D
givenby(Weidenschilling1977;Whipple1972) Figure10.Thetorqueperunitplanetesimalmassasafunctionoftimefor
runsinwhichthemassoftheplanetesimalisrespectively1/10(bolddashed
1
FD = 2CDπa2ρu2, (5) line),1(boldsolidline),10(thinsolidline)and100(thindashedline)times
themassofagasparticle.
where aisthe radiusof theplanetesimal, uisthe velocity of the
planetesimalrelativetothegas,ρisthegasdensityandC isthe
D
dragcoefficient. 35
normal run
Thusthetimescaleonwhichrelativemotionbetweentheplan- run with planets of 1/10th the mass
etesimalandthediscgasisdampedbygasdragis 30 rurunn w w/ /p plalanneetst so of f1 1000 ttiimmeess tthhee mmaassss
t = Mpu. (6) 25
e F
D
20
Inthe case that weareconsidering here, therelativevelocity be- AU]
tweentheplanetesimalandthediscgasdoesnotderive-asinthe R [ 15
caseusuallyconsidered(e.g.Weidenschilling1977)-fromthedif-
ferencebetweenKeplerianplanetesimalmotionandsub-Keplerian 10
gasmotionsubjecttooutwardpressureforces.Instead,thisrelative
5
velocity,u,derivesfromtheeccentricityoftheplanetesimalorbits
inducedbythespiralstructureinthedisc,sothatwehaveu=eVK, 0
whereV isthelocalKeplerianvelocity. 0 20 40 60 80 100 120 140
K
t [yr]
Assumingasphericalplanetesimalofdensityρ ,wecanthen
p
writeteas Figure11. Thedistancetothestarisplottedagainsttimeforoneplanetes-
8ρ R imalinthenormalrun(solidline),therunwithatenthplanetesimalmass
t = p (7) (longdashedline),therunwith10timeshigherplanetesimalmasses(short
e 3C ρeV
D K dashedline)andtherunwith100timeshigherplanetesimalmasses(very
or,equivalently(andadoptingC intherange0.5−1,whichis shortdashedline).
D
appropriatetothecaseconsideredherewheretheReynoldsnumber
exceedsafewhundred)
tothetorquesontheplanetesimalsonlyariseswhentheplanetesi-
t ∼ 8 ρp aΩ−1 (8) malmassisraisedabovetentimestheoneusedinourstandardrun,
e 1.32 ρ R e corresponding toabout 8×1027g.Wereachthesameconclusion
frominspectionofFig.11,whereagainweseeamarkeddeviation
whereRistheorbitalradiusoftheplanetesimal.Ifwethenrequire
that t > 104 years (∼ 300Ω−1), and adopting typical values of inorbitalevolutionforthelargestmassplanetesimal.
e∼0e.1,ρ∼afew×10−11 gcm−3 andρ ∼afewgcm−3,wethen Sosummarising,wehaveasaregionofvalidityofoursimu-
p lations:
obtaintherequirementthatamustbeatleastoforderakilometre
orsobeforeonecanneglecttheeffectofgasdrag. 2·10−8M⊕ ∼1016g<Mp <8·1027g=1M⊕ (9)
Wethusconcludethattheeffectofgasdragisnegligiblefor
orinplanetesimalradii
planetesimalsofkilometrescaleorabove(correspondingtobodies
ofmasses∼1016g.) 1km<a<5000km. (10)
Oursimulationsarevalidforallhighermassesaslongasthe
accelerationsexperiencedbytheplanetesimalsremainindependent
of mass. But theseaccelerations change as soon asthe planetesi-
4 DISCUSSION
malshaveenoughmasstogravitationallyinfluencethegasdensity.
Thusweget anupper limitfor thevalidityof our simulations by Wecansummariseournumericalfindingsasfollows.Planetesimals
checkingatwhatmassscalethetorqueperplanetesimalmassT/M inthesize/mass rangetowhichour calculationsapply(seeequa-
p
changeswithmass.InFigure10weplottheevolutionofT/M for tion10)aresubjecttolargescalefluctuationsinorbitalradiusdue
p
asingleplanetesimal inthecase that ithas themassindicated in tointeractionwithregenerativespiralstructureinthedisc.Wecan
eachofthecurves.ItcanbeseeninFig.10thatsignificantchanges discernnonet signtothemeanmigrationrateover thetimescale
8 M. Britsch,C. J. ClarkeandG. Lodato
ofthesimulationandarethusunabletoansweroneofthegoalsof ρ2/3M1/3
t = p (11)
thepresentstudy,i.e.whetherweexpectplanetesimalstomigrate grow Σ Ω(1+4GM /(R σ2))
p p p
intothestarduringtheself-gravitatingphaseofthedisc.Sincethe
whereΣ isthesurfacedensityofplanetesimals,σistheplanetes-
disc’slifetimeintheself-gravitatingphaseexceedsthedurationof p
imal velocity dispersion and Ω is the local Keplerian orbital fre-
thepresentexperimentbyaboutafactor30itwillbecomputation-
quency.Foraplanetesimalofdensity∼3g/cm3locatedat10A.U.
allychallengingtoanswerthisquestion,especiallysince(giventhe
inadiscwithsurfacedensityofplanetesimalsequalto0.4g/cm2,
stochastic nature of the orbital evolution) it will be necessary to
thegrowthtimescaleviatwobodycollisionsis:
assessthefractionofplanetesimalsretainedinthediscthroughfol-
lowingtheevolutionofalargeensembleofparticles. M˜1/3
t =5×109 yrs (12)
grow 1+M˜2/3σ˜−2
We can however make some comments about how the dy- where M˜ = (Mp/610−4M⊕)andσ˜ = (σ/1kms−1).Evidentlythis
namical evolution of our particles is likely to affect their col-
timescaleisfartoolongtobeofinterestunlessthesecondtermon
lisonalgrowthduringtheself-gravitatingphase.AsnotedbyNel-
thedenominator(duetoenhancement ofthecollisionalcrosssec-
son(2005)inthecontextoftheevolutionofplanetesimalsinadisc
tionbygravitationalfocusing)islarge.Ifwesimplysetσtobethe
subject to MRI, the fluctuations inorbital radius and eccentricity
productoftheorbitaleccentricityandlocalKeplerianvelocity(an
couldinprincipleeitheraidorhinderplanetesimalgrowth.Onthe
assumptionwerevisitbelow),itfollowsthatatagivenlocationin
positive side, orbital migration provides a mechanism for avoid-
thediscandforagivenmassscaleofplanetesimal,thenthegrowth
ing the self-limitation of planetesimal growth through the clear-
timescaleinthegravitationallyfocusedregimescalesasthesquare
ing out of a proto-planet’s “feeding zone”. In conventional mod-
oftheeccentricitydividedbythesurfacedensitynormalisation.In
els,thisimposesan“isolationmass”limitwhichmakesitdifficult
otherwords,foraplanetesimalatagivenmassscaletobeableto
togrowbeyondaterrestrialmassscaleintheinnerregionsofthe
growsignificantlyoveradiscphaseofdurationt werequire
life
disc(Ida&Lin2004).Althoughthislimitationisby-passedinour
t
self-gravitatingmodels,sincethefeedingzonewouldbecontinu- grow <1 (13)
ouslyreplenished by planetesimalsundergoing an orbitalrandom tlife
walk,thisfactor isnotverysignificantsinceinanycasethehigh whichtranslatesintoalowerlimitontheproductt Σe−2(whereΣ
life
surfacedensityduringtheself-gravitatingphasemeansthatisola- isthesurfacedensitynormalisationofthedisc).
tion masses are high (higher, that is, than the upper limit for the Wecan now return tothe issue raisedin theIntroduction of
applicabilityof our calculations). Another potentiallypositiveas- whetherornotplanetesimalgrowthisfavouredinthehighdensity
pectofevolutioninaself-gravitatingdiscisthepossibilityofcon- butshortlivedphaseduringwhichthediscisself-gravitating.As-
centratingplanetesimalsinhighdensityspiralfeatures,wherethey suming for simplicitythat the self-gravitating phase lasts ∼ 10%
mightbemorepronetoundergocollisionalgrowth.Inthecaseof ofthesubsequent lifetimeoftheprimordialgas/dust discbutthat
the more massive – 0.5M⊙ – disc (see Section 4), there is some its surface density (mass) is typically ten times greater, we see
evidenceforthedominanceattimesofspiralmodeswhichpersist that these two factors roughly cancel. Thus whether the growth
overanumberofdynamicaltimes,allowingparticlestosettleinto of planetesimals is favoured during the self-gravitating phase ac-
orbits which conserve a Jacobi integral in the co-rotating frame. tually boils down to the relative values of the eccentricity in the
Some of these orbits involve a preference for particles to spend twophases.
timeinthespiralarms,animpressionstrengthened byinspection Our simulations have shown that this factor is decisive in
ofanimationswhichshowparticleswhichco-rotateforawhilein disfavouring the self-gravitating regime. In conventional mod-
the vicinity of an arm and which then respond to changes in the els for the dynamical evolution of planetesimal swarms (e.g.
spiral structure by passing rapidly to another spiral feature. This Kokubo&Ida2000)theequilibriumeccentricity(attainedthrough
tendency will not lead to a significantly enhanced collision rate, thebalancebetweenstirringbylargerbodiesanddampingbygas
however, unlessthereisastrongenhancement inthesolidtogas drag)isverylow(oforder0.01).Wehavehoweverseenthatinthe
ratiointhespiralarms:sincetheamplitudeofthegaseousarmsis self-gravitatingphase,thetypicaleccentricityis0.1ormore.Since
aroundafactortwoatmost,thenevenaparticlethatspentallits thegrowthtimescalesasthesquareoftheeccentricity,thismeans
timeinthearmswouldonly,atfixedsolidtodustratio,experience that the self-gravitating phase is highly unfavourable (compared
an enhancement in the local planetesimal density of a factor two withthelater morequiescent phase) for thecollisional growthof
atmost.Sincewehaveonlystudiedtheevolutionoftenplanetesi- planetesimals. We note that the above estimates for conventional
mals,wecannot,onthebasisofthepresentstudy,commentonthe proto-planetarydiscsassumealaminardiscstructure;eccentricity
degreetowhichplanetesimalsaredifferentiallyconcentratedinthe drivingbyinteractionwithturbulentstructuresinadiscexhibiting
arms.Nevertheless,theresultsofRiceetal.(2004)suggestthatthe MRI would further hinder collisional growth, whether or not the
differentialaccumulationofsolidsinthearmsisweakatsizescales discisself-gravitating.
wheregasdragisunimportant. There however remain two possible qualifiers of the above
conclusion.Firstly,itisonlyvalidtosetthelocalvelocitydisper-
Wenowturntothenegativeimplicationofstochasticmigra- sion equal to the product of the eccentricity and local Keplerian
tionnotedbyNelson(2005),i.e.thegrowthofplanetesimalveloc- velocity inthe case that the kinematics iswell described as aset
ity dispersion. This hinders planetesimal growth in two respects: ofellipticalorbitswhichareelongatedinarandomdirection.Such
highercollisionalvelocitiesareassociatedwithdisruptionofparent adescriptiondoes not of course allow for thepossibilityof local
bodiesinsteadofcollisionalgrowth(e.g.Leinhardtetal2007)and velocitycoherenceinaplanetesimalswarminwhichvelocitiesare
alsosuppresstheroleofgravitationalfocusingineffectingphysi- significantly non-circular. In this sense, therefore, the above esti-
calcollisions.Quantifyingthesecondeffect,thetimescaleforcol- matesarepessimistic.Furthercalculations,involvingalargeswarm
lisionalgrowthofaplanetesimalofmassM ,densityρandradius of planetesimals from which one could derive the local velocity
p
R isgivenby: dispersion, are required in order to settle this issue. We however
p
Eccentricitygrowthof planetesimalsin a self-gravitatingprotoplanetarydisc 9
notethatourdifficultyinidentifyingcoherentorbitalfamilies,in- IdaS.,LinD.N.C.,2004,ApJ,604,388
teractingwithalonglivedspiralmode,probablyimpliesthatthis KenyonS.J.,BromleyB.C.,2004,AJ,127,513
effectisunlikelytoimprovetheprospectsforparticlegrowthsig- KokuboE.,IdaS.,2000,Icarus,143,15
nificantly.Thesecondpossiblequalifieristhat,asnotedinSection Leinhardt,Z.M.,Stewart,S.T.,Schultz,P.H.,2007,astroph0705.3943
LodatoG.,RiceW.K.M.,2004,MNRAS,351,630
2,oursimulationsemployacoolingtimescalewhichisclosetothe
LodatoG.,RiceW.K.M.,2005,MNRAS,358,1489
critical one for disc fragmentation and therefore represent condi-
LucyL.B.,1977,AJ,82,1013
tions where spiral modes are of relatively high amplitude. In the
MejiaA.C.,DurisenR.H.,PickettM.K.,CaiK.,2005,ApJ,619,1098
innerregionsofprotoplanetarydiscs,however,coolingtimescales
MonaghanJ.J.,1992,ARA&A,30,543
are longer than employed here and thus one would expect lower NelsonR.P.,2005,A&A,443,1067
amplitude spiral structures, with a correspondingly lower ampli- PollackJ.B.,etal.,1996,Icarus,124,62
tude of stochastic driving of planetesimal orbits. This possibility RiceW.K.M.,LodatoG.,ArmitageP.J.,2005,MNRAS,364,L56
needstobequantifiedthroughfurthernumericalinvestigation. RiceW.K.M.,LodatoG.,PringleJ.E.,ArmitageP.J.,BonnellI.A.,
Ifoneleavesasidethesepossibilities,thenwehaveshownthat 2004,MNRAS,355,543
theself-gravitatingphaseisunfavourableforcollisionalgrowthin RiceW.K.M.,LodatoG.,PringleJ.E.,ArmitageP.J.,BonnellI.A.,
2006,MNRAS,372,L9
the regime of planetesimal scale studied. We have not ruled out
ThommesE.W.,DuncanM.J.,LevisonH.F.,2003,Icarus,161,431
the possibility of collisional growth at smaller scales (where gas
ToomreA.,1964,ApJ,139,1217
dragcaneffectasignificantenhancement ofthelocalsolidtogas
VorobyovE.I.,ShchekinovY.A.,2006,NewAstronomy,11,240
ratioinspiralfeatures; Riceetal.2004) nor that (wheresuch en-
WeidenschillingS.,1977,MNRAS,180,57
hancement exists) it may not be possible to produce much larger WhippleF.L.,1972,inElviusA.,ed.,FromPlasmatoPlanetsWiley,
bodies through the gravitational fragmentation of the solid com- London
ponent (e.g., Riceetal. 2006). We nevertheless consider it likely WyattM.C.,GreavesJ.S.,DentW.R.F.,CoulsonI.M.,2005,ApJ,620,
that whatever the sizescaleattained inthisway, there willbe no 492
further collisional growth of planetesimals while the disc is self-
gravitating.
Thisobviouslylimitstheextenttowhichplanetformationcan
getunderwayintheself-gravitatingphase.Eitheramassiveplanet
isproducedviagasphaseJeansinstability(Boss2000),thusobvi-
atingtheneedforcollisionalgrowth,orelse,ifself-gravityaidsthe
initialaccumulationofsolidbodies(Riceetal.2004,2006),thenit
isunlikelythatsuchbodiescanundergofurthercollisionalgrowth
oncetheybecomelargeenoughforgasdragtobeunimportant.In
this latter case, further growth of such bodies could only resume
oncethedischadevolvedtothepoint whereitsself-gravity(and
theassociatedpumpingupoftheplanetesimalvelocitydispersion)
hasbecomeunimportant.Thevitalunansweredquestion,therefore,
iswhetherorbitalmigrationwouldpreventtheplanetesimalsfrom
remaining in the disc throughout the self-gravitating phase. This
question can only be answered through further, computationally
challenging,calculations.
5 ACKNOWLEDGMENTS
M.BritschgratefullyacknowledgesthehospitalityoftheIOAand
financialsupportviatheMarieCurieEARAEarlyStageTraining
program.Wethankthereferee,PhilArmitage,forconstructivecrit-
icism.
REFERENCES
AndrewsS.M.,WilliamsJ.P.,2005,ApJ,631,1134
ArmitageP.J.,ClarkeC.J.,PallaF.,2003,MNRAS,342,1139
BateM.R.,BonnellI.A.,PriceN.M.,1995,MNRAS,277,362
BateM.R.,LubowS.H.,OgilvieG.I.,MillerK.A.,2003,MNRAS,213,
341
BenzW.,1990,inBuchlerJ.,ed.,TheNumericalModelingofNonlinear
StellarPulsationsKluwer,Dordrecht
BossA.P.,2000,ApJ,536,L101
deVal-BorroM.,etal.,2006,MNRAS,370,529
EisnerJ.,CarpenterJ.,2006,ApJ641,1162
GammieC.F.,2001,ApJ,553,174
HaischK.,LadaE.,LadaC.,2001,ApJ,553,L153