Hybrid methods in planetesimal dynamics: Formation of protoplanetary systems and the mill condition

0.84 MB·English

by  Pau Amaro-Seoane

Mon.Not.R.Astron.Soc.000,1–15(2013) Printed18January2013 (MNLATEXstylefilev2.2)
Hybrid methods in planetesimal dynamics:
Formation of protoplanetary systems and the mill
condition
Pau Amaro-Seoane1 (cid:63), Patrick Glaschke2 & Rainer Spurzem3,4,2
3 1Max Planck Institut fu¨r Gravitationsphysik (Albert-Einstein-Institut), D-14476 Potsdam, Germany
1 2Astronomisches Rechen-Institut, M¨onchhofstraße 12-14, Zentrum fu¨r Astronomie, Universit¨at Heidelberg, Germany
0 3National Astronomical Observatories of China, Chinese Academy of Sciences, 20A Datun Lu, Chaoyang District, 100012, Beijing, China
2 4Kavli Institute for Astronomy and Astrophysics, Peking University, China

n
a
J draft18January2013

6
1
ABSTRACT

] The formation and evolution of protoplanetary discs remains a challenge from both
P
a theoretical and numerical standpoint. In this work we first perform a series of
E
tests of our new hybrid algorithm presentedin Glaschke, Amaro-Seoane andSpurzem
h. 2011 (henceforth Paper I) that combines the advantages of high accuracy of direct-
p summation N−body methods with a statistical description for the planetesimal disc
- based on Fokker-Planck techniques. We then address the formation of planets, with a
o
focusontheformationofprotoplanetsoutofplanetesimals.Wefindthattheevolution
r
t of the system is driven by encounters as well as direct collisions and requires a careful
s
modelling of the evolution of the velocity dispersion and the size distribution over a
a
[ largerangeofsizes.Thesimulationsshownoterminationoftheprotoplanetaryaccre-
tionduetogapformation,sincethedistributionoftheplanetesimalsisonlysubjected

1 to small fluctuations. We also show that these features are weakly correlated with the
v
positions of the protoplanets. The exploration of different impact strengths indicates
0
that fragmentation mainly controls the overall mass loss, which is less pronounced
1
during the early runaway growth. We prove that the fragmentation in combination
9
with the effective removal of collisional fragments by gas drag sets an universal upper
3
. limit of the protoplanetary mass as a function of the distance to the host star, which
1
we refer to as the mill condition.
0
3 Key words: protoplanetary discs, planets and satellites: dynamical evolution and
1 stability, methods: numerical, methods: N-body, methods: statistical
:
v
i
X
r
1 INTRODUCTION Understanding planet formation comprises many chal-
a
lenges, such as hydrodynamics of the protoplanetary disc,
The origin of our solar system remains to be one of the
chemical evolution of the embedded dust grains, migration
most exciting problems of today’s astronomy. For a long
ofplanetsandplanetesimalsandevenstar-starinteractions
time it has been the only known planetary system. While
indenseyoungstarclusters(seeArmitage2010,forareview
it is still the only planetary system that can be studied in
andreferencestherein,andalsotheintroductionofPaperI,
detail,progressinobservationtechniqueshasledtothedis-
for a brief summary). All these components constitute the
covery of extrasolar planets and even some extrasolar plan-
framefortheessentialprocessofplanetformation:Anenor-
etary systems. The wealth of observational data raised the
mous growth from dust-sized particles to the final planets,
question of how a planetary system forms in general. As of
accompanied by a steady decrease of the number of parti-
writing these lines, 859 planets and 676 planetary systems
cles which contain most of the mass over many orders of
are known1. Most of these planetary systems are very dif-
magnitude. The particle number changes over many orders
ferent compared to our solar system.
of magnitude as planetary growth proceeds. There is active
researchoneachofthedifferentaspectsofplanetformation,
(cid:63) E-mail: [email protected] (PAS); butthecurrenteffortsarefarfromaunifiedmodelofplanet
[email protected] (PG); [email protected] formation (Goldreich et al. 2004; Lissauer 1993).
heidelberg.de(RS)
1 http://exoplanet.eu/catalog-all.php We address the study of this many–to–few transition
(cid:13)c 2013RAS
2 P. Amaro-Seoane, P. Glaschke & R. Spurzem
fromplanetesimalstofewprotoplanets.Thisstageisofpar- 5
N-body 100%
ticular interest, as it links the early planetesimal formation N-body 50%
Statistic 50%
to the final planet formation. Collisions still play a major 4 Statistic 100%
role in the evolution of the system, and the close interplay
between the change of the size distribution and the evolu- 1/2>
tionoftherandomvelocitiesrequiresacarefultreatmentof 2h 3
the complete size range. 2<i/
Small N−body simulations (i.e. with less than few 104 1/2> ,
particles) have been useful in exploring the basic growth 2h 2
mechanisms at the price of a modified timescale and an ar- 2<e/

tificiallyreducedparticlenumber(e.g.Kokubo1995,1996).
1
Statisticalcodesexploredthelimitoflargeparticlenumbers
in the early phases and are now tentatively applied to the
fullplanetformationprocess.Anefficientsolutionwouldbe 0
0 50 100 150 200
thecombinationofthesetwoapproachesinonehybridcode
T[yrs]
to unify the advantages of both methods.
Inthiswork,wepresentaseriesoftestsandfirstresults Figure 1. Test simulations T1a–T1c (uniform mass, see ta-
ofournewhybridcode,whichwaspresentedinPaperI.As ble1).WeshowtheresultsfromtheN−bodycalculation(100%
describedinthefirstpaper,itcombinestheNbody6code(a N−body),thestatisticalcalculation(100%Statistic)andthehy-
descendant of the widespread N−body family (see Aarseth bridcalculation(50%Statisticreferstothestatisticalcomponent,
whereas50%N−bodyistheN−bodypart).Theredcurveisthe
1999; Spurzem 1999; Aarseth 2003) with a new statistical
eccentricitydata,andthegreencurvetheinclination.
codewhichusesrecentworksonthestatisticaldescriptionof
planetesimal systems. The new hybrid code includes a con-
sistent modelling of the velocity distribution and the mass
spectrum over the whole range of relevant sizes, which al-
lows us to apply a detailed collision model rather than the
perfect-merger assumption used in previous N−body simu-
lations. We then apply this new code to follow the forma- Weuseahomogeneousringofplanetesimalsasourfirst
tion of protoplanets out of 1–10 km sized planetesimals. In test. The reason is that we can analyse the evolution with
section 2 we present a series of tests that check for the ru- three different setups – a pure N−body calculation, a pure
bustness of the code. In section 3.1 we explain the initial statistical calculation and a mixed hybrid calculation. All
conditions we use for our numerical experiments, which we threeapproachesshouldinprinciplereproducethesamere-
show in section 3.2. In section 4 we derive a useful relation sult.HenceweprepareasmallN−bodytestrun(T1a)and
that allows us to introduce an universal upper limit of the let the system evolve (see figure 1). As a second test run,
protoplanetarymassasafunctionofthedistancetothehost we shift one half of the bodies to the statistical model and
star.Finally,insection5wediscussourprogressandresults conducttheintegrationagain(T1b).Whilethisusageofthe
and potential future applications. hybrid code is somewhat artificial, it provides an excellent
setup to examine the interplay between N−body and sta-
tistical part, since neither component dominates the result.
2 VALIDATING THE CODE Finally, we run a complete statistical calculation (T1c).
We can see that all different approaches are in good
Thenewhybridcoderequirestheimplementationofrather agreement.AlthoughtheaccordancebetweenN−bodyand
different methods within a single framework. We have here statisticalcalculationisnotanewfinding–itmerelyshows
two possible sources of problems. First, the method is new that the stirring terms provide a proper description of a
andthereforeitmustbecarefullyassessedwithotherwork; planetesimal system (this was already shown by Ohtsuki
ontheotherhand,theimplementationmustalsobechecked etal.(2002))–wedeemthetesttobenecessarytodemon-
meticulously, since it combines two rather different ap- strate that the agreement holds in our approach and, in
proaches. We hence present in this section a number of particular, that the accuracy in the integration of the sta-
tests to check all code components, namely the evolution tistical model is robust. A more stringent test is posed by
of the velocity dispersion, the accuracy of the solver of the the hybrid run, which proves that the pseudo-force method
coagulation equation, the proper joining of statistical and linksbothcodecomponentsinaconsistentwaywithoutspu-
N−body component and an overall comparison of statisti- rious energy transfer. Inthis respect, figure1 includes both
cal, N−body and hybrid calculations. Table 1 summarises componentsofthehybridcalculationseparately,butthedif-
the selected test runs with the respective initial conditions. ference is so small that they are hardly distinguishable.
We also run a second test runs that follow the same
approach but with a bimodal mass distribution, with the
2.1 Energy Balance
sametotalmassinbothcomponents.Thefirstcase,T2a,is
The first test run is dedicated to a careful check of the in- a pure N−body calculation, whereas the second one treats
terplaybetweenstatisticalcomponentandN−bodycompo- the smaller particles with the statistical model. This test is
nentwithrespecttotheevolutionofthevelocitydispersion. particularlyinterestingbecauseitisclosetotherealpurpose
We therefore exclude from this first test collisions and ac- of our hybrid code. In figure 2 we can see that there is a
cretion. satisfactory agreement between the two test runs.
(cid:13)c 2013RAS,MNRAS000,1–15
Hybrid methods: Formation of protoplanetary systems 3
No. Σ ∆a N N e2/h2 i2/h2 m Type
rad
T1a 1.1251×10−6 0.02 1000 – 0.04 0.01 1.41×10−10 N−body
T1b 1.1251×10−6 0.02 500 10 0.04 0.01 1.41×10−10 Hybrid
T1c 1.1251×10−6 0.02 – 10 0.04 0.01 1.41×10−10 Statistic
T2a 0.5626×10−6 0.08 800 – 4 1 5×10−10 N−body
0.5626×10−6 200 – 4 1 2×10−9
T2b 0.5626×10−6 0.08 – 10 4 1 5×10−10 Hybrid
0.5626×10−6 200 – 4 1 2×10−9
T3 Safronov – – – – – – Statistic
T4a 1.1251×10−6 0.02 10.000 – 4 1 1.41×10−11 N−body
T4b 1.1251×10−6 0.02 – 10 4 1 1.41×10−11 Hybrid
T4c 1.1251×10−6 0.02 – 10 4 1 1.41×10−11 Statistic
T5 1.8789×10−6 – – – 620 155 2.4×10−15 Statistic
Table 1. Parameters of all test simulations (and hence preceded by a “T”). The transition mass in T4b is mtrans =3.1×10−10 Only
simulationsT3,T4a–T4candT5includecollisions.AllvaluesuseMc=G=r0=1.
0.01 Safronov Test Case
N-body 100% m
N-body 50% m 2 0.6
N-body 100% m21 Analyτ=ti0c
0.008 Statistic 50% m1 τ=2
0.5 τ=4
τ=6
τ=8
1/2> 0.006 0.4
21/22e> , <i 0.004 2mn/n0 0.3
<
0.2
0.002
0.1
0 0
0 500 1000 1500 2000 2500 3000 0 5 10 15 20 25 30 35 40 45
T[yrs] Bin Number
Figure2. ComparisonoftheN−bodycalculationT2awiththe Figure3. Testofthesolutionofthecoagulationequation(T3).
hybridcalculationT2b.Thecodingisthesameasinfigure1. TheanalyticalsolutionispresentedinPaperI.
2.2 Coagulation Equation 2.3 Testing the complete code
In thissection we verify the numerical solutionof the coag- Themostrobusttestofourhybridcode(orthestand-alone
ulationequationbyrunningacomparisonwiththeanalytic statisticalcode)isacomparisonwithapureN−bodysimu-
solution of the Safronov problem, as presented in Paper I. lationwiththesameinitialconditions.Whilealargeparticle
The collisional cross-section is assumed to be propor- numberisdesirabletocoveralargerangeinmasses,weare
tionaltothesumofthemassesofthecollidingbodies.Thus, limited in the number of particles to be used in the direct
the coagulation kernel is known and an additional integra- N−body techniques to a few 104.
tion of the velocity dispersions is not necessary. Figure 3 Wethereforechooseasingle-masssystemwithinitially
summarises the numerical solution, simulation T3, of the 10,000 particles. We enlarge the radii of the planetesimal
Safronov test. by a factor f =5, which speeds up the calculation without
Themassbinsarespacedbyafactorδ=2.Whilesome modifying the growth mode. The transition mass is twenty
slight differences emerge near the maximum of the density times larger than the initial planetesimal mass, keeping the
distribution, the overall shape is well conserved throughout particlenumbercoveredbythestatisticalcomponentlarger
the integration. This proves that a spacing within a factor than a few thousands.
two still guarantees a reliable solution of the coagulation We compare a full N−body run with a hybrid calcu-
equation without a modified timescale for the growth. lation and a pure statistical calculation. Though the stand-
(cid:13)c 2013RAS,MNRAS000,1–15
4 P. Amaro-Seoane, P. Glaschke & R. Spurzem
Nbody Nbody
Σ/Σ* Tr
1E-06 1E-04
1E-05
1E-07 1E-06
1E-07
1E-08
1E-08
1E-09 1E-09
20000 20000
15000 15000
1E-11 10000T[yrs] 1E-11 10000T[yrs]
1Em-1/0M* 1E-9 1E-8 1E-7 5000 1Em-/1M0* 1E-9 1E-8 1E-7 5000
Figure 4. SurfacedensityandradialvelocitydispersionoftheN−bodymodel(T4a).
Hybrid Hybrid
Σ/Σ* Tr
1E-06 1E-04
1E-05
1E-07 1E-06
1E-07
1E-08
1E-08
1E-09 1E-09
20000 20000
15000 15000
1E-11 1E-10 1E-9 500010000T[yrs] 1E-11 1E-10 1E-9 500010000T[yrs]
m/M* 1E-8 1E-7 m/M* 1E-8 1E-7
Figure 5. Surfacedensityandradialvelocitydispersionofthehybridmodel(T4b).
Statistic Statistic
Σ/Σ* Tr
1E-06 1E-04
1E-05
1E-07 1E-06
1E-07
1E-08
1E-08
1E-09 1E-09
20000 20000
15000 15000
1E-11 10000T[yrs] 1E-11 10000T[yrs]
1Em-1/0M* 1E-9 1E-8 1E-7 5000 1Em-1/M0* 1E-9 1E-8 1E-7 5000
Figure 6. Surfacedensityandradialvelocitydispersionofthestatisticalmodel(T4c).
(cid:13)c 2013RAS,MNRAS000,1–15
Hybrid methods: Formation of protoplanetary systems 5
T = 20,000 yrs 1e+09 1000
A: 500 yr
Nbody 1e+08 B: 10000 yr
C: 25000 yr
Hybrid 1e+07 D: 50000 yr E
Statistic E: 100000 yr
1000 1e+06 F: 150000 yr 100 D
F
m) 100 N(>m) 1 1000000000 σ [m/s]h 10 B C
> 1000
N( E A
100
F
10 A B D 1
10
1 C
1e+20 1e+22 1e+24 1e+26 1e+20 1e+22 1e+24 1e+26
m[g] m[g]
1
1e-11 1e-10 1e-09 1e-08 1e-07 Figure10. ComparativecalculationT5whichadoptstheinitial
m/M* conditionsofInabaetal.2001(theirfigure9,bottom).
Figure 7. Cumulativesizedistributionofthecomparativeruns
T4a–T4catT=20,000yr.
alonestatisticalcalculationincludesthepropertreatmentof
the runaway bodies via the gravitational range method, if
T = 20,000 yrs onlyfewparticlesresideinonemassbinwedonottakeinto
accountsuppressionofself-accretionandself-stirring.While
Nbody
Hybrid thehybridapproachdescribesthisregimeinmuchmorede-
Statistic
tail, we include the full statistical calculation nevertheless
for completeness.
0.01
Infigures4–6wehaveanoverviewofthetimeevolution
21/2e> owfhtohlee ssyysstteemm., wAhllercealaclullaqtuioanntsitsieeesmaretoinatgergereatreadthoevrerwtehlle,
<
although the statistical noise in the N−body calculation
andthehybridcalculationisquitestrongduetotheparticle
number. Runaway growth leads to the fast formation of a
few protoplanets on a timescale of a few thousand years,
with a good agreement of the fast initial growth phase in
allthreetestruns.Theboundarybetweensmoothevolution
1e-11 1e-10 1e-09 1e-08 1e-07
m/M andnoisydatamarksthelocationofthetransitionmassin
*
the hybrid calculation.
Figure 8. Mean square eccentricities of the comparative runs
We compare the size distribution and the velocity dis-
T4a–T4c at T=20,000 yr. Error bars indicate the spread due to
persion at the end of the integration, which is 20,000 yr, in
therayleighdistributionoftheeccentricity.
more detail in figures 7–9. Both the N−body data and the
hybrid data are projected on to the same grid as the full
T = 20,000 yrs statistical calculation to allow a convenient comparison.
0.01
Nbody The agreement of the size distribution N(> m) is ex-
Hybrid cellent; the small deviations are within the statistical error.
Statistic
We note that the strong variations in the size distributions
of figures 4–6 are located at the high mass end, where only
few particles dominate the surface density. In addition, the
1/2> growth in the statistical model seems to be faster than the
2<i N−body reference calculation. However, the density at the
highest masses refers to less than one particle. As we noted
before, this is due to the poor treatment of the few–body
limit.
The comparison of the velocity dispersions yields good
results,inparticularinthelow–massregime,wherethesta-
0.001
1e-11 1e-10 1e-09 1e-08 1e-07 tistical error is small. The high mass regime does not only
m/M suffer from bad statistics, but also from a pronounced time
*
variability, as we can see by comparing the fluctuations in
Figure 9. The same as figure 8 for the inclination. The strong
figures 4–6. Taking these variations into account, all three
deviationatm=3×10−9 isduetoasingleparticle.
calculationsareingoodagreement.Asbefore,thedeviation
at m=3×10−9 is due to a single particle.
2.4 The statistical code
Inaba et al. (2001) presented a high accuracy statistical
(cid:13)c 2013RAS,MNRAS000,1–15
6 P. Amaro-Seoane, P. Glaschke & R. Spurzem
code. In this section we run our last test calculation by
comparing with this work, in particular with their figure ηDisc 0.01
9, bottom. They included in their code approximately the ηreg 0.002
same physics and interpolation formulae, with only minor
differences to our approach. ηirr 0.001
While our approach allows us to set a spacing up to Rmin 0.95 AU
δ = 2, their solution of the coagulation requires a smaller
Rmax 1.05 AU
spacing,ofδ=1.1,toguaranteeareliablesolution.Thefew-
bodylimitishandledproperly,withanadditionaltreatment m¯ 3×1018 g
oftheprotoplanetsviathegravitational-rangeapproach.In
ρ 2.7 g/cm3
figure10weshowourcomparison,simulationT5,withtheir
runs. Again, we find a good agreement but for minor devi- δ 2
ations.Thesearelikelyrelatedtothedifferentimplementa-
∆vg 60 m/s
tion of the collisional probability.
Table2.Generalparameterscommontoallsimulationslistedin
table3.
3 SIMULATIONS: PROTOPLANETARY
GROWTH
3.1 Initial conditions
r =1 M =1 G=1 (5)
c c
Weapplyourhybridcodenowtoawelldefinedinitialsetup
of a planetesimal disc. All simulations use a homogenous In table 2 we summarise the main parameters of the
ring of planetesimals extending from an inner boundary simulations, fixed to the same values for all of them.
R to an outer boundary R . Since radial migration
min max
is not included, all planetesimals are bounded to this vol-
3.2 Main objectives of the analysis
umethroughoutthesimulation.Thecentralstarhasamass
of one solar mass. Each simulation starts with no N−body Theschemewehavedevelopedisinprinciplereadytosolve
particles, so that we need only to specify the setup for the the complete planetesimal problem, at least concerning the
statistical part of the calculation. The differential surface large range of sizes. However, in practise we are limited by
density as a function of mass is the computational power. A small ring with a width of 0.1
AU centred at 1 AU with a moderate size for the lower
dΣ m cut-off requires some days of integration, with the largest
=Σ exp (−m/m¯) (1)
dm 0m¯2 fraction of time spentin the statistical model. While we fo-
Var(m)=m¯2, (2) cus on these initial conditions for our simulations, we also
present some more refined models that required larger cal-
where Σ0 is the total surface density and m¯ is the mean culations. We adapt a surface density Σ=10 g/cm2 in the
mass. Equation 1 provides a smooth variation over a few simulations,whichcanbeenvisagedasanominalvalueused
massbins,whichavoidsnumericalproblemsatthebeginning in the related literature. In the remaining of this work, we
ofthesimulation.Theinitialvelocitydispersionisrelatedto focus on the following aspects of protoplanetary growth:
the mean escape velocity v of the initial size distribution
∞
(i) Different collision models: This represents a fun-
defined by Eq. 1
damental uncertainty, since the impact physics of planetes-
imals is not well established yet. In order to do realistic
1
v2 =T +T +T , (3) modelsofplanetesimalcollisionsweneedtounderstandthe
100 ∞ r φ z
internal structure of the bodies taking part in the colli-
with the ratio of the velocity dispersions sion. Planetesimals emerge as fragile dust aggregates and
evolveintosolidbodies,sothattheirinternalstructureand
strength is time-dependent.
T =4T =4T (4)
r φ z (ii) Spatial (radial) density structure (e.g. gap for-
We adopt a rather small initial velocity dispersion to mation) This is related to the slowly evolving inhomo-
avoid strong spurious fragmentation due to an overestima- geneities introduced by the growing protoplanets. It has
tion of the velocity dispersion. Furthermore, strong relax- beenarguedthatgapopeningintheplanetesimaldisccould
ation in the initial phase of the calculation quickly estab- stop the accretion well before the isolation mass is reached
lishesanequilibriumvelocitydispersion.Thetimestepcon- (Rafikov2001).Ourhybridcodeincludesanaccuratetreat-
trolparametersarechosensuchthattheenergyerror∆E/E mentofspatialstructuring,sothatweareinthepositionof
of the N−body component remains always smaller than ascertainingtheroleofgapformationintheprotoplanetary
10−8 throughoutthesimulation.Likewise,ourchoiceofthe growth process.
parametersofthestatisticalcomponentassuresthatthesta- (iii) Resolution effects hinge on the limitation of com-
tisticalmodelissolvedaccuratelyandremainsstable,asin- putingpower.Sincethesolutionofthecoagulationequation
dicated by the set of comparative runs. All runs simulate scales with the third power of the number of grid cells, the
only a narrow ring centred at a distance r and we choose choiceofarealisticcut-offmassmaybeprohibitivelyexpen-
c
the following units: sive.
(cid:13)c 2013RAS,MNRAS000,1–15
Hybrid methods: Formation of protoplanetary systems 7
Code Strength NM NR Σ[g/cm2] mmin/Mc mtrans/Mc ρg[g/cm3]
S1FB B&A1999 24 50 10 3.48×10−18 3.89×10−11 10−9
S2FH H&H1990 24 50 10 3.48×10−18 3.89×10−11 10−9
S3FN PerfectMerger 24 50 10 3.48×10−18 3.89×10−11 10−9
S4FBN B&A1999 24 50 10 3.48×10−18 3.89×10−11 0
S5FBL B&A1999 24 5 10 3.48×10−18 3.89×10−11 10−9
S6FBH B&A1999 24 100 10 3.48×10−18 3.89×10−11 10−9
S7FB2 B&A1999 40 50 10 5.31×10−23 3.89×10−11 10−9
S8 S2 B&A1999 15 50 2 3.48×10−18 4.87×10−12 2×10−10
S9 S100 B&A1999 27 50 100 3.48×10−18 3.11×10−10 10−8
Table3.Completelistofallsimulations(thenamesofthemodelsarethereforeprecededwithan“S”).Thefirstgroupexaminesdifferent
collisionalmodels,thesecondgroupresumesthenominalsimulationS1FBwithdifferentresolutionsandthethirdgroupexploresdifferent
surfacedensities.
We present the results of the simulations of the four
No. A B C D E F
different approaches in a figure with four panels: Figure 11
T[yr] 0 1,000 10,000 20,000 50,000 100,000 shows model S3FN, figure 12 model S2FH, figure 13 model
S1FBandfigure14modelS4FBN.Inthesefigureswedepict
intheupper,leftpanelthecumulativesizedistributionN(>
Table4.IntegrationtimesfromtheevolutionarystagesAtoF.
m), which allows us to see the distribution of particles as a
function of the range of masses at different moments of the
(iv) Different surface densities: To address this, we integration(T =0,103,104,2×104,5×104 and 105yrs,and
conduct a small set of different surface densities with our we follow in the figures the notation of table 4).
reference fragmentation model (Benz & Asphaug 1999,im- On the upper right panel we display the evolution of
pact strength, referred to as B&A 1999 hereafter). the surface density per bin Σ . Since we are using a loga-
∆
rithmically equal spacing of the mass grid, Σ is related to
In table 3 we summarise the various parameters of our ∆
the differential surface density
simulations. In the following subsections we discuss each
simulation in more detail.
We project the N−body data on to an extended mass 2 ∂Σ
Σ ≈ , (6)
gridderivedfromthestatisticalmodeltogenerateaunified ∆ 3∂ln(m)
representationofahybridrun.Thisissobecausethehybrid
codeusesbothastatisticalrepresentationandN−bodydata where we assume δ=2.
to integrate the planetesimal disc. The lower left and right panels show the radial (T )
r
and vertical (T ) velocity dispersion of the system at the
z
different times of table 4.
3.3 Fragmentation models
Oneconclusionthatwecanderiveimmediatelyinview
Thetreatmentofcollisionsisakeyelementinanysimulation of these figures is that in spite of the rather different ini-
of planetesimal growth. In this section we explore different tialapproachesofthemodels,theirtimeevolutionisrather
collisionalmodelswithfourdifferentsetups,soastoanalyse similar. The runaway growth sets in after some 104 years,
its influence on the final results. i.e. around stage C in the figures. This is relatively easy to
Theperfectmergerassumption(S3FNhereafter,seeta- see because of the pronounced peak at the high mass end.
ble 3) is the simplest approach for mutual collisions among The onset of runaway growth roughly coincides with the
smaller planetesimals. This rather simplistic approach can creation of the first N−body particles. Contrary to previ-
beenvisagedasawaytoderiveanupperlimitforthegrowth ous work done with statistical calculations (Wetherill 1989,
speed in our models. The second and third model use our 1993), we find in our models no gap in the size distribu-
detailedcollisionalmodel(seesection“Collisionalandfrag- tion, but a smooth transition from the slowly growing field
mentation model” of Paper I) with the B&A 1999 impact planetesimals (peak around 1019 g) to the rapidly growing
strength (S1FB) and the approach of Housen & Holsapple protoplanets.
(1990) for the impact strength (S2FH from now onwards). The initiation of runaway growth is associated with a
These two approaches roughly delimit the range of possible qualitative change in the velocity dispersion. While the ini-
values (see e.g. the overview in Benz & Asphaug 1999). tialchoiceofthevelocitydispersionquicklyrelaxestoacon-
The fourth model (S4FBN) assumes that the gaseous stantvalueatsmallersizes(transitionstageA→B),dynami-
disc has dispersed early, so that we have a gas-free system. calfrictionestablishesenergyequipartitionamongthelarger
Thismodelprovidesuswithadifferentevolutionfortheran- masses (see e.g. Khalisi et al. 2007,in the context of stellar
domvelocities,whichleadstoadifferentroleofthecollisions dynamics). The turnover point between these two regimes
All other simulations neglect the dispersion of the gaseous refers to a balance between the stirring due to larger bod-
disc,sincethesimulationtimeisstillshortcomparedtothe ies and damping due to encounters with smaller planetesi-
disc lifetime. mals (Rafikov 2003). In addition, the smaller planetesimals
(cid:13)c 2013RAS,MNRAS000,1–15
8 P. Amaro-Seoane, P. Glaschke & R. Spurzem
10 10
1e+09 1e+09
1e+08 1e+08
1e+07 1e+07
1e+06 1 1e+06 1
N(>m) 1 0100000000 2Σ[g/cm] ∆ EF N(>m) 1 1000000000 2Σ[g/cm] ∆ E F
1000 0.1 1000 0.1
100 A B C AB 100 A B C AB
E F
10 D F C D 10 DE C D
1 0.01 1 0.01
1e+16 1e+18 1e+20 1e+22 1e+24 1e+26 1e+16 1e+18 1e+20 1e+22 1e+24 1e+26 1e+16 1e+18 1e+20 1e+22 1e+24 1e+26 1e+16 1e+18 1e+20 1e+22 1e+24 1e+26
m[g] m[g] m[g] m[g]
100 100 100 100
F F F F
E E
D E D E
1/2T[m/s]r 10 B C 1/2T[m/s]z 10 B C D 1/2T[m/s]r 10 B C 1/2T[m/s]z 10 B C D
1 1 1 1
A A
A A
0.1 0.1
0.1 0.1 1e+16 1e+18 1e+20 1e+22 1e+24 1e+26 1e+16 1e+18 1e+20 1e+22 1e+24 1e+26
1e+16 1e+18 1e+20 1e+22 1e+24 1e+26 1e+16 1e+18 1e+20 1e+22 1e+24 1e+26
m[g] m[g]
m[g] m[g]
Figure 12. SummaryofsimulationS2FH,whichusestheH&H
Figure11. SummaryofsimulationS3FN,whichassumesperfect
1990strength.Table4givesthetimecodingofthelabelsA–F.
mergers.Table4givesthetimecodingofthelabelsA–F.
10
1e+09
are subjected to damping by the gaseous disc, which sig-
nificantly reduces the velocity dispersion at smaller sizes. 1e+08
Hence this damping is absent in the gas-free case, which 1e+07
can be seen by comparing the flat distribution of S4FBN, 1e+06 1
figurWe 1e4embopthtoamsis,ewtihthattahlelositmheurlamtioodneslsd.o not generate any N(>m) 1 1000000000 2Σ[g/cm] ∆ F
artifacts which could be attributed to an improper joining
1000 0.1
of the statistical and the N−body component. Some non-
100 A B C AB E
smooth structure is visible at the high mass end (i.e., it is F
related to data from the N−body component), but these 10 DE C D
variationsdonotexceedthefluctuationsthatwecanexpect 1 0.01
1e+16 1e+18 1e+20 1e+22 1e+24 1e+26 1e+16 1e+18 1e+20 1e+22 1e+24 1e+26
from small number statistics. m[g] m[g]
All simulations with destructive collisions exhibit the
evolutionofafragmenttail.Theexpectedequilibriumslope
100 100
isroughlyk≈2(seesection“Collisionalcascades”ofPaper F F
I),whichreferstoasteepsizedistributionandaratherflat
E
density distribution: E
N(>m)∝m−1 1/2T[m/s]r 10 B CD 1/2T[m/s]z 10 B C D
1 1
Σ ≈constant (7)
∆
A A
Simulation S1FB (B&A 1999 strength, figure 13) and
S4FBN (gas–free, figure 14) show a clear plateau in the 0.1 0.1
densitydistributionaround0.1,inaccordancewiththepre- 1e+16 1e+18 1e+20 1e+22 1e+24 1e+26 1e+16 1e+18 1e+20 1e+22 1e+24 1e+26
m[g] m[g]
vious estimate, Eq. 7. In contrast, simulation S2FH (H&H
1990 strength, figure 12) evolves a second maximum at the Figure 13. SummaryofsimulationS1FB,whichusestheB&A
lower boundary of the mass grid. Although this structure 1999strength.Table4givesthetimecodingofthelabelsA–F.
is partly due to the lower grid boundary, the main cause is
the reduced H&H 1990 impact strength at sizes of a few 10
kilometres (as compared to the B&A 1999 strength), which
leads to the quick destruction of the remaining field plan-
etesimals at masses around 1018 g.
(cid:13)c 2013RAS,MNRAS000,1–15
Hybrid methods: Formation of protoplanetary systems 9
10 2.5e+25
1e+09 B&A 1999
H&H 1990
1e+08 Perfect Merger
1e+07 2e+25 GLoaws-eFrr eCeut-Off
1
1e+06
N(>m) 1 01 001000000000 2Σ[g/cm] ∆ 0.1 E F m[g]max 1.5e+25
1e+25
100 A B C AB
F
10 DE C
D
5e+24
1 0.01
1e+16 1e+18 1e+20 1e+22 1e+24 1e+26 1e+16 1e+18 1e+20 1e+22 1e+24 1e+26
m[g] m[g]
0
0 20000 40000 60000 80000 100000
100 100
T[yr]
F F
E
Figure15. Largestbodyinthesimulationasafunctionoftime
D E forthedifferentcollisionmodelsS1FB(B&A1999),S2FH(H&H
1/2T[m/s]r 10 C 1/2T[m/s]z 10 C D 1w9e9a0l)s,oSi3nFcNlud(peesrifmecutlamtieorngeSr7)FaBnd2Sw4itFhBaNlo(gwaesr-fcrueet-)o.ffInmaadsdsi.tion,
B B
1 1
A A 0.12
B&A 1999
H&H 1990
Perfect Merger
0.1 0.1 0.1 Gas-Free
1e+16 1e+18 1e+20 1e+22 1e+24 1e+26 1e+16 1e+18 1e+20 1e+22 1e+24 1e+26 Lower Cut-Off
m[g] m[g]
0.08
Figure14. SummaryofsimulationS4FBN,whichusestheB&A Mtot
1o9f9t9hestlraebnegltshAa–nFd.agas–freesystem.Table4givesthetimecoding /Nbody 0.06
M
0.04
Theoverallagreementofthedifferentsimulationsisre-
0.02
flectedbythegrowthofthelargestmassinthesystem,aswe
depictinfigure15.Upto2×104 years,allsimulationsagree
0
well. Later on simulation S3FN (which follows the approxi-
0 20000 40000 60000 80000 100000
mationofperfectmergers)exhibitsthelargestgrowthrate, T[yr]
as one would naturally expect. Although simulation S1FB
(which uses the B&A 1999 strength prescription) seems to Figure 16. The same as figure 15 for the total mass in the
N−bodycomponent.
show a slower growth than simulation S2FH (which follows
the recipe of H&H 1990 for strength), this is only due to a
different sequence of major impacts. In fact, the B&A 1999
3.4 Spatial distribution
strength simulation makes possible a much faster growth,
inaccordancewiththetotalmasscontainedintheN−body Howwellthecodecantreatspatialinhomogeneitiesdepends
component,whichisdisplayedinfigure16.Thegas-freesim- onthechoiceofthespatialresolution.Wehencecomparea
ulation S4FBN exhibits the slowest growth among the four lowresolutionmodel,modelS5FBLoftable3,whichvirtu-
test cases. allyinhibitsanyspatialstructuring,withamodelthatuses
A further examination of the mass loss – which we de- ourfiducialresolution,modelS1FB,aswellaswithamodel
fine as the mass in planetesimals which cross the lower grid that has a finer resolution, S6FBH. We adjust the fiducial
boundary– reveals the cause of this different behaviour: A resolutiontothewidthoftheheatingzoneofaplanetesimal
pronounced mass loss in simulation S2FH slows down the at the transition mass.
protoplanetary growth reducing the surface density. In the In the left panel of figure 17 we have the spatial struc-
gas-free case, the accretion rate is mainly reduced because ture at T = 30,000 yr, i.e. shortly after stage D (nominal
of a larger velocity dispersion, although we can still notice model S1FB). While the protoplanets are already massive
some enhanced mass loss by comparing the lower panels of enough after a few 104 years (stage C) to open gaps in the
figures 13 and 14. planetesimal component, there is only a weak correlation
We find no accelerated growth due to the inclusion between the radial structures and the location of the most
of fragmentation events, contrary to the work of Wetherill massiveprotoplanets.Acloserexaminationofthetimeevo-
(1989).Wefindthatalowerimpactstrengthortheabsence lution of the radial structure reveals that most features are
ofgasdampingslowsdownthegrowthbyanincreasedmass “fossils” from the first emerged N−body particles, which
loss.ThetotalmassintheN−bodycomponentisstillsmall areslowlywashedoutbythediffusionofthefieldplanetesi-
at the end of the simulations, of about ≈ 10% of the total mals.Intherightpanelofthesamefigure17weconfirmthe
mass, as shown in table 5 of next section. furthersmoothingoftheradialfeatures.Whilemajormerg-
(cid:13)c 2013RAS,MNRAS000,1–15
10 P. Amaro-Seoane, P. Glaschke & R. Spurzem
T=30,000 yr T=100,000 yr
1.3 1e-08 1
Statistic Statistic
N-Body N-Body 1e-08
1.2 0.9
1e-09
0.8
1.1 1e-09
Σ0 M*Σ0 M*
Σ/ m/Σ/ 0.7 m/
1
1e-10
1e-10
0.6
0.9
Grid Resolution: 0.5 Grid Resolution:
0.8 1e-11 1e-11
0.94 0.96 0.98 1 1.02 1.04 1.06 0.94 0.96 0.98 1 1.02 1.04 1.06
r r
Figure 17. Left y-axis and solid, red curve of the left panel: Radial density structure of the statistical component of model S1FB at
T =3×104 yr.Righty-axisandbluedotsoftheleftpanel:SemimajoraxisandmassesoftheN−bodyparticlesinthesimulationafter
thesameamountoftime.Theerrorbarsare10HillradiiwideandrefertotheheatingzoneofeachN−bodyparticle.Wealsodisplay
thegridresolutionasareferencepoint.IntherightpanelwedepictthesameafterT =105 yrs.
2e+25 N =5 isalreadydynamicallytoohottoallowthesystemtodevelop
R
1.8e+25 NR=50 radialstructures.Theeccentricitiesofthefieldplanetesimals
N =100
R are comparable to the width of the heating zone (compare
1.6e+25
figure13,bottom),andhenceanyplanetesimalthatisscat-
1.4e+25
tered to larger (or smaller) radii immediately encounters a
1.2e+25 neighbouring protoplanet.
[g]max 1e+25 In summary, the protoplanets (or rather their precur-
m sors) are too abundant when the system is dynamically
8e+24
cool enough, but when a group of mature protoplanets has
6e+24 evolved, the system is already too hot. Thus, we expect an
4e+24 evenlesseffectiveradialstructuringforlargersurfacedensi-
ties.Whilesystemswithalowersurfacedensitymayleadto
2e+24
theformationofgap-likestructures,theyevolveslowlythat
0 planet formation may never reach the final growth phases.
0 20000 40000 60000 80000 100000
T[yr]
Figure18. Largestbodyinthesimulationasafunctionoftime 3.5 Resolution
for the different resolutions S5FBL (NR =5), S1FB (NR =50) We can further evaluate the (minor) role of gap formation
andS6FBH(NR=100).
by comparing the growth process for the three different ra-
dial resolutions N = 5, N = 50 and N = 100. Besides
R R R
somevariationsduetoadifferentsequenceofmajorimpacts
ers among the protoplanets still lead to distinct features in
(see figure 18), all three simulations are in excellent agree-
the surface density even after a few 104 years, any further
mentwithrespecttothemasslossandthetotalmassinthe
structuring ceases at the end of the simulation.
N−body component.
The absence of any prominent gap formation (fluctua-
Accordingly,wefindnodifferencesbetweenthevarious
tions are smaller than 20 %) is related to the evolution of
fragmentationmodels(S1FB,S2FH,S3FN,S4FBN)withre-
the overall size distribution. Though the gap opening crite-
specttopossibleemerginggaps,exceptanearlierhomogeni-
rion (see section “Protoplanet growth” of Paper I, and we
sationinthegas-freecaseS4FBNduetothestrongerheat-
reproduce here the relevant equation for convenience)
ing of the smaller planetesimals.
We conduct one additional simulation, named S7FB2,

Mgap ≈  ΣMac2 (cid:16)Mmc(cid:17)1/3 if v(cid:46)ΩrHill (8) sbeoeunfidgaurrye b2y0,ainfacwtohrich105w.eArltehdouucgehththeeloswtaenrdamrdasschgoricide
Mc  ΣMac2 (cid:16)Mmc(cid:17)1/3(cid:16)ΩrvHill(cid:17)2 if v(cid:29)ΩrHill, wmhmeirne=m6ig.9ra×tio1n01w5ogulids irnemacocvoerdthaencsemwailtlehrtfhraegsmizeenrtes,gitmhee
is formally satisfied by all protoplanets during the runaway actual mass cut-off is less sharp as we estimated in section
phase,thedenseoverlappingoftheassociatedheatingzones “Collisional cascades” of Paper I. A reduced lower cut-off
(seefigure17)inhibitstheevolutionofanygap-likefeature. increases the dwell time of collisional fragments in the sys-
Astheprotoplanetsgrow,theyexertagrowinginfluenceon tem, thus increasing the mass fraction which could be ac-
the dynamics of the planetesimal system. While this domi- cretedbytheprotoplanets.Asaresult,masslossisreduced
nancecouldinprincipleenhancegapformation,thesystem by30%comparedtoourfiducialcaseS1FB,aswecanseein
(cid:13)c 2013RAS,MNRAS000,1–15