Hybrid methods in planetesimal dynamics (I) : Description of a new composite algorithm
MMon. Not. R. Astron. Soc. , 1–42 (2011) Printed 29 October 2018 (MN L A TEX style file v2.2)
Hybrid methods in planetesimal dynamics (I) :Description of a new composite algorithm
P. Glaschke , P. Amaro-Seoane (cid:63) , , & R. Spurzem , , Astronomisches Rechen-Institut, M¨onchhofstraße 12-14, Zentrum f¨ur Astronomie, Universit¨at Heidelberg, Germany Max Planck Institut f¨ur Gravitationsphysik (Albert-Einstein-Institut), D-14476 Potsdam, Germany National Astronomical Observatories of China, Chinese Academy of Sciences, 20A Datun Lu, Chaoyang District, 100012, Beijing, China Institut de Ci`encies de l’Espai (CSIC-IEEC), Campus UAB, Torre C-5, parells, na planta, ES-08193, Bellaterra, Barcelona, Spain Kavli Institute for Astronomy and Astrophysics, Peking University, China draft 29 October 2018
ABSTRACT
The formation and evolution of protoplanetary systems, the breeding grounds of planetformation, is a complex dynamical problem that involves many orders of magnitudes.To serve this purpose, we present a new hybrid algorithm that combines a Fokker-Planck approach with the advantages of a pure direct-summation N − body scheme,with a very accurate integration of close encounters for the orbital evolution of thelarger bodies with a statistical model, envisaged to simulate the very large number ofsmaller planetesimals in the disc. Direct-summation techniques have been historicallydevelopped for the study of dense stellar systems such as open and globular clus-ters and, within some limits imposed by the number of stars, of galactic nuclei. Thenumber of modifications to adapt direct-summation N − body techniques to planetarydynamics is not undemanding and requires modifications. These include the way closeencounters are treated, as well as the selection process for the “neighbour radius” of theparticles and the extended Hermite scheme, used for the very first time in this work,as well as the implementation of a central potential, drag forces and the adjustment ofthe regularisation treatment. For the statistical description of the planetesimal disc weemploy a Fokker-Planck approach. We include dynamical friction, high- and low-speedencounters, the role of distant encounters as well as gas and collisional damping andthen generalise the model to inhomogenous discs. We then describe the combinationof the two techniques to address the whole problem of planetesimal dynamics in arealistic way via a transition mass to integrate the evolution of the particles accordingto their masses. Key words: protoplanetary discs, planets and satellites: dynamical evolution andstability, methods: numerical, methods: N-body, methods: statistical
The formation of a planetary system is closely related tothe formation of the host star itself. Cool molecular cloudscollapse and fragment into smaller substructures which arethe seeds for subsequent star formation. Angular momen-tum conservation in the forming clumps forces the infallingmatter into a disc-like structure. The subsequent viscousevolution of the disc leads to a transport of angular momen-tum which channels gas to the protostar in the centre. Theseprotoplanetary discs are the birth place of planets (for a de-tailed review see Armitage 2010). Embedded dust grains are (cid:63)
E-mail: [email protected] (PG); [email protected] (PAS, corresponding author);[email protected] (RS) the seed for the enormous growth to bodies of planetary size.The first hint to the structure of protoplanetary discs hasbeen provided by our own solar system. Through “smearingout” all planets and adding the missing fraction of volatileelements, one can estimate the structure and mass of theprotoplanetary disc. Since the efficiency of planet formationis unknown, this yields only a lower limit – the minimummass solar nebula (Hayashi 1981). The inferred surface den-sity decreases with the distance from the sun as ∝ r − / .Further insight has been obtained by the detection of aninfrared excess of some stars, i. e. additional infrared radia-tion that does not originate from the star but an unresolveddisc. Advancements in observation led in the mid-90s to thedirect imaging of nearby star forming regions which opened anew flourishing field in astronomy (see O’dell et al. 1993,for c (cid:13) a r X i v : . [ a s t r o - ph . E P ] M a y P. Glaschke, P. Amaro-Seoane & R. Spurzem an example with Hubble Space Telescope images)). Sincethen a big amount of observations of protoplanetary discsacross the electromagnetic range has been obtained and in-terpreted, both space-based (ISO, Spitzer and now Herschel)and in the ground (VLT, VLT-I, Subaru, Keck). We refer thereader to “The Star Formation Newsletter” URL for anoverview of the recent papers in star- and planet formationand the book review “Protostars and Planets V” (Reipurthet al. 2007).Protoplanetary discs masses cover a range from 10 − to 0 . M (cid:12) with a peak around 0 . M (cid:12) (data from Tau-rus/Ophiuchus Beckwith 1996), in accordance with mass es-timates deduced from the minimum mass solar nebula. Sincethe disc lifetime can not be measured directly, it is derivedindirectly from the age of young, naked (i. e. discless) starswhich sets an upper limit. Pre-main sequence evolutionarytracks are used to gauge the stellar ages, providing lifetimesof a few 10 years. The subsequent evolution of the discproceeds in several stages.Two different scenarios have been proposed to explainthe formation of kilometre-sized planetesimals.(i) One process is the gravitational instability of the dustcomponent in a protoplanetary disc that leads to the directformation of larger bodies. Goldreich & Ward (1973) pro-pose that an initial growth phase of dust grains leads to athin dust disc that undergoes a gravitational collapse. Asthe dense dust layer decouples from the gas, it rotates withthe local Keplerian velocity, whereas the gas component ro-tates slower as it is partially pressure supported. This givesrise to a velocity shear at the boundary, which may exciteturbulence through the Kelvin-Helmholtz-instability. Sincethe motion of small dust grains in the boundary layer is cou-pled to the gas, the turbulent velocity field could suppressthe formation of a stratified dust layer, which is a necessaryprerequisite for the gravitational instability (Weidenschilling1977).(ii) The collisional agglomeration of dust particles is anopposed formation mechanism. Relative velocities are dom-inated by the Brownian motion in the early phases of thegrowth process. This mechanism becomes increasingly in-efficient with growing mass, but successively the particlesdecouple from the gas and settle to the midplane – a pro-cess that yields even larger velocities with increasing mass.The sedimentation initiates a growth mode that is similarto the processes in rain clouds: Larger grains drop faster,thus accreting smaller grains on their way to the midplane.Turbulence may modify this basic growth scenario by forc-ing the dust grains in a convection-like motion. Dust grainsstill grow during the settling process, but the turbulent ve-locity field could mix up dust from the midplane, and a newcycle begins. Each cycle adds a new layer to the dust grains– a mechanism that also operates in hail clouds – until thegrains are large enough to decouple from the turbulent mo-tion. Again, turbulence plays an important role in determin-ing the growth mode and the relative velocities. While therelative velocities are high enough to allow for a fast growth,it is not clear a priory that collisions are sticky enough to allow for a net growth. High speed encounters lead to frag-mentation, which counteracts agglomeration (e.g. Blum &Wurm 2000,and references therein). An important bottle-neck in the agglomeration process is the fast orbital decayof 1 m sized boulders. Their orbital lifetime is as short as100 years, and a quick increase in size – at least over oneorder of magnitude – is needed to reduce the radial driftsignificantly.To overcome the difficulties associated with each of thesescenarios, modifications have been proposed. Magneto-hydrodynamic simulations include electro-magnetic interac-tions in hydrodynamical calculations. See the reviews of Bal-bus & Hawley 1998 and Balbus. 2003. The MHD simulationsby Johansen et al. (2006) show that trapping of larger parti-cles in turbulent vortices helps in increasing the orbital life-time, but could also trigger local instabilities that may leadto the direct formation of planetesimals (Inaba et al. 2005).Johansen et al. (2007) describe a gravoturbulence mecha-nism as a feasible pathway to planetesimal formation in ac-creting circumstellar discs.The details of agglomeration have drawn a lot of at-tention and are still under question (see Kempf et al. 1999;Paszun & Dominik 2009; Wada et al. 2009), but the succes-sive agglomeration of planetesimals is commonly accepted. The further growth of planetesimals proceeds through mu-tual collisions, where the initial phase involves a large num-ber of particles and is well described by a coagulation equa-tion (Safronov 1969). While earlier works (e.g. Nakagawaet al. 1983) focussed mainly on the evolution of the size dis-tribution, subsequent refinements of the evolution of the ran-dom velocities showed that it is important to evolve the sizedistribution and the velocity dispersion in a consistent way.A fixed velocity dispersion is an oversimplification, whichchanges the growth mode and increases the growth timescaleas well (Wetherill 1989).The initial growth is quite democratic. All planetesimalsgrow roughly at the same rate and the maximum of the sizedistribution is shifted gently towards larger sizes. As soonas gravitational focusing and dynamical friction become im-portant, the growth mode changes to a qualitatively differ-ent mechanism. Efficient gravitational focusing leads to agrowth timescale (which we denote as M/ ˙ M ) that decreaseswith mass. Hence larger bodies grow faster than the smallerplanetesimals, a trend that is further supported by energyequipartition due to planetesimal–planetesimal encounters.This dynamical friction keeps the largest bodies on nearlycircular orbit, thus the relative velocities are small and gravi-tational focusing remains efficient. Smaller planetesimals arestirred up into eccentricity orbits, which slows down theirgrowth rate compared to the largest bodies. This acceler-ated growth, denoted as runaway growth of planetesimals(see e.g Greenberg et al. 1978; Wetherill 1989, 1993), short-ens the growth timescale to a few 10 years. The term run-away growth stresses that the growth timescale of a particledecreases with mass, hence the largest body “runs away” tothe high mass end of the distribution (see Kokubo 1995).While energy equipartition increases the velocity dis-persion with decreasing mass of the planetesimals, addi- c (cid:13) , 1–42 ybrid methods: A new composite algorithm tional damping due to the gaseous disc leads to a turn-overat smaller sizes. However, higher relative velocities may leadnevertheless to destructive encounters, but these fragmen-tation events could even speed up the growth (Chambers2006; Bromley & Kenyon 2006, 2010). Since smaller bodiesare more subjected to gas drag, their velocity dispersion issmaller which allows an efficient accretion by the runawaybodies. Moreover, smaller particles damp the velocity dis-persion of the largest bodies more efficiently.As runaway proceeds, the system becomes more andmore dominated by few big bodies – the protoplanets . Dueto the dominance of few, very large bodies one can not usestatistical methods anymore to simulate the problem.Kobayashi et al. (2010) derived analytical expressionsfor the final masses of these planetary “embryos”, as theycall them, including the role of planetesimal depletion due tocollisional disruption. They conclude that the final mass inthe minimum-mass solar nebula at several AU can achieve ∼ . years. The runaway growth of large bodies (i. e. protoplanets)ceases to be efficient as soon as the protoplanets start tocontrol the velocity dispersion of the remaining planetesi-mals in their vicinity. Gravitational focusing becomes lesseffective, therefore the growth timescale increases with sizeand the growth mode changes to oligarchic growth. The pro-toplanets still grow faster than the field planetesimals , butthe masses of the protoplanets remain comparable.A combination of dynamical friction due to the fieldplanetesimals and perturbations from the neighbouring pro-toplanets conserves a separation of five to ten Hill radiibetween neighbouring bodies. Therefore only planetesimalsfrom a limited area, the feeding zone, are accreted by a givenprotoplanet. If this zone is emptied, they have reached theirfinal isolation mass (Kokubo et al. 1998; Kokubo & Ida 2000,2002).As the damping of the remaining field planetesimals isweak enough, further growth is dominated by mutual per-turbations among the protoplanets, which leads to giant im-pacts. Protoplanets beyond the “snow line” (or ice conden-sation point) can grow larger than 5–15 earth masses andinitiate the formation of giant planets. If the protoplanetarydisc is very massive in the inner planetary system, this maylead to an in-situ formation of hot Jupiters (see Bodenheimeret al. 2000).Kokubo et al. (2006) used the oligarchic growth modelof protoplanets to statistically quantify the giant impactstage (i.e. collisions of protoplanets to form planets) with N − body simulations. They found that for steeper surfacedensity profiles, large planets usually form closer to the star. The proposed three-stage scenario of planet formation cov-ers the dominant growth processes, but a major mechanism The term “field planetesimals” denotes in the following thesmooth component of smaller planetesimals. is still missing – the migration of bodies in the system. Mi-gration is a generic term that summarises a set of differentmechanisms that lead to secular radial drift of bodies (seee.g. the review of Papaloizou 2006; Armitage 2010):(i) The dissipation due to the remaining gas disc leadsto an orbital decay of the planetesimals. While this poses asevere problem for 1 m–sized objects, larger bodies drift veryslowly inward. One denotes this process as type 0 migration.(ii) Planets which are embedded in a gaseous disc launchspiral density waves at the inner and outer Lindblad reso-nances, which leads to an exchange of angular momentumwith the resonant excited waves. This type I migration leadsto a robust inward migration independent of the density pro-file. The perturbation from the planet is small, hence linearperturbation theory is in principle applicable (Ward 1986),but recent work proves that non-linear and radiative effectscan be quite important (Paardekooper et al. 2011).(iii) If the protoplanet is massive enough, it opens a gapin the gaseous disc and excites waves through tidal interac-tion with the gaseous disc. An imbalance of the exchangeof angular momentum with the inner and outer part of thedisc leads to type II migration. The strong interaction be-tween planet and disc leaves the linear regime and requiresa numerical solution of the hydrodynamic equations (Lin &Papaloizou 1979).(iv) Even with an opened gap, the planet still channelsgas between the inner and outer part of the disc. While thiscorotational flow is to some extent already present duringtype II migration, it dominates the angular momentum bal-ance in the case of type III migration due to an asymmetryin the leading and trailing part of the flow. An imbalanceallows for an efficient exchange of angular momentum withthis corotational gas stream, which gives rise to a remark-ably fast migration (Masset & Papaloizou 2003). Type IIImigration is subjected to a positive feedback: A faster mi-gration increases the asymmetry in the corotational flow,which speeds up the migration. Both inward and outwardmigration are possible.The four migration types modify the three-stage scenario indifferent ways.
Type 0 migration is most efficient for small planetesimals(i. e. 1 km or smaller). It poses a severe problem during theearly stages of planet formation, as it may induce a signifi-cant loss of solid material, accompanied by a global changein the initial surface density (Kornet et al. 2001). The im-portance of this process diminishes as planetesimal growthproceeds, but it still leads to the loss of collisional fragmentsduring the final disc clearing.
Type I is most important for protoplanets (i. e. 0.1 earthmasses or larger). It leads to an orbital decay of protoplan-ets, but this does not only imply a loss of protoplanets, butalso breaks the conditions for isolation. Migrating bodies canaccrete along their way through the disc and are thus notconstrained to a fixed feeding zone, which may increase theisolation mass.
Type II and type III mainly influences giant planets (largerthan 10 earth masses) causing an inward or outward migra-tion, depending on the angular momentum exchange bal-ance. It can explain the large number of giant planets closeto their host star (hot Jupiters) found in extrasolar systems.An important issue is the timescale of the migration process. c (cid:13) , 1–42 P. Glaschke, P. Amaro-Seoane & R. Spurzem
If the migration is too fast, virtually all planets spiral inwardand leave an empty system behind.Migration is a powerful process with the capability toreshape an entire planetary system, at least for gas giants,since terrestrial planets are likely not affected by migrationin gaseous discs (Bromley & Kenyon 2011). However, it alsorequires some “parking mechanism” which terminates mi-gration before all planets (or protoplanets) are lost. Inho-mogeneities in the gaseous disc may change the crucial mo-mentum balance of the inner and outer part of the disc, thusstalling or even reversing the drift of a planet.The migration processes end after the dissipation of thegaseous disc due to photoevaporation or star-star encounters(i. e. after a few 10 to 10 years).The formation of a planetary system is a vital pro-cess that is driven by the interplay between the differentgrowth phases and the migration of planets and protoplan-etary cores (i. e. the precursors of giant planets). While thepreceding sections only summarised the main evolutionaryprocesses, even more processes could influence the forma-tion of planetary systems. A fast accretion of giant planetsin the outer parts of a planetary system could introduce fur-ther perturbations on the inner part and may even triggerthe formation of terrestrial planets. Moreover, the stellar en-vironment in dense star clusters and multiple stellar systemsalso perturbs planet formation, which would require an evenbroader view on the problem. Last but not least, the exci-tation due to gravitational coupling to gas surface densityfluctuations plays an important role (see e.g. Ida et al. 2008).Any approach to planet formation can hardly includethis wealth of different phenomena, thus it is important tofocus on a well-defined subproblem. In this work we focuson the formation of protoplanets for the following reasons:(i) The size, growth timescale and spacing of the proto-planets is a key element in the planet formation process.(ii) The protoplanet growth is well-defined by differentgrowth modes. It starts with the already formed planetes-imals ( ≈ M , the Hill’s problem and the protoplanetgrowth from a theoretical point of view, respectively. Subsequently, in section 6 and its corresponding subsec-tions, we ascertain the direct-summation part of the generalalgorithm, as well as the required modifications and addi-tions to tailor it to the specific planetesimal problem. Insection 9 we introduce our collisional treatment, includingfragmentation, collisional cascades, migration and coagula-tion. The statistical part of the algorithm is described indetail in section 16. Finally, in section 17 we explain howto bring together both schemes into a single numerical tooland in section 18 we give our conclusions and compare ourscheme to other numerical approaches. The basis for all planet formation models is the structure ofthe protoplanetary disc. We summarise the pioneering workof Hayashi (1981) to have a robust initial model at hand.Subsequent evolution of the disc may change this simpleapproach, but it is still a valuable guideline.A basic estimate of the minimum surface density of solidmaterial in the disc can be deduced from the mass and lo-cation of the present planets in the solar system:Σ solid ( r ) = (cid:40) . r/ − / g / cm . (cid:54) r (cid:54) . . r/ − / g / cm . (cid:54) r (cid:54) . gas ( r ) = 1700 (cid:16) r (cid:17) − / gcm (1)Since the dust content is rather low, the gaseous componentis transparent to the visible solar radiation. Thus the gastemperature follows from the radiation balance: T = T (cid:16) r (cid:17) − / (cid:18) LL (cid:12) (cid:19) / T = 280 K (2) L is the solar luminosity during the early stages, normalisedby the present value L (cid:12) . The three-dimensional densitystructure is given by an isothermal profile ρ gas ( r, z ) = Σ gas √ πh exp (cid:18) − z h (cid:19) (3)with a radially changing scale height h : h = c s Ω c s = (cid:18) k B Tµm H (cid:19) / µ = 2 . c (0) s Ω (cid:16) r (cid:17) / (cid:18) LL (cid:12) (cid:19) / c (0) s = 993 .
56 ms (4) c s is the sound velocity of an ideal gas with a mean molec-ular weight µ in units of the hydrogen mass m H . Since thedensity profile is related to a radially varying pressure, thegas velocity deviates from the local Keplerian velocity. Thebalance of forces relates the angular velocity Ω g to the pres-sure gradient: r Ω g = r Ω + 1 ρ dPdr (5) c (cid:13) , 1–42 ybrid methods: A new composite algorithm Thus the angular velocity Ω g of the gas is (see e.g. Adachiet al. 1976): Ω g = Ω (cid:112) − η g ( r ) η g = − d ln( ρ gas c s ) d ln( r ) (cid:18) c s v K (cid:19) (6)It is more appropriate to formulate the rotation of thegaseous disc in terms of a velocity lag ∆ v g normalised tothe local Keplerian velocity v K :∆ v g = r (Ω g − Ω) ≈ − η g v K (7)A typical value of ∆ v g for the minimum mass solar nebulaat 1 AU is ∆ v g = −
60 m / s.This simple model provides a brief description of theinitial disc. However, the subsequent evolution further mod-ifies the structure of the protoplanetary disc. Since embed-ded dust grains are coupled to the gas, it is likely that aglobal migration of solid material changes the surface den-sity. Moreover, the dust grains are chemically processed, de-pending on the local temperature and composition which in-troduces additional spatial inhomogeneities. When the grow-ing particles pass the critical size of ∼ Planetesimals in a protoplanetary disc are subjected to var-ious perturbations: Close encounters change their orbits, asmall but steady gas drag gives rise to a radial drift andaccretion changes the mass of the planetesimals. While alltheses processes drive the disc evolution on a timescale of atleast a few thousand years, each planetesimal moves most ofthe time on an orbit close to an unperturbed Kepler ellipse.Though the protoplanetary disc introduces additional per-turbations, the central potential dominates for typical discmasses around 0 . M (cid:12) . Therefore the classical orbital ele-ments still provide a proper framework to study planetesimaldynamics.The orbital elements of a test particle moving around amass M are: E = v − GMra = − GM Ee = (cid:114) − L GMa cos( i ) = L z L (8) L = r × v e cos( φ E ) = rv GM − e sin( φ E ) = r · v √ GMa (9) a is the semimajor axis, e is the eccentricity and i is theinclination of the orbit. As long as no dominant body isstructuring the protoplanetary disc, it is justified to assumeaxisymmetry. Hence the argument of the perihelion ω , thelongitude of the ascending node Ω and the eccentric anomaly φ E are omitted in the statistical description.The deviation of planetesimal orbits from a circle isquite small. Thus it is appropriate to expand the above setof equations. A planetesimal at a distance r , in a distance z above the midplane and with a velocity ( v r , v φ , v z ) withrespect to the local circular velocity v K has orbital elements(leading order only): a ≈ r + 2 r v φ v K (10) e ≈ v r + 4 v φ v K (11) i ≈ z r + v z v K (12)These expressions allow us later a convenient transforma-tion between the statistical representation through orbitalelements and the utilisation of a velocity distribution func-tion. When two planetesimals pass close by each other, they ex-change energy and angular momentum and separate withmodified orbital elements. Successive encounters transfer en-ergy between planetesimals with different masses, driving anevolution of the overall velocity distribution.It seems that an encounter is a two-body problem, asthere are only two planetesimals involved, but the centralmass has also a major influence turning the problem intoa three-body encounter The complexity of the problem isconsiderably simplified by reducing it to Hill’s problem (Hill1878).Consider two masses m and m that orbit a muchlarger mass M c , where both masses are small compared tothe central mass M c . The mass ratio m : m could be arbi-trary. This special type of a three body problem is denotedas Hill’s problem, originally devised to calculate the orbitof the moon. It provides a convenient framework to exam-ine planetesimal encounters in the potential of a star. Theequations of motion of the two planetesimals including the The term three-body encounter does not imply a close passageof all involved bodies, but emphasises the strong influence of athird one. The following derivation is quite common to the literature (seee.g. Ida 1990; H´enon 1986). The later work uses a slightly differentscaling.c (cid:13) , 1–42
P. Glaschke, P. Amaro-Seoane & R. Spurzem central potential and their mutual interaction are: ¨r = − r GM c r − ( r − r ) Gm r ¨r = − r GM c r − ( r − r ) Gm r (13)We now introduce the relative vector r and the centre-of-mass R : r = r − r R = m r + m r m + m (14)Furthermore, the equations of motion are transformed to acorotating set of coordinates which are scaled by the mutualHill radius r Hill and the local Kepler frequency Ω x = r (cid:48) − a r Hill y = a ( φ − Ω t ) r Hill z = z (cid:48) r Hill r Hill = a (cid:114) m + m M c (15)where ( r (cid:48) , φ, z (cid:48) ) are heliocentric cylindrical coordinates. a is the radius of a properly chosen reference orbit, while themutual Hill radius r Hill is an intrinsic length scale of theproblem : Since the two orbiting planetesimals are smallcompared to the central star, the Hill radius is much smallerthan the size a of the reference orbit. Hence it is possibleto expand the central potential about the reference orbit.This yields approximate equations for the relative motionand the centre-of-mass:¨ x = 2 ˙ y + 3 x − x/r ¨ y = − x − y/r ¨ z = − z − z/r (16)¨ X = 2 ˙ Y + 3 X ¨ Y = − X ¨ Z = − Z (17)The equations 16 and 17 have some interesting proper-ties: Firstly, the centre-of-mass motion separates from theinteraction of the two bodies. Secondly, the scaled relativemotion is independent of the masses m and m , imply-ing a fundamental similarity of planetesimal encounters .As Eq. 17 is a simple linear differential equation, one read-ily obtains the solution of the centre-of-mass motion X = b − e cos( t − τ ) Y = − bt + ψ + 2 e sin( t − τ ) Z = i sin( t − ω ) (18)which is equivalent to a first-order expansion of a Keplerellipse. ω and τ are the longitudes of the ascending nodeand the pericentre, while e and i are the eccentricity and section 16 makes extensive use of this property. -2-1.5-1-0.5 0 0.5 1 1.5 2 2.5 -3 -2 -1 0 1 2 3 x yL L U * Hill sphere
Figure 1.
Equipotential lines for the effective potential U at z = 0 (see Eq. 20). U = U ∗ refers to the largest allowed volume,which is enclosed by the Hill sphere and the two Lagrange points L and L . inclination scaled by the reduced (i. e. dimensionless) Hillradius r Hill /a . The value of b depends on the choice of thereference orbit, but it is natural to set b = 0 which impliesthat the centre-of-mass defines the reference orbit.While the nonlinear nature of the relative motion (seeEq. 16) prevents any general analytical solution, Eq. 18 pro-vides at least an asymptotic solution for a large separation ofthe planetesimals, where b can be interpreted as an impactparameter. Nevertheless, small b do not necessarily implyclose encounters, as opposed to the standard definition ofthe impact parameter. However, the expression b = a − a provides a measure of the distance of the two colliding bod-ies without invoking the complicated encounter geometry.A special solution to Eq. 16 are the unstable equilibriumpoints L and L at ( x, y, z ) = ( ± , , . In addition, an inspection of Eq. 16 revealsthat the Jacobi energy E J is conserved: E J = 12 (cid:0) ˙ x + ˙ y + ˙ z + z − x (cid:1) − r (19)Since the kinetic energy is always a positive quantity, thefollowing inequality holds: E J (cid:62) U = 12 (cid:0) z − x (cid:1) − r (20)Thus the allowed domain of the particle motion is enclosedby the equipotential surfaces of the effective potential U . Asubset of these equipotential surfaces restricts the alloweddomain to the vicinity of the origin (see Fig. 1). The largestof these surfaces passes through the Lagrange points L and L . Hence we identify the Hill radius (or Hill sphere) as themaximum separation which allows the bound motion of twoplanetesimals . The additional Lagrange points L – L are missing due to theHill approximation. The same argument applies to the tidal boundary in clusterc (cid:13) , 1–42 ybrid methods: A new composite algorithm Beside the numerical solution of the equations of mo-tion, it is useful to define osculating (or instantaneous) or-bital elements b = 4 x + 2 ˙ yi = z + ˙ z e = ˙ x + (3 x + 2 ˙ y ) (21) ψ = − x + y + 32 btω = t − arctan( z, ˙ z ) τ = t − arctan( ˙ x, x + 2 ˙ y ) (22) E J = 12 ( e + i ) − b − r (23)which provide a convenient description of the initial relativeorbit and the modified orbit after the encounter. Our work follows the evolution of a planetesimal disc intofew protoplanets, including the full set of interaction pro-cesses. Hence we summarise the main aspects of protoplanetgrowth first to provide a robust framework.Although the sizes of the planetesimal cover a widerange, they virtually form two main groups: The smallerfield planetesimals and the embedded protoplanets (or theirprecursors). This two-group approximation (e.g. Wetherill1989; Ida 1993) allows one to have a clearer insight into thegrowth process.During the initial phase all planetesimals share the samevelocity dispersion independently of their mass. The initialrandom velocities are low enough for an efficient gravita-tional focusing. Hence, the growth rate of a protoplanet withmass M radius R can be estimated as (e.g. Ida et al. 1993):˙ M ≈ v rel Σ H πR (cid:18) GMRv (cid:19) (24)Eq. 24 is the two-body accretion rate, which should be mod-ified due to the gravity of the central star. Nevertheless thisapproximation gives an appropriate description to discussthe basic properties of the growth mode. The scale height H (Eq. 215) and the relative velocity v rel are related to themean eccentricity e m = (cid:112) (cid:104) e (cid:105) of the field planetesimals: H ≈ v rel / Ω v rel ≈ e m a Ω (25)Thus the accretion rate in the limit of strong gravitationalfocusing (2
GM/R (cid:29) v ) is:˙ M ≈ πR Σ GMa Ω e m (26) ∝ M / (27)If the protoplanets are massive enough, they start to controlthe velocity dispersion of the planetesimals in their vicinity. dynamics or the Roche lobe in stellar dynamics, which are equiv-alent to the Hill sphere. The width ∆ a of this sphere of influence, the heating zone,is related to the Hill radius R Hill of the protoplanet (Ida1993): ∆ a = ∆˜ aR Hill = 4 R Hill (cid:114)
43 (˜ e m + ˜ i m ) + 12 h = (cid:114) M M c (28)˜ e m and ˜ i are eccentricity and inclination of the field plan-etesimals, scaled by the reduced Hill radius h of the pro-toplanet. M c is the mass of the central star. The conditionthat the protoplanet controls the velocity dispersion of thefield planetesimals reads (Ida et al. 1993):2 M πa ∆ a > Σ m (29)This condition is equivalent to a lower limit of the proto-planetary mass: Mm > (cid:18) π ∆˜ a √ (cid:19) / (cid:18) Σ a M c (cid:19) / (cid:18) mM c (cid:19) / (30) M/m depends on several parameters, but reasonable valuesyield
M/m ≈ v ≈ R Hill
Ω (31)which gives an interesting relation to the condition thatleads to gap formation. A protoplanet can open a gap inthe planetesimal component if it is larger than a criticalmass M gap (Rafikov 2001) M gap M c ≈ Σ a M c (cid:16) mM c (cid:17) / if v (cid:46) Ω r HillΣ a M c (cid:16) mM c (cid:17) / (cid:16) Ω r Hill v (cid:17) if v (cid:29) Ω r Hill (32)where r Hill is the Hill radius of the field planetesimals. If thevelocity dispersion v is controlled by the protoplanet, Eq. 32together with Eq. 31 demonstrate that the condition for gapformation is equivalent to Eq. 30, i. e. the efficient heating ofthe field planetesimals implies gap formation and vice versa.The higher velocity dispersion of the field planetesimals (seeEq. 31) reduces the growth rate given by equation 26 to˙ M ≈ π ΣΩ RR Hill ˜ e m (33) ∝ M / (34)Different mass accretion rates imply different growth mode.If two protoplanets have different masses M and M , theirmass ratio evolves as: ddt M M = M M (cid:18) ˙ M M − ˙ M M (cid:19) (35)When the growth timescale M/ ˙ M decreases with mass, asmall mass difference increases with time. This is the casefor Eq. 26, which gives rise to runaway accretion. As soon asthe protoplanets control the velocity dispersion of the fieldplanetesimals, the growth timescale increases with mass andtherefore the protoplanet masses become more similar asthey grow oligarchically.The field planetesimals also damp the excitation due c (cid:13) , 1–42 P. Glaschke, P. Amaro-Seoane & R. Spurzem
Σ[g / cm ] M iso /M (cid:12) M iso /M ⊕ . × − . × − . × − Table 1.
Isolation mass for different surface densities at r = 1AU and M c = 1 M (cid:12) . to protoplanet–protoplanet interactions and keep them onnearly circular orbits. The balance between these scatteringsand the dynamical friction due to smaller bodies establishesa roughly constant orbital separation b (Kokubo 1997): b = R Hill (cid:114) e m M π Σ aR (36)˜ b = b/R Hill (37) R is the radius of the protoplanet, M is its mass and ˜ e is the reduced eccentricity of the field planetesimals. Thestabilised spacing prevents collisions between protoplanets,but it also restricts the feeding zone – the area from which aprotoplanet accretes. If all matter in the feeding zone is ac-creted by the protoplanet, it reaches its final isolation mass(Kokubo & Ida 2000): M iso = 2 πba Σ (38)Inserting Eq. 36 yields the isolation mass in units of themass of the host star M c : M iso /M c = (112 π ) / (cid:18) (cid:19) / (cid:18) π (cid:19) / (cid:0) ˜ e m (cid:1) / (cid:18) a Σ M c (cid:19) / (cid:18) a ρM c (cid:19) / ≈ . × (cid:0) ˜ e m (cid:1) / (cid:18) a Σ M c (cid:19) / (cid:18) a ρM c (cid:19) / (39)There is a weak dependence on the density ρ of the pro-toplanet, but the most important parameter is the surfacedensity Σ. If we take the minimum mass solar nebula as anexample, the radial dependence of the surface density im-plies an isolation mass that grows with increasing distanceto the host start. Hence protoplanets beyond some criticalradial distance are massive enough (larger than ≈ M ⊕ Bodenheimer 1986) to initiate gas accretion from the proto-planetary disc.As the protoplanets approach the isolation mass, inter-actions with the gaseous disc and neighbouring protoplanetsbecome increasingly important. We estimate the onset of or-bit crossing by a comparison of the perturbation timescale τ pert of protoplanet–protoplanet interactions with the damp-ing timescale τ damp due to planetesimal–protoplanet scatter-ings. Since the protoplanets are well separated (˜ b ≈ . . . τ pert ≈ ˜ b h Ω (40)We anticipate section 16 (see Eq. 221) to derive the damping timescale τ damp ≈ T / r √ πG ln(Λ)( M + m ) n m (41)where T r and T z are the radial and vertical velocity disper-sion of the field planetesimals. Hence the criterion for theonset of orbital crossing is: τ pert < τ damp (42)As the protoplanets control the velocity dispersion of thefield planetesimals (see Eq. 31), this condition reduces to:Σ M > Σ m ln(Λ) 727 π (cid:18) ˜ b ˜ e (cid:19) > Σ m × f (43)Thus orbital crossing sets in when the mean surface den-sity Σ M of the protoplanets exceeds some fraction f of thefield planetesimal density Σ m . While the factor f dependsstrongly on the separation ˜ b of the protoplanets, a fiducialvalue is f ≈
1, in agreement with the estimates of Goldreichet al. (2004).The onset of migration and the resonant interactionof protoplanets with the disc and other protoplanets ter-minates the local nature of the protoplanet accretion pro-cess and requires a global evolution of the planetary system.While the final stage deserves a careful analysis, further re-search is beyond the scope of this work.
The protoplanet formation is essentially an N –body prob-lem. Although we seek for a more elaborated solution to thisproblem which benefits from statistical methods, the pure N –body approach is a logical starting point. Direct calcu-lations with a few thousand bodies have provided us witha valuable insight into the growth mode (Ida 1992; Kokubo1996,see e.g.), but they are also powerful guidelines that helpdeveloping other techniques. Statistical calculations rely ona number of approximations and “exact” N –body calcula-tions provide the necessary, unbiased validation of the de-rived formula.The choice of the integrator is a key element in thenumerical solution of the equations of motion. Our require-ments are the stable long-term integration of a few ten thou-sand planetesimals with the capability of treating close en-counters, collisions and the perspective to evolve it into animproved hybrid code. Approximative methods like the FastMultipole Method or Tree codes have a scaling of the com-putational time close to N , but the accuracy in this regimeis too poor to guarantee the stable integration of Keple-rian orbits (compare the discussion in Hernquist et al. 1993;Spurzem 1999).The class of exact methods (all scale asymptoticallywith N ) roughly divides in two parts:(i) Symplectic methods (see e.g. Wisdom 1991 or the Symba code, Duncan et al. 1998) rely on a careful expan-sion of the Hamiltonian which guarantees that the numericalintegration follows a perturbed Hamiltonian. While there is c (cid:13)000
The protoplanet formation is essentially an N –body prob-lem. Although we seek for a more elaborated solution to thisproblem which benefits from statistical methods, the pure N –body approach is a logical starting point. Direct calcu-lations with a few thousand bodies have provided us witha valuable insight into the growth mode (Ida 1992; Kokubo1996,see e.g.), but they are also powerful guidelines that helpdeveloping other techniques. Statistical calculations rely ona number of approximations and “exact” N –body calcula-tions provide the necessary, unbiased validation of the de-rived formula.The choice of the integrator is a key element in thenumerical solution of the equations of motion. Our require-ments are the stable long-term integration of a few ten thou-sand planetesimals with the capability of treating close en-counters, collisions and the perspective to evolve it into animproved hybrid code. Approximative methods like the FastMultipole Method or Tree codes have a scaling of the com-putational time close to N , but the accuracy in this regimeis too poor to guarantee the stable integration of Keple-rian orbits (compare the discussion in Hernquist et al. 1993;Spurzem 1999).The class of exact methods (all scale asymptoticallywith N ) roughly divides in two parts:(i) Symplectic methods (see e.g. Wisdom 1991 or the Symba code, Duncan et al. 1998) rely on a careful expan-sion of the Hamiltonian which guarantees that the numericalintegration follows a perturbed Hamiltonian. While there is c (cid:13)000 , 1–42 ybrid methods: A new composite algorithm still an integration error, all properties of a Hamiltonian sys-tem like conservation of phase-space volume are conservedby the numerical integration. The drawback of these very el-egant methods is that the symplecticity is immediately bro-ken by adaptive time steps, collisions or complicated exter-nal forces if no special precaution is taken. Even the numeri-cal truncation error breaks the symplecticity to some extend(Skeel 1999). Moore & Quillen (2010) recently developed aparallel integrator for graphics processing unit (GPU) thatuses symplectic and Hermite algorithms according to theresolution needed. The algorithm is similar to Symba , butless accurate.(ii) The second group represents the “classical” methodsthat are based on Taylor expansions of the solution. Theycome in different flavours like implicit methods, predictor–corrector integrators or iterated schemes. Time symmetricmethods stand out among these different approaches, asthey show no secular drift in the energy error. These in-tegrators are so well-behaved that one may even call them“nearly symplectic”.Taking all the requirements into account we have chosen
Nbody6++ , an integrator which is the most recent descen-dant from the famous N –body code family from S. Aarseths’“factory”. This version was parallelised by Spurzem (1999),which opened the use of current supercomputers.This parallel version, named Nbody6++ , offers manyversatile features that were included over the past years andmore and more refined as time passed by. While all theseelaborated tools deserve attention, we restrict ourselves forbrevity to the components which contribute to the planetes-imal problem. The main components of the code are:(i) Individual time steps and a block time step scheme.(ii) Ahmad–Cohen neighbour scheme.(iii) Hermite scheme.(iv) KS–Regularisation for close two-body encounters.We will explain each of these features and highlight theiradvantages for the main goal of this work. This is impor-tant to understand the immediate first level of modifica-tions required to adapt the numerical scheme to planetarydynamics. We present later, in section 8, the next level ofcomplexity in the modifications to be carried out to moldthe numerical tool to our purposes.
The choice of the time step controls the accuracy as wellas the efficiency of any given integrator. Too small timesteps slow down the integration without necessity, whereastoo large values increase the error. An efficient solution tothis dilemma is an adaptive time step that is adjusted af-ter each integration step according to a specified accuracylimit. While the idea is quite clear, there is no unique receipthow to choose the proper time step. A common approach for N –body systems is to use the parameters at hand (like par-ticle velocity, force, etc.) to derive a timescale of the particlemotion. This procedure leaves enough space for a wealth of Aarseth (1999) gives a nice review on the remarkable historyof the
Nbody -codes, more details are given in Aarseth (2003). different time step criteria.
Nbody6++ uses the standardAarseth expression (Aarseth 1985)∆ t = (cid:115) η F F (2) + ( F (1) ) F (1) F (3) + ( F (2) ) (44)which makes use of the force and the time derivatives up tothird order.The time step choice is not unique in multi-particle sys-tems. One solution is to take the minimum of all these valuesas a shared (or global) time step, but this is not recom-mended unless all individual steps are of the same size.The second option is to evolve each particle track withits own, individual time step. This method abandons theconvenient synchronisation provided by a global step, there-fore each integration of a particle demands a synchronisationof all particles through predictions. Since the prediction ofall particles is O ( N ), it is counter-balanced by the savingof force computations. Nevertheless the overhead is still sig-nificant, so a further optimisation might be desired. Thebasic idea of the block time step method is to force parti-cles in groups that are integrated together, which reducesthe number of necessary predictions by a factor comparableto the mean group size. These groups are enforced throughtwo constraints. The first condition is a discretisation of thesteps in powers of two:∆ t = 2 − k k (cid:62) T i ≡ t i (46)i. e. the particle time T i is an integer multiple of the actualtime step ∆ t i so that all particles with the same step sharethe same block. Note that a step can not be increased atany time T i , but only when the second condition (Eq. 46)for the larger step is met. The first N –body codes calculated always the full force (i. e.summation over all particles) to integrate a particle. Butnot all particles contribute with the same weight to the to-tal force. Distant particles provide a smooth, slowly vary-ing regular force, whereas the neighbouring particles form arapidly changing environment which gives rise to an irregu-lar force. The Ahmad-Cohen neighbour scheme (Ahmad &Cohen 1973) takes advantage of this spatial hierarchy by di-viding the surrounding particles in the already mentionedtwo groups according to the neighbour sphere radius R s .Both partial forces fluctuate on different timescales, whichare calculated according to Eq. 44. The key to the efficiencyof the method is the inequality∆ t irr (cid:28) ∆ t reg (47)Regular forces are extrapolated between two full force cal-culations F reg = F (0)reg + ∆ t F (1)reg + 12 (∆ t ) F (2)reg + 16 (∆ t ) F (3)reg (48) c (cid:13) , 1–42 P. Glaschke, P. Amaro-Seoane & R. Spurzem Rs Figure 2.
All particles inside the neighbour sphere R s are se-lected as neighbours (filled circles). As the neighbour list is fixedduring regular steps (marked by black lines), a second shell in-cludes possible intruders. with a third order accurate expression, whereas the plain Nbody6++ uses only linear extrapolation (see the next sec-tions for a detailed discussion).
The Hermite scheme is a fourth-order accurate integratorthat was applied first by Makino (1992) to the integra-tion of a planetesimal system. This scheme is used to in-tegrate single particles and CM-bodies in the main part of Nbody6++ . It is accomplished by a predictor and a correc-tor step. The prediction is second order accurate: x p = x + v ∆ t + 12 a (∆ t ) + 16 ˙a (∆ t ) v p = v + a ∆ t + 12 ˙a (∆ t ) (49)Now the acceleration is evaluated at the predicted positionto derive the second and third order time derivatives of theforce: a (2)0 = − ˙a − ˙a p ∆ t + − a + 6 a p ∆ t (50) a (3)0 = 6 ˙a + 6 ˙a p ∆ t + 12 a − a p ∆ t (51)Using the additional derivatives one can improve the predic-tion: x c = x p + 124 a (2)0 (∆ t ) + 1120 a (3)0 (∆ t ) + O (∆ t )= x p + (cid:18) − a + 320 a p (cid:19) (∆ t ) − (cid:18) ˙a + 130 ˙a p (cid:19) (∆ t ) + O (∆ t ) (52) v c = v p + 16 a (2)0 (∆ t ) + 124 a (3)0 (∆ t ) + O (∆ t )= v p + 12 ( − a + a p )∆ t + (cid:18) − ˙a − ˙a p (cid:19) (∆ t ) + O (∆ t ) (53) CM denotes centre–of–mass, see the section on KS-Regularisation for more details.
The corrected positions are fourth order accurate. While theHermite scheme is robust and stable, even in combinationwith the neighbour scheme, it is not accurate enough tointegrate planetesimal orbits efficiently.
The Hermite scheme bears the potential for a powerful ex-tension, since it is a predictor–corrector scheme. An essentialpart of this scheme is the calculation of the new forces withthe predicted positions, but it should improve the accuracyif they are recalculated using the corrected positions. Thenew forces are readily used to improve the corrected values,which closes the scheme to an iteration loop – the Hermiteiteration (Kokubo et al. 1998).There are only few modifications necessary to obtainthe iterated scheme. The particle prediction remains second-order accurate: x p = x + v ∆ t + 12 a (∆ t ) + 16 ˙a (∆ t ) v p = v + a ∆ t + 12 ˙a (∆ t ) (54)The force and its first time derivative are calculated to derivehigher derivatives according to Eq. 50 and 51: a p = f ( x p , v p ) (55) ˙a p = f ( x p , v p ) (56)The new corrector omits the highest order term in the posi-tion, making the velocity and the position to the same orderaccurate: x c = x p + 124 a (2)0 (∆ t ) + O (∆ t )= x p + 14 ( a p − a )(∆ t ) + (cid:18) − ˙a − ˙a p (cid:19) (∆ t ) + O (∆ t ) (57) v c = v p + 12 ( − a + a p )∆ t + (cid:18) − ˙a − ˙a p (cid:19) (∆ t ) + O (∆ t ) (58)It seems unreasonable to drop one order in the position, buta reformulation of the predictor–corrector step reveals thatthis slight change yields a time symmetric scheme: v c = v + 12 ( a p + a )∆ t −
112 ( ˙a p − ˙a )(∆ t ) x c = x + 12 ( v c + v )∆ t −
112 ( a p − a )(∆ t ) (59)The iteration is achieved by returning to Eq. 55–56 with thepredicted positions replaced by the more accurate values x c , v c . Although the integration does not need the secondand third time derivative of the forces explicitly, it is usefulto provide them at the end of the iteration for the calculationof the new time step: a (2) ( t + ∆ t ) = 2 ˙a + 4 ˙a p ∆ t + 6 a − a p ∆ t (60) a (3) ( t + ∆ t ) = 6 ˙a + 6 ˙a p ∆ t + 12 a − a p ∆ t (61)Numerical tests show that convergence is reached after oneor two iterations, making this approach very efficient. Itneeds less force evaluations than the Hermite scheme for c (cid:13)000
112 ( a p − a )(∆ t ) (59)The iteration is achieved by returning to Eq. 55–56 with thepredicted positions replaced by the more accurate values x c , v c . Although the integration does not need the secondand third time derivative of the forces explicitly, it is usefulto provide them at the end of the iteration for the calculationof the new time step: a (2) ( t + ∆ t ) = 2 ˙a + 4 ˙a p ∆ t + 6 a − a p ∆ t (60) a (3) ( t + ∆ t ) = 6 ˙a + 6 ˙a p ∆ t + 12 a − a p ∆ t (61)Numerical tests show that convergence is reached after oneor two iterations, making this approach very efficient. Itneeds less force evaluations than the Hermite scheme for c (cid:13)000 , 1–42 ybrid methods: A new composite algorithm the same accuracy. Moreover, its time symmetry suppressesa secular drift of the energy error. The Hermite scheme is an integral part of
Nbody6++ and proved its value in many applications. It would havebeen natural to improve the performance with an additionaliteration, but our first tentative implementations showedrather negative results: The iterated scheme was more un-stable, slower and even less accurate than the plain Hermitescheme. An inspection of the code structure revealed thatthe Ahmad–Cohen neighbour scheme is the cause of this.Each particle integration is composed of two parts –frequent neighbour force evaluations and less frequent totalforce evaluations including derivative corrections. Every reg-ular correction leads to an additional change in the positionof a particle, which introduces a spurious discontinuity inthe neighbour force and its derivatives. The Hermite itera-tion reacts to this artificial jump in two ways: It increases theregular correction, and – what is more important – it am-plifies any spurious error during the iteration which leads toan extreme unstable behaviour.Since the Hermite iteration is a key element in the ef-ficient integration of planetesimal orbits, we sought for amodification of the Hermite scheme that circumvents the de-picted instability. The problem gives already an indication ofa possible solution. A scheme with much smaller correctionswould not suffer from the feedback of spurious errors.
Nbody6++ stores already the second and third timederivative of the forces for the time step calculation. A man-ifest application of these derivatives at hand is the improve-ment of the predictions to fourth order: x p = x + v ∆ t + 12 a (∆ t ) + 16 ˙a (∆ t ) + 124 a (2)0 (∆ t ) + 1120 a (3)0 (∆ t ) (62) v p = v + a ∆ t + 12 ˙a (∆ t ) + 16 a (2)0 (∆ t ) + 124 a (3)0 (∆ t ) (63)The prediction to fourth order was used in the iterativeschemes of Kokubo et al. (1998); Mikkola & Aarseth (1998).Again, the new forces a p and ˙a p are calculated to improve x p and v p – but with a modified corrector: a p = f ( x p , v p ) (64) ˙a p = f ( x p , v p ) (65) a (2) n ( t ) = − ˙a − ˙a p ∆ t + − a + 6 a p ∆ t (66) a (3) n ( t ) = 6 ˙a + 6 ˙a p ∆ t + 12 a − a p ∆ t (67) x c = x p + 124 ( a (2) n − a )(∆ t ) +1120 ( a (3) n − a )(∆ t ) (68) v c = v p + 16 ( a (2) n − a (2) )(∆ t ) +124 ( a (3) n − a (3) )(∆ t ) (69) Finally, the derivatives are updated: a ( t + ∆ t ) = a p (70) ˙a ( t + ∆ t ) = ˙a p (71) a (2) ( t + ∆ t ) = a (2) n + ∆ t a (3) n (72) a (3) ( t + ∆ t ) = a (3) n (73)The new scheme has an appealing property, which is relatedto the usage of the higher force derivatives. As the predictoris fourth-order accurate, it is equivalent to one full Hermitestep. Since the corrector uses new forces to improve the twohighest orders, it is equivalent to a first iteration step. Thuswe obtained a one-fold iterated Hermite scheme at no ex-tra cost. This extended Hermite scheme reduces the energyerror by three orders of magnitude, compared to the plain Nbody6++ with the same number of force evaluations.
Two bodies undergoing a close encounter are integrated ina special set of regular coordinates that separates the rela-tive motion from the motion of the centre-of-mass. A closeencounter poses a numerical problem due to the singular-ity of the gravitational forces at zero separation. While thegrowing force amplifies roundoff errors as the two bodiesapproach each other closely, the collision is only an appar-ent singularity since the analytic solution stays well-defined.This opens the possibility of a proper coordinate transfor-mation which removes the singularity from the equations ofmotion. The Kustaanheimo–Stiefel regularisation takes ad-vantage of a four-dimensional set of variables to transformthe Kepler problem into a harmonic oscillator (Kustaan-heimo & Stiefel 1965). Perturbations are readily includedin the new set of equations of motion.The centre-of-mass is added as a pseudo-particle, theCM-body, which is integrated as a normal particle plus aperturbation force due to the deviation from a point mass.See Mikkola (1997) or Mikkola & Aarseth (1998) for moredetails.
Nbody6++ only includes the gravitational interaction of allparticles, therefore additional forces have to be added “perhand”. A planetesimal disc requires two new forces: Thepresence of a central star introduces an additional centralpotential, while the gaseous component of the protoplane-tary disc is the source of a friction force. It is importantthat the new forces are properly included in the neighbourscheme to assure that regular steps remain larger than ir-regular steps. Since a dissipative force breaks the energyconservation, one has to integrate the energy loss as well tomaintain a valid energy error control. In the next subsectionswe describe how we have done this. c (cid:13) , 1–42 P. Glaschke, P. Amaro-Seoane & R. Spurzem
A star is much heavier than a planetesimal. Thus, the centralstar is introduced as a spatially fixed Keplerian potential:Φ c = − GM c r F = − GmM c x | x | (74)Since the orbital motion of the planetesimals sets the domi-nant (and largest) dynamical timescale in the system, we in-cluded the central force as a component of the regular force.Moreover, the central potential also introduces a strong syn-chronisation, since planetesimals in a narrow ring share vir-tually the same regular block time step. As the whole planetesimal system is embedded in a dilutegaseous disc, each planetesimal is subjected to a small, butnoticeable drag force. The drag regime depends on the gasdensity and the size of the planetesimals. Kilometer-sizedplanetesimals are subjected to the deceleration d v dt = − πC D m ρ gas R | v − v g | ( v − v g ) (75) C D = 0 . R ( m ) of theplanetesimal in this drag regime. v g is the rotational velocityof the gaseous disc, which rotates slower than the planetes-imal system as it is partially pressure supported. The dragforce leads to an orbital decay ˙ a of the semimajor axis of aplanetesimals: ˙ a = − C D ρ gas ρ Body (cid:104) (∆ v ) (cid:105) R ( m )Ω (77)Thus smaller particles migrate faster, with a maximum at R ≈ W drag = F drag · v ˙ W drag = ˙F drag · v + F drag · F tot (78)We integrate the dissipation rate W drag to maintain a validenergy error:∆ E = (cid:90) t t W drag dt (79)∆ t = t − t (80)∆ E = 12 ( W drag , + W drag , ) ∆ t + 112 (cid:16) ˙ W drag , − ˙ W drag , (cid:17) ∆ t + O (∆ t ) (81)The expression is fourth order accurate in accordance withthe order of the extended Hermite scheme. The main drag regimes are Stokes (laminar flow), Epstein(mean free molecular path larger than object size) and Newton’sdrag law (turbulent flow). Weidenschilling (1977) provides a nicereview on the different drag regimes.
Both new forces also demand a modification of the regu-larisation treatment. They perturb the relative motion of aKS–pair and modify the orbit of the centre-of-mass. Whilethe modification of the equations of motion is rather clear,the neighbour scheme requires some additional work.Let r , r be the positions of the two regularised particles.The equations of motion read ( G = 1) r = r − r ¨r = − M c r r + m r r + F drag , ¨r = − M c r r − m r r + F drag , (82)where the perturbations by other particles have been omit-ted for clarity. Centre-of-mass motion and the orbital motionare separated: ¨r = − M r r − M c r r + M c r r + F drag , − F drag , (83) ¨R = − M c m M r r − M c m M r r + m M F drag , + m M F drag , (84) r = r − r M = m + m (85) R = 1 M ( m r + m r ) (86)Two new contributions show up due to the externalforces: The KS–pair is tidally perturbed by the central starand influenced by the gaseous disc. While the aerodynamicproperties of a single particle are well understood, two bod-ies revolving about each other may induce complex gas flowsin their vicinity, which could invalidate the linear combina-tion of the drag forces on each component. Therefore we dropthe drag force term to avoid spurious dissipation. Since thedynamic environment allows virtually no stable binaries in a planetesimal disc, the influence of the drag force on theencounter dynamics is negligible.We further decompose the additional acceleration of thecentre-of-mass motion, since the neighbour scheme benefitsfrom a clear separation of the timescales. Therefore, the tidalperturbation is split into a smooth mean force and a pertur-bation force: ¨R = F mean + F pert F mean = − M c R R F pert = M c R R − M c m M r r − M c m M r r (87)= 0 + O ( r ) (88)The mean forces varies on the orbital timescale and is hence Tidal capturing of moons starts in the late stages of planetformation, but is limited to the planets or their precursors. How-ever, the quiescent conditions in an early Kuiper belt lead to amore prominent role of binaries. See the summary of Astakhovet al. (2005). c (cid:13)000
Both new forces also demand a modification of the regu-larisation treatment. They perturb the relative motion of aKS–pair and modify the orbit of the centre-of-mass. Whilethe modification of the equations of motion is rather clear,the neighbour scheme requires some additional work.Let r , r be the positions of the two regularised particles.The equations of motion read ( G = 1) r = r − r ¨r = − M c r r + m r r + F drag , ¨r = − M c r r − m r r + F drag , (82)where the perturbations by other particles have been omit-ted for clarity. Centre-of-mass motion and the orbital motionare separated: ¨r = − M r r − M c r r + M c r r + F drag , − F drag , (83) ¨R = − M c m M r r − M c m M r r + m M F drag , + m M F drag , (84) r = r − r M = m + m (85) R = 1 M ( m r + m r ) (86)Two new contributions show up due to the externalforces: The KS–pair is tidally perturbed by the central starand influenced by the gaseous disc. While the aerodynamicproperties of a single particle are well understood, two bod-ies revolving about each other may induce complex gas flowsin their vicinity, which could invalidate the linear combina-tion of the drag forces on each component. Therefore we dropthe drag force term to avoid spurious dissipation. Since thedynamic environment allows virtually no stable binaries in a planetesimal disc, the influence of the drag force on theencounter dynamics is negligible.We further decompose the additional acceleration of thecentre-of-mass motion, since the neighbour scheme benefitsfrom a clear separation of the timescales. Therefore, the tidalperturbation is split into a smooth mean force and a pertur-bation force: ¨R = F mean + F pert F mean = − M c R R F pert = M c R R − M c m M r r − M c m M r r (87)= 0 + O ( r ) (88)The mean forces varies on the orbital timescale and is hence Tidal capturing of moons starts in the late stages of planetformation, but is limited to the planets or their precursors. How-ever, the quiescent conditions in an early Kuiper belt lead to amore prominent role of binaries. See the summary of Astakhovet al. (2005). c (cid:13)000 , 1–42 ybrid methods: A new composite algorithm γ ~ t r-/r Hill
Time Step Ratio σ / Ω r Hill =0 σ / Ω r Hill =2.2 σ / Ω r Hill =10 Analytic
Figure 3.
Time step ratio for N nb = 100. Curves are plotted fordifferent values of σ v / ( r Hill
Ω). The dotted line is approximationEq. 95. included as a regular force component, while the perturba-tion is treated as an irregular force as it changes with theinternal orbital period of the pair. N − BODY TO THE DISC
An astrophysical simulation is a tool to analyse problemsand predict dynamical systems which are not accessibleto experiments. The design of a new simulation tool doesnot only require the careful implementation of the invokedphysics, but also an analysis of the code performance tomake best use of the available hardware.We want to apply
Nbody6++ for the first time toplanet formation, a subject that is quite different to stellarclusters. The central star forces the planetesimals on regu-lar orbits which need higher accuracy than the motion ofstars in a cluster. In addition, the orbital motion also in-troduces a strong synchronisation among the planetesimals,thus allowing a more efficient integration.We examine the differences due to the integration of adisc system in the following sections. In particular, we willaddress in detail the role of the geometry of the problemand the neighbour scheme, the prediction of the number ofneighbouring particles, the communication, the block sizedistribution and the optimal neighbour particle number forthe direct-summation of the massive particles in the proto-planetary system
The introduction of the neighbour scheme by Ahmad & Co-hen (1973) has provided us with a technique to save a con-siderable amount of computational time in star cluster simu-lations. Since the average ratio of the regular to the irregulartime step γ t is of the order of 10, the integration is speededup by the same factor. One may expect a similar speedupfor planetesimal systems, but in this case the time step ratiois roughly three. The time step is calculated with the stan-dard Aarseth time step criterion (it should be mentioned, however, that the relation of regular and irregular costs ismore complicated with GPU technology)∆ t = (cid:115) η F F (2) + ( F (1) ) F (1) F (3) + ( F (2) ) (89)where F ( i ) are the force and its time derivatives. It is appliedto the calculation of the regular step using the regular forceand accordingly to the irregular step based on the irregularforce. The regular time step of a particle orbiting the centralmass M c at a distance r is:∆ t r = √ η r
1Ω Ω = (cid:115) GM c r (90)For simplicity, we introduce the scaled timestep ratio ˜ γ t = γ t (cid:112) η irr /η reg . The free parameters of the problem are themean particle distance ¯ r , the velocity dispersion σ v (addi-tional to the Keplerian shear), the particle mass m i and theneighbour number N nb . We employ Hill’s approximation forthe central potential and obtain:∆ t i ≈ √ η i Ω f (Ω , r Hill , ¯ r, σ v , N nb ) (91) r Hill = r (cid:114) m i M c (92) f is a yet unknown function. Dimensional analysis leads to∆ t i ≈ √ η i Ω f (cid:18) σ v r Hill Ω , ¯ rr Hill , N nb (cid:19) (93)˜ γ t ≈ f (cid:18) σ v r Hill Ω , ¯ rr Hill , N nb (cid:19) (94)which shows that the time step ratio is essentially con-trolled by the interparticle distance and the velocity disper-sion. We generated different random realisations of planetes-imals discs with different densities and velocity dispersionsto cover the range of possible values. The neighbour numberis fixed to N nb = 100 to reduce the noise due to small num-ber statistics, but γ t converges to a value independently ofthe neighbour number already for N nb >
10. Fig. 3 showsthe numerical calculation of the time step ratio for variousvalues of ¯ r and σ v . A good approximation to the calculatedvalues of ˜ γ t is:˜ γ t ≈ . × (cid:114) . σ v ¯ r Ω + 0 . r ¯ r (95)Planetesimal discs have usually a small velocity dispersion(compared to the orbital velocity) and a low density in termsof the Hill radius, which leaves a major influence to the Ke-plerian shear. Since the shear motion is directly linked tothe local Keplerian frequency, this synchronisation reduces γ t to values smaller than ten. The numerical calculationsshow larger time step ratios with increasing velocity disper-sion and for high densities , but planetesimal discs are farfrom these extreme parameter values. c (cid:13) , 1–42 P. Glaschke, P. Amaro-Seoane & R. Spurzem N r e g / ∆ T Exponent(1) V-Exponent=0(2) M-Exponent=1/6
Figure 4.
Regular steps per particle and per 1 N –body timein the inner core ( R i < .
5) of a 5000 particle plummer model.Plotted are (1) different mass exponents with velocity exponent0 and (2) different velocity exponents with mass exponent 1 / The standard neighbour criterion uses the geometrical dis-tance: Particles are neighbours if their distance to the refer-ence particle is smaller than a limit R s . This criterion is sim-ple and probably the best choice for an equal mass system.However, a multi-mass system may require a different crite-rion, since a massive particle outside the neighbour spherecould have a stronger influence than lighter particles insidethe neighbour sphere. Also, relative effects are smaller atlarge distances. A more appropriate selection should rely onsome “perturbation strength” of a particle.It turned out that a better criterion is the magnitudeof the fourth time derivative of the pairwise force F (4) ij ,i. e. those particles are selected as neighbours which pro-duce the largest integration error in accordance with theHermite scheme. F (4) ij is a complicated expression (compareAppendix A), but the leading term can be estimated viadimensional analysis: F (4) ij ∝ m j v ij r ij (96)We use this expression to define a new apparent distancebetween the integrated particle i and a neighbour j : r app = r ij (cid:18) m i m j (cid:19) / (cid:18) v s v ij (cid:19) / (97) v s is an arbitrary scaling velocity to obtain a distance withdimension length. This new distance definition moves mas-sive or fast neighbours to an apparently smaller distance,thus enforcing that these particles are preferentially includedin the neighbour list. In addition, the modified distance isreadily included in the conventional neighbour scheme. Wetested different mass and velocity exponents to verify thatEq. 97 is the optimal choice. Figure 4 shows that these ex-ponents are indeed the optimal choice for a Plummer model(Plummer 1911) with mass spectrum. The new scheme saves25% of the force evaluations in the core, but the impact on ¯ r/r Hill < a planetesimals system is smaller, as it is the case for theneighbour scheme. While a velocity dependent distance re-duces the number of necessary full force evaluations, it in-troduces a distance changing with time which destabilisesthe integration. The result is a much larger energy errorcompared to the achieved speedup. Therefore we only rec-ommend the mass modification of the apparent distance. The rate at which the neighbours of a given particle changehas a noticeable influence on the accuracy of the code. Dur-ing the course of an integration the second and third timederivative of the regular and irregular force are calculatedfrom an interpolation formula (see Eq. 50 and Eq. 51).Whenever a particle leaves (or enters) the neighbour sphere,these derivatives are corrected by analytic expressions .Hence many neighbour changes lead to a pronounced spuri-ous difference.We estimated the rate at which particles cross theneighbour sphere boundary to quantify this effect. Neigh-bour changes are due to the Keplerian shear and the super-imposed random velocities of the particles. The two effectslead to: N + / − = ∆ t r N nb T orb Shear N + / − = 32 ∆ t r N nb σ v R nb = ∆ t r N nb T orb πσ v R nb Ω Dispersion (98) R nb and N nb are the neighbour sphere radius and theneighbour number, respectively. In practice, the neighbourchanges due to the shear account for up to 80% of thetotal neighbour changes. The standard regular time step∆ t r = 2 − and 50 neighbours yield a change of one par-ticle per regular step, which is fairly safe. Each integration step is preceded by the prediction of allneighbours of the particles that are due. A regular step re-quires the full prediction of all particles, so there is no possi-bility to save computing time. In contrast, an irregular stepcalculates only neighbour forces, which requires the predic-tion of less particles. Thus the prediction of all particlesto prepare an irregular step is a simple, but, depending onthe block size, computational costly solution. It seems to bemore efficient to predict only the required particles, but ran-dom access to the particle data and the complete check ofall neighbour list entries introduces an additional overhead.Therefore large block sizes should favour the first approach,whereas the second approach is more suitable for small blocksizes.Both regimes are separated by a critical block size N ∗ irr .If N irr particles with (cid:104) N nb (cid:105) neighbours are due, then only N merge particles need to be predicted: Appendix A gives a complete set of the force derivatives up tothird order. c (cid:13) , 1–42 ybrid methods: A new composite algorithm Process 0 1 2 3 4 5 6 7Send to 1 2 3 4 5 6 7 0Receive from 7 0 1 2 3 4 5 6
Table 2.
Ring Communication. Communication partners arefixed, while the exchanged data varies. n p − Table 3.
Hierarchical Communication. Communication partnerschange after every cycle. The exchanged data amount doubleswith every new cycle, hence only ln ( n p ) cycles are needed. N merge ≈ N tot (cid:18) − exp (cid:18) − N irr (cid:104) N nb (cid:105) N tot (cid:19)(cid:19) (cid:54) N irr (cid:104) N nb (cid:105) (99)The size N merge of the merged neighbour lists is smallerthan the total number of neighbour list entries, since someparticles are by chance members of more than one neighbourlist. Performance measurements show that the prediction ofthe merged neighbour lists is 10 % more costly (per particle)than the full prediction, mainly due to additional sorting anda random memory access. Thus N ∗ irr satisfies: N tot = 1 . × N merge (100)Inserting Eq. 99 yields the critical block size: N ∗ irr ≈ . N tot (cid:104) N nb (cid:105) (101)The prediction mode is chosen according to the actual blocksize. Nbody6++ is parallelised using a copy algorithm. A com-plete copy of the particle data is located on each node, sothe integration step of one particle does not need any com-munication. Therefore a block of N bl particles is dividedin n p parts ( n p is the processor number), which are inte-grated by different processors in parallel. The integrationstep is completed by an all–to–all communication of thedifferent subblocks to synchronise the particle data on allnodes. Hence the amount of communicated data is propor-tional to N bl × n p . A communication in a ring-like fashion(see table 2) needs n p − ( n p ) communication cycles. The differ-ence between the two approaches remains small, as long asthe communication is bandwidth limited, i. e. the blocks arelarge. Small block sizes shift the bottleneck to the latency,which is significantly reduced by the second scheme – espe-cially if the code runs on many processors. Block n p α [ µ s] τ l [ µ s] A BIrregular 10 0.35 51 145 4.5Regular 10 0.22 113 512 40Irregular 20 0.35 308 877 8.8Regular 20 0.22 368 1668 75 Table 4.
Timings on a Beowulf cluster (Hydra, see table 5). Seetext for an explanation of the variables. Timings are obtainedfor a maximal neighbour number LMAX=64. In practice, B istwice as large due to storage rearrangements in Nbody6++ . SeeAppendix 5 for details on the computers.Block n p α [ µ s] τ l [ µ s] A BIrregular 8 0.29 255 837 1.7Regular 8 0.60 981 1763 4.1Irregular 16 0.28 188 700 1.7Regular 16 0.60 306 561 6.7Irregular 64 0.27 241 887 7.7Regular 64 0.46 401 871 21.7 Table 6.
Timings on the IBM. More than 32 processors requiremore than one node.
A hierarchical scheme reduces the latency, but never-theless it is possible that the parallel integration is actuallyslower than a single CPU integration. We estimated boththe runtime on one CPU and on a parallel machine to ex-plore the transition between these two regimes. The latencytime τ l per communication is included in the wallclock timeexpressions for one regular/irregular step: τ l = αAt single = αN bl N nb (102) t par = α (cid:16) N bl N nb n p (cid:124) (cid:123)(cid:122) (cid:125) Arithmetic + A ln ( n p ) (cid:124) (cid:123)(cid:122) (cid:125) Latency + BN bl (cid:124) (cid:123)(cid:122) (cid:125) Communication (cid:17) (103)If t single (runtime on a single CPU) is equal to t par (parallelcomputation), one can deduce the critical block size N min which gives the minimal block size for efficient parallelisa-tion: t single = t par N min = A ln ( n p ) n p N nb ( n p − − Bn p (104)The hierarchical communication gives a minimal block sizethat increases logarithmically with the processor number.Eq. 103 gives immediately the speedup S and the optimalprocessor number for a certain block size N bl : S = n p An p ln ( n p ) N bl N nb + B n p N nb n p, opt ( N bl ) = ln(2) N bl × N nb A − Hierarchical (105)A comparison to the optimal processor number for a ring c (cid:13) , 1–42 P. Glaschke, P. Amaro-Seoane & R. Spurzem
Name Beowulf (Hydra) Titan JUMPInstitute ARI/ZAH ARI/ZAH Forschungszentrum J¨ulichLocation Heidelberg Heidelberg J¨ulichProcessors 20 64 1248Speed 2.2GHz 3.2GHz 1.7GHzProcessors/Node 2 2 32Network Myrinet Infiniband Gigabit–EthernetBandwidth 2Gbit/sec 20 Gbit/sec 10 Gbit/sec
Table 5.
Specs of the different supercomputers used in running the algorithm communication n p, opt ( N bl ) = (cid:114) N bl × N nb A − Ring t par = α (cid:16) N bl N nb n p + An p + BN bl (cid:17) (106)stresses the efficiency of the hierarchical communication,since it allows a much larger processor number for a givenproblem size. Equation 102 and 103 are also useful to derivethe total wallclock time, since the total runtime scales withthe number of regular and irregular blocks: N reg ≈ T N / N / √ η reg (107) N irr ≈ T N / √ η irr (108)These equations are only approximate expressions, but theygive the right order of magnitude without detailed calcula-tions that need a precise knowledge of the N –body model.Table 4 and table 6 summarise the timing parameters drawnfrom our experience with the Hydra and JUMP parallel su-percomputers. The preceding section showed that the block size is closelyrelated to the efficiency of the parallelisation. Small blocksare dominated by the latency and the parallelisation couldbe even slower than a single CPU calculation. Therefore wederive the block size distribution for the block time stepscheme to asses its influence on the efficiency.Suppose that the time steps h of all N particles in themodel are distributed according to some known function f : dp = f ( N, h ) dh (109) f is in most cases a complicated function. It involves spa-tial averaging and integration over the velocity distribution,which could be quite complicated even for simple time stepformulæ. Nevertheless there is a constraint on the time stepdistribution, simply because every particle has a neighbourwithin a finite distance: There is some upper limit h max that We use h instead of ∆ t in this section to avoid unclear nota-tion. hdpdh f h max Figure 5.
Timestep distribution f = dp/dh . The short-dashedline on the left indicates approximation Eq. 110, whereas thedashed line on the right defines a reasonable upper limit h max . restricts the major fraction of the time steps to a finite in-terval. Thus it is possible to capture the main features ofthe time step distribution with an expansion around h = 0(Fig. 5 sketches this approximation): f ≈ C ( N ) h a h (cid:54) h max (110) a is the lowest non-vanishing order of the expansion. Nowwe consider a block level with the largest possible time step h k . The number of particles N bl in this block is: N bl = N (cid:90) h k fdh ≈ C ( N ) a + 1 ( h k ) a +1 (111)According to the block time step scheme the number ofblocks per time with the largest possible time step h k is pro-portional to ( h k ) − . Therefore the probability that a blocksize is in the range [ N bl , N bl + dN bl ] is dp ∝ (cid:88) k δ (cid:18) N bl − C ( N ) a + 1 h a +1 k (cid:19) h k dN bl (112)where δ is Dirac’s delta function. The sum over the loga-rithmically equidistant time steps h k is approximated by an c (cid:13)000
Timestep distribution f = dp/dh . The short-dashedline on the left indicates approximation Eq. 110, whereas thedashed line on the right defines a reasonable upper limit h max . restricts the major fraction of the time steps to a finite in-terval. Thus it is possible to capture the main features ofthe time step distribution with an expansion around h = 0(Fig. 5 sketches this approximation): f ≈ C ( N ) h a h (cid:54) h max (110) a is the lowest non-vanishing order of the expansion. Nowwe consider a block level with the largest possible time step h k . The number of particles N bl in this block is: N bl = N (cid:90) h k fdh ≈ C ( N ) a + 1 ( h k ) a +1 (111)According to the block time step scheme the number ofblocks per time with the largest possible time step h k is pro-portional to ( h k ) − . Therefore the probability that a blocksize is in the range [ N bl , N bl + dN bl ] is dp ∝ (cid:88) k δ (cid:18) N bl − C ( N ) a + 1 h a +1 k (cid:19) h k dN bl (112)where δ is Dirac’s delta function. The sum over the loga-rithmically equidistant time steps h k is approximated by an c (cid:13)000 , 1–42 ybrid methods: A new composite algorithm p ( > N ) N/N tot
Plummer Model N=2500Theory a=2
Figure 6.
Cumulative irregular block size distribution for a N =2500 particle Plummer model. integral: dp ∝ (cid:90) ∞ δ (cid:18) N bl − C ( N ) a + 1 h a +1 (cid:19) d ln( h ) h dN bl ≈ a + 1 N − ( a +2) / ( a +1)bl dN bl (113)Thus the average block size and the median of the block sizedistribution are: (cid:104) N bl (cid:105) ≈ a N a/ ( a +1) median( N bl ) ≈ a +1 (114)Special expressions for the average block size were alreadyderived by Makino (1988), but the general relation of thetime step distribution to the block size distribution is a newresult. The median is surprisingly independent of the parti-cle number, i. e. 50 % of all blocks are always smaller thana fixed value. It seems that this is a threat to the efficiencyof the method, but the median of the wallclock timemedian( N ) ≈ N ( a +1) /a (115)demonstrates that these small blocks account only for asmall fraction of the total CPU time. We confirmed the de-rived block size distribution (Eq. 113) by numerical calcula-tions (see Fig. 6). The order parameter a is roughly two in(at least locally) homogenous systems, while an additionalKeplerian potential reduces the order to a = 1. A planetesi-mal disc – or more precisely, a narrow ring of planetesimals– has a very narrow distribution of time steps since all par-ticles share nearly the same orbital period. Thus the regularblock size is always equal to the total particle number mak-ing the parallelisation very efficient. We treated the mean neighbour number N nb so far as somefixed value. But it is also a mean to optimise the speed ofthe integration. Large neighbour spheres reduce fluctuationsin the regular forces allowing larger regular steps, whichreduces the total number of force evaluations. But largerneighbour lists also imply a larger communication overhead, N nb N N nb,opt N N Figure 7.
Optimal neighbour number as function of particlenumber N . The plot includes the numerical solution of Eq. 124and the two asymptotic solutions. Timing constants are takenfrom a Beowulf cluster. as all the neighbour lists have to be broadcast to synchro-nise the different nodes. The best choice balances these twoextremes, thus maximising the speed.Before we derive the optimal neighbour number on aparallel machine, we briefly summarise the known solutionfor a single CPU run (Makino 1988) for an extensive deriva-tion. The computational effort of the irregular steps is pro-portional to the neighbour number, while the number offorce evaluations for the regular steps is proportional to thetotal number of particles, reduced by the time step ratio γ t : γ t := ∆ t reg ∆ t irr T CPU = f ( N ) (cid:18) N nb + Nγ t ( N nb ) (cid:19) (116) γ t ( N nb ) ≈ N / (117) f ( N ) collects all factors depending only on the total num-ber of particles. Optimisation with respect to the neighbournumber N nb yields the well known result:0 = ddN nb T CPU N nb , opt ∝ N / (118)The calculation of the elapsed time for Nbody6++ on aPC cluster includes more terms. For clearness, we restrictourselves to a rather simple model that involves only thedominant terms to show how parallelisation influences theoptimal neighbour number. We make the following approx-imations:(i) We only take the force calculation and communicationinto account.(ii) We use the same time constants for regular and irreg-ular expressions.(iii) We neglect all numerical factors that are comparableto unity.The total CPU time is an extension of Eq. 106, which isapplied to the regular and the irregular step. A new constant B n includes the neighbour list communication separately, c (cid:13) , 1–42 P. Glaschke, P. Amaro-Seoane & R. Spurzem while all factors depending on N are represented by f ( N ): N bl ≈ N / γ t ≈ N / T irr = f ( N ) (cid:18) N bl N nb p + Ap + BN bl (cid:19) T reg = 1 γ t f ( N ) (cid:18) N bl Np + Ap + ( B + B n N nb ) N bl (cid:19) T tot = T irr + T reg (119)Optimisation with respect to the processor number p leadsto: 0 = ∂∂p T tot p opt = (cid:118)(cid:117)(cid:117)(cid:116) N / N / + N / A ( N / + 1) ≈ N / √ AN / (120)Further optimisation with respect to the neighbour numbergives the expression:0 = ∂∂N nb T tot (121)0 = N / − N − AN − / p − Bp +23 N nb B n p (122)For a fixed p or B n = 0 (very fast neighbour list communi-cation), we recover for large N : N nb , opt ∝ N / (123)In general, one can not neglect the neighbour list commu-nication. Therefore we seek for the optimal choice of p and N nb , thus combining Eq. 120 and Eq. 122: AN / + (cid:18) N nb − B B n (cid:19) B n √ AN / N / =13 ( N / + 1) NA (124)Since this equation has no closed solution, we identifythe dominant terms in Eq. 124 to calculate the asymptoticsolution for large N : N nb , opt ≈ (cid:18) A B n (cid:19) / N / N (cid:29) (cid:18) A B n (cid:19) / (125)For small N we get the approximated solution: N nb , opt ≈ (cid:18) N (cid:19) / N < (cid:18) A B n (cid:19) / (126)Fig. 7 compares the approximate expressions with the nu-merical solution of equation 124. In spite of the complicatedstructure of Eq. 124, both approximate expressions are re-liable solutions. The example uses timing constants derivedform our local Beowulf cluster: A ≈ B ≈ B n ≈ . N nb increases so slowly with the particlenumber that a neighbour number around 100 is a safe choice. The growth of planetesimals proceeds through collisionsamong planetesimals which form (at least in a sufficient frac-tion of incidents) larger bodies with a net gain of accretedmatter. But some collisions are mere destructive events thatshatter and disperse the colliding planetesimals. Small bod-ies are more susceptible to destruction, but they are alsodriven to high relative velocities due to the global energyequipartition making them even more vulnerable. A modelthat attempts to cover the full size range from objects rang-ing between a kilometer and the size of Mars needs a realisticcollision algorithm that covers both fragmentation and ac-cretion. Some examples in the literature are Cazenave et al.(1982); Beauge & Aarseth (1990). In our algorithm we usethe approach of Glaschke (2003), which was applied to as-teroid families. The following section describes the approachto fragmentation that is implemented in the code.
Two colliding bodies are equal in the sense that their intrin-sic properties are not different. Only the comparison of twobodies defines the larger body – usually denoted as target– and the smaller one denoted as projectile . The two termsstem from laboratory experiments where they indicate muchmore than different sizes. A small projectile is shot on atarget at rest to study the various parameters related tofragmentation. In the following, projectile and target onlyindicate the relative size of the two bodies.The collision of two bodies initiates a sequence of com-plex phenomena. Shock waves run through the material,flaws start to grow rapidly breaking the bodies apart inmany pieces. Some kinetic energy is transferred to the frag-ments, which leads to the ejection of fragments at differ-ent velocities in various directions. If the fragment cloud ismassive enough, some of the larger fragments may capturedebris. This post-collisional accretion is denoted as reaccu-mulation.Although the depicted scenario is quite complex, thereare a few measures that describe the most important as-pects:(i) Mass of the largest fragment M L , or dimensionless f l = M L /M where M is the combined mass of the twocolliding bodies.(ii) f l < refers to fragmentation, whereas f l > isdenoted as cratering.(iii) Energy per volume S that yields f l = is denotedas impact strength.(iv) f KE := 2 E fragkin /E kin : Fraction of the impact energythat is converted into kinetic energy of the fragments.Different fragment sizes and velocities are summarised byappropriate distribution functions. m i , D i and v i are mass,diameter and modulus of the velocity of a given fragment,respectively. c (cid:13)000
Two colliding bodies are equal in the sense that their intrin-sic properties are not different. Only the comparison of twobodies defines the larger body – usually denoted as target– and the smaller one denoted as projectile . The two termsstem from laboratory experiments where they indicate muchmore than different sizes. A small projectile is shot on atarget at rest to study the various parameters related tofragmentation. In the following, projectile and target onlyindicate the relative size of the two bodies.The collision of two bodies initiates a sequence of com-plex phenomena. Shock waves run through the material,flaws start to grow rapidly breaking the bodies apart inmany pieces. Some kinetic energy is transferred to the frag-ments, which leads to the ejection of fragments at differ-ent velocities in various directions. If the fragment cloud ismassive enough, some of the larger fragments may capturedebris. This post-collisional accretion is denoted as reaccu-mulation.Although the depicted scenario is quite complex, thereare a few measures that describe the most important as-pects:(i) Mass of the largest fragment M L , or dimensionless f l = M L /M where M is the combined mass of the twocolliding bodies.(ii) f l < refers to fragmentation, whereas f l > isdenoted as cratering.(iii) Energy per volume S that yields f l = is denotedas impact strength.(iv) f KE := 2 E fragkin /E kin : Fraction of the impact energythat is converted into kinetic energy of the fragments.Different fragment sizes and velocities are summarised byappropriate distribution functions. m i , D i and v i are mass,diameter and modulus of the velocity of a given fragment,respectively. c (cid:13)000 , 1–42 ybrid methods: A new composite algorithm Figure 8.
Section for f l = 0 .
04 and n = 3. The largest fragmentis coloured in dark-grey. In this calculation 60 × ×
60 grid cellsare used. Note the decomposition in grid cells and the Voronoypolyhedra which form the fragments. (i) Fragment size distribution:(a) N m ( m ) : Number of all fragments with a mass m i (cid:62) m ,(b) M ( m ) : Mass of all fragments with a mass m i (cid:62) m ,(c) N D ( D ) : Number of all fragments with a diameter D i (cid:62) D .The distribution functions are related to each other: N m ( m ) = N D ( D ( m )) M ( m ) = (cid:90) ∞ m x (cid:12)(cid:12)(cid:12)(cid:12) d N m ( x )d x (cid:12)(cid:12)(cid:12)(cid:12) dxN m ( m ) = (cid:90) ∞ m x (cid:12)(cid:12)(cid:12)(cid:12) d M ( x )d x (cid:12)(cid:12)(cid:12)(cid:12) dxD ( m ) is the size–mass relation.(ii) Velocity distribution:(a) ¯ v ( m ): mean velocity as a function of mass. Any theoretical or empirical prescription of a collision has torelate the afore mentioned parameters, namely the impactenergy, to the sizes and velocities of the produced fragments.The central quantity is the impact strength, which is a mea-sure for the overall stability of a body. Objects smaller than1 m are accessible to laboratory experiments, while collisionsof larger bodies up to asteroid size have to be analysed bycomplex computer simulations. Asteroid families, which areremnants of giant collisions in the asteroid belt, provide in-dependent insight, although the data is difficult to interpret. S [ J / m ] v [ m / s ] R[km] Housen 1990Benz 1999
Figure 9.
Impact strength according to Eq. 128 and Eq. 130.The right axis gives the corresponding impact velocity accordingto S = 1 / ρv with ρ = 2 . / cm . We selected two different impact strength models asreference for our work. The first was obtained by Housen(1990) through the combination of asteroid family data andlaboratory experiments via scaling laws: S = S (cid:18) R (cid:19) − . (cid:34) . × − (cid:18) R (cid:19) . (cid:35) f KE = 0 . S = 1 . × Jm − = 1 . × erg cm − (129)Later, Benz (1999) obtained another result throughSPH simulations (for basalt, v =3 km / s): S = S (cid:18) R (cid:19) − . (cid:34) . × − (cid:18) R (cid:19) . (cid:35) f KE ≈ .
01 (130) S = 6 . × Jm − (131) ρ = 2 . (132) f KE is a measure of the kinetic energy that is transferred tothe fragments: E fragkin = f KE E kin (133)We introduce a dimensionless measure γ of the relative im-portance of gravity for the result of a collision. It is definedas the ratio of the energy per volume S G that is necessaryto disperse the fragments to the impact strength S : S G = 2 π − √ f KE GR ρ γ := S G /S (134)The first step towards the prediction of a collisional out-come is to relate the impact energy and the impact strengthto ascertain the size of the largest fragment f l . Laboratory c (cid:13) , 1–42 P. Glaschke, P. Amaro-Seoane & R. Spurzem experiments and simulations indicate the functional form (cid:15) ( f l ) = (cid:26) − f l ) for f l > (2 f l ) − K otherwise. (cid:15) = E kin ρ SM (135)which is both valid in the fragmentation regime and the cra-tering limit. The size of the largest fragment is used to derivethe full size distribution. To accomplish the decomposition“seed fragments” are distributed inside the target accord-ing to the largest desired fragment. The full set of fragmentis derived from a Voronoy tessellation using these seedpoints. Fig. 8 depicts the result of such a decomposition.The fragment velocities are calculated from the total kineticenergy after the collision to initiate a post-collisional N –body calculation to treat reaccumulation.We conducted a large set of such calculations to cover asufficient range in f il (i. e. impact energy) and γ (i. e. bodysize). Table 7 summarises the derived values of the largestand second largest fragment including reaccumulation.
10 COLLISIONAL CASCADES
A first well-defined application of the fragmentation modelis a collisional cascade. The term cascade denotes that frag-ments of one collision in a many-body system may hit otherbodies, whose fragments further shatter even more bodies.Thus the particle number increases exponentially with everysubsequent collision.Although the formation of planets requires a net growthdue to collisions, this destructive process plays a role in theformation of larger bodies as the overall size distributioncontrols the accretion rate of the protoplanets. Therefore itis worth to have a closer look into this mechanism.
A system of colliding bodies is usually embedded in abroader context, like stars moving in a galaxy or asteroidsorbiting in our own solar system. First, we simplify this dy-namical background as well as some aspects of the collisionsto make the problem tractable.The first step is to decompose an inhomogeneous sys-tem into smaller subvolumes which are locally homogenous.Furthermore, it is assumed that these subvolumes hardlyinteract with each other. Hence it is possible to apply the particle–in–a–box–method (Safronov 1969) to analyse colli-sions within the small subvolumes:(i) All particles are contained in a constant volume.(ii) The particle sizes are described by a distribution func-tion n ( m ), i. e. the particle number per volume and massinterval.(iii) For convenience, we assume a constant (or typical)relative velocity for a given pair of colliding bodies. The Voronoy tessellation assigns every volume element to theclosest seed point. First applications date back to the 17th cen-tury, but the Russian mathematician of Ukrainian origin GeorgyFeodosevich Voronoy put it on a general base in 1908.
The distribution function is evolved by the coagulation equa-tion. We modified the equation given by Tanaka et al. (1996)by introducing a new function M red to arrive at a more con-cise expression:0 = ∂∂t mn ( t, m ) + ∂∂m F m ( t, m ) (136)The mass flux F m is given by: F m = − (cid:90) (cid:90) n ( t, m ) n ( t, m ) ξdm dm ξ ≡ σ ( m , m ) v rel M red ( m, m , m ) (137) M tot = (cid:90) n ( t, m ) mdm (138) ∂∂t M tot = F m ( m min ) , (139)where ξ is the coagulation kernel, n is the already intro-duced size distribution, v rel is the mean relative velocity, σ isthe cross section for colliding bodies ( m , m ) and M red isthe newly introduced fragment redistribution function. M red contains all information on the fragments arising from thebreakup of body m due to the impact of body m . Its defi-nition avoids double counting of collisions in the above inte-gral. The redistribution function is related to the differentialnumber distribution function n coll ( m , m , m ), i. e. the num-ber of fragments produced by a collision per mass interval.Since the target m formally disappears, it is included as anegative contribution: M red ( m, m , m ) := (cid:90) m (cid:0) n coll ( m , m , ˜ m ) − δ ( ˜ m − m ) (cid:1) ˜ md ˜ m (140)Mass conservation in each collision is reflected by M red (0 , m , m ) = M red ( ∞ , m , m ) = 0. The cross section σ depends on the velocities and radii R i of the particles. Asimple approach is the geometric cross section: σ ( m , m ) = π ( R + R ) (141)If gravity plays an important role during encounters, twocolliding bodies move on hyperbolic orbits with a pericentredistance that is smaller than the impact parameter. Thisleads to an additional enlargement of the cross section, de-noted as gravitational focusing: σ ( m , m ) = π ( R + R ) (cid:18) G ( m + m ) v ( R + R ) (cid:19) (142)A special class of collisional models are self-similar collisions.Self-similarity implies an invariance of the collisional out-come with respect to the scale of the colliding bodies. If thetarget mass as well as the projectile mass are enlarged by afactor of two, then only the masses of all fragments doubleswithout further changes in the collisional outcome. They al-low us to introduce the convenient dimensionless fragmentredistribution function f m : M red ( m, m , m ) = mf m ( m /m, m /m ) (143)We follow Tanaka et al. (1996) and employ the substitu-tion m = mx , m = mx to simplify Eq. 137: A similar approach to the solution of the coagulation equationis the Zakharov transformation, see Connaughton et al. (2004).c (cid:13) , 1–42 ybrid methods: A new composite algorithm Largest Fragment γ f il γ f il Table 7.
Data compilation of the fragmentation calculations. F m = − (cid:90) (cid:90) n ( t, mx ) n ( t, mx ) m / σ ( x , x ) v rel f m ( x , x ) dx dx (144)A simple solution is a steady-state cascade with F m = const.The loss of bodies of a given size is balanced by the frag-ment supply from larger bodies, hence the system maintainsa steady-state ∂∂t n ( t, m ) = 0. Eq. 144 inspires the ansatz n ( m ) ∝ m − k , which yields k = 11 /
6. This is the well knownequilibrium slope in self-similar collisional cascades, whichwas already found by Dohnanyi (1969). Strong gravitationalfocusing changes the exponent to k = 13 /
6. Both steady-state solutions seem to be rather artificial, as they containan infinite amount of mass and require a steady mass influxfrom infinity. However, they provide an appropriate descrip-tion for the relaxed fragment tail of a size distribution, aslong as the largest bodies provide a sufficient flux of newfragments. Once the largest bodies start to decay, the finiteamount of mass in the system leads to an overall decay of thecollisional cascade. Thus we seek for a more general solution Tanaka et al. (1996) state that k < to Eq. 136 using the ansatz n ( t, m ) = a ( t ) n ( m ): ∂∂t a ( t ) = − Ca ( t ) (145) mn ( m ) = 1 C ∂∂m F m (146) C is determined by fixing n at an arbitrary value m ∗ . a ( t )is independent of the collision model: a ( t ) = 11 + Ct (147) C ∝ n ( m ∗ ) (148)A power law solution is n ( m ) ∝ − Cm − k +1 which is onlyvalid for C <
C >
0, we perturb the already known equilibrium solution: n ( m ) = N m − k − CN m − k +2 + O ( C ) (149)1 N = (2 − k ) (cid:90) (cid:90) x − k x − k +22 σ ( x , x ) v rel ( f m ( x , x ) + f m ( x , x )) dx dx (150) N is small if the integral on the right hand side is large. Thisis the case for a sufficiently large impact strength. Eq. 149has the interesting property that n ( m (cid:48) ) = 0 for some mass m (cid:48) , given that k <
2. This mass m (cid:48) represents the largestbody in the system, e. g. the largest asteroid in a fictitiousasteroid belt. c (cid:13) , 1–42 P. Glaschke, P. Amaro-Seoane & R. Spurzem
11 SIZE–DEPENDENT STRENGTH
Self-similarity is an enormous help in analysing the coagu-lation equation. It releases completely the need to know anyspecific details of the collisional process and provides valu-able insight at the same time. But self-similarity is also astrong limitation on the underlying collisional physics.A major component of a fragmentation model is theknowledge of the impact strength as a function of size. Sim-ulations as well as asteroid families establish that it is notsome fixed value, but changes with size which immediatelybreaks the self-similarity. Larger bodies are weaker due toan increasing number of flaws (there are no big monocrys-tals), but then gravity leads to a turnover and increases thestrength.We model the size dependent strength S with a powerlaw to examine the influence on the equilibrium solution.The velocity dispersion v and the collisional cross section σ are also modelled with power laws to account for relaxationprocesses: v = v (cid:18) mm (cid:19) w (151) σ = σ (cid:18) mm (cid:19) s (152) S = S (cid:18) mm (cid:19) α (153)The subscript “0” denotes values for an arbitrarily chosenscaling mass. Since smaller bodies are more abundant thanlarger ones, we safely assume that most collisions involve alarge mass ratio. In addition, we assume w <
0, since weexpect energy equipartition to some degree in most cases.These restrictions lead to the following simplifications ( m >m ): σ ( m , m ) ≈ σ ( m ) (154) v rel ≈ v ( m ) (155) (cid:15) ≈ m ρv m S (156)Therefore the smaller body m enters only through the spe-cific energy (cid:15) : F m ≈ − (cid:90) (cid:90) n ( t, m ) n ( t, m ) σ ( m ) v rel ( m ) m f m ( m /m, (cid:15) ) dm dm (157)We introduce new dimensionless quantities with the help ofEq. 154–156 to simplify the integral: m = mx m = m (cid:18) m m (cid:19) α w (cid:18) S ρv (cid:19) w (cid:15) w (158)Again we assume a power law for the density n ∝ m − k andchange the integration parameters to ( x , (cid:15) ). Applying theconstant–flux condition yields the equilibrium exponent k ≈ s + 3 + α + w (2 s + α + 5)2 + α + 2 w (159) and the scaling exponent k (cid:48) of the total mass loss: k (cid:48) ≈ s − w + 12 + α + 2 w∂∂t M tot ∝ − ˜ S − k (cid:48) ˜ S = 2 S ρv (160)The exponent k (cid:48) in Eq. 160 is close to unity for realisticvalues of the free parameters. Thus the mass loss is roughlyinversely proportional to the strength of the bodies. The gen-eral formula Eq. 159 contains the special solution of O’Brien(2003), who concentrated on the parameters s = 2 / w = 0and a special collisional model. In fact, the derivation appliesto a much wider class of collisional models that we denote asscalable collisional models. Scalable indicates that the modelis self-similar except a scaling of the impactor mass.
12 PERTURBATION OF EQUILIBRIUM
The derived scaling relations provide insight into the over-all properties of a collisional cascade, which is in (or closeto) equilibrium. However, they do not provide informationon how the equilibrium is attained or how the system re-sponds to various external perturbations. A rigorous ap-proach would be the approximate solution of the coagulationequation , which is by no means simple since it requires acareful analysis of the collision model.Hence we turn to perturbations of the equilibrium sizedistribution, as it is easier to asses the quality of the derivedexpression for a variety of collision models. In addition, allequations are linear in the perturbation, allowing the de-tailed analysis of the solution.If the equilibrium solution n ( m ) = n ( m/m ) − k is per-turbed with a small deviation ∆ n ( m ), we get to first order:0 = ∂∂t m ∆ n ( m ) + ∂∂m F p ( t, m ) (161) F p = − (cid:90) (cid:90) ∆ n ( m ) n ( t, m ) σ ( m , m ) v rel × ( M red ( m, m , m ) + M red ( m, m , m )) dm dm Despite of the expansion in ∆ n , Eq. 161 is still a compli-cated integro-differential equation. Thus it is not possibleto obtain a solution without further information about theproblem. While there is no general solution, we restrict ourattention to self-similar collisional processes. In virtue of thisassumption it is possible to simplify Eq. 161, as we can see inEq. 162 and 162. In those expressions σ and v are velocityand cross section of an arbitrarily chosen scaling mass m . F ( x ) contains all information about the collisional process.If collisions do not result in extreme outcomes, like crateringor a complete destruction of the target, most of the fragmentmass is contained in bodies with similar size as the parentbody. Hence we expect that F ( x ) peaks around x ≈ x gets larger (or smaller). We introducethe dimensionless relative perturbation g ( m ): g ( m ) = ∆ n ( m ) n ( m ) = ∆ n ( m ) m k n m k (164) Appendix C highlights a possible approach.c (cid:13)000
The derived scaling relations provide insight into the over-all properties of a collisional cascade, which is in (or closeto) equilibrium. However, they do not provide informationon how the equilibrium is attained or how the system re-sponds to various external perturbations. A rigorous ap-proach would be the approximate solution of the coagulationequation , which is by no means simple since it requires acareful analysis of the collision model.Hence we turn to perturbations of the equilibrium sizedistribution, as it is easier to asses the quality of the derivedexpression for a variety of collision models. In addition, allequations are linear in the perturbation, allowing the de-tailed analysis of the solution.If the equilibrium solution n ( m ) = n ( m/m ) − k is per-turbed with a small deviation ∆ n ( m ), we get to first order:0 = ∂∂t m ∆ n ( m ) + ∂∂m F p ( t, m ) (161) F p = − (cid:90) (cid:90) ∆ n ( m ) n ( t, m ) σ ( m , m ) v rel × ( M red ( m, m , m ) + M red ( m, m , m )) dm dm Despite of the expansion in ∆ n , Eq. 161 is still a compli-cated integro-differential equation. Thus it is not possibleto obtain a solution without further information about theproblem. While there is no general solution, we restrict ourattention to self-similar collisional processes. In virtue of thisassumption it is possible to simplify Eq. 161, as we can see inEq. 162 and 162. In those expressions σ and v are velocityand cross section of an arbitrarily chosen scaling mass m . F ( x ) contains all information about the collisional process.If collisions do not result in extreme outcomes, like crateringor a complete destruction of the target, most of the fragmentmass is contained in bodies with similar size as the parentbody. Hence we expect that F ( x ) peaks around x ≈ x gets larger (or smaller). We introducethe dimensionless relative perturbation g ( m ): g ( m ) = ∆ n ( m ) n ( m ) = ∆ n ( m ) m k n m k (164) Appendix C highlights a possible approach.c (cid:13)000 , 1–42 ybrid methods: A new composite algorithm ∂∂t m ∆ n ( m ) − n m σ v ∂∂m (cid:90) ∆ n ( t, mx ) F ( x )( mx /m ) k dx (162) F ( x ) = (cid:90) m k − x − k x − k σ ( x , x ) σ v rel v ( f m ( x , x ) + f m ( x , x )) dx (163) G ( u ) S ~ k ’ u S~=1000S~=10S~=0.1S~=0.001 Figure 10.
Scaled fragmentation kernel G ( u ) for a simplefragmentation model (see Eq. 200) and different scaled impactstrength ˜ S . Thus the new differential equation reads:0 = ∂∂t ( m/m ) − k ∆ g ( m ) − n m σ v ∂∂m (cid:90) ∆ g ( t, mx ) F ( x ) dx (165)We change to logarithmic coordinates to arrive at a convo-lution integral: u = ln( m/m ) u = ln( x ) (166)Furthermore we define a collisional timescale τ τ = ( n m σ v ) − (167)to obtain a more concise expression. The transformed equa-tion is:0 = ∂∂t g ( t, u ) e u (2 − k ) − τ ∂∂u (cid:90) g ( t, u + u ) G ( u ) du (168) G ( u ) = F ( e u ) e u (169)If g ( u ) is varying on a scale larger than the width of thekernel G ( u ) (compare Fig. 10), it is justified to expand g ( u )under the integral. We retain the first two moments of G ( u ):0 = ∂∂t g ( t, u ) e u (2 − k ) − G τ ∂∂u g ( t, u ) − G τ ∂ ∂u g ( t, u )(170) G k = (cid:90) u k G ( u ) du (171)The first order moment G , which introduces a diffusiveterm, is omitted in the following for clarity . We introduce The study of wave-like structures in the size distribution(Bagatin et al. 1994,see e.g.) requires even the second order mo-ment G . a fragmentation time τ frag ( u ) and transform Eq. 170 backto m : 0 = ∂∂t g ( t, m ) − mτ frag ( m ) ∂∂m g ( t, m ) (172) τ frag = τ G e u (2 − k ) = τ G ( m/m ) − k (173)Eq. 172 is a modified advection equation, which conservesthe total mass. It is possible to derive equations similar toEq. 172 for any collisional model. However, the general ap-proach is less fruitful, as it lacks a robust frame of a knownequilibrium solution and reliable scaling relations. Thereforewe provide only the extension to scalable collisional modelsin Appendix B. We readily obtain the general solution: g ( t, m ) = f (cid:18) t + τ ( m/m ) (2 − k ) G (2 − k ) (cid:19) ∆ n ( t, m ) = n ( m ) f (cid:18) t + τ ( m/m ) (2 − k ) G (2 − k ) (cid:19) (174)The function f is determined by the initial value g (0 , m )of the perturbation. As the collisional cascade evolves, theinitial perturbation function is shifted as a whole to smallermasses. This evolution becomes clearer if we attach labels M (0) to the initial perturbation function and follow the timeevolution of these tags. The functions M ( t ) are the charac-teristics of the differential equation 172: M ( t ) = m (cid:104) ( M (0) /m ) (2 − k ) − t/τ G (2 − k ) (cid:105) / (2 − k ) (175)The meaning of the fragmentation time τ frag becomes clearby the relation M ˙ M = − τ frag (176)which is the time until a body has lost a significant fractionof its mass due to destructive collisions. A comparison of theperturbation equation 172 with the scaling relations fromthe previous section gives the scaling of the zeroth ordermoment G with respect to the impact strength: G = G (cid:48) ˜ S − k (cid:48) (177) G (cid:48) should only depend on the fragmentation model (i. e.fragment size distribution as a function of the largest frag-ment f l ) within the limits of this approximation. Fig. 10shows that the scaling with the impact strength works quitewell, except slight variations which are small compared tothe covered range of impact strengths. Likewise, it is pos-sible to restate the total equilibrium flux F eq in terms of G (cid:48) : F eq ( m ) ≈ − G (cid:48) n ( m ) σ ( m ) m v rel ˜ S − k (cid:48) (178) In general, characteristics of a partial differential equation arepaths along which the solution is constant.c (cid:13) , 1–42 P. Glaschke, P. Amaro-Seoane & R. Spurzem
The fragmentation timescale τ frag ( m ) allows a more intuitiveexpression: F eq ( m ) ≈ − n ( m ) m τ frag ( m ) (179)Our simple collisional model (see Fig. 10 and Eq. 200) refersto: F eq ( m ) = − (1 . . . × n ( m ) σ ( m ) m v rel ˜ S − k (cid:48) (180)
13 MIGRATION AND COLLISIONS
The local perturbation analysis is only applicable to a plan-etesimal disc, if the migration velocity of the planetesimalsis negligible small. This assures that collisional cascades atdifferent radial distances do not couple to each other, so thatthe whole disc is composed of many local cascades. Whilethis assumption is justified for larger bodies, migration isstrongly influencing bodies below 1 km in size. Hence weextend our analysis to examine the influence of migrationon the (no longer) local collisional processes.We assume that the collisional evolution of the systemleads to an equilibrium planetesimal distribution everywherein the disc: Σ ( r, m ) = Σ r, ( r ) C ( m ) (181)Σ r ( r ) is the total surface density at a given distance r , while C ( m ) is the universal equilibrium distribution. Though theplanetesimal distribution at larger sizes is likely different atdifferent locations in the disc, we only demand a universalfunction at smaller sizes, where migration is important. Thepower law exponent k depends on the details of the invokedphysics, but numerical simulations show that k ≈ r, m ) = g ( r, m )Σ ( r, m ) (182)where the dimensionless function g contains the changes dueto migration. The collisional evolution is governed by thecontinuity equation with an additional collisional term ∂ Σ( r, m ) ∂t − r ∂∂r ( v ( r, m ) r Σ( r, m )) = ˙Σ coll (183)where v ( r, m ) is the migration velocity (see Eq. 77), definedsuch that positive v imply an inward migration. We expressthe collisional term with the help of Eq. 172 and seek for asteady-state solution ˙Σ = 0:1 τ frag ( m, r ) ∂g∂m m Σ r, ( r ) + 1 r ∂∂r ( gvr Σ r, ( r )) = 0 (184) τ frag ( m, r ) is the fragmentation timescale of a mass m ata distance r . Since the surface density Σ and the variouscontributions to the drag force are well described by a powerlaw (with respect to radius), Eq. 184 further simplifies to:1 τ frag ( m, r ) ∂g∂m m + ∂g∂r v − br gv = 0 (185) b is a combination of the various invoked power lawexponents. As the surface density Σ and the gas densitydrop with increasing radius in any realistic disc model, it issafe to assume b >
0. We choose a self-similar ansatz for g : g ( r, m ) = g ( ζ ) , ζ = mg m ( r ) (186) g ζ k r =-1k r =0k r =1 Figure 11.
Cut-off function g according to Eq. 192. The massexponent is k m = 1 /
3, while the mass influx exponent is b = 1 . The new differential equation is1 τ frag ( m, r ) dgdζ mg m ( r ) + m dgdζ dg m dr v − br gv = 0 (187)which is equivalent to the more concise expression: d ln( g ) d ln( ζ ) (cid:18) rvτ frag + d ln( g m ) d ln( r ) (cid:19) = b (188)We assume a power-law dependence for the timescale ratio τ mig /τ frag : rvτ frag = τ mig τ frag = ( m/m ) k m ( r/r ) k r (189)The cut-off mass m at a distance r has a timescale ratio τ mig /τ frag = 1, which defines a proper lower cut-off withinthis context. Hence the solution is: b = d ln( g ) d ln( ζ ) (cid:18) ζ k m + k r k m (cid:19) (190) g m ( r ) = ( r/r ) k r /k m m (191) g ( ζ ) = (cid:18) k r k m ζ k m (cid:19) − b/k r (192)Though the analytical solution Eq. 192 provides a completedescription of the lower cut-off of the size distribution, itis more appropriate within the frame of this discussion totranslate the equilibrium solution to an equilibrium massloss due to migration:˙Σ mig ( r, m ) = − bvr Σ + bvr Σ k r /k m ζ k m + k r /k m = − b Σ τ mig + k r /k m τ frag (193)An inspection of the timescale ratio shows that the massexponent k m should be positive, whereas simple estimationsof k r on the basis of the minimum mass solar nebula aresomewhat inconclusive. The value of k r is so close to zerothat any change in the assumed equilibrium slope or theimpact strength scaling gives easily both positive and neg-ative values. Moreover, Eq. 193 requires a globally relaxedplanetesimal disc, but the huge spread in the various in-volved timescales at different radii inhibits any significantrelaxation in the early stages. c (cid:13) , 1–42 ybrid methods: A new composite algorithm However, it is possible to gain valuable information fromthe two limiting cases k r > k r <
0. Both values of k r give the proper limit g → k r reduces the effective mass lossdue to migration, as fragments from the outer part of thedisc replenish the local mass loss. Hence the fragmentationtimescale controls the net loss of smaller planetesimals. Incontrast, a negative exponent k r leads to a pronounced cut-off in the size distribution, since only larger planetesimalsare replenished through inward migration. Though the massloss rate is singular at some mass m (cid:48) , this sharp cut-off isan artifact due to the perturbation approximation.Our analysis is subjected to several restrictions. We ap-plied the perturbation equation to values of g that exceedthe limit for a safe application (i. e. g (cid:54)≈
1) of the perturba-tion expansion. Furthermore, the steady-state solution re-quires a global relaxation of the collisional processes, whichis practically never obtained during the disc evolution. De-spite of these restrictions, we gained insight on a more qual-itative level. Numerical calculations indicate that the per-turbation approximation is inappropriate close to the lowercut-off of the size distribution. However, a comparison of dif-ferent exponents k r (see Fig. 11) attributes only a minor roleto the replenishment of fragments due to inward migration.Only unrealistic small slopes b of the migrational mass influxwould strengthen the importance of this process. Thoughtemporally non-equilibrium phenomena are not ruled outby the previous derivation, their study would require theglobal simulation of the system.
14 COAGULATION
While most coagulation kernels are only restricted to a lim-ited analytical analysis (e. g. scaling relations), there existsome special kernels that allow the closed solution of thecoagulation equation. All rely on the assumption of perfectmergers, which allows the reformulation of the general equa-tion 136 to ∂n ( m, t ) ∂t = 12 (cid:90) m K ( m − m (cid:48) , m (cid:48) ) n ( m − m (cid:48) , t ) n ( m (cid:48) , t ) dm (cid:48) − n ( m, t ) (cid:90) ∞ K ( m, m (cid:48) ) n ( m (cid:48) , t ) dm (cid:48) where K is the coagulation kernel. One of these particularkernels was introduced by Safronov (1969): K ( m , m ) = A ( m + m ) (194)This coagulation kernel implies perfect mergers, where thecoalescence rate of two particles m and m is assumed to beproportional to the sum of their masses. It seems that thisis an artificial choice, devised to allow an analytic solution.However, the Safronov cross section provides an intermedi-ate case between a geometric cross section ( σ ∝ m / ) andstrong gravitational focusing ( σ ∝ m / ). A special solutionto the initial condition n ( m,
0) = n m exp ( − m/m ) (195) ρ k / f KE K Table 8.
Main parameters of the collisional model. is the function (Ohtsuki et al. 1990,see e.g.) n ( m, τ ) = n ˜ gm √ − ˜ g exp( − (2 − ˜ g ) m/m ) I (2 m/m (cid:112) − ˜ g )(196)˜ g = exp( − τ ) τ = A ρtρ = (cid:90) ∞ mn ( m ) dm = n m (197)where τ is the dimensionless time and I is a modified Besselfunction of the first kind.
15 MODELS FOR M RED
Though we already obtained insight into the nature of col-lisional cascades without a detailed specification of the co-agulation kernel, any detailed study of a collisional systemrequires the specification of a realistic collisional model.First, we restate the well-known perfect accretionmodel. While it is a gross oversimplification for collisionsamong kilometre–sized planetesimals, its simplicity allows areliable code testing and eases the comparison with otherworks: M red ( m, m , m ) = − m Θ( m − m ) − m Θ( m − m )+( m + m )Θ( m − m − m ) (198)Although our fragmentation model (see section 9) pro-vides a very detailed description of the outcome of a colli-sion, we abandon most of the details for the following rea-sons. The computational effort of the numerical solution ofthe coagulation equation scales with the third power of thenumber of mass bins. Hence we chose a mass grid whose res-olution is by far smaller than the information provided bythe detailed collisional model. As a mismatch of the massresolution could produce undesired artifacts, a lower resolu-tion of the collisional model is needed for consistence. Thusonly the largest fragment f l ( f il , γ ) and the second fragment f (2) l ( f il , γ )(which contains information on reaccumulation)enter the fragment size distribution: M red ( xM ) M = x (cid:62) f l − f l if f l > x (cid:62) f (2) l (1 − f l − f (2) l )( x/f (2) l ) f l otherwise (199)Both values f l and f (2) l are interpolated from table 7, wherethe initial fragment size f il is calculated from the dimen-sionless impact energy (cid:15) . We used a reduced fragmentationmodel for test purposes: M red ( xM ) /M = (cid:26) x (cid:62) f l (1 − f l )( x/f l ) f l otherwise (200) c (cid:13) , 1–42 P. Glaschke, P. Amaro-Seoane & R. Spurzem
Table 8 summarises the most important model parameters.
16 STATISTICAL MODEL
The direct approach to the integration of an N –body systemis, in principle, possible for any particle number. While thisprocedure becomes computationally too expensive for verylarge particle numbers, a by far more efficient approach isapplicable in this regime. Instead of tracking all particle or-bits, a distribution function f (also phase-space density ),which gives the probability to find a particle at a position x with a velocity v , contains the state of the system: dp = f ( x , v ) d xd v (201)As long as only dynamical interactions are taken into ac-count, the number of all particles (e. g. stars, planetesimals)is conserved. The continuity equation reads:0 = ∂f∂t + v · ∇ f − ∇ Φ · ∂f∂ v (202)This is the collisionless Boltzmann equation. Collisions leadto an additional term (cid:18) ∂f∂t (cid:19) coll = ∂f∂t + v · ∇ f − ∇ Φ · ∂f∂ v (203)which will be discussed later. f is a function of six variables,so an exact solution is usually very complicated or even im-possible. However, it is possible to gain valuable insight intothe problem by taking the moments of the distribution func-tion: (cid:104) x ni v mj (cid:105) = (cid:90) f ( x , v ) x ni v mj d xd v n, m > ν ( x ) = (cid:90) f ( x , v ) d v (205)Integration of Eq. 203 over all velocities yields the corre-sponding continuity equation: ∂ν∂t + ∂ν ¯ v i ∂x i = (cid:18) ∂ν∂t (cid:19) coll (206)The first order moment with respect to velocity gives thetime evolution of the mean velocity ¯ vν ∂ ¯ v j ∂t + ν ¯ v i ∂ ¯ v j ∂x i = − ν ∂ Φ ∂x j − ∂ ( νσ ij ) ∂x i + ν (cid:18) ∂ ¯ v j ∂t (cid:19) coll (207)¯ v i = 1 ν (cid:90) f ( x , v ) v i d vσ ij = v i v j − ¯ v i ¯ v j (208)where σ ij is the anisotropic velocity dispersion and the con-tinuity equation was used to arrive at a more concise for-mulation. Equations 206 and 207 are the Jeans equations.While the structure of the moment equations is already fa-miliar from hydrodynamics, they do not provide a closed setof differential equations, since each differential equation ofa given moment is related to (yet unknown) higher ordermoments. Hence any finite set of momenta needs a closurerelation – additional constraints that relate the highest ordermoments to known quantities. The choice of this relation isa key element in the validity of the equations, but it is not unique and depends well on the problem at hand (Larson1970,compare e.g.).Owing to the geometry of a planetesimal disc, it is usefulto express the Boltzmann equation in cylindrical coordinates ∂f∂t + v r ∂f∂r + v z ∂f∂z + (cid:32) v φ r − ∂ Φ ∂r (cid:33) ∂f∂v r − v r v φ r ∂f∂v φ − ∂ Φ ∂z ∂f∂v z = 0 (209)where all derivatives with respect to φ have been droppeddue to the assumed axisymmetry of the disc. Any statistical description of a planetesimal disc requiresthe knowledge of the distribution function. Since the fullproblem including collisions, encounters and gas drag hasno analytic solution, a collisionless planetesimal disc (i. e. noperturbations) is a natural basis for further investigations.The distribution function that describes such a simplifiedsystem is a solution of the Boltzmann equation. A specialsolution to Eq. 209 is a thin homogenous planetesimal disc f ( z, v ) = ΩΣ2 π T r T z m exp (cid:32) − v r + 4 v φ T r − v z + Ω z T z (cid:33) (210)provided that the radial velocity dispersion T r and the ver-tical dispersion T z are small compared the mean orbital ve-locity v K . The azimuthal velocity dispersion T φ is locked to T r by the local epicyclic frequency κ in a central potential,where the ratio 1 : 4 is a special solution of (Binney 1994,seee.g) κ T r = 4Ω T φ (211)All velocities v r , v φ and v z refer to the local Keplerian ve-locity. The normalisation is the same as in Stewart (2000): (cid:90) d vdzf ( z, v ) = Σ m (212)A planetesimal disc is a slowly evolving system compared tothe orbital time, hence it is reasonable to use Eq. 210 as ageneral solution of the perturbed problem. Σ, T z and T r arenow functions of time and of the radial distance to the star.All information on the system is contained in these threemomenta of the distribution function, where higher ordermoments can be deduced from Eq. 210. Thus the functionalform of the distribution function represents an implicit clo-sure relation.The validity of this approximation can be further as-sessed by a closer examination of the Boltzmann equation.We summarise all perturbations in an evolution timescale T evol and reduce the radial structure to some typical lengthscale ∆ r to estimate the deviation from the functional formEq. 210. A comparison with Eq. 209 shows that the differ-ence is small if the migration timescale and the evolutiontimescale are large compared to the orbital time T : T (cid:28) ∆ r/ (cid:104) v r (cid:105) (213) T (cid:28) T evol (214) c (cid:13)000
The direct approach to the integration of an N –body systemis, in principle, possible for any particle number. While thisprocedure becomes computationally too expensive for verylarge particle numbers, a by far more efficient approach isapplicable in this regime. Instead of tracking all particle or-bits, a distribution function f (also phase-space density ),which gives the probability to find a particle at a position x with a velocity v , contains the state of the system: dp = f ( x , v ) d xd v (201)As long as only dynamical interactions are taken into ac-count, the number of all particles (e. g. stars, planetesimals)is conserved. The continuity equation reads:0 = ∂f∂t + v · ∇ f − ∇ Φ · ∂f∂ v (202)This is the collisionless Boltzmann equation. Collisions leadto an additional term (cid:18) ∂f∂t (cid:19) coll = ∂f∂t + v · ∇ f − ∇ Φ · ∂f∂ v (203)which will be discussed later. f is a function of six variables,so an exact solution is usually very complicated or even im-possible. However, it is possible to gain valuable insight intothe problem by taking the moments of the distribution func-tion: (cid:104) x ni v mj (cid:105) = (cid:90) f ( x , v ) x ni v mj d xd v n, m > ν ( x ) = (cid:90) f ( x , v ) d v (205)Integration of Eq. 203 over all velocities yields the corre-sponding continuity equation: ∂ν∂t + ∂ν ¯ v i ∂x i = (cid:18) ∂ν∂t (cid:19) coll (206)The first order moment with respect to velocity gives thetime evolution of the mean velocity ¯ vν ∂ ¯ v j ∂t + ν ¯ v i ∂ ¯ v j ∂x i = − ν ∂ Φ ∂x j − ∂ ( νσ ij ) ∂x i + ν (cid:18) ∂ ¯ v j ∂t (cid:19) coll (207)¯ v i = 1 ν (cid:90) f ( x , v ) v i d vσ ij = v i v j − ¯ v i ¯ v j (208)where σ ij is the anisotropic velocity dispersion and the con-tinuity equation was used to arrive at a more concise for-mulation. Equations 206 and 207 are the Jeans equations.While the structure of the moment equations is already fa-miliar from hydrodynamics, they do not provide a closed setof differential equations, since each differential equation ofa given moment is related to (yet unknown) higher ordermoments. Hence any finite set of momenta needs a closurerelation – additional constraints that relate the highest ordermoments to known quantities. The choice of this relation isa key element in the validity of the equations, but it is not unique and depends well on the problem at hand (Larson1970,compare e.g.).Owing to the geometry of a planetesimal disc, it is usefulto express the Boltzmann equation in cylindrical coordinates ∂f∂t + v r ∂f∂r + v z ∂f∂z + (cid:32) v φ r − ∂ Φ ∂r (cid:33) ∂f∂v r − v r v φ r ∂f∂v φ − ∂ Φ ∂z ∂f∂v z = 0 (209)where all derivatives with respect to φ have been droppeddue to the assumed axisymmetry of the disc. Any statistical description of a planetesimal disc requiresthe knowledge of the distribution function. Since the fullproblem including collisions, encounters and gas drag hasno analytic solution, a collisionless planetesimal disc (i. e. noperturbations) is a natural basis for further investigations.The distribution function that describes such a simplifiedsystem is a solution of the Boltzmann equation. A specialsolution to Eq. 209 is a thin homogenous planetesimal disc f ( z, v ) = ΩΣ2 π T r T z m exp (cid:32) − v r + 4 v φ T r − v z + Ω z T z (cid:33) (210)provided that the radial velocity dispersion T r and the ver-tical dispersion T z are small compared the mean orbital ve-locity v K . The azimuthal velocity dispersion T φ is locked to T r by the local epicyclic frequency κ in a central potential,where the ratio 1 : 4 is a special solution of (Binney 1994,seee.g) κ T r = 4Ω T φ (211)All velocities v r , v φ and v z refer to the local Keplerian ve-locity. The normalisation is the same as in Stewart (2000): (cid:90) d vdzf ( z, v ) = Σ m (212)A planetesimal disc is a slowly evolving system compared tothe orbital time, hence it is reasonable to use Eq. 210 as ageneral solution of the perturbed problem. Σ, T z and T r arenow functions of time and of the radial distance to the star.All information on the system is contained in these threemomenta of the distribution function, where higher ordermoments can be deduced from Eq. 210. Thus the functionalform of the distribution function represents an implicit clo-sure relation.The validity of this approximation can be further as-sessed by a closer examination of the Boltzmann equation.We summarise all perturbations in an evolution timescale T evol and reduce the radial structure to some typical lengthscale ∆ r to estimate the deviation from the functional formEq. 210. A comparison with Eq. 209 shows that the differ-ence is small if the migration timescale and the evolutiontimescale are large compared to the orbital time T : T (cid:28) ∆ r/ (cid:104) v r (cid:105) (213) T (cid:28) T evol (214) c (cid:13)000 , 1–42 ybrid methods: A new composite algorithm An order–of–magnitude estimate of the evolution time sup-ports condition 213 and 214. Furthermore, numerical calcu-lations confirm that the velocity distribution stays triaxialGaussian (see Ida 1992).The distribution function is equivalent to an isothermalvertical density structure with scale height h : h = (cid:114) T z Ω (215) ρ ( z ) = ρ exp (cid:18) − z h (cid:19) (216)Thus the central density ρ and the mean density (cid:104) ρ (cid:105) arerelated to the surface density in a simple way: ρ = Σ √ πh (cid:104) ρ (cid:105) = ρ √ e and i : dn ( e , i ) = 1 (cid:104) e (cid:105)(cid:104) i (cid:105) exp (cid:18) − e (cid:104) e (cid:105) − i (cid:104) i (cid:105) (cid:19) de di (cid:104) e (cid:105) = 2 T r (Ω r ) (cid:104) i (cid:105) = 2 T z (Ω r ) (218)Planetesimal encounters couple the time evolution of eccen-tricity and inclination, so that the ratio i /e tends to anequilibrium value after a few relaxation times. It is closeto 1 / Planetesimal–planetesimal scatterings change the velocitydistribution through two different processes. Firstly, it isunlikely that two planetesimals scatter each other on circu-lar orbits. Thus we expect a steady increase of the velocitydispersion due to this viscous stirring. Secondly, encountersbetween unequal masses lead successively to energy equipar-tition, slowing down the larger bodies through dynamicalfriction. The later mechanism is not related to the disc ge-ometry at all, but operates in any multi-mass system. Aspecial case is the systematic deceleration of a massive body M in a homogeneous sea of lighter particles m with density n , which is given by the Chandrasekhar dynamical frictionformula (Chandrasekhar 1942) d v M dt = − v M π ln Λ G ( M + m ) n mv M (cid:18) erf( X ) − X √ π e − X (cid:19) X = v M √ σ v (219)where σ v is the velocity dispersion of the lighter particles.The Coulomb logarithm Λ arises from an integration overall impact parameters smaller than an upper limit l max andis given by Λ ≈ σ v l max G ( m + M ) (220) Eq. 10–12 provide the coordinate transformation.
Although encounters in the gravitational field of the sun de-viate from pure two-body scatterings, it is safe to neglectthe presence of the sun if the encounter velocity is largecompared to the Hill velocity Ω R Hill . Thus the classicaldynamical friction formula is also applicable to planetesimalencounters in the high velocity regime, though a generalisa-tion to triaxial velocity distributions σ i is necessary (Binney1977,see e.g.): dv M,i dt = − v M,i √ πG ln(Λ)( M + m ) n mB i (221) B i = (cid:90) ∞ exp (cid:18) − (cid:88) v j σ j + u (cid:19) du (cid:112) ( σ + u )( σ + u )( σ + u )( σ i + u ) (222)An additional complication is the choice of l max (i. e. thechoice of the Coulomb logarithm). There are several scalelengths, which could determine the largest impact parame-ter l max : The scale height of the planetesimal disc, the radialexcursion due to the excentric motion of the planetesimalsand the Hill radius of the planetesimals. As it is not possibleto derive a unique expression for l max from first principles, aproper formula is often fitted to N –body calculations (com-pare Eq. 235). The velocity dispersion of a planetesimal discis triaxial with T φ /T r = 1 / T z /T r ≈ /
4. We takethese values and expand Eq. 221 for small velocities v M : dv M,r dt ≈ − . v M,r √ πG ln(Λ)( M + m ) n mT / r dv M,φ dt ≈ − . v M,φ √ πG ln(Λ)( M + m ) n mT / r dv M,z dt ≈ − . v M,z √ πG ln(Λ)( M + m ) n mT / r (223)The derived expressions provide a compact tool to analysedynamical friction in disc systems. However, the involvedapproximations are too severe compared to the needs ofan accurate description. While these concise expressions arevaluable for basic estimations, the following sections deriveviscous stirring and dynamical friction formulæ for a plan-etesimal system in a rigorous way. We return to the Boltzmann equation as a starting point forthe derivation of the scattering coefficients: (cid:18) ∂f∂t (cid:19) coll = ∂f∂t + v · ∇ f − ∇ Φ · ∂f∂ v (224)In virtue of the ansatz for the distribution function (seeEq. 210), it is sufficient to derive the time derivative of thesecond order velocity moments T r and T z . Since the distri-bution function is time independent in the absence of en-counters, only the collisional term contributes to the timederivative of the velocity dispersions T k ( k ∈ ( r, z, φ ) in the Whenever relative velocities are classified as “high” or “low”in the following sections, a comparison with the Hill velocity isimplied.c (cid:13) , 1–42 P. Glaschke, P. Amaro-Seoane & R. Spurzem following): dρT k dt = (cid:90) d vmv k (cid:18) ∂f∂t (cid:19) coll (225)The collisional term invokes the averaging over many differ-ent scattering trajectories and is, given that the underlyingencounter model is analytically solvable, still too complexto derive an exact expression. If most of the encounters areweak – a realistic assumption in a planetesimal disc – it ispossible to expand the collisional contribution in terms ofthe velocity change ∆ v i . This is the Fokker-Planck approx-imation (see e.g. Binney 1994)) (cid:18) ∂f∂t (cid:19) coll = − (cid:88) i ∂∂v i [ fD (∆ v i )]+12 (cid:88) i,j ∂ ∂v i ∂v j [ D (∆ v i , ∆ v j )] (226)where the diffusion coefficients D contain all information onthe underlying scattering process. Next we consider two in-teracting planetesimal populations m, m ∗ with distributionfunctions f = ΩΣ2 π T r T z m exp (cid:32) − v r + 4 v φ T r − v z + Ω z T z (cid:33) f ∗ = ΩΣ ∗ π T ∗ r T ∗ z m ∗ exp (cid:32) − v r + 4 v φ T ∗ r − v z + Ω z T ∗ z (cid:33) (227)to evaluate the terms in equation 226. We follow Stewart(2000) except some minor changes in the notation. The col-lisional term requires an averaging over the velocities of thetwo interacting planetesimals m and m ∗ : d (cid:104) ρv k (cid:105) dt = 2 πG mm ∗ (cid:90) d v (cid:90) d v ∗ ff ∗ × (cid:20) − A ( m + m ∗ ) u k v k m ∗ u + Bu + (2 C − B )3 u k ) u (cid:21) u k = v k − v ∗ k A = ln(Λ + 1) C = Λ Λ + 1 B = A − C (228)A coordinate transformation to the relative velocity u andthe modified centre-of-mass velocity ww k = V k + ( m ∗ T ∗ r − mT r ) u k ( m + m ∗ )( T r + T ∗ r ) for k ∈ { r, φ } V k + ( m ∗ T ∗ z − mT z ) u k ( m + m ∗ )( T z + T ∗ z ) for k = z V = m v + m ∗ v ∗ m + m ∗ (229)further simplifies the double integral. Thus the integrationseparates in a simple integral over w and a more demanding u –integration: d (cid:104) ρv r (cid:105) dt = 2 πG mm ∗ (cid:90) d w (cid:90) d uff ∗ × (cid:20) A ( m ∗ T ∗ r − mT r ) u r m ∗ ( T r + T ∗ r ) u + B ( u − u r ) u (cid:21) d (cid:104) ρv φ (cid:105) dt = 2 πG mm ∗ (cid:90) d w (cid:90) d uff ∗ × (cid:34) A ( m ∗ T ∗ r − mT r ) u φ m ∗ ( T r + T ∗ r ) u + B ( u − u φ ) u (cid:35) d (cid:104) ρv z (cid:105) dt = 2 πG mm ∗ (cid:90) d w (cid:90) d uff ∗ × (cid:20) A ( m ∗ T ∗ z − mT z ) u z m ∗ ( T z + T ∗ z ) u + B ( u − u z ) u (cid:21) All integrals are solvable and give the result d (cid:104) v r,φ (cid:105) dt = G ρ ∗ √ T r + T ∗ r ) / [ B ( T ∗ r + T r ) m ∗ J r,φ ( β ) + 2 A ( T ∗ r m ∗ − mT r ) H r,φ ( β )] d (cid:104) v z (cid:105) dt = G ρ ∗ √ T r + T ∗ r ) / ( T z + T ∗ z )[ B ( T ∗ z + T z ) m ∗ J z ( β ) + 2 A ( T ∗ z m ∗ − mT z ) H z ( β )] β := T z + T ∗ z T r + T ∗ r (230)where six auxiliary functions are introduced to arrive at amore compact notation: a = (cid:112) − x b = (cid:112) − (1 − β ) x H r := 8 √ π (cid:90) x ab dxH φ := 8 √ π (cid:90) − x a ( βa + b ) dxH z := 8 √ π (cid:90) β (1 − x ) b ( βa + b ) dxJ r := − H r + H φ + H z J φ := H r − H φ + H z J z := H r + H φ − H z (231)Since these are non-trivial functions, we apply a standardChebyshev approximation for β ∈ [0 , f ( x ) ≈ (cid:88) k =0 c k T k ( x ) − c (232)Table 9 summarises the Chebyshev coefficients. A final z –averaging yields the expressions: d (cid:104) v r,φ (cid:105) dt = G ΩΣ ∗ √ π ( T r + T ∗ r ) / ( T z + T ∗ z ) / × (233)[ B ( T ∗ r + T r ) m ∗ J r,φ ( β ) + 2 A ( T ∗ r m ∗ − mT r ) H r,φ ( β )] d (cid:104) v z (cid:105) dt = G ΩΣ ∗ √ π ( T r + T ∗ r ) / ( T z + T ∗ z ) / × (234)[ B ( T ∗ z + T z ) m ∗ J z ( β ) + 2 A ( T ∗ z m ∗ − mT z ) H z ( β )]The determination of a proper Coulomb logarithm Λ leavesroom for further optimisation. A careful comparison with N –body models gives rise to the empirical choice (Ohtsuki c (cid:13) , 1–42 ybrid methods: A new composite algorithm f J r J φ J z H r H φ H z c -10.34660733 1.81674741 8.52985992 11.00434580 6.94989422 4.71219005 c c -1.25533220 -1.18724874 2.44258094 0.60969641 0.58700192 -0.62294130 c c -0.07040537 -0.11070339 0.18110876 0.03112047 0.04455314 -0.05271757 c Table 9.
Chebyshev coefficients of the auxiliary functions J k and H k . et al. 2002): Λ = 112 ( (cid:104) ˜ e (cid:105) + (cid:104) ˜ i (cid:105) ) (cid:104) ˜ i (cid:105) / (235)˜ e = √ T r Ω R Hill ˜ i = √ T z Ω R Hill (236)Ohtsuki et al. 2002 also report a further improvement bysetting B ≡ A . Encounters in the low velocity regime exhibit a wealth of dif-ferent orbits, as the solar gravity field perturbs the two-bodyscattering. Only a small subset of the trajectories representssimple, regular orbits like Tadpole or Horseshoe orbits .Hence an examination of this velocity regime is done bestwith a numerical study of the parameter space by integrat-ing the equations of motions numerically (see Eq. 16).Ohtsuki et al. (2002) integrated a large set of planetesi-mal encounters and extracted fitting formulæ that cover thelow velocity regime. Their expressions for viscous stirringare: dT r dt = Gr Ω h Σ ∗ m + m ∗ ) 73 C m ∗ dT z dt = Gr Ω h Σ ∗ m + m ∗ ) C m ∗ (cid:18) (cid:104) ˜ i (cid:105) + 0 . (cid:104) ˜ e (cid:105) ) / (cid:113) (cid:104) ˜ i (cid:105) (cid:19) ˜ e = e/h ˜ i = i/hC := ln(10Λ / ˜ e + 1)10Λ / ˜ e C := ln(10Λ √ ˜ e + 1)10Λ √ ˜ e (237)The stirring rate of the radial velocity dispersion approachesa finite value for very low velocity dispersions, while the stir-ring rate for the vertical velocity dispersion drops to zero asthe velocity dispersion decreases. This different behaviour ofthe two limits is due to the encounter geometry: If two plan-etesimals have zero inclination, they may still excite highereccentricities during an encounter, but they remain confined The most famous example of such a regular orbit are the twosaturnian moons Janus and Epimetheus which share nearly thesame orbit. to the initial orbital plane preventing any excitation of in-clinations.The respective expressions for the dynamical frictionrates are: dT r dt = Gr Ω h Σ ∗ m + m ∗ )( T r + T ∗ r ) 10 C (cid:104) ˜ e (cid:105) ( T ∗ r m ∗ − T r m ) dT z dt = Gr Ω h Σ ∗ m + m ∗ )( T z + T ∗ z ) 10 C (cid:104) ˜ i (cid:105) ( T ∗ z m ∗ − T z m ) C := ln(10Λ + 1)10Λ (238)As the stirring rates are only valid in the low velocity regime,Ohtsuki et al. (2002) introduced special interpolation coef-ficients C i . These coefficients tend to unity for very smallvelocity dispersions, and drop to zero in the high veloc-ity regime. Thus the interpolation formulæ are properly“switched off” in the high velocity regime, so they do notinterfere with the known high velocity stirring rates. All formulæ include only the stirring rates due to close en-counters, but non-crossing orbits also contribute to the over-all change of the velocity distribution. As these distant en-counters lead to small changes of the orbital elements, theproblem is accessible to perturbation theory; see Hasegawa(1990) for a detailed treatment. Stewart (2000) integratedthe perturbation solution over all impact parameters to de-rive the collective effect of all distant encounters: d (cid:104) e (cid:105) dt = Ω m ∗ Σ ∗ r ( m + m ∗ ) (cid:104) P VS , dist (cid:105)(cid:104) P VS , dist (cid:105) = 7 . α ( m + m ∗ ) M c EXINT (cid:16) α h ( (cid:104) e (cid:105) + (cid:104) e ∗ (cid:105) ) (cid:17) − EXINT (cid:16) α h ( (cid:104) i (cid:105) + (cid:104) i ∗ (cid:105) ) (cid:17) (cid:104) e (cid:105) + (cid:104) e ∗ (cid:105) − (cid:104) i (cid:105) − (cid:104) i ∗ (cid:105) EXINT( x ) := exp( x )Γ(0 , x ) h = (cid:114) m + m ∗ M c α ≈ α accounts for the uncertainty in the smallest impact param-eter that is regarded as a distant encounter. While distantencounters are already included in the interpolation formula c (cid:13) , 1–42 P. Glaschke, P. Amaro-Seoane & R. Spurzem of the low–velocity regime, we use the modified expression: (cid:18) dT r dt (cid:19) dist = 12 (Ω r ) (cid:18) d (cid:104) e (cid:105) dt (cid:19) dist (1 − C )= GM c r Ω m ∗ Σ ∗ m + m ∗ ) (cid:104) P VS , dist (cid:105) (1 − C ) (240)Stewart (2000) omitted the change in the inclination, as itis small due to the encounter geometry. Nevertheless we de-rived the integrated stirring rate for completeness, which wegive in 241 A close inspection of the integrated perturbationshows that the above formula is roughly a factor (cid:104) i (cid:105) + (cid:104) i ∗ (cid:105) smaller than the corresponding changes in the eccentricity. The presence of a gaseous disc damps the velocity dispersionof the planetesimals and introduces a slow inward migration.Adachi et al. (1976) used the drag law Eq. 75 to approxi-mate the average change of the orbital elements: τ = 2 mπC D ρ g R v K η g = | v K − v g | v K ddt e ≈ − e τ (cid:18) . e + 0 . i + 32 η g (cid:19) (242) ddt i ≈ − i τ (cid:18) . e + 0 . i + 12 η g (cid:19) ddt a ≈ − aτ η g (0 . e + 0 . i + η g ) (243) η g is the dimensionless velocity lag of the sub–Keplerianrotating gaseous disc. All expressions for the different velocity regimes are con-structed such that a smooth transition between the differentregimes is assured. Thus, a simple addition of all contribu-tions yields already the unified expressions dT r dt = (cid:18) dT r dt (cid:19) high + (cid:18) dT r dt (cid:19) low + (cid:18) dT r dt (cid:19) gas + (cid:18) dT r dt (cid:19) dist (244) dT z dt = (cid:18) dT z dt (cid:19) high + (cid:18) dT z dt (cid:19) low + (cid:18) dT z dt (cid:19) gas + (cid:18) dT z dt (cid:19) dist (245)which cover the full range of relative velocities. Althoughonly two populations m and m ∗ were assumed, Eq. 244 andEq. 245 are readily generalised to a multi-mass system byadding a summation over all masses. The preceding derivations assumed a homogeneous disc,which simplified the calculation, since the integration overall impact parameters needed no special precaution. A more A formal expansion at e = 0, i = 0, η g = 0 is not possible,since the drag law involves the modulus of the relative velocity.Kary et al. (1993) corrected a missing factor 3 / sophisticated consequence is that the spatial density and thedensity in semimajor axis space are equal:Σ( r ) = Σ( a ) = Σ (246)Density inhomogeneities break this simple relation, as par-ticles at the same radial distance could have different semi-major axes, and particles with the same semimajor axis arelocated at different positions at a given time. While bothrepresentations are equivalent (i. e. describe the same sys-tem in different ways), we chose the density in semimajoraxis space as the primary density . The spatial density isderived as:Σ( r ) = (cid:90) (cid:112) πa (cid:104) e ( a ) (cid:105) exp (cid:18) − ( a − r ) a (cid:104) e ( a ) (cid:105) (cid:19) Σ( a ) da (247)Likewise, T r and T z are also functions of the semimajor axis.Furthermore, an inhomogeneous surface density invali-dates the averaging over all impact parameters. Planetesimalencounter are most efficient for impact parameters smallerthan a few Hill radii, so the derivation is still valid if the sur-face density is roughly constant on that length scale. How-ever, a planetesimal that is large enough will “feel” the spa-tial inhomogeneities or even generates density fluctuations.Hence it is essential to extend the validity of the averagedexpressions to inhomogeneous systems. We use the averagedexpressions (cid:28) dT r,z dt (cid:29) = Σ( a ) (cid:90) ∞−∞ d ˜ T r,z ( b ) dt db (248)as a starting point ( d ˜ T r,z /dt excludes the surface density, asopposed to the averaged expressions). The (yet unknown)scattering contribution d ˜ T r,z /dt as a function of the impactparameter b is our starting point for a general expression fora varying surface density: dT ( a ) r,z dt = (cid:90) ∞−∞ Σ( a + b ) d ˜ T r,z ( b ) dt db (249)We restate Eq. 249 in terms of a weight function w ( b ): dT ( a ) r,z dt = (cid:28) dT r,z dt (cid:29) a ) (cid:90) ∞−∞ Σ( a + b ) w ( b ) db (250) w ( b ) = Σ( a ) d ˜ T r,z ( b ) dt (cid:28) dT r,z dt (cid:29) − (251)The numerical solution of the Hill problem (see Eq. 16)gives some insight into how the weight function w ( b ) changeswith the impact parameter. Fig. 12 illustrates the change in e of the relative motion during an encounter of two plan-etesimals that were initially on circular orbits. While thedetails depend on the initial inclination and eccentricity aswell as on the selected orbital element, all result share somebasic features. Small (compared to the Hill radius) impactparameters allow for a horseshoe orbit and the change inthe orbital elements is small except a change in the semima-jor axis. Intermediate impact parameters which lead to closeencounters provide the strongest perturbation, but they arealso more susceptible to complicated dynamics (compare the We denote Σ( a ) also as “surface density” and refer to a as aradial coordinate. However, all formulæ are precise in discrimi-nating both representations in r and a .c (cid:13)000
Chebyshev coefficients of the auxiliary functions J k and H k . et al. 2002): Λ = 112 ( (cid:104) ˜ e (cid:105) + (cid:104) ˜ i (cid:105) ) (cid:104) ˜ i (cid:105) / (235)˜ e = √ T r Ω R Hill ˜ i = √ T z Ω R Hill (236)Ohtsuki et al. 2002 also report a further improvement bysetting B ≡ A . Encounters in the low velocity regime exhibit a wealth of dif-ferent orbits, as the solar gravity field perturbs the two-bodyscattering. Only a small subset of the trajectories representssimple, regular orbits like Tadpole or Horseshoe orbits .Hence an examination of this velocity regime is done bestwith a numerical study of the parameter space by integrat-ing the equations of motions numerically (see Eq. 16).Ohtsuki et al. (2002) integrated a large set of planetesi-mal encounters and extracted fitting formulæ that cover thelow velocity regime. Their expressions for viscous stirringare: dT r dt = Gr Ω h Σ ∗ m + m ∗ ) 73 C m ∗ dT z dt = Gr Ω h Σ ∗ m + m ∗ ) C m ∗ (cid:18) (cid:104) ˜ i (cid:105) + 0 . (cid:104) ˜ e (cid:105) ) / (cid:113) (cid:104) ˜ i (cid:105) (cid:19) ˜ e = e/h ˜ i = i/hC := ln(10Λ / ˜ e + 1)10Λ / ˜ e C := ln(10Λ √ ˜ e + 1)10Λ √ ˜ e (237)The stirring rate of the radial velocity dispersion approachesa finite value for very low velocity dispersions, while the stir-ring rate for the vertical velocity dispersion drops to zero asthe velocity dispersion decreases. This different behaviour ofthe two limits is due to the encounter geometry: If two plan-etesimals have zero inclination, they may still excite highereccentricities during an encounter, but they remain confined The most famous example of such a regular orbit are the twosaturnian moons Janus and Epimetheus which share nearly thesame orbit. to the initial orbital plane preventing any excitation of in-clinations.The respective expressions for the dynamical frictionrates are: dT r dt = Gr Ω h Σ ∗ m + m ∗ )( T r + T ∗ r ) 10 C (cid:104) ˜ e (cid:105) ( T ∗ r m ∗ − T r m ) dT z dt = Gr Ω h Σ ∗ m + m ∗ )( T z + T ∗ z ) 10 C (cid:104) ˜ i (cid:105) ( T ∗ z m ∗ − T z m ) C := ln(10Λ + 1)10Λ (238)As the stirring rates are only valid in the low velocity regime,Ohtsuki et al. (2002) introduced special interpolation coef-ficients C i . These coefficients tend to unity for very smallvelocity dispersions, and drop to zero in the high veloc-ity regime. Thus the interpolation formulæ are properly“switched off” in the high velocity regime, so they do notinterfere with the known high velocity stirring rates. All formulæ include only the stirring rates due to close en-counters, but non-crossing orbits also contribute to the over-all change of the velocity distribution. As these distant en-counters lead to small changes of the orbital elements, theproblem is accessible to perturbation theory; see Hasegawa(1990) for a detailed treatment. Stewart (2000) integratedthe perturbation solution over all impact parameters to de-rive the collective effect of all distant encounters: d (cid:104) e (cid:105) dt = Ω m ∗ Σ ∗ r ( m + m ∗ ) (cid:104) P VS , dist (cid:105)(cid:104) P VS , dist (cid:105) = 7 . α ( m + m ∗ ) M c EXINT (cid:16) α h ( (cid:104) e (cid:105) + (cid:104) e ∗ (cid:105) ) (cid:17) − EXINT (cid:16) α h ( (cid:104) i (cid:105) + (cid:104) i ∗ (cid:105) ) (cid:17) (cid:104) e (cid:105) + (cid:104) e ∗ (cid:105) − (cid:104) i (cid:105) − (cid:104) i ∗ (cid:105) EXINT( x ) := exp( x )Γ(0 , x ) h = (cid:114) m + m ∗ M c α ≈ α accounts for the uncertainty in the smallest impact param-eter that is regarded as a distant encounter. While distantencounters are already included in the interpolation formula c (cid:13) , 1–42 P. Glaschke, P. Amaro-Seoane & R. Spurzem of the low–velocity regime, we use the modified expression: (cid:18) dT r dt (cid:19) dist = 12 (Ω r ) (cid:18) d (cid:104) e (cid:105) dt (cid:19) dist (1 − C )= GM c r Ω m ∗ Σ ∗ m + m ∗ ) (cid:104) P VS , dist (cid:105) (1 − C ) (240)Stewart (2000) omitted the change in the inclination, as itis small due to the encounter geometry. Nevertheless we de-rived the integrated stirring rate for completeness, which wegive in 241 A close inspection of the integrated perturbationshows that the above formula is roughly a factor (cid:104) i (cid:105) + (cid:104) i ∗ (cid:105) smaller than the corresponding changes in the eccentricity. The presence of a gaseous disc damps the velocity dispersionof the planetesimals and introduces a slow inward migration.Adachi et al. (1976) used the drag law Eq. 75 to approxi-mate the average change of the orbital elements: τ = 2 mπC D ρ g R v K η g = | v K − v g | v K ddt e ≈ − e τ (cid:18) . e + 0 . i + 32 η g (cid:19) (242) ddt i ≈ − i τ (cid:18) . e + 0 . i + 12 η g (cid:19) ddt a ≈ − aτ η g (0 . e + 0 . i + η g ) (243) η g is the dimensionless velocity lag of the sub–Keplerianrotating gaseous disc. All expressions for the different velocity regimes are con-structed such that a smooth transition between the differentregimes is assured. Thus, a simple addition of all contribu-tions yields already the unified expressions dT r dt = (cid:18) dT r dt (cid:19) high + (cid:18) dT r dt (cid:19) low + (cid:18) dT r dt (cid:19) gas + (cid:18) dT r dt (cid:19) dist (244) dT z dt = (cid:18) dT z dt (cid:19) high + (cid:18) dT z dt (cid:19) low + (cid:18) dT z dt (cid:19) gas + (cid:18) dT z dt (cid:19) dist (245)which cover the full range of relative velocities. Althoughonly two populations m and m ∗ were assumed, Eq. 244 andEq. 245 are readily generalised to a multi-mass system byadding a summation over all masses. The preceding derivations assumed a homogeneous disc,which simplified the calculation, since the integration overall impact parameters needed no special precaution. A more A formal expansion at e = 0, i = 0, η g = 0 is not possible,since the drag law involves the modulus of the relative velocity.Kary et al. (1993) corrected a missing factor 3 / sophisticated consequence is that the spatial density and thedensity in semimajor axis space are equal:Σ( r ) = Σ( a ) = Σ (246)Density inhomogeneities break this simple relation, as par-ticles at the same radial distance could have different semi-major axes, and particles with the same semimajor axis arelocated at different positions at a given time. While bothrepresentations are equivalent (i. e. describe the same sys-tem in different ways), we chose the density in semimajoraxis space as the primary density . The spatial density isderived as:Σ( r ) = (cid:90) (cid:112) πa (cid:104) e ( a ) (cid:105) exp (cid:18) − ( a − r ) a (cid:104) e ( a ) (cid:105) (cid:19) Σ( a ) da (247)Likewise, T r and T z are also functions of the semimajor axis.Furthermore, an inhomogeneous surface density invali-dates the averaging over all impact parameters. Planetesimalencounter are most efficient for impact parameters smallerthan a few Hill radii, so the derivation is still valid if the sur-face density is roughly constant on that length scale. How-ever, a planetesimal that is large enough will “feel” the spa-tial inhomogeneities or even generates density fluctuations.Hence it is essential to extend the validity of the averagedexpressions to inhomogeneous systems. We use the averagedexpressions (cid:28) dT r,z dt (cid:29) = Σ( a ) (cid:90) ∞−∞ d ˜ T r,z ( b ) dt db (248)as a starting point ( d ˜ T r,z /dt excludes the surface density, asopposed to the averaged expressions). The (yet unknown)scattering contribution d ˜ T r,z /dt as a function of the impactparameter b is our starting point for a general expression fora varying surface density: dT ( a ) r,z dt = (cid:90) ∞−∞ Σ( a + b ) d ˜ T r,z ( b ) dt db (249)We restate Eq. 249 in terms of a weight function w ( b ): dT ( a ) r,z dt = (cid:28) dT r,z dt (cid:29) a ) (cid:90) ∞−∞ Σ( a + b ) w ( b ) db (250) w ( b ) = Σ( a ) d ˜ T r,z ( b ) dt (cid:28) dT r,z dt (cid:29) − (251)The numerical solution of the Hill problem (see Eq. 16)gives some insight into how the weight function w ( b ) changeswith the impact parameter. Fig. 12 illustrates the change in e of the relative motion during an encounter of two plan-etesimals that were initially on circular orbits. While thedetails depend on the initial inclination and eccentricity aswell as on the selected orbital element, all result share somebasic features. Small (compared to the Hill radius) impactparameters allow for a horseshoe orbit and the change inthe orbital elements is small except a change in the semima-jor axis. Intermediate impact parameters which lead to closeencounters provide the strongest perturbation, but they arealso more susceptible to complicated dynamics (compare the We denote Σ( a ) also as “surface density” and refer to a as aradial coordinate. However, all formulæ are precise in discrimi-nating both representations in r and a .c (cid:13)000 , 1–42 ybrid methods: A new composite algorithm d (cid:104) i (cid:105) dt = Ω m ∗ Σ ∗ r ( m + m ∗ ) (cid:104) Q VS , dist (cid:105)(cid:104) Q VS , dist (cid:105) = 0 . α ( m + m ∗ ) M c × (cid:104) e (cid:105) + (cid:104) e ∗ (cid:105) − (cid:104) i (cid:105) − (cid:104) i ∗ (cid:105) × (cid:34) − αh (cid:104) i (cid:105) + (cid:104) i ∗ (cid:105) EXINT (cid:18) α h ( (cid:104) i (cid:105) + (cid:104) i ∗ (cid:105) ) (cid:19) − ( (cid:104) i (cid:105) + (cid:104) i ∗ (cid:105) ) EXINT (cid:16) α h ( (cid:104) e (cid:105) + (cid:104) e ∗ (cid:105) ) (cid:17) − EXINT (cid:16) α h ( (cid:104) i (cid:105) + (cid:104) i ∗ (cid:105) ) (cid:17) (cid:104) e (cid:105) + (cid:104) e ∗ (cid:105) − (cid:104) i (cid:105) − (cid:104) i ∗ (cid:105) (cid:35) (241) ∆ ( e / h ) b/H Figure 12.
Change of the relative eccentricity e due to an en-counter of two bodies initially on circular orbits. b/H is the impactparameter in units of the Hill radius. The plot was obtained byintegrating Eq. 16. resonant structures in Fig. 12). As the gravitational attrac-tion drops with increasing distance, non-crossing orbits yieldever smaller perturbations with increasing impact parame-ter. Aside from this qualitative behaviour, it is very difficultto derive precise expressions. While the limit of high ve-locities is accessible through the two-body approximation,any general formula involves some empiric interpolation tocover the full parameter space (Rafikov 2003b,a,see the ap-proximations of). Therefore we decided to approximate theweight function such that the main features of the trueweight function w ( b ) are reproduced. While this approachis less accurate, it provides better insight into the involvedapproximations. We expand the surface density under theintegral in Eq. 250 and compare the expansion coefficientsfor w ( b ) and the approximation ˜ w ( b ) to derive constraintson the choice of ˜ w ( b ). The lowest non-vanishing order is: l = (cid:90) ∞−∞ b w ( b ) db = (cid:90) ∞−∞ b ˜ w ( b ) db (252) l can be interpreted as the width of the heating zone. Con-dition 252 inspires our choice of the weight function ˜ w ( b )˜ w ( b ) = 1 √ πl exp (cid:18) − b l (cid:19) l = 1Ω ( T ( i ) r + T ( j ) r ) + R (253)where T ( i ) r and T ( j ) r are the radial velocity dispersions ofthe interacting radial bins and l is adjusted to the findingsof Ida (1993). The advantage of the bell curve is that it hasa discrete counterpart˜ w ( b ) ≈ a N (cid:32) Nb/ ∆ a + N/ (cid:33) N = 4 l / (∆ a ) (254) which makes the weight function readily applicable to thesummation on an equidistant radial grid with spacing ∆ a . We concentrated on the evolution of the velocity dispersionso far, but scatterings among planetesimals also change thesemimajor axis of the disc particles, inducing a diffusive evo-lution of the surface density: ∂ Σ ∂t = ∆ a ( D Σ) (255)The diffusion coefficient D is related to the typical changein semimajor axis ∆ a and the timescale T on whichplanetesimal encounters operate: D ≈ (∆ a ) T (256)If we neglect the radial displacement during an encounter,the change in semimajor axis is solely due to the change ofthe velocity: − GM a = − GMr + 12 v ∆ a ≈ a GM v · ∆ v (257)An average over all orientations of the velocity v and thevelocity change ∆ v yields the mean square change in semi-major axis: (cid:104) (∆ a ) (cid:105) ≈ a GM (cid:104) (∆ v ) (cid:105) = 43Ω (∆ T r + ∆ T φ + ∆ T z ) (258)This yields the mean diffusion coefficient D ≈ (cid:18) ddt T r + ddt T z (cid:19) (259)where the time derivatives of the velocity dispersions T r and T z are taken with respect to encounters. We already stated the coagulation equation for a multi-masssystem: 0 = ∂∂t mn ( t, m ) + ∂∂m F m ( t, m ) (260) F m = − (cid:90) (cid:90) n ( t, m ) n ( t, m ) σ ( m , m ) v rel M red ( m, m , m ) dm dm (261)Since the vertical density profile of a planetesimal disc isspecified by the known distribution function, we insert the c (cid:13) , 1–42 P. Glaschke, P. Amaro-Seoane & R. Spurzem isothermal density profile (see Eq. 216) in the coagulationequation 260 and integrate over z :0 = ∂∂t Σ( t, m ) + ∂∂m F m ( t, m ) (262) F m = − (cid:90) (cid:90) (cid:112) π ( h ( m ) + h ( m ) ) Σ ( t, m ) m Σ ( t, m ) m ˜ σ ( m , m ) v rel M red ( m, m , m ) dm dm (263)Σ( m ) is a short-hand notation for the differential surfacedensity d Σ dm . Further integration over all masses gives thetotal mass balance:Σ tot = (cid:90) m max m min Σ( t, m ) dm∂∂t Σ tot = F m ( m min ) − F m ( m max ) (264)The calculation of collisional cross sections is closely relatedto the underlying encounter dynamics. A homogenous sys-tem introduces no systematic perturbation, hence an en-counter is a pure two-body problem which is analyticallysolvable. Thus it is possible to derive the cross section with-out any approximation. Since encounters in the field of acentral star deviate noticeably from the pure Kepler solu-tion, the cross sections are also modified. While the crosssection in the high velocity regime reduces to the two-bodyformula (except minor corrections), the low velocity regimeis explored best by numerical calculations. It is not appropri-ate to disentangle the different contributions in Eq. 263, butto combine the various terms to the collisional probability P coll = ˜ σ ( m , m ) v rel (cid:112) π ( h ( m ) + h ( m ) ) (265)which can be easily deduced from the fraction of collidingorbits in Monte–Carlo simulations. An accurate expressionfor the collisional probability should include the two-bodycross section in the limit of high velocities and the numericaldata for the low velocity regime as well. We use numericalcalculations from Greenberg (1991); Greenzweig (1992) as a basis for a unified fitting formula˜ σ = σ × .
572 (1 + 3 . v Hill /v rel ) (cid:18) . σ Ω T z (cid:19) − / (266) σ = σ geom (cid:18) v ∞ v + 1 . v (cid:19) (267) v = 12 ( T r + T φ + T z ) v Hill = Ω r Hill (268)which gives an effective cross section ˜ σ for planetesimal–planetesimal encounters. Eq. 266 reduces to the well-knowngravitational focusing formula in the limit of high velocities:˜ σ ∝ σ geom (cid:18) v ∞ v (cid:19) (269) v ∞ = 2 G ( m + m ) R + R (270)If the vertical velocity dispersion is small, the disc becomestwo-dimensional and the cross section is proportional to R . Their work includes an averaging over the Rayleigh distributedinclinations and eccentricities of the colliding planetesimals.
The main differences to the two-body cross section 269 is afinite gravitational focusing factor, since the Keplerian shearinhibits a zero relative velocity, and a finite collisional prob-ability for very small velocities, again due to the shear whichprovides a finite influx of particles.The precise calculation of the coagulation kernel shouldinclude an integration over all semimajor axes with a properweighting kernel. As collisions among particles in the statis-tical model play only a major role when the system is stillhomogenous, we omitted this contribution. In addition, thishelps saving computational time, since the solution of thecoagulation equation is very costly. However, interactionsbetween N –body particles and the statistical model includespatial inhomogeneities properly (see section 17). Collisions are a dissipative process that removes kinetic en-ergy from the planetesimal system and damps the eccentric-ity. Low speed encounters leave the colliding bodies intactand damp the relative velocities through inelastic collisions.In contrast, high velocity encounters disrupt the collidingbodies and turn them into an expanding cloud of fragments.As a major part of the initial kinetic energy is convertedinto heat, the fragments disperse with rather low velocitiesthus reducing the overall velocity dispersion. We formulatethe dissipation due to collisions analogue to Eq. 262:0 = ∂∂t T k Σ( t, m ) + ∂∂m F Q,k ( t, m ) F Q,k = − (cid:90) (cid:90) (cid:112) π ( h ( m ) + h ( m ) ) Σ ( t, m ) m Σ ( t, m ) m ˜ σ ( m , m ) v rel Q red ,k ( m, m , m ) dm dm k ∈ { r, z } (271) Q red ,k is the kinetic energy redistribution function and F Q,k is the associated flux across the mass distribution. k indi-cates the two velocity dispersions. Q red ,k is a complex func-tion, since the disruption of a planetesimal produces frag-ments with a large scatter in velocities and a complicatedvelocity field. The velocity of a fragment consists of two con-tributions: The ejection velocity relative to the target andthe velocity ¯ v of the target within the corotating coordinatesystem. Owing to the strong dissipation, fragment velocitiesare dominated by the motion of the centre-of-mass of thetwo colliding bodies. Thus we neglect the ejection velocitiesand estimate the centre-of-mass motion. The initial kineticenergy of two colliding bodies is: E ini = 12 m v + 12 m v (272)We separate E ini into centre-of-mass motion and relativemotion and average over v , v : E ini = 12 M (cid:16) m v + m v M (cid:17) + 12 µ ( v + v ) (cid:104) E ini (cid:105) = 12 M (cid:0) m T + m T (cid:1) + 12 µ ( T + T ) (273)Most of the relative kinetic energy is dissipated during thecollision, so we neglect the relative motion after the incident.The final energy is (cid:104) E final (cid:105)(cid:104) E ini (cid:105) = m T + m T M ( m T + m T ) (274) c (cid:13) , 1–42 ybrid methods: A new composite algorithm which gives the drift motion ¯ v of the expanding fragmentcloud: ¯ v = 2 M (cid:104) E final (cid:105) (275) Q red is therefore coupled to the fragment redistribution func-tion M red Q red = ¯ v M red + Q diss (276)where the additional function Q diss removes the dissipatedenergy. The statistical model of a planetesimal disc does not only re-quire a large particle number to assure a proper descriptionof the system by a distribution function, but also the un-correlated motion of the planetesimals. Each of the formulæderived before involves the averaging over different impactparameters to some extend, in combination with the vitalassumption that all distances are equally probable. As longas all particles are subjected to perturbations by surround-ing bodies, strong correlations are suppressed. This appliesto the early stages, but the formation of protoplanets intro-duces a few dominant bodies that are not susceptible to theperturbation of the field planetesimals. Orbit repulsion givesrise to a regular spacing of the protoplanets, which preventsmutual collisions. Therefore not all impact parameters areequally probable due to this strong correlation. Hence a sta-tistical model is inherently not applicable to the late stagesof protoplanetary growth.Since statistical models are superior to N –body calcula-tions with respect to speed and (effective) particle number,modifications have been proposed to remedy this problem.The statistical model by Wetherill (1993) uses the fol-lowing solution: A gravitational range ∆ a (or buffer zone)is attached to each planetesimal, which represents the min-imal spacing that allows for stable orbits. They propose theexpression ∆ a = f ∆ R Hill + (cid:112) T r / Ω (277)where f ∆ is the minimal spacing in terms of the Hill radius.The value of f ∆ is adopted from Birn (1973), who derivedthe minimal spacing that allows for stable orbits: f ∆ = 2 √ m sep by the as-sumption that all bodies larger than this critical mass main-tain a clear buffer zone: f = (cid:90) ∞ m sep d Σ dm πa ∆ am dmf = 1 (279) f is the area covered by the buffer zones (overlapping is nottaken into account, therefore f > m sep can not en-force a minimum distance to their neighbours, as the wholedisc surface is already covered by the buffer zones of thelargest bodies. Owing to the regular spacing introduced bythe buffer zones, planetesimals larger than the critical massare not allowed to collide with each other. This approachhas also been employed by Inaba et al. (2001), who adopted f C f f ∆ =10f ∆ =2 √ Figure 13.
Covered fraction f C as a function of f for simulationS1FB at T = 10 years. The protoplanets are already formed andgrow oligarchically. f ∆ = 10, which is the mean distance of protoplanets accord-ing to the orbit repulsion mechanism.To shed light on the proper gravitational range and thevalidity of this approach, we defined an additional quantity f C which is the true area (i. e. overlapping is handled prop-erly) covered by the buffer zones in terms of the total area: f C = (cid:90) ∞ m sep d Σ dm πa (∆ a ) C m dmf C (cid:54) f C ( m ) as a functionof the integrated buffer zones f ( m ) for one of our hybrid cal-culations. Both values for f ∆ are included as well as the twolimiting cases random placing and perfect ordering. Thoughwe tested also other values of f ∆ , a spacing of ten Hill radiiproved to be the best choice.Our own experience with this method indicates that itworks reasonably well and agrees with Inaba et al. (2001)who used the same technique. However, this modification in-cludes the regular spacing of the protoplanets in a prescribedway, so any exploration of later stages, like the initiationof orbit crossing, is not accessible through this approach.Therefore we use it for comparison purposes only, since thehybrid code (see section 17) provides a much more generalframework. All involved quantities are only functions of a and m . There-fore we introduce a two dimensional grid, where Σ, T r and T z are cell centred quantities. Fig. 14 summarises the defini-tion of the two-dimensional grid. Since the full planetesimalsize range covers several orders of magnitude in mass, wechose a logarithmically equidistant discretisation in mass tocover the necessary mass range in a reliable way. The ra-dial spacing of the grid cells is equidistant. Thus the gridsetup for the mass discretisation reads ( N grid cells from c (cid:13) , 1–42 P. Glaschke, P. Amaro-Seoane & R. Spurzem mm i m i+1 aa i a i+1 F i F i+1 Σ i T r T z C o lli s i on s Encounters
Figure 14.
Numerical grid. The arrows indicate transport of ki-netic energy (red), spatial transport of mass(green) and accretion(black). Non-neighbouring cells are coupled by the coagulationkernel and the radial interpolation kernel. m min . . . m max ): m i = m min δ − i (1 / δ/ i = 1 , . . . , N ∆ m i = m min δ − i (1 − δ )Σ i = d Σ dm ∆ m i δ = (cid:18) m min m max (cid:19) /N (281)The grid spacing δ controls the number of cells which arenecessary to cover a specified mass range. As the evaluationof the coagulation equation scales with the third power ofthe number of grid cells, δ should be as large as possible. Ifthe flux integral (see Eq. 260) is approximated in a standardway F i = − N (cid:88) j =1 N (cid:88) k =1 F ( jk ) i F ( jk ) i = 1 (cid:113) π ( h j + h k ) Σ j m j Σ k m k × σ ( m j , m k ) v rel M red ( m i − ∆ m i / , m j , m k ) (282)a spacing δ much smaller than 2 is required to guarantee asufficient accuracy . However, it is possible to use a spacingof 2 if special precaution is taken. Spaute et al. (1991) ap-proximated the surface density with a power law, thus tak-ing the gradient with respect to mass into account. Whilethey reached only a sufficient accuracy with further specialadaptations, we use a more rigorous approach. A large spac-ing δ reduces the accuracy, since the partial flux (Eq. 282) isstrongly varying even inside one grid cell. Thus we rearrange Ohtsuki et al. (1990) give a thorough analysis of the impor-tance of the resolution. this expression to identify the most important terms, as wecan see in Eq. 283.The strongest varying contribution F V is now approxi-mated by a power law with respect to m j : F V ( m ) ≈ F V ( m j ) (cid:18) mm j (cid:19) q (284)Thus Eq. 284 is used to provide improved partial fluxes F ( jk ) i and thus we come to Eq. 285.Since the fragment redistribution function is a piecewisepower law, an analytical solution of the integral is possible.Eq. 285 gives reliable results even with a spacing δ = 2. Thetime derivative of the surface density reads˙Σ im = − F i +1 + F i (286)which assures the conservation of mass within numerical ac-curacy. All contributions to the evolution of the surface density Σand the velocity dispersions T r and T z are summarised bythe following set of differential equations: D = 43Ω (cid:18)
54 ˙ T r, enc + ˙ T z, enc (cid:19) d Σ dt = ∆ a ( D Σ) + ˙Σ coll d Σ T r dt = ∆ a ( D Σ T r ) + Σ ˙ T r + ddt (Σ T r ) coll d Σ T z dt = ∆ a ( D Σ T z ) + Σ ˙ T z + ddt (Σ T z ) coll (287)The Laplace operator is approximated in accordance withthe equidistant radial grid:∆ a f = 1 a ∂∂a (cid:18) a ∂∂a f (cid:19) ∆ a f ≈ f i +1 (1 + ∆ a/ (2 a i )) − f i + (1 − ∆ a/ (2 a i )) f i − (∆ a ) (288)We chose the Heun method as the basic integrator for thestatistical model. It is a second order accurate predictor–corrector scheme ( X is a vector containing all the abovequantities): dXdt = f ( X ) X p = X n + ∆ tf ( X n ) X n +1 = X p + 12 ∆ t ( f ( X p ) − f ( X n )) + O (∆ t ) (289)The Heun method is readily extended to an iterate scheme,which is equivalent to the implicit expression: X n +1 = X n + 12 ∆ t ( f ( X n +1 ) + f ( X n )) (290)This adds stability to the method and allows the secure in-tegration of stiff configurations that may appear during the The name of this method is not unique. Some texts denote itas the modified Euler method. The Heun method is a special caseof the Runge–Kutta methods. c (cid:13) , 1–42 ybrid methods: A new composite algorithm F ( jk ) i = v rel (cid:113) π ( h j + h k ) Σ k m k Σ j m j × σ ( m j , m k ) m j (cid:124) (cid:123)(cid:122) (cid:125) F V ( m j ) m j M red ( m i − ∆ m i / , m j , m k ) (283) j (cid:62) kF ( jk ) i = v rel (cid:113) π ( h j + h k ) Σ k m k (cid:90) m j +∆ m j / m j − ∆ m j / F V ( m ) m ∆ m j M red ( m i − ∆ m i / , m, m k ) dm (285) j (cid:62) k Pseudo-Force
Disc Particlesx i m i v i Σ (m,a)T r (m,a)T z (m,a) m trans xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx m min Accretion New ParticlesAccretion/Erosion CollisionsHeat Disc Scatter Disc ? Fragments
Figure 15.
Interplay between the N –body component and thestatistical component of the hybrid code. Black arrows indicatemass transfer, red arrows exchange of kinetic energy and greenarrows indicate spatial structuring, respectively. runaway accretion phase. In practice, three iterations aresufficient to guarantee a stable integration. As the diffu-sive part is discretised with a first order accurate formula(see Eq. 288), the whole iterated scheme is equivalent to theCrank–Nicolsen method. We choose a global time step forthe statistical model according to the expression∆ t = min (cid:18) η Disc X ˙ X , X ∈ { Σ , T r , T z } (cid:19) (291)where the hybrid code (see next section) applies an addi-tional discretisation in powers of two to achieve a bettersynchronisation with the N –body code component.
17 BRINGING THE TWO SCHEMESTOGETHER: THE HYBRID CODE
We introduced two different methods to solve the plan-etesimal growth problem. On the one hand, we modified
Nbody6++ , which has been used so far mainly for the sim-ulation of stellar clusters, to adapt it to the special require-ments of a long-term integration of planetesimal orbits. Onthe other hand, we developed a new statistical code with aconsistent evolution of the velocity dispersion, the capabil-ity to treat spatial inhomogeneities and a thoroughly con-structed collision treatment. Neither of the two approaches is powerful enough to provide a complete and accurate de-scription of the planetesimal problem, since each method isconfined to a certain range of the particle number. However,these restrictions are complementary in the sense that eachmethod covers a regime where the other method fails. Thisintriguing relation stimulated the construction of a hybridcode which combines the benefits of both methods.The basic idea is to introduce a transition mass m trans ,which separates the two mass regimes. Particles with a lowermass are treated by the statistical model, whereas larger par-ticles belong to the N –body model. Though both parts areclearly divided in different mass ranges, they are connectedby various interdependencies:(i) Direct collisions between particles lead to a mass ex-change. One process is the accretion of small particles by N –body particles, but agglomeration within the statisticalmodel can also produce particles larger than the transitionmass. This requires the generation of new N –body parti-cles. Energetic impacts may erode larger particles, so a cor-responding particle removal is also required for consistency.(ii) Mutual scatterings among N –body particles andsmaller planetesimals transfer kinetic energy. While energyequipartition leads to a systematic heating of the smallerfield planetesimals, a consistent treatment has to includeboth transfer directions.(iii) Accretion and scattering by the N –body particles in-duce spatial inhomogeneities or even gaps in the planetesi-mal component, if the particles have grown massive enough.Likewise, the small particles could induce some structure inthe distribution of the N –body particles. Since the spatialstructure is dominated by the stirring from few protoplanets,we neglect the latter process.Fig. 15 summarises this brief overview of the interactionsbetween the two code components in a schematic diagram.The following sections explain the implementation of eachinteraction term in more detail. An N –body particle accretes smaller particles in its vicinity.We already derived expressions which describe agglomera-tion within the statistical model, so it is manifest to applythese formulæ to derive the accretion rate of an N –bodyparticle.Most of the material is accreted within the cross-sectional area σ (see Eq. 267), but the finite eccentricityof an orbit extends the accessible radial feeding zone. Thus c (cid:13) , 1–42 P. Glaschke, P. Amaro-Seoane & R. Spurzem we assign the following surface density to each particleΣ( a ) = M πa √ πl exp (cid:18) − ( a − a ) l (cid:19) (292) l = σ/π + 12 a e + T r / Ω (293)by smearing it out over its feeding zone. T r is the radialvelocity dispersion of particles in the statistical model withsemimajor axis a . This density distribution is projected ontothe radial grid to calculate the accretion rate. As the timestep of the statistical model is much larger than the regu-lar step of an N –body particle, the particle mass update issynchronised with the statistical integration. The projectiontechnique allows the calculation of the accretion rates in asimple way, which gives the right size of the feeding zoneand the proper total accretion rate.Particle generation is included in the following way: A“virtual” mass bin is introduced as the boundary betweenthe statistical grid (denoted by the dashed area in Fig. 15)and the N –body component. Its sole task is to store massand kinetic energy that drives the statistical model towardshigher masses. If the mass content exceeds one m trans , a newparticle is created with inclination and eccentricity accord-ing to the stored velocity dispersions.The masses of the N –body particles are regularlychecked to detect any particle which dropped below the tran-sition mass. While this procedure would remove the particleand transfer the associated quantities back to the grid, wenever observed such a particle erosion. The projection of an N –body particle onto the grid withthe help of a proper weight function is also useful for thecalculation of the disc excitation due to stirring by the largerparticles. Since the Hill radius sets the proper length scalefor planetesimal encounters, the weight function is modifiedto Σ( a ) = M πa √ πl exp (cid:18) − ( a − a ) l (cid:19) l = R + 12 a e + T r / Ω (294)where T r is the radial velocity dispersion of the heated plan-etesimal component. The velocity dispersion of the stirring N –body particle is (in accordance with Eq. 11 and Eq. 12): T r, = 12 (Ω a ) e T z, = 12 (Ω a ) i (295)We employ the orbital elements as mediators between thefast varying instantaneous position and velocity of a parti-cle and the slow evolution of the statistical model, whichoperates on a longer relaxation timescale. In virtue of theprojection of the particle, we readily apply the standard in-teraction terms (see section 16) to evaluate the additionalheating due to the presence of N –body particles. While an N –body particle is moving through the disc, it alsointeracts gravitationally with the particles in the statisticalmodel. The collective effect of all these encounters leads to a change in the orbital elements of the N –body particle.Again, we project the N –body particle onto the grid andevaluate the stirring rates ˙ T r and ˙ T z , which correspond to achange in the orbital elements: ddt e = 2 ˙ T r (Ω a ) ddt i = 2 ˙ T z (Ω a ) (296)These time derivatives of eccentricity and inclination aretranslated to a pseudo-force that effects the desired changeof the orbital elements. We chose the ansatz F x,y = C r ( v x,y − ( v K ) x,y ) F z = C z v z (297)where v K is the local Keplerian velocity. In addition, wetried a simpler expression F x,y = 2 C r r x,y r · v r F z = C z v z (298)without any significant differences in the accuracy or thesimulation outcome. The proper friction coefficients are: C r = ˙ T r T r C z = ˙ T z T z (299)Since the relevant quantities are the time derivatives of theorbital elements, any other pseudo-force is also applicable.Though this approach yields the right mean change of theorbital elements, it lacks the statistical fluctuations fromthe particle disc. Hence the distribution of the orbital ele-ments of the N –body particles is artificially narrowed, whichis especially important when the N –body particles and thestatistical particles have a comparable mass. As the masscontrast between the two code parts is quite significant inplanet formation simulations, it is safe to neglect the fluctu-ating part without major restrictions on the realism of thesimulations.The friction coefficients C i are kept constant betweentwo integration steps of the statistical model. While a morefrequent update of the coefficients would be easily possible, aregular update on the basis of the statistical time step is ac-curate enough. Moreover, each update poses a considerablecomputational effort (roughly equivalent to 1000 force eval-uations), so our approach also saves valuable computationaltime. The first insight into planetesimal formation was obtainedby the particle–in–a–box method, which invokes the under-lying assumption that the planetesimal disc stays homoge-neous throughout the protoplanet growth (Greenberg et al.1978,see e.g.). While few large bodies introduce some coarse-graininess of the surface density, all smaller bodies are as-sumed to be evenly spread in the disc. Research on the in-teraction of protoplanets showed that this is an oversimpli-fication, as bodies that are massive enough could open gapsin their vicinity (Lin & Papaloizou 1979; Rafikov 2001,see c (cid:13) , 1–42 ybrid methods: A new composite algorithm No Σ ∆ a N N rad e /h i /h m TypeG1 1 . × − × − N –bodyPerturber — 1 – e = 6 . × − i = 3 . × − × − G2 1 . × − × − StatisticPerturber — 1 – e = 6 . × − i = 3 . × − × − Table 10.
Parameters of the statistical and the N –body gap simulation. The perturber is placed at the centre of the ring. Σ / Σ rGap at T=2000 yr N-BodyStatistic Figure 16.
Gap opening in a planetesimal disc. The gap is fullydeveloped after 2000 years. Table 10 summarises the initial con-ditions for the comparative runs. e.g.). Gap formation does not only change the overall surfacedensity, but also controls the accretion onto the protoplanetthrough the amount of planetesimals in the feeding zone. Ifgap formation is too effective, the growth of the protoplanetmay well stop before the isolation mass is reached. Henceany hybrid code should provide a framework that allowsthis mechanism to operate. A necessary condition is a radialdensity grid with a sufficient resolution to describe possiblyemerging gaps. A too low resolution suppresses local pertur-bations from the protoplanets by a simple averaging, thusinhibiting the formation of any spatial inhomogeneities. Asecond requirement is that the interaction terms relatingstatistical model and N –body model include the local inter-action between particles and the statistical component in aproper way.Our hybrid approach includes gap formation implicitlythrough the diffusive terms. A protoplanet heats only theplanetesimals in its vicinity (defined by the heating kernel),thus also increasing locally the diffusion coefficient. Hencethe surface density drops due to outward diffusion of theplanetesimals, given that the protoplanet is massive enough.The minimum gap opening mass is set by the condition thatthe protoplanet controls the random velocities of the fieldplanetesimals in its heating zone (see Eq. 30), which is equiv-alent to the independently derived gap formation criterion(compare Eq. 32).Although our algorithm invokes a simplified picture ofthe protoplanet–planetesimal interaction, it is surprisinglyaccurate with respect to the width of the forming gap andthe opening criterion. Fig. 16 shows a simulation which ex- < e / H > / , < i / H > / rT=2000 yr N-Body (e)Statistic (e)N-Body (i)Statistic (i) Figure 17.
Mean square eccentricity and inclination of thesmaller planetesimals in terms of the reduced Hill radius H ofthe protoplanet according to simulation G1 ( N –body) and G2(Statistic). amines the accuracy of our approach. The overall perfor-mance of the statistical code is quite remarkable, except asignificant overestimation of the surface density at the gapboundary compared to the N –body model. This deviationis due to the improper treatment of strong planetesimal–protoplanet encounters, which exceed the diffusive approxi-mation. Moreover, the higher concentration of planetesimalsnear the gap boundary leads to an additional overestimationof the velocity dispersion of the smaller planetesimals in thestatistical calculation (see Fig. 17). While the comparisonwith the N –body calculation clearly indicates a necessaryimprovement of the treatment of spatial inhomogeneities,our approach catches the main features of gap formation. Since the inventory of the new hybrid code is now completed,we turn to the specification of the transition mass m trans .The mass boundary between statistical and N –body parthas a major influence on the realism and the speed of thesimulation. On the one hand, optimisation with respect tospeed favours a large transition mass, whereas a reasonableresolution of the transition between the two components in-troduces some upper limit.Hence we identify first the set of large masses, whichcontrols the velocity dispersion of the disc, since these ob-jects are also possible candidates for gap opening. The in-spection of all involved stirring terms gives approximately c (cid:13) , 1–42 P. Glaschke, P. Amaro-Seoane & R. Spurzem the inequality: (cid:90) m trans d Σ dm mdm < (cid:90) ∞ m trans d Σ dm mdm (300)While this is a necessary condition to select all potential ma-jor perturbers, criterion 300 does not imply that all particlesin the selected mass range exert indeed a strong influenceon the disc. The number of possible gaps – and therefore thenumber of perturbers associated with them – is ultimatelylimited by the available space. Thus we integrate the areaof all potential gaps (width ≈ f ∆ R Hill ) and normalise it tothe total disc area: f C ≈ (cid:90) ∞ m trans f ∆ d Σ dm πaR Hill m dm (301)If the covered fraction f C is much larger than one, it ispossible to increase the transition mass until the condition f C (cid:47) m ).Therefore we chose a priori a fiducial value of the transitionmass, run the simulation and conduct an a posteriori check,whether the initial choice matches our requirements at anyevolutionary stage of the disc. A reliable value for a solarsystem analogue at 1 AU is m trans ≈ × − M (cid:12) (303)which restricts the number of N –body particles to atractable amount. Later stages would allow an even largertransition mass, but the current hybrid code does not in-clude any dynamical adjustment of the transition mass atruntime. Any numerical simulation is limited to a finite simulationvolume and a finite time interval. Therefore it is mandatoryto introduce proper boundary conditions which provide areasonable closure of the simulation volume.While boundary conditions with respect to time are thefamiliar initial conditions, the choice of the spatial bound-ary conditions for the various involved quantities depends onthe problem at hand and the type of the boundary. A sim-ulation boundary can be due to physical reasons (like wallsof a concert hall, surface of a terrestrial planet) or simplydue to a limitation in computational power that inhibits thecomplete numerical coverage of the problem.The current capability of the hybrid code sets limits onthe radial range as well as on the covered mass range, whicha simulation can handle in a reasonable time. Hence we haveto introduce artificial boundaries in radius, and a lower limitfor the mass grid.Any migration process couples the evolution of a lo-cal ring area in the planetesimal disc to the evolution ofthe whole disc. Inward (or outward) migrating material alsotransports information on the radial zone where the materialoriginated from. As this information is not available within the frame of a local simulation, any choice of the boundarycondition alters the evolution to some extend.However, we focus on a formation stage where migra-tion is not a dominant process, but provides only removalof the smaller collisional fragments. Thus we apply closedboundary conditions for the outer and inner radius of thering area (i. e. all fluxes vanish at the boundary), and anopen boundary for the lower end of the mass range. Whilethese conditions exclude the study of migrational processes,we gain clearer insight into the protoplanet growth.
18 DISCUSSION AND CONCLUSIONS
The formation of planetary systems represents a challengefrom a numerical standpoint. The dynamical problem spansover many orders of magnitudes in length and demands thecombination of different techniques. We have presented acomposite algorithm that brings together the advantagesof direct-summation tools and statistics for the descriptionof the planetesimal disc. Direct-summation N − body tech-niques have been around for some decades and have proventheir accuracy in a very large number of studies of stellarclusters such as galactic nuclei and globular and open clus-ters. We deem it to be the numerical tool to integrate themotion of the bodies for the very precise integration of theorbits and treatment of close encounters. Typically, in a sim-ulation of a stellar system, the energy is conserved in eachtimestep by E/ ∆ E ∼ − (where E is the total energyand ∆ E the difference between the former and current to-tal energy for a specific time), so that even if we integratefor a long time the cluster, the accumulated energy erroris negligible. Nevertheless, porting the numerical tool to theproblem of planetary dynamics is not straighforward and re-quires important modifications and additions. In this workwe present them in detail: the neighbour radius selection forthe protoplanets, the Hermite iteration and we introducefor the very first time the new extended Hermite scheme,since the usual Hermite scheme is not sufficient to integrateplanetesimal orbits accurately enough. Then we bring in newforces to the problem, namely the introduction of the centralpotential of the star, as well as the drag forces, which dependon the gas density and size of the planetesimals. Hence, theregularisation scheme, crucial to exactly integrate the closeencounters, has to be accordingly modified. We then intro-duce the disc geometry and discuss the required changes tothe neighbour scheme and prediction, as well as the commu-nication algorithm and block size distribution.For the statistical description of the planetesimal discwe employ a Fokker-Planck approach. We include dynamicalfriction, high- and low-speed encounters, the role of distantencounters as well as gas and collisional damping and thengeneralise the model to inhomogenous discs. We then de-scribe the combination of the two techniques to address thewole problem of planetesimal dynamics in a realistic way viaa transition mass to integrate the evolution of the particlesaccording to their masses.In particular, we introduce and describe the extendedHermite scheme, which reduces the the energy error by threeorders of magnitude with the same number of force evalua-tions, compared to the standard version of Nbody6++ .While the implementation and some code details are c (cid:13) , 1–42 ybrid methods: A new composite algorithm No. Σ ∆ a N N rad e /h i /h m TypeT1a 1 . × − . × − N –bodyT1b 1 . × − . × − HybridT1c 1 . × − . × − StatisticT2a 0 . × − × − N –body0 . × −
200 – 4 1 2 × − T2b 0 . × − × − Hybrid0 . × −
200 – 4 1 2 × − T3 Safronov – – – – – – StatisticT4a 1 . × − . × − N –bodyT4b 1 . × − . × − HybridT4c 1 . × − . × − StatisticT5 1 . × − – – – 620 155 2 . × − Statistic
Table 11.
Parameters of all test simulations. The transition mass in T4b is m trans = 3 . × − Only simulations T3, T4a–T4c andT5 include collisions. All values are scaled to M c = G = r = 1. newly introduced to the field of planet formation simula-tions, the first hybrid approach was developed in the early90’s. Spaute et al. (1991) (further improved in Weiden-schilling et al. 1997) constructed a hybrid code with a statis-tical component to treat the smaller particles and a specialtreatment for the larger particles. A statistical model cov-ers the field planetesimals with the help of a distributionfunction (similar to Wetherill (1989)), whereas the largerparticles are individually stored and characterised by mass,semimajor axis, eccentricity and inclination. While the inter-action between these single particles and the statistical com-ponent is expressed by standard viscous stirring and dynam-ical friction terms, perturbations among the single particlesare equated in a different way. First, the probability of anencounter of two neighbouring particles is calculated. Thisprobability is used in a second step to decide whether a (nu-merically integrated) two-body encounter of the neighbour-ing particles is carried out to derive the change in the orbitalelements. Though these two well-defined code componentsjustify to speak about a hybrid approach, the Monte–Carlolike integration of the largest particles is still closely relatedto a statistical treatment.A modified N –body approach is used in the work of Lev-ison & Morbidelli (2007). Their method covers the largestparticles by a direct N –body code, which includes thesmaller particles as “tracer” particles. The term “tracer”indicates that each particle represents a whole ensemble ofplanetesimals. In a similar line of approach and inspired bythis idea, Levison et al. (2010) modified a symplectic algo-rithm, Symba , to study the formation of giant planet cores.However, they made some assumptions in order to calcu-late the gravitational interaction between the planetesimals.In particular, they ignored totally close encounters betweenplanetesimals.Ormel & Spaans (2008) present in their work a schemebased on Monte Carlo techniques to cover the vast range of sizes. For this, they assign more resolution to those parti-cles that are more relevant to the interactions, typically thelargest bodies. Smaller particles are grouped and treatedcollectively, which means that they all share the same massand structural parameters. This classification is done in ac-cordance to the “zoom factor”, a free paramenter. Later,Ormel et al. (2010a) presented an detailed comparison oftheir Monte Carlo code with other techniques, in particularwith pure direct-summation N − body results and other sta-tistical studies and found that system leaves the runaway ata larger radius, in particular at the outer disc. With theirsimulations, the authors propose a new criterion for the run-away growth-oligarchy transition: from several hundreds ofkm in the inner disk regions up to a thousand km for theouter disc (Ormel et al. 2010b).Bromley & Kenyon (2006) published a description of ahybrid method with a basic approach similar to our work.They employ two velocity dispersions and the surface den-sity of the planetesimals to describe the planetesimal system.The statistical component includes migration of the plan-etesimals and dust particles due to gas drag and Pointing–Robertson drag. In contrast to our approach, they did notinclude mass transport due to the diffusion of the planetesi-mals, which precludes the study of spatial structures inducedby the protoplanets. One must note also that their methoduses the standard discretisation of the collisional flux (seeEq. 282) and thus restrict the spacing factor to δ (cid:46) . c (cid:13) , 1–42 P. Glaschke, P. Amaro-Seoane & R. Spurzem cretion on to massive cores, as well as accretion of small par-ticles in planetary atmospheres (Bromley & Kenyon 2010).While a variety of hybrid approaches emerged over thepast years, this technique is still far from a routinely appli-cation and is still challenged by many open issues. Hybridcodes bear the potential to address the dynamical evolutionof a whole planetary system, the later stages of protoplanetformation initiate a strong interaction with the gaseous disc,which may require more diligence than the inclusion of a fewadditional interaction terms. However, the development ispicking up speed, which places our work in a good positionfor further research.
APPENDIX A: CENTRAL FORCE –DERIVATIVES
Central force F per mass (i. e. acceleration) and its timederivatives are: F = − x Mx F (1) = − v Mx − A FF (2) = − a Mx − A ˙F − B FF (3) = − ˙a Mx − A F (2) − B F (1) − C Fa = ˙v A = x · v x B = v x + x · a x + A = ˙ A + 3 A C = 3 v · a x + x · ˙a x + A (3 B − A ) (A1)The F ( i ) denote the central force and its time deriva-tives, whereas a and ˙a refer to the total acceleration of theparticle. The assumption that x , v , a and ˙a are independentof each other allows the derivation of averaged expressionsfor particle–particle interactions: (cid:104) ( F ) (cid:105) = m x (cid:104) ( F (1) ) (cid:105) = m v x (cid:104) ( F (2) ) (cid:105) = m (cid:18) v x + 2 a x (cid:19) (cid:104) ( F (3) ) (cid:105) = m (cid:18) v x + 126 a v x + 2 ˙ a x (cid:19) (A2)We combine these expressions with Aarseth’s time step for-mula to derive the regular time step as a function of theneighbour sphere radius R s :∆ t reg ≈ √ η reg R s ¯ v
11 + (cid:112) R s /R R = 4 ¯ v a ≈ v ¯ r Gm = 4 ¯ r r close (A3)¯ r is the average particle distance and r close is the impactparameter for a 90–degree deflection. APPENDIX B: SCALABLE COLLISIONS FLUX
The mass flux according to the perturbation equation 161is: F p = − (cid:90) (cid:90) ( n ( m )∆ n ( m ) + n ( m )∆ n ( m )) σ ( m ) v ( m ) m f m ( m /m, (cid:15) ) dm dm = F (1) + F (2) (B1)Firstly, we employ the substitution m = mx m = m (cid:18) m m (cid:19) α w (cid:16) ˜ S (cid:17) w (cid:15) w (B2)to solve for the partial flux F (1) : F (1) = − n m σ v (cid:90) g ( mx ) F ( x ) dx F ( x ) = ˜ S − k (cid:48) (cid:90) (cid:15) − w + s +3+ α α +2 w f m ( x , (cid:15) ) x (1 + 2 w ) d(cid:15) (B3)The second contribution F (2) requires a slightly differenttransformation: m = mx (cid:15) − / (1+ α ) m = m (cid:18) mx m (cid:19) α w (cid:16) ˜ S (cid:17) w (B4)Thus the partial flux F (2) is: F (2) = − n m σ v (cid:90) g ( m ) F ( x ) dx F ( x ) = ˜ S − k (cid:48) (cid:90) (cid:15) − w + s +3+ α α +2 w f m ( x (cid:15) − / (1+ α ) , (cid:15) ) x (1 + 2 w ) d(cid:15) (B5)We change to a new set of logarithmic coordinates u = ln( m/m ) u = ln( x ) ˜ s = ln( ˜ S )1 + α (B6)which transforms the total flux F m to a convolution integral: F p = − n m σ v (cid:90) [ g ( u + u ) G ( u )+ g ( p ( u + u + s )) G ( u )] du (B7) p = 1 + α w (B8) p = 1 refers to the already derived solution for self-similarcollisions. Hence we expand Eq. B7 at p = 1 and retain onlythe zeroth-order moment of the fragmentation kernel: F p = − n m σ v (cid:104) g ( u ) G , + ( g ( u )+ u ( p − ∂g∂u ) G , (cid:105) (B9)This expression is equivalent to F p = − n m σ v ( g ( u ) G , +[ g ( u ) + ( p − g ( u ) − g (0))] G , ) (B10)where higher derivatives of g ( u ) are neglected. Hence werecover the same functional form of the perturbed mass flux F p as for self-similar collisions: F p = − n m σ v g ( u ) ( G , + p G , ) + const. ∝ ˜ S − k (cid:48) (B11) c (cid:13) , 1–42 ybrid methods: A new composite algorithm APPENDIX C: COAGULATION EQUATION
While the success of a general approximation of the coag-ulation equation depends heavily on the used coagulationkernel, we nevertheless provide a more general approach toembed section 10 in a broader context. The standard coag-ulation equation is:0 = ∂∂t mn ( t, m ) + ∂∂m F m ( t, m ) F m = − (cid:90) (cid:90) n ( t, m ) n ( t, m ) σ ( m , m ) v rel M red ( m, m , m ) dm dm (C1)In virtue of our experience drawn from the perturbationexpansion, we transform the coagulation equation to loga-rithmic coordinates u = ln( m ) (C2)and employ the size distribution g ( u ) relative to the steady-state solution n eq ( m ):0 = ∂∂t g ( u, t ) n eq ( u ) e u + ∂∂u F u ( t, m ) F u = − (cid:90) (cid:90) g ( t, u ) g ( t, u ) K ( u, u , u ) du du (C3) K ( u, u , u ) is the properly transformed new coagulationkernel. g ( u ) is expanded under the integral to arrive at amoment expansion of the flux F u : F u = − K ( u ) g ( u ) − ( K ( u ) + K ( u )) g ( u ) ∂g∂u + . . .K ij = (cid:90) (cid:90) K ( u, u , u ) u i u j du du (C4)Retaining only the leading order terms, we recover an ap-proximate coagulation equation which is similar to the in-viscid Burgers’ Equation :0 = ∂∂t g ( u, t ) n eq ( u ) e u − ∂∂u (cid:0) K ( u ) g ( u ) (cid:1) (C5) ACKNOWLEDGMENTS
It is a pleasure to thank Sverre Aarseth, Cornelis Dullemondand Phil Armitage for comments on the manuscript. PASthanks the National Astronomical Observatories of China,the Chinese Academy of Sciences and the Kavli Institutefor Astronomy and Astrophysics in Beijing, for an extendedvisit, as well as the Aspen Center of Physics and the orga-nizers of the summer meeting, where this work was finished.PAS expresses his utmost gratitude to Hong Qi, WenhuaJu and Xian Chen for their hospitality during his stay inBeijing. He is also somehow marginally indebted with the2011 winter strain of German H1N2, which allowed him toskip all possible duties and focus for an extended periodof time on this work at home. RS acknowledges supportby the Chinese Academy of Sciences Visiting Professorshipfor Senior International Scientists, Grant Number 2009S1-5(The Silk Road Project). The special supercomputer Laohu This notion goes back to Burgers (1948), but the equation wasalready introduced by Bateman (1915). at the High Performance Computing Center at National As-tronomical Observatories, funded by Ministry of Finance un-der the grant ZDYZ2008-2, has been used. Simulations werealso performed on the GRACE supercomputer (grants I/80041-043 and I/84 678-680 of the Volkswagen Foundation and823.219-439/30 and /36 of the Ministry of Science, Researchand the Arts of Baden-urttemberg). Computing time on theIBM Jump Supercomputer at FZ J¨ulich is acknowledged.
REFERENCES
Aarseth S., 1985, Academic Press Orlando, 378—, 1999, PASP, 111, 1333—, 2003, Cambridge University PressAdachi I., Hayashi C., Nakazawa K., 1976, Progress of The-oretical Physics, 56, 1756Ahmad A., Cohen L., 1973, Journal of ComputationalPhysics, 12, 389Armitage P. J., 2010, ArXiv e-printsAstakhov S., Lee E., Farrelly D., 2005, MNRAS, 360, 401Bagatin A., Cellino A., Davis D., Farinella P., Paolicchi P.,1994, Planet. Space. Sci., 42, 1079Balbus. S., 2003, Ann. Rev. Astron. Astrophys., 41, 555Balbus S., Hawley J., 1998, Reviews of Modern Physics,70, 1Bateman H., 1915, Monthly Weather Review, 43, 163Beauge C., Aarseth S. J., 1990, MNRAS, 245, 30Beckwith S., 1996, Nature, 383, 139Benz W., 1999, Icarus, 142, 5Binney J., 1977, MNRAS, 181, 735—, 1994, Princeton University Press, Princeton, New Jer-seyBirn J., 1973, Astronomy and Astrophysics, 24, 283Blum J., Wurm G., 2000, Icarus, 143, 138Bodenheimer P., 1986, Icarus, 67, 391Bodenheimer P., Hubickyj O., Lissauer J., 2000, Icarus,143, 2Bromley B. C., Kenyon S. J., 2006, AJ, 131, 2737—, 2010, ArXiv e-prints—, 2011, ArXiv e-printsBurgers J., 1948, Adv. Appl. Mech., 1, 171Cameron A. G. W., 1973, Space Science Reviews, 15, 121Cazenave A., Lago B., Dominh K., 1982, Icarus, 51, 133Chambers J., 2006, Icarus, 180, 496Chandrasekhar S., 1942, Astrophysical Journal, 97, 255Connaughton C., Rajesh R., Zaboronski O., 2004, PhysicalReview E, 69, 061114Dohnanyi J., 1969, Journal of Geophys. Research, 74, 2531Duncan M. J., Levison H. F., Lee M. H., 1998, 116, 2067Glaschke P., 2003, in ESA Special Publication, Vol. 539,Earths: DARWIN/TPF and the Search for ExtrasolarTerrestrial Planets, M. Fridlund, T. Henning, & H. La-coste, ed., pp. 425–428Goldreich P., Lithwick Y., Sari R., 2004, AstrophysicalJournal, 614, 497Goldreich P., Ward W. R., 1973, ApJ, 183, 1051Greenberg R., 1991, Icarus, 94, 98Greenberg R., Hartmann W., Chapman C., Wacker J.,1978, Icarus, 35, 1Greenzweig Y., 1992, Icarus, 100, 440Hasegawa M., 1990, Astronomy and Astrophysics, 227, 619 c (cid:13) , 1–42 P. Glaschke, P. Amaro-Seoane & R. Spurzem
Hayashi C., 1981, Progress of Theoretical Physics Supple-ment, 70, 35H´enon M., 1986, Celestial Mechanics, 38, 67Hernquist L., Hut P., Makino J., 1993, ApJ Lett., 402,L85+Hill G., 1878, American J. Math., 1, 5Housen K., 1990, Icarus, 84, 226Ida S., 1990, Icarus, 88, 129—, 1992, Icarus, 96, 107—, 1993, Icarus, 106, 210Ida S., Guillot T., Morbidelli A., 2008, ApJ, 686, 1292Ida S., Kokubo E., Makino J., 1993, MNRAS, 263, 875Inaba S., Barge P., Daniel E., Guillard H., 2005, Astronomyand Astrophysics, 431, 365Inaba S., Tanaka H., Nakazawa K., Wetherill G., KokuboE., 2001, Icarus, 149, 235Johansen A., Klahr H., Henning T., 2006, ApJ, 636, 1121Johansen A., Oishi J. S., Mac Low M.-M., Klahr H., Hen-ning T., Youdin A., 2007, Nat, 448, 1022Kary D., Lissauer J., Greenzweig Y., 1993, Icarus, 106, 288Kempf S., Pfalzner S., Henning T. K., 1999, Icarus, 141,388Kenyon S., 1998, Astronomical Journal, 115, 2136Kobayashi H., Tanaka H., Krivov A. V., Inaba S., 2010,Icarus, 209, 836Kokubo E., 1995, Icarus, 114, 247—, 1996, Icarus, 123, 180—, 1997, Icarus, 131, 171Kokubo E., Ida S., 2000, Icarus, 143, 15—, 2002, ApJ, 581, 666Kokubo E., Kominami J., Ida S., 2006, ApJ, 642, 1131Kokubo E., Yoshinaga K., Makino J., 1998, MNRAS, 297,1067Kornet K., Stepinski T., R´o˙zyczka M., 2001, Astronomyand Astrophysics, 378, 180Kustaanheimo P. E., Stiefel E. L., 1965, J. Reine Angew.Math.Larson R., 1970, MNRAS, 147, 323Levison H., Morbidelli A., 2007, Icarus, 189, 196Levison H. F., Thommes E., Duncan M. J., 2010, 139, 1297Lin D., Papaloizou J., 1979, MNRAS, 188, 191Makino J., 1988, Astrophysical Journal Supplement, 68,833—, 1992, PASJ, 44, 141Masset F. S., Papaloizou J. C. B., 2003, ApJ, 588, 494Mikkola S., 1997, Celestial Mechanics and Dynamical As-tronomy, 67, 145Mikkola S., Aarseth S. J., 1998, New Astronomy, 3, 309Moore A., Quillen A. C., 2010, ArXiv e-printsNakagawa Y., H., N., 1983, Icarus, 54, 361O’Brien D., 2003, Icarus, 164, 334O’dell C., Wen Z., Hu X., 1993, Astrophysical Journal, 410,696Ohtsuki K., Nakagawa Y., Nakazawa K., 1990, Icarus, 83,205Ohtsuki K., Stewart G., Ida S., 2002, Icarus, 155, 436Ormel C. W., Dullemond C. P., Spaans M., 2010a, ApJLett., 714, L103—, 2010b, Icarus, 210, 507Ormel C. W., Spaans M., 2008, ApJ, 684, 1291Paardekooper S.-J., Baruteau C., Kley W., 2011, MNRAS,410, 293 Papaloizou J., 2006, Rep. Prog. Phys., 69, 119Paszun D., Dominik C., 2009, A&A, 507, 1023Petit J.-M., 1986, Icarus, 66, 536Plummer H. C., 1911, MNRAS, 71, 460Rafikov R., 2001, Astronomical Journal, 122, 2713—, 2003a, Astronomical Journal, 125, 922—, 2003b, Astronomical Journal, 125, 906Reipurth B., Jewitt D., Keil K. E., 2007, Lunar and Plan-etary Information Bulletin, issue 110, p. 25 (May 2007),110, 25Safronov V., 1969, NASA-TTFSkeel R., 1999, Applied Numerical Mathematics, 29, 3Spaute D., Weidenschilling S., Davis D., Marzani F., 1991,Icarus, 92, 147Spurzem R., 1999, The Journal of Computational and Ap-plied Mathematics, Special Volume Computational Astro-physics, 109, 407Stewart G., 2000, Icarus, 143, 28Tanaka H., Inaba S., Nakazawa K., 1996, Icarus, 123, 450Wada K., Tanaka H., Suyama T., Kimura H., YamamotoT., 2009, ApJ, 702, 1490Ward W., 1986, Icarus, 67, 164Weidenschilling S., Spaute D., Davis D., Marzari F., Oht-suki K., 1997, Icarus, 128, 429Weidenschilling S. J., 1977, MNRAS, 180, 57Wetherill G., 1989, Icarus, 77, 330—, 1993, Icarus, 106, 190Wisdom J., 1991, Astronomical Journal, 102, 1528 c (cid:13)000