The Multifaceted Planetesimal Formation Process
Anders Johansen
Lund UniversityJürgen Blum
Technische Universität Braunschweig

Hidekazu Tanaka
Hokkaido University

Chris Ormel
University of California, Berkeley

Martin Bizzarro
Copenhagen University

Hans Rickman
Uppsala University
Polish Academy of Sciences Space Research Center, Warsaw
Accumulation of dust and ice particles into planetesimals is an important step in the planet
formation process. Planetesimals are the seeds of both terrestrial planets and the solid cores ofgas and ice giants forming by core accretion. Left-over planetesimals in the form of asteroids,trans-Neptunian objects and comets provide a unique record of the physical conditions inthe solar nebula. Debris from planetesimal collisions around other stars signposts that theplanetesimal formation process, and hence planet formation, is ubiquitous in the Galaxy. The
planetesimal formation stage extends from micrometer-sized dust and ice to bodies which canundergo run-away accretion. The latter ranges in size from 1 km to 1000 km, dependent on
the planetesimal eccentricity excited by turbulent gas density fluctuations. Particles face many
barriers during this growth, arising mainly from inefficient sticking, fragmentation and radial
drift. Two promising growth pathways are mass transfer, where small aggregates transfer up
to 50% of their mass in high-speed collisions with much larger targets, and fluffy growth,where aggregate cross sections and sticking probabilities are enhanced by a low internal
density. A wide range of particle sizes, from mm to 10 m, concentrate in the turbulent gas
flow. Overdense filaments fragment gravitationally into bound particle clumps, with mostmass entering planetesimals of contracted radii from 100 to 500 km, depending on local discproperties. We propose a hybrid model for planetesimal formation where particle growth starts
unaided by self-gravity but later proceeds inside gravitationally collapsing pebble clumps toform planetesimals with a wide range of sizes.

1.

INTRODUCTION

causing particles of sizes from centimeters to meters to drift
towards the star (Whipple, 1972; Weidenschilling, 1977a).
Drift speeds depend on the particle size and hence differentsized particles experience high-speed collisions.
The growth from dust and ice grains to planetesimals –the seeds of both terrestrial planets as well as the cores ofgas giants and ice giants – is an important step towards theassembly of a planetary system. A planetesimal can be defined as a body which is held together by self-gravity rather
than material strength, corresponding to minimum sizes oforder 100–1,000 meters (Benz, 2000). The planetesimalformation stage, on the other hand, must extend to sizes

Most stars are born surrounded by a thin protoplanetary
disc with a characteristic mass between 0.01% and 10% of
the mass of the central star (Andrews and Williams, 2005).
Planetesimal formation takes place as the embedded dust
and ice particles collide and grow to ever larger bodies.
Tiny particles collide gently due to Brownian motion, whilelarger aggregates achieve higher and higher collision speeds
as they gradually decouple from the smallest eddies of the
turbulent gas flow and thus no longer inherit the incompressibility of the gas flow (Voelk et al., 1980). The gas
disc is partially pressure-supported in the radial direction,1

where the escape speed of the largest bodies exceeds therandom motion of the planetesimals, to enter the next stageof gravity-driven collisions. The random speed of planetesimals and preplanetesimals (the latter can be roughly definedas bodies larger than 10 meters in size) is excited mainly by
the stochastic gravitational pull from density fluctuations in
the turbulent gas; under nominal turbulence conditions thelargest planetesimals need to reach sizes as large as 1,000
km to start run-away growth (Ida et al., 2008). In disc regions with less vigorous turbulence, planetesimal sizes between 1 and 10 kilometers may suffice to achieve significant
gravitational focusing (Gressel et al., 2012).
Planetesimals must grow to run-away sizes despite
bouncing and disruptive collisions (Blum and Wurm, 2000).
Planetesimals must also grow rapidly – radial drift timescales of cm-m-sized particles are as short as a few hundred
orbits. Despite these difficulties there is ample evidence
from cosmochemistry that large planetesimals formed in
the solar nebula within a few million years (Kleine et al.,
2004; Baker et al., 2005), early enough to melt and differentiate by decay of short-lived radionuclides.
The time since the last Protostars & Planets review onplanetesimal formation (Dominik et al., 2007) has seen alarge expansion in the complexity and realism of planetesimal formation studies. In this review we focus on threeareas in which major progress has been obtained, namely(i) the identification of a bouncing barrier at mm sizes andthe possibility of growth by mass transfer in high-speed collisions (Wurm et al., 2005; Johansen et al., 2008; Güttler
et al., 2010; Zsom et al., 2010; Windmark et al., 2012b), (ii)numerical simulations of collisions between highly porous
ice aggregates which can grow past the bouncing and radial drift barriers due to efficient sticking and cross sections that are greatly enhanced by a low internal density(Wada et al., 2008, 2009; Seizinger and Kley, 2013), and(iii) hydrodynamical and magnetohydrodynamical simulations which have identified a number of mechanisms forconcentrating particles in the turbulent gas flow of protoplanetary discs, followed by a gravitational fragmentationof the overdense filaments to form planetesimals with characteristic radii larger than 100 km (Rice et al., 2004; Johansen et al., 2007; Lyra et al., 2008a, 2009; Johansen
et al., 2009a; Bai and Stone, 2010a; Raettig et al., 2013).
The review is laid out as follows. In section 2 we givean overview of the observed properties of planetesimal belts
around the Sun and other stars. In section 3 we review thephysics of planetesimal formation and derive the planetesimal sizes necessary to progress towards planetary systems.
Section 4 concerns particle concentration and the density
environment in which planetesimals form, while the following section 5 describes laboratory experiments to determine
the outcome of collisions and methods for solving the coagulation equation. In section 6 we discuss the bouncing
barrier for silicate dust and the possibility of mass transfer
in high-speed collisions. Section 7 discusses computer simulations of highly porous ice aggregates whose low internal
density can lead to high growth rates. Finally in section 8

we integrate the knowledge gained over recent years to propose a unified model for planetesimal formation involving
both coagulation, particle concentration and self-gravity.
2.

PLANETESIMAL BELTS

The properties of observed planetesimal belts provide
important constraints on the planetesimal formation process. In this section we review the main properties of theasteroid belt and trans-Neptunians population as well as debris discs around other stars.
2.1.

Asteroids

The asteroid belt is a collection of mainly rocky bodies
residing between the orbits of Mars and Jupiter. Asteroid
orbits have high eccentricities (e∼0.1) and inclinations, excited by resonances with Jupiter and also a speculated population of embedded super-Ceres-sized planetary embryos
which were later dynamically depleted from the asteroid
belt (Wetherill, 1992). The high relative speeds imply thatthe asteroid belt is in a highly erosive regime where collisions lead to fragmentation rather than to growth.
Asteroids range in diameters from D ≈ 1000 km (Ceres)
down to the detection limit at sub-km sizes (Gladman et al.,2009). The asteroid size distribution can be parameterised
in terms of the cumulative size distribution N> (D) or thedifferential size distribution dN> /dD. The cumulative
size distribution of the largest asteroids resembles a broken power law N> ∝ D−p , with a steep power law indexp ≈ 3.5 for asteroids above 100 km diameter and a shallower power law index p ≈ 1.8 below this knee (Bottke
et al., 2005).
Asteroids are divided into a number of classes basedon their spectral reflectivity. The main classes are themoderate-albedo S-types which dominate the region from
2.1 to 2.5 AU and the low-albedo C-types which dominate
regions from 2.5 AU to 3.2 AU (Gradie and Tedesco, 1982;
Bus and Binzel, 2002). If these classes represent distinct
formation events, then their spatial separation can be used
to constrain the degree of radial mixing by torques fromturbulent density fluctuations in the solar nebula (Nelson
and Gressel, 2010). Another important class of asteroids
is the M-type of which some are believed to be the metallic
cores remaining from differentiated asteroids (Rivkin et al.,2000). Two of the largest asteroids in the asteroid belt,
Ceres and Vesta, are known to be differentiated (from theshape and measured gravitational moments, respectively,
Thomas et al., 2005; Russell et al., 2012).
Asteroids are remnants of solar system planetesimals
that have undergone substantial dynamical depletion andcollisional erosion. Dynamical evolution models can beused to link their current size distribution to the primordial
birth size distribution. Bottke et al. (2005) suggested, basedon the observed knee in the size distribution at sizes around
100 km, that asteroids with diameters above D ≃ 120 km
are primordial and that their steep size distribution reflects
their formation sizes, while smaller asteroids are fragments2

of collisions between their larger counterparts. Morbidelli
et al. (2009) tested a number of birth size distributions forthe asteroid belt, based on either classical coagulation models starting with km-sized planetesimals or turbulent concentration models predicting birth sizes in the range between 100 and 1000 km (Johansen et al., 2007; Cuzzi et al.,2008), and confirmed that the best match to the current size
distribution is that the asteroids formed big, in that the asteroids above 120 km in diameter are depleted dynamically but
have maintained the shape of the primordial size distribution. On the other hand, Weidenschilling (2011) found thatan initial population of 100-meter-sized planetesimals canreproduce the current observed size distribution of the asteroids, including the knee at 100 km. However, gap formation
around large planetesimals and trapping of small planetesimals in resonances is not included in this particle-in-a-box
approach (Levison et al., 2010).
2.2.

Fig. 1.— X-ray elemental maps composed of Mg (red), Ca(green), and Al (blue) of a fined-grained CAI from the Efremovka
CV3 chondrite (A) and two barred-olivine chondrules (B, C) fromthe NWA 3118 CV3 chondrite.

Meteorites

Direct information about the earliest stages of planet formation in the solar nebula can be obtained from meteorites
that record the interior structure of planetesimals in the asteroid belt. The reader is referred to the chapter by Davis etal. for more details on the connection between cosmochemistry and planet formation. Meteorites are broadly characterized as either primitive or differentiated (Krot et al.,2003). Primitive meteorites (chondrites) are fragments ofparent bodies that did not undergo melting and differentiation and, therefore, contain pristine samples of the earlysolar system. In contrast, differentiated meteorites are fragments of parent bodies that underwent planetesimal-scalemelting and differentiation to form a core, mantle and crust.
The oldest material to crystallise in the solar nebula
is represented in chondrites as mm- to cm-sized calciumaluminium-rich inclusions (CAIs) and ferromagnesian silicate spherules (chondrules) of typical sizes from 0.1 mm to
1 mm (Fig. 1). The chondrules originated in an unidentified thermal processing mechanism (or several mechanisms) which melted pre-existing nebula solids (e.g. Desch and Connolly, 2002; Ciesla et al., 2004); alternatively
chondrules can arise from the crystallisation of splash ejecta
from planetesimal collisions (Sanders and Taylor, 2005;
Krot et al., 2005).
The majority of CAIs formed as fine-grained condensates from an 16 O-rich gas of approximately solar composition in a region with high ambient temperature (>1300K) and low total pressures (∼10−4 bar). This environment
may have existed in the innermost part of the solar nebula during the early stage of its evolution characterized byhigh mass accretion rates (∼10−5 M⊙ yr−1 ) to the protoSun (D’Alessio et al., 2005). Formation of CAIs near the
proto-Sun is also indicated by the presence in these objectsof short-lived radionuclide 10 Be formed by solar energetic
particle irradiation (McKeegan et al., 2000). In contrast,
most chondrules represent coalesced 16 O-poor dust aggregates that were subsequently rapidly melted and cooled inlower-temperature regions (<1000 K) and higher ambient
vapor pressures (≥ 10−3 bar) than CAIs, resulting in igneous porphyritic textures we observe today (Scott, 2007).
Judging by their sheer abundance in chondrite meteorites,
the formation of chondrules must reflect a common process that operated in the early solar system. Using theassumption-free uranium-corrected Pb-Pb dating method,
Connelly et al. (2012) recently showed that CAIs define abrief formation interval corresponding to an age of 4567.30± 0.16 Myr, whereas chondrule ages range from 4567.32 ±0.42 to 4564.71 ± 0.30 Myr. These data indicate that chondrule formation started contemporaneously with CAIs andlasted ∼3 Myr.
A consequence of accretion of planetesimals on Myr
time-scales or less is the incorporation of heat-generating
short-lived radioisotopes such as 26 Al, resulting in widescale melting, differentiation and extensive volcanic activity. Both long-lived and short-lived radioisotope chronometers have been applied to study the timescale of asteroidal
differentiation. Of particular interest are the 26 Al-26 Mg and
182
Hf-182 W decay systems, with half-lives of 0.73 Myr and
9 Myr, respectively. Al and Mg are refractory and lithophileelements and thus remain together in the mantle after coreformation. In contrast, the different geochemical behavior of Hf and W during melting results in W being preferentially partitioned in the core and Hf being partitioned
into the silicate mantle and crust. Therefore, the short-lived
182
Hf-182 W system is useful to study the timescales of coreformation in asteroids, as well as in planets, which can beused to constrain planetesimal formation models.
Eucrite and angrite meteorites are the two most common groups of basaltic meteorites, believed to be derived
from the mantles of differentiated parent bodies. The HED(howardite-eucrite-diogenite) meteorite clan provides ourbest samples of any differentiated asteroid and could comefrom the 500-km-diameter asteroid 4 Vesta (Binzel and Xu,
3

after solar system formation. Therefore, chondrites may
represent samples of asteroidal bodies that formed after theaccretion of differentiated asteroids, at a time when the levels of 26 Al were low enough to prevent significant heating
and melting.
2.3.

Trans-Neptunian objects

The trans-Neptunian objects are a collection of rocky/icy
objects beyond the orbit of Neptune (Luu and Jewitt, 2002).
Trans-Neptunian objects are categorised into several dynamical classes, the most important for planetesimal formation being the classical Kuiper belt objects, the scattereddisc objects and the related centaurs which have orbits thatcross the orbits of one or more of the giant planets. The
classical Kuiper belt objects do not approach Neptune atany point in their orbits and could represent the pristine population of icy planetesimals in the outer solar nebula. The
so-called cold component of the classical Kuiper belt has alarge fraction (at least 30%) of binaries of similar size (Nollet al., 2008) which strongly limits the amount of collisional
grinding that these bodies can have undergone (Nesvorný
et al., 2011). The scattered disc objects on the other hand
can have large semi-major axes, but they all have perihelia
close to Neptune’s orbit. The centaurs move between the giant planets and are believed to be the source of short period
comets (see below).
The largest trans-Neptunian objects are much larger thanthe largest asteroids, with Pluto, Haumea, Makemake andEris defined as dwarf planets of 1.5-3 times the diameter of Ceres (Brown, 2008). Nevertheless the size distribution of trans-Neptunian objects shows similarities withthe asteroid belt, with a steep power law above a knee ataround D∼100 km (Fuentes and Holman, 2008). The turnover at the knee implies that there are fewer intermediatemass planetesimals than expected from an extrapolationfrom larger sizes – this has been dubbed the “Missing
Intermediate-Sized Planetesimals” problem (Sheppard andTrujillo, 2010; Shankman et al., 2013) and suggests that thecharacteristic planetesimal birth size was ∼100 km. The
accretion ages of Kuiper belt objects are not known, in
contrast to the differentiated asteroids where the inclusionof large amounts of 26 Al requires early accretion. Whilethe largest Kuiper belt objects are likely differentiated intoa rocky core and an icy mantle – this is clear e.g. for
Haumea which is the dense remnant of a differentiated body
(Ragozzine and Brown, 2009) – differentiation could be
due to long-lived radionuclides such as U, Th and 40 K andhence happen over much longer time-scales (∼100 Myr)
than the ∼Myr time-scale characteristic of asteroid differentiation by 26 Al decay (McKinnon et al., 2008).

Fig. 2.— Surface and interior structure of the 500-km-diameter
asteroid 4 Vesta, obtained by the NASA Dawn satellite. The ironcore is 110 km in radius, surrounded by a silicate mantle (green)
and a basaltic crust (gray). Image credit: NASA/JPL-Caltech1993), although see Schiller et al. (2011) for a different
view. The interior structure of Vesta was recently determined by the Dawn mission (Russell et al., 2012), indicating an iron core of approximately 110 km in radius (see Fig.
2). The 26 Al-26 Mg systematics of the angrite and HED meteorites indicate that silicate differentiation on these parent
bodies occurred within the first few Myr of solar system formation (Schiller et al., 2011; Bizzarro et al., 2005; Schiller
et al., 2010; Spivak-Birndorf et al., 2009). Similarly, theMg isotope composition of olivines within pallasite meteorites - a type of stony-iron meteorites composed of cmsized olivine crystals set in a iron-nickel matrix - suggests
silicate differentiation within 1.5 Myr of solar system formation (Baker et al., 2012). Chemical and isotopic diversityof iron meteorites show that these have sampled approximately 75 distinct parent bodies (Goldstein et al., 2009).
The 182 Hf-182 W systematics of some magmatic iron meteorites require the accretion and differentiation of their parent bodies to have occurred within less than 1 Myr of solar
system formation (Kleine et al., 2009).
An important implication of the revised absolute chronology of chondrule formation of Connelly et al. (2012), which
is not based on short-lived radionuclides for which homogeneity in the solar nebula has to be assumed, is that the accretion and differentiation of planetesimals occurred during
the epoch of chondrule formation. This suggests that, similarly to chondrite meteorites, the main original constituentsof early-accreted asteroids may have been chondrules. The
timing of accretion of chondritic parent bodies can be constrained by the ages of their youngest chondrules, which
requires that most chondrite parent bodies accreted >2 Myr

2.4.

CometsComets are typically km-sized volatile-rich bodies
which enter the inner solar system. They bring with them awealth of information about the conditions during the planetesimal formation epoch in the outer solar system. Out-

4

gassing provides knowledge about the compositions andheating histories of icy planetesimals, and the volumes andmasses of comet nuclei can be used to derive their density
which can be compared to the density expected from planetesimal formation models. Comet nuclei are (or are relatedto) icy planetesimals, i.e. planetesimals that formed beyond
the snow line in a formation zone extending from roughly
15 to 30 AU from the Sun – at least in the framework of
the so-called Nice Model where the giant planets form in
a compact configuration between 5 and 12 AU and later
migrate by planetesimal scattering to their current positions(Levison et al., 2011). This zone is wide enough to allow
for the appearance of chemical zoning, but such a zoning
would likely not be reflected in any separation between the
short- and long-period comet source regions.
Comets have to be considered together with the Centaurs and trans-Neptunians, because these are dynamicallyrelated to the Jupiter Family comets, and most of the data
we have on comets come from observations and modelingof Jupiter Family comets. The emergent picture is a very
wide size spectrum going from sub-km to 103 km in diameter. However, the slope of the size distribution is difficult to establish, since small objects are unobservable in theouter populations, and large ones are lacking in the Jupiter
Family. Generally, the measurements of masses and sizesof comet nuclei are consistent with a rather narrow range
of densities at about 0.5 g/cm3 (Weissman et al., 2004;
Davidsson et al., 2007). Most of the determinations use
the non-gravitational force due to asymmetric outgassing.
Assuming a composition that is roughly a 50-50 mix of ice
and refractories (Greenberg and Hage, 1990) the porosity
comes out as roughly 2/3. But the mass determinations only
concern the bulk mass, so one cannot distinguish between
meso-scale-porosity (rubble piles) and small-scale porosity(pebble piles). Estimates of comet tensile strengths havegenerally been extremely low, as expected for porous bodies. For comet Shoemaker-Levy 9, modeling of its tidal
breakup led to (non-zero) values so low that the object was
described as a strengthless rubble pile (Asphaug and Benz,
1996). The non-tidal splittings often observed for othercomets appear to be so gentle that, again, an essentially zerostrength has been inferred (Sekanina, 1997).
The current results on volatile composition of comets
have not indicated any difference between the Jupiter Family comets – thought to probe the scattered disc – and theHalley-type and long-period comets, which should comefrom the Oort Cloud (A’Hearn et al., 2012). From a dynamical point of view, this result is rather expected, sinceboth these source populations have likely emerged from thesame ultimate source, a disc of icy planetesimals extending beyond the giant planet zone in the early Solar System(Brasser and Morbidelli, 2013).
The origin of comets is closely related to the issue of thechemical pristineness of comets. If the total mass of theplanetesimal disc was dominated by the largest planetesimals, several hundreds of km in size, then the km-sized
comet nuclei that we are familiar with could arise fromcollisional grinding of the large planetesimals. Nevertheless there are severe problems with this idea. One is thatcomet nuclei contain extremely volatile species, like theS2 molecule (Bockelée-Morvan et al., 2004), which would
hardly survive the heating caused by a disruptive impact on
the parent body. The alternative that tidal disruptions oflarge planetesimals at mutual close encounters may lead tolarge numbers of smaller objects seems more viable, but thistoo suffers from the second pristineness problem, namely,
the geologic evolution expected within a large planetesimaldue to short-lived radio nuclei (Prialnik et al., 2004). If thelarge planetesimal thus becomes chemically differentiated,
there seems to be no way to form the observed comet nuclei
by breaking it up – no matter which mechanism we invoke.
The pristine nature of comets is thus consistent with theformation of small, km-scale cometesimals in the outer partof the solar nebula, avoiding due to their small size melting and differentiation due to release of short-lived radionuclides. Another possibility is that comets accreted frommaterial comprising an early formed 26 Al-free component(Olsen et al., 2013), which would have prevented differentiation of their parent bodies. The recent chronologyof 26 Al-poor inclusions found in primitive meteorites indicate that this class of objects formed coevally with CAIs
(which record the canonical 26 Al/27 Al value of 5 × 10−5 ,Holst et al., 2013). This provides evidence for the existence of 26 Al-free material during the earliest stages of the
protoplanetry disc. The discovery of refractory material
akin to CAIs in the Stardust samples collected from Comet
81P/Wild 2 show efficient outward transport of material
from the hot inner disc regions to cooler environments far
from the Sun during the epoch of CAI formation (Brownlee
et al., 2006). Efficient outward transport during the earliest stages of solar system evolution would have resulted in
the delivery of a significant fraction of the 26 Al-poor material to the accretion region of cometary bodies. Analysis of
the Coci refractory particle returned from comet 81P/Wild
2 did not show detectable 26 Al at the time of its crystallization (Matzel et al., 2010) which, although speculative,
is consistent with the presence an early-formed 26 Al-freecomponent in comets.
2.5.

Debris discs

Planetesimal belts around other stars show their presence
through the infrared emission of the dust produced in collisions (Wyatt, 2008). The spectral energy distribution of
the dust emission reveals the orbital distance of the planetesimal belt. Warm debris discs resembling the asteroid
belt in the Solar System are common around young starsof around 50 Myr age (at least 30%), but their occurrence
falls to a few percent within 100–1,000 Myr (Siegler et al.,2007). Thus planetesimal formation appears to be ubiquitous around the Sun as well as around other stars. This isin agreement with results from the Kepler mission that flat
planetary systems exist around a high fraction of solar-typestars (Lissauer et al., 2011; Tremaine and Dong, 2012; Jo-5ρo = 1 g cm−3101

ρo = 10−5 g cm−3

105

St = 10

St = 1

St = 10

St = 0.1

St = 1

10−1

104

103

a [m]

a [m]

100

St = 0.1

St = 0.01

10−2

102
To Stokes dragTo non−linear dragTo quadratic drag

10−3
0.1

St = 0.01

101
1.0

10.0

0.1

r [AU]

1.0

10.0
r [AU]

Fig. 3.— The particle size corresponding to Stokes numbers from 0.01 to 10, the left plot for compact particles with material densityρ• = 1 g cm−3 and the right for extremely porous aggregates with ρ• = 10−5 g cm−3 . Such fluffy particles may be the result ofcollisions between ice aggregates (discussed in section 7). Transitions between drag force regimes are indicated with large dots.

gas sound speed cs and the gas density ρg . Particles withsizes above 9/4 times the mean free path of the molecules
enter the Stokes regime (see Whipple, 1972, and references
therein), withRρ• 4 R.
(3)
τf =cs ρg 9 λhansen et al., 2012b; Fang and Margot, 2012).
3.

THE PLANETESIMAL FORMATION STAGE

The dust and ice particles embedded in the gas in protoplanetary discs collide and merge, first by contact forces
and later by gravity. This process leads eventually to theformation of the terrestrial planets and the cores of gas giants and ice giants forming by core accretion. The planetesimal formation stage can broadly be defined as the growthfrom dust grains to particle sizes where gravity contributes
significantly to the collision cross section of two colliding
bodies.
3.1.

Here the friction time is proportional to the squared radius
and independent of gas density, since λ is inversely proportional to the gas density. The flow Reynolds number
past the particle, Re = (2Rδv)/ν, determines the furthertransition to drag regimes that are non-linear in the relative
speed, with the kinematic viscosity given by ν = (1/2)cs λ.
At unity flow Reynolds number the drag transitions to anintermediate regime with the friction time proportional to(δv)−0.4 . Above Re = 800 the drag force finally becomes
quadratic in the relative velocity, with friction time

Drag force

Small particles are coupled to the gas via drag force. The
acceleration by the drag force can be written asv̇ = −

1(v − u) ,
τf

(1)

τf =

where v is the particle velocity, u is the gas velocity at theposition of the particle and τf is the friction time which contains all the physics of the interaction of the particle with thegas flow (Whipple, 1972; Weidenschilling, 1977a). The friction time can be divided into different regimes, dependingon the mean free path of the gas molecules, λ, and the speedof the particle relative to the gas, δv = |v − u|. The Epstein regime is valid when the particle size is smaller thanthe mean free path. The flux of impinging molecules is set
in this regime by their thermal motion and the friction timeis independent of the relative speed,
τf =

Rρ•
.
cs ρg

6Rρ•
.
(δv)ρg(4)

Following the descriptions above, the step-wise transition
from Epstein drag to fully quadratic drag happens in the optically thin minimum mass solar nebula (MMSN, Hayashi,1981), with power-law index −1.5 for the surface density(Weidenschilling, 1977b) and −0.5 for the temperature, atparticle sizes
9λ
4
ν
R2 =
2(δv)
800ν
R3 =
2(δv)
R1 =

(2)

=

3.2 cm

≈6.6 cm≈ r 2.75
AU r 2.5,AU r 2.5.
52.8 m
AU

,

(5)
(6)
(7)

Here R1 denotes the Epstein-to-Stokes transition, R2 theStokes-to-non-linear transition and R3 the non-linear-to-Here R is the radius of the particle, assumed to be spherical. The other parameters are the material density ρ• , the6The azimuthal drift peaks at vφ = vK −∆v for the smallest Stokes numbers where the particles are carried passively
with the sub-Keplerian gas. The radial drift peaks at unity
Stokes number where particles spiral in towards the star at
vr = −∆v. The radial drift of smaller particles is slowed
down by friction with the gas, while particles larger than
Stokes number unity react to the perturbing gas drag by entering mildly eccentric orbits with low radial drift.
Radial drift puts requirements on the particle growth at
Stokes number around unity (generally from St = 0.1 to
St = 10) to occur within time-scale tdrift ∼r/∆v∼100–
1000 orbits (Brauer et al., 2007, 2008a), depending on thelocation in the protoplanetary disc. However, three considerations soften the at first glance very negative impact of theradial drift. Firstly, the ultimate fate of drifting particles is
not to fall into the star, but rather to sublimate at evaporation
fronts (or snow lines). This can lead to pile up of material
around evaporation fronts (Cuzzi and Zahnle, 2004) and toparticle growth by condensation of water vapour onto existing ice particles (Stevenson and Lunine, 1988; Ros andJohansen, 2013). Secondly, the radial drift flow of particles is linearly unstable to streaming instabilities (Youdin
and Goodman, 2005), which can lead to particle concentration in dense filaments and planetesimal formation by gravitational fragmentation of the filaments (Johansen et al.,2009a; Bai and Stone, 2010a). Thirdly, very fluffy particles with low internal density reach unity Stokes number,where the radial drift is highest, in the Stokes drag regime(Okuzumi et al., 2012). In this regime the Stokes number,which determines radial drift, increases as the square of theparticle size and hence growth to “safe” Stokes numberswith low radial drift is much faster than for compact particles. These possibilities are discussed more in the following
sections.

quadratic transition. The latter two equations are only approximate because of the dependence of the friction time on
the relative speed. Here we used the sub-Keplerian speed∆v as the relative speed between particle and gas, an approximation which is only valid for large Stokes numbers(see equation 8 below).
A natural dimensionless parameter to construct from thefriction time is the Stokes number St = Ωτf , with Ω denoting the Keplerian frequency at the given orbital distance.
The inverse Keplerian frequency is the natural referencetime-scale for a range of range of physical effects in protoplanetary discs, hence the Stokes number determines (i)turbulent collision speeds, (ii) sedimentation, (iii) radial
and azimuthal particle drift, (iv) concentration in pressure
bumps and vortices and (v) concentration by streaming instabilities. Fig. 3 shows the particle size corresponding toa range of Stokes numbers, for both the nominal densityof 1 g/cm3 and extremely fluffy particles with an internal
density of 10−5 g/cm3 (which could be reached when ice
aggregates collide, see section 7).
3.2.

Radial drift

Protoplanetary discs are slightly pressure-supported in
the radial direction, due to gradients in both mid-plane density and temperature. This leads to sub-Keplerian motion of
the gas, vgas = vK − ∆v with vK denoting the Keplerian
speed and the sub-Keplerian velocity difference ∆v definedas (Nakagawa et al., 1986)∆v ≡ ηvK = −

1
2



H
r

2

∂ ln PvK .
∂ ln r(8)

In the MMSN the aspect ratio H/r rises proportional to
r1/4 and the logarithmic pressure gradient is ∂ ln P/∂ ln r =−3.25 in the mid-plane. This gives a sub-Keplerian speed
which is constant ∆v = 53 m/s, independent of orbital distance. Radiative transfer models with temperature and density dependent dust opacities yield disc aspect ratios H/rwith complicated dependency on r and thus a sub-Keplerianmotion which depends on r (e.g. Bell et al., 1997). Nevertheless, a sub-Keplerian speed of ∼50 m/s can be used as
the nominal value for a wide range of protoplanetary discmodels.
The drag force on the embedded particles leads to particle drift in the radial and azimuthal directions (Whipple,
1972; Weidenschilling, 1977a)
vr

=

vφ

=

2∆v
,
St + St−1∆v
vK −.
1 + St2−

3.3.

Collision speeds

Drag from the turbulent gas excites both large-scale random motion of particles as well as collisions (small-scalerandom motion). The two are distinct because particles may
have very different velocity vectors at large separations, butthese will become increasingly aligned as the particles approach each other due to the incompressible nature of thegas flow. Particles decouple from turbulent eddies with
turn-over times shorter than the friction time, and the gradual decoupling from the smallest eddies results in crossing
particle trajectories. This decoupling can be interpreted as
singularities in the particle dynamics (so-called “caustics”,
Gustavsson and Mehlig, 2011). Caustics in turn give rise
to collisions as small particles enter regions of intense clustering (clustering is discussed further in section 4.1.1). The
random motion is similar to the turbulent speed of the gas
for small particles, but particles experience decreased random motion as they grow to Stokes numbers above unity.
The problem when calculating the turbulent collision
speed is that at close separations (when the particles areabout to collide) they interact with the same eddies, which

(9)
(10)

These equations give the drift speed directly in the Epsteinand Stokes regimes. In the non-linear and quadratic dragregimes, where the Stokes number depends on the relative
speed, the equations can be solved using an iterative method
to find consistent vr , vφ and St.

7

200.0

30.0

104

3.0

α = 10−2, γ = 0

10.0

100.0

200.0

104

102

100.0

α = 10−2, γ = 6×10−4

102

100.0

a2 [m]

100.0

α = 10−4, γ = 0

102

10−6
10−6

104

104

10−4

10−2

200.0

100.0

3.0

100
a1 [m]

100
a1 [m]

102

α = 10−4, γ = 6×10−5

104

10.0

30.0

10−2

30.0
10.0
3.0

10−4

1.0

104

100.0

10−4

100.0

3.0

10−4

10−6
10−6

100.0

10−2

10.0
30.0

10−2

100

10.0
30.0

a2 [m]

30.0

100

0.1

102

102
10.0
30.0
10.0

30.0

a2 [m]

10−2

100

30.0

10−2
0.1
1.0
3.0
10.0
30.0

0.1

10−4

10−6
10−6

10.0

30.0
10.0

10−4

10−2

100
a1 [m]

1.0
3.0
10.0
30.0

a2 [m]

3.0

100

10−4

102

10−6
10−6

104

10−4

10−2

100
a1 [m]

102

104

Fig. 4.— The collision speed, in meters per second, of two particles of size a1 and a2 , with contributions from Brownian motion,differential radial and azimuthal drift, and gas turbulence. The upper panels show collision speeds for α = 10−2 and the lower panels
show collision speeds for α = 10−4 . The gravitational pull from turbulent gas density fluctuations is included in the right panels. The
red line marks the transition from dominant excitation by direct drag to dominant excitation by turbulent density fluctuations. The “oasis”
of low collision speeds for particles above 10 meters vanishes when including eccentricity pumping by turbulent density fluctuations.
causes their motions to become highly correlated. The
framework set out by Voelk et al. (1980) which employsa Langevin approach, is still widely used, and Ormel andCuzzi (2007) provided closed-form analytical approximations to their results (but see Pan and Padoan, 2010, fora criticism of the simplifications made in the Völk model).
The closed-form expressions of Ormel and Cuzzi (2007) require numerical solution of a single algebraic equation foreach colliding particle pair (defined by their friction times).
With knowledge of the properties of the turbulence, particularly the turbulent rms speed and the frequency of the smallest and the largest eddies, the collision speeds can then be
calculated at all locations in the disc.
Another important contribution to turbulent collision
speeds is the gravitational pull from turbulent gas density
fluctuations. The eccentricity of a preplanetesimal increases
as a random walk due to uncorrelated gravitational kicks
from the turbulent density field (Laughlin et al., 2004; Nelson and Papaloizou, 2004). The eccentricity would growunbounded with time as e ∝ t1/2 in absence of dissipation. Equating the eccentricity excitation time-scale withthe time-scale for damping by tidal interaction with the gas

disc (from Tanaka and Ward, 2004), aerodynamic gas drag,and inelastic collisions with other particles, Ida et al. (2008)
provide parameterisations for the equilibrium eccentricity
as a function of particle mass and protoplanetary disc properties. The resulting collision speeds dominate over thecontributions from the direct drag from the turbulent gas at
sizes above approximately 10 meters.
Ida et al. (2008) adopt the nomenclature of Ogihara
et al. (2007) for the eccentricity evolution, where a dimensionless parameter γ determines the proportionality between the eccentricity and t1/2 . The parameter γ is expected to scale with the density fluctuations δρ/ρ but canbe directly calibrated with turbulence simulations. The
shearing box simulations by Yang et al. (2009) of turbulence caused by the magnetorotational √instability (Balbus
and Hawley, 1991) suggest that δρ/ρ ∝ α, where α is thedimensionless measure of the turbulent viscosity (Shakura
and Sunyaev, 1973). In their nominal ideal-MHD turbulence model with α ≈ 0.01, Yang et al. (2012) find γ ≈6 × 10−4 . This leads to an approximate expression for √γ asa function of the strength of the turbulence, γ ≈ 0.006 α.
The resulting eccentricities are in broad agreement with
8

103

10−1

1 AU, α =1 × 10−2 , γ =10−3 (ideal)sc

5 AU, α =2 × 10−4 , γ =10−5 (dead)
102

10−2

Ideal MRI
10−3

101

100

mm

10−4

Dead zone

cm

10 cm

m
km
10 m
100 m
Dust particle or planetesimal size

Eccentricity (5 AU curves)

5 AU, α =1 × 10−2 , γ =10−3 (ideal)ve

Turbulent rms-velocity ∆v [m s−1 ]1 AU, α =2 × 10−4 , γ =10−5 (dead)

10 km

100 km

1 000 kmFig. 5.— The collision speed of equal-sized particles as a function of their size, from dust to planetesimals, based on the stirring modelof Ormel and Okuzumi (2013). The collision speeds of small particles (below approximately 10 m) are excited mainly by direct gas drag,while the gravitational pull from turbulent gas density fluctuations dominates for larger particles. The transition from relative-speeddominated to escape-speed-dominated can be used to define the end of the planetesimal formation stage (dashed line). The planetesimal
sizes that must be reached for run-away accretion range from near 1000 km in highly turbulent discs (α = 10−2 ) to around 1–10 km indisc regions with an extended dead zone that is stirred by sheared density waves from the active regions (α = 2 × 10−4 ).

the resistive magnetohydrodynamics simulations of Gressel
et al. (2012) where the turbulent viscosity from Reynolds
stresses fall below α = 10−4 in the mid-plane.
In Fig. 4 we show the collision speeds of two particlesof sizes from µm to 100 km, in a figure similar to Fig. 3of the classical Protostars and Planets III review by Weidenschilling and Cuzzi (1993). We take into account the
Brownian motion, the differential drift, and the gas turbulence (from the closed-form expressions of Ormel andCuzzi, 2007). We consider an MMSN model at r = 1 AU
and α = 10−2 in the upper panels and α = 10−4 in thelower panels. The gravitational pull from turbulent density
fluctuations driven by the magnetorotational instability is
included in the right panels. The collision speed approaches
100 m/s for meter-sized boulders when α = 10−2 and 30
m/s when α = 10−4 , due to the combined effect of thedrag from the turbulent gas and differential drift. The collision speed of large, equal-sized particles (with St>1) would
drop asrα(11)δv = cs
St

The collision speed of equal-sized particles is exploredfurther in Fig. 5. Here we show results both of a modelwith fully developed MRI turbulence and a more advanced
turbulent stirring model which includes the effect of a deadzone in the mid-plane, where the ionisation degree is too
low for the gas to couple with the magnetic field (Fleming and Stone, 2003; Oishi et al., 2007), with weak stirringby density waves travelling from the active surface layers(Okuzumi and Ormel, 2013; Ormel and Okuzumi, 2013).
The collision speeds rise with size until peaking at√unity
Stokes number with collision speed approximately αcs .
The subsequent decoupling from the gas leads to a decline
by a factor ten, followed by an increase due to the turbulent density fluctuations starting at sizes of approximately
10 – 100 meters. The planetesimal formation stage can be
defined as the growth until the gravitational cross section of
the largest bodies is significantly larger than their geometriccross section, a transition which happens when the escapespeed of the largest bodies approach their random speed. In
Fig. 5 we indicate also the escape speed as a function of
the size of the preplanetesimal. The transition to run-away
accretion can only start at 1000 km bodies in a disc withnominal turbulence α = 10−2 (Ida et al., 2008). Lower
values of α, e.g. in regions of the disc where the ionisation
degree is too low for the magnetorotational instability, leadto smaller values for the planetesimal size needed for runaway accretion, namely 1–10 km.

in absence of turbulent density fluctuations (for a discussionof this high-St regime, see Ormel and Cuzzi, 2007). This region has been considered a safe haven for preplanetesimals
after crossing the ridges around unity Stokes number (Weidenschilling and Cuzzi, 1993). However, the oases vanishes
when including the gravitational pull from turbulent gas
density fluctuations, and instead the collision speeds continue to rise towards larger bodies, as they are damped lessand less by gas drag (right panels of Fig. 4).

9

3.4.

Sedimentation

Mid-plane solids-to-gas ratios above unity are reached
for St = 1 when δ < 10−4 . Weaker stirring allows successively smaller particles to reach unity solids-to-gas ratio in
the mid-plane. This marks an important transition to whereparticles exert a significant drag on the gas and become concentrated by streaming instabilities (Youdin and Goodman,
2005; Johansen et al., 2009a; Bai and Stone, 2010a,c). Thiseffect will be discussed further in the next section.

While small micrometer-sized grains follow the gas density tightly, larger particles gradually decouple from the gas
flow and sediment towards the mid-plane. An equilibrium
mid-plane layer is formed when the turbulent diffusion of
the particles balances the sedimentation. The Stokes number introduced in section 3.1 also controls sedimentation(because the gravity towards the mid-plane is proportionalto Ω 2 ). Balance between sedimentation and turbulent diffusion (Dubrulle et al., 1995) yields a mid-plane layer thickness Hp , relative to the gas scale-height H, of
r
δ
Hp
=
.
(12)
H
St + δ

4.

High local particle densities can lead to the formation ofplanetesimals by gravitational instability in the sedimented
mid-plane layer. Particle densities above the Roche density
ρR =Here δ is a measure of the diffusion coefficient D = δcs H,
similar to the standard definition of α for turbulent viscosity(Shakura and Sunyaev, 1973). In general δ ≈ α in turbulence driven by the magnetorotational instability (Johansen
and Klahr, 2005; Turner et al., 2006), but this equality may
be invalid if accretion is driven for example by disc winds(Blandford and Payne, 1982; Bai and Stone, 2013).
Equation (12) was derived assuming particles to fall
towards the mid-plane at their terminal velocity vz =−τf Ω 2 z. Particles with Stokes number larger than unity
do not reach terminal velocity before arriving at the midplane and hence undergo oscillations which are damped bydrag and excited by turbulence. Nevertheless Carballido
et al. (2006) showed that equation (12) is in fact valid forall values of the Stokes number. Youdin and Lithwick (2007)
interpreted this as a cancellation between the increased sedimentation time of oscillating particles and their decreased
reaction to the turbulent motion of the gas.
The particle density in the equilibrium mid-plane layerisr
St + δ
H
= Zρg
.
(13)
ρp = Zρg
Hp
δ

9Ω 2
4πG

=
=

4.3 × 10−7 g cm−3
 r −1/4
315 ρgAU r −3AU(15)

are bound against the tidal force from the central star andcan contract towards solid densities. The scaling with gas
density assumes a distribution of gas according to the minimum mass solar nebula. Sedimentation increases the particle density in the mid-plane and can trigger bulk gravitational instabilities in the mid-plane layer (Goldreich andWard, 1973). Sedimentation is nevertheless counteracted
by turbulent diffusion (see section 3.4 and equation 12) andparticle densities are prevented from reaching the Roche
density by global turbulence or by mid-plane turbulence
induced by the friction of the particle on the gas (Weidenschilling, 1980; Johansen et al., 2009a; Bai and Stone,
2010a).
High enough densities for gravitational collapse can nevertheless be reached when particles concentrate in the gas
turbulence to reach the Roche density in local regions.
4.1.

Passive concentration

The turbulent motion of gas in protoplanetary discs canlead to trapping of solid particles in the flow. The size ofan optimally trapped particle depends on the length-scale
and turn-over time-scale of the turbulent eddies. We describe here particle trapping progressively from the smallest scales unaffected by disc rotation (eddies) to the largest
scales in near-perfect geostrophic balance between the pressure gradient and Coriolis accelerations (vortices and pressure bumps).
We consider rotating turbulent structures with lengthscale ℓ, rotation speed ve and turn-over time te = ℓ/ve .
These quantities are approximate and all factors of orderunity are ignored in the following. Dust is considered tobe passive particles accelerating towards the gas velocity
in the friction time τf , independent of the relative velocity.
Particle trapping by streaming instabilities, where the particles play an active role in the concentration, is described in
section 4.2. The three main mechanisms for concentrating
particles in the turbulent gas flow in protoplanetary discs
are sketched in Fig. 6.

Here Z is the ratio of the particle column density to the gas
column density. In the limit of large particle sizes the midplane layer thickness is well approximated by Hp = cp /Ω,
where cp is the random particle motion. This expression is
obtained by associating the collision speed of large particles
in equation (11) with their random speed (with α = δ) andinserting this into equation (12) in the limit St ≫ δ. The
mass flux density of particles can then be written as
F = cp ρp = Zρg HΩ ,PARTICLE CONCENTRATION(14)

independent of the collision speed as well as the degreeof sedimentation, since the increased mid-plane density oflarger particles is cancelled by the decrease in collision
speeds. Hence the transition from Stokes numbers above
unity to planetesimal sizes with significant gravitationalcross sections is characterised by high collision speeds andslow, ordered growth (Ṙ is independent of size when F
is constant) which does not benefit from increased particle
sizes nor from decreased turbulent diffusion.
10

Eddies

Pressure bumps / vortices

Streaming instabilitiesl ~ η ~ 1 km, St ~ 10−5−10−4

l ~ 1−10 H, St ~ 0.1−10

l ~ 0.1 H, St ~ 0.01−1

Fig. 6.— The three main ways to concentrate particles in protoplanetary discs. Left panel: turbulent eddies near the smallest scalesof the turbulence, η, expel tiny particles to high-pressure regions between the eddies. Middle panel: the zonal flow associated withlarge-scale pressure bumps and vortices, of sizes from one scale height up to the global scale of the disc, trap particles of Stokes number
from 0.1 to 10. Right panel: streaming instabilities on intermediate scales trap particles of Stokes number from 0.01 to 1 by accelerating
the pressure-supported gas to near the Keplerian speed, which slows down the radial drift of particles in the concentration region.

4.1.1.

Isotropic turbulence

1991; Fessler et al., 1994). In an astrophysics context, suchturbulent concentration of sub-mm-sized particles between
small-scale eddies has been put forward to explain the narrow size ranges of chondrules found in primitive meteorites(Cuzzi et al., 2001), as well as the formation of asteroids
by gravitational contraction of rare, extreme concentration
events of such particles (Cuzzi et al., 2008). This model wasnevertheless criticised by Pan et al. (2011) who found that
efficiently concentrated particles have a narrow size range
and that concentration of masses sufficiently large to form
the primordial population of asteroids is unlikely.

On the smallest scales of the gas flow, where the Coriolis
force is negligible over the turn-over time-scale of the eddies, the equation governing the structure of a rotating eddy
is
1 ∂P
dvr
(16)
=−
≡ fP .
dtρ ∂rHere fP is the gas acceleration caused by the radial pressure
gradient of the eddy. We use r as the radial coordinate in aframe centred on the eddy. The pressure must rise outwards,∂P/∂r > 0, to work as a centripetal force. In such lowpressure eddies the rotation speed is set by
fP = −ve2ℓ.

4.1.2.

On larger scales of protoplanetary discs, gas and particle motion is dominated by Coriolis forces and shear. We
first expand our particle-trapping framework to flows dominated by Coriolis forces and then generalise the expression
to include shear.
In a gas rotating rigidly at a frequency Ω, the equilibriumof the eddies is now given by(17)

Very small particles with τf ≪ te reach their terminal velocity
vp = −τf fP(18)
on a time-scale much shorter than the eddy turn-over timescale. This givesv2
τf
vp = −τf fP = τf e = ve .
ℓteTurbulence in rigid rotation2Ωve −

(19)

v2
1 ∂P
=− e .
ρ ∂r
ℓ

(20)

For slowly rotating eddies with ve /ℓ ≪ Ω we can ignore
the centripetal term and get

The largest particles to reach their terminal velocity in theeddy turn-over time-scale have τf ∼ te . This is the optimal particle size to be expelled from small-scale eddies
and cluster in regions of high pressure between the eddies.
Larger particles do not reach their terminal velocity before
the eddy structure breaks down and reforms with a newphase, and thus their concentration is weaker.
Numerical simulations and laboratory experiments haveshown that particles coupling at the turn-over time-scale ofeddies at the Kolmogorov scale of isotropic turbulence experience the strongest concentrations (Squires and Eaton,

ve = −

fP
.
2Ω

(21)

High pressure regions have ve < 0 (clockwise rotation),while low pressure regions have ve > 0 (counter-clockwise
rotation).
The terminal velocity of inertial particles can be found11

by solving the equation system
dvx
dt
dvy
dt

1
vx ,
τf
1
−2Ωvx − (vy − ve ) ,τf=+2Ωvy −=

Haghighipour and Boss, 2003). Rice et al. (2004) identified this mechanism as the cause of particle concentrationin spiral arms in simulations of self-gravitating protoplanetary discs. Shearing box simulations of turbulence causedby the magnetorotational instability (Balbus and Hawley,1991) show the emergence of long-lived pressure bumps
surrounded by a super-Keplerian/sub-Keplerian zonal flow
(Johansen et al., 2009b; Simon et al., 2012). Similarly
strong pressure bumps have been observed in global simulations (Fromang and Nelson, 2005; Lyra et al., 2008b).
High-pressure anticyclonic vortices concentrate particles inthe same way as pressure bumps (Barge and Sommeria,1995). Vortices may arise naturally through the baroclinic
instability thriving in the global entropy gradient (Klahr and
Bodenheimer, 2003), although the expression of the baroclinic instability for long cooling times and in the presenceof other sources of turbulence is not yet clear (Lesur and Papaloizou, 2010; Lyra and Klahr, 2011; Raettig et al., 2013).
Pressure bumps can also be excited by a sudden jumpin the turbulent viscosity or by the tidal force of an embedded planet or star. Lyra et al. (2008a) showed that the
inner and outer edges of the dead zone, where the ionisation degree is too low for coupling the gas and the magnetic field, extending broadly from 0.5 to 30 AU in nominal
disc models (Dzyurkevich et al., 2013, chapter by Turner
et al.), develop steep pressure gradients as the gas piles up
in the low-viscosity region. The inner edge is associatedwith a pressure maximum and can directly trap particles.
Across the outer edge of the dead zone the pressure transitions from high to low, hence there is no local maximum
which can trap particles. The hydrostatic equilibrium nevertheless breaks into large-scale particle-trapping vortices
through the Rossby wave instability (Li et al., 2001) in thevariable-α simulations of Lyra et al. (2008a). Dzyurkevich
et al. (2010) identified the inner pressure bump in resistive
simulations of turbulence driven by the magnetorotational
instability. The jump in particle density, and hence in ionisation degree, at the water snow line has been proposed tocause a jump in the surface density and act as a particle trap
(Kretke and Lin, 2007; Brauer et al., 2008b; Dra̧żkowska
et al., 2013). However, the snow line pressure bump remains to be validated in magnetohydrodynamical simulations.
Lyra et al. (2009) found that the edge of the gap carved
by a Jupiter-mass planet develops a pressure bump which
undergoes Rossby wave instability. The vortices concentrate particles very efficiently. Observational evidence for alarge-scale vortex structure overdense in mm-sized pebbles
was found by van der Marel et al. (2013). These particles
have approximately Stokes number unity at the radial distance of the vortex, in the transitional disc IRS 48. This
discovery marks the first confirmation that dust traps exist
in nature. The approximate locations of the dust-trapping
mechanisms which have been identified in the literature areshown in Fig. 7.

(22)
(23)

for vx ≡ vp . Here we have fixed a coordinate system in thecentre of an eddy with x pointing along the radial direction
and y along the rotation direction at x = ℓ. The terminal
velocity is
ve
vp =
.
(24)
(2Ωτf )−1 + 2Ωτf
Thus high-pressure regions, with ve < 0, trap particles,
while low-pressure regions, with ve > 0, expel particles.
The optimally trapped particle has 2Ωτf = 1. Since we are
now on scales with te ≫ Ω −1 , the optimally trapped particle has τf ≪ te and thus these particles have ample time toreach their terminal velocity before the eddies turn over. An
important feature of rotating turbulence is that the optimallytrapped particle has a friction time that is independent of theeddy turn-over time-scale, and thus all eddies trap particlesof τf ∼ 1/(2Ω) most efficiently.
4.1.3.

Turbulence with rotation and shear

Including both Keplerian shear and rotation the terminal
velocity changes only slightly compared to equation (24),
vp =

2ve
.
−1
(Ωτf ) + Ωτf(25)

Equilibrium structures in a rotating and shearing frame
are axisymmetric pressure bumps surrounded by superKeplerian/sub-Keplerian zonal flows (Fig. 6). The optimally trapped particle now has friction time Ωτf = 1. This
is in stark contrast to isotropic turbulence where each scaleof the turbulence traps a small range of particle sizes withfriction times similar to the turn-over time-scale of the eddy.
The eddy speed is related to the pressure gradient throughve = −

1 1 ∂P2Ω ρ ∂r(26)

Associating ve with −∆v defined in equation (8), we recover the radial drift speed
vr = −

2∆v.
St−1 + St(27)

However, this expression now has the meaning that particles
drift radially proportional to the local value of ∆v. Particles
pile up where ∆v vanishes, i.e. in a pressure bump.
4.1.4.

Origin of pressure bumps

Trapping of St∼1 particles (corresponding to cm-sizedpebbles to m-sized rocks and boulders, depending on theorbital distance, see Fig. 3) in pressure bumps in protoplanetary discs has been put forward as a possible way to crossthe meter-barrier of planetesimal formation (Whipple, 1972;
12

MRI = magnetorotational instability
DZE = dead zone edge
SI = streaming instability
BI = baroclinic instability
EF = evaporation front (snow line)
PGE = planet gap edge
GI = gravitational instability

MRI

GI
MRI
DZE
EF
SI

PGE

SI

BI

DZE

0.1

1

10

100
AU

Fig. 7.— Sketch of the particle concentration regions in a wedge of a protoplanetary disc seen from above. Regions where themagnetorotational instability is expected to operate are marked with red, while the extent of the dead zone in a nominal protoplanetary
disc model is marked with blue. The particle trapping mechanisms are described in the main text.

4.2.

Streaming instability

different wavenumbers and friction times. These modes
have subsequently been used in other papers (e.g. Bai andStone, 2010b) to test the robustness and convergence of numerical algorithms for the coupled dynamics of gas andsolid particles.
The non-linear evolution of the streaming instability canbe studied either with or without particle stratification. The
case without particle stratification is closest to the linear stability analysis. In this case the mean particle mass-loadingin the simulation domain must be specified, as well as thefriction time of the particles. The initial gas and particle
velocities are set according to drag force equilibrium (Nakagawa et al., 1986), with particle drifting in towards the starand gas drifting out.
In simulations not including the component of the stellar gravity towards the mid-plane, i.e. non-stratified simulations, high particle densities are reached mainly for particles
with St = 1 (Johansen and Youdin, 2007; Bai and Stone,
2010b). At a particle mass-loading of 3 (times the mean gas
density), the particle density reaches almost 1000 times thegas density. The overdense regions appear nearly axisymmetric in both local and semi-global simulations (Kowalik
et al., 2013). Smaller particles with St = 0.1 reach only
between 20 and 60 times the gas density when the particle
mass-loading is unity of higher. Little or no concentration
is found at a dust-to-gas ratio of 0.2 for St = 0.1 particles.
Bai and Stone (2010b) presented convergence tests of nonstratified simulations in 2-D and found convergence in theparticle density distribution functions at 10242 grid cells.

The above considerations of passive concentration ofparticles in pressure bumps ignore the back-reaction friction force exerted by the particles onto the gas. The radial
drift of particles leads to outwards motion of the gas in the
mid-plane, because of the azimuthal frictional pull of theparticles on the gas.
Youdin and Goodman (2005) showed that the equilibrium streaming motion of gas and particles is linearly unstable to small perturbations, a result also seen in the simplified mid-plane layer model of Goodman and Pindor (2000).
The eight dynamical equations (six for the gas and particle
velocity fields and two for the density fields) yield eightlinear modes, one of which is unstable and grows exponentially with time. The growth rate depends on both thefriction time and the particle mass-loading. Generally thegrowth rate increases proportional to the friction time (upto St∼1), as the particles enjoy increasing freedom to move
relative to the gas. The dependence on the particle massloading is more complicated; below a mass-loading of unity
the growth rate increases slowly (much more slowly thanlinearly) with mass-loading, but after reaching unity in dustto-gas ratio the growth rate jumps by one or more orders ofmagnitude. The e-folding time-scale of the unstable modeis as low as a few orbits in the regime of high mass-loading.
The linear mode of the streaming instability is an exact
solution to the coupled equations of motion of gas and particles, valid for very small amplitudes. This property canbe exploited to test numerical algorithm for solving the full
non-linear equations. Youdin and Johansen (2007) tested
their numerical algorithm for two-way drag forces against
two modes of the streaming instability, the modes having

13

105
643
1283
2563

max(ρp/ρg)10410ρp = 516.2 ρgρp = 2599.9 ρgρp = 12509.0 ρg

3

102
101
Z=0.01

100
0

10

20

Z=0.0230
40t/Torb50

60

70

Fig. 8.— Particle concentration by the streaming instability is mainly elongated along the Keplerian flow direction, as seen in the left
panel showing the column density of particles with St = 0.3 in a local frame where the horizontal axis represents the radial direction andthe vertical axis the orbital direction. The substructure of the overdense region appears fractal with filamentary structure on many scales.
The right panel shows the measured maximum particle density as a function of time, for three different grid resolutions. The metallicity
is gradually increased from the initial Z = 0.01 to Z = 0.02 between 30 and 40 orbits, triggering strong particle concentration. Higher
resolution resolves smaller embedded substructures and hence higher densities.

4.2.1.

Stratified simulations

sand times the local gas density even for relatively small
particles of St = 0.3 (Johansen et al., 2012a). Measurements of the maximum particle density as a function of time
are shown in Fig. 8, together with a column density plot of
the overdense filaments.
Particles of St = 0.3 reach only modest concentration in non-stratified simulations. The explanation for thehigher concentration seen in stratified simulations may bethat the mid-plane layer is very thin, on the order of 1%
of a gas scale height, so that slowly drifting clumps are
more likely to merge in the stratified simulations. The maximum particle concentration increases when the resolution
is increased, by a factor approximately four each time thenumber of grid cells is doubled in each direction. Thisscaling may arise from thin filamentary structures which
resolve into thinner and thinner filaments at higher resolution. The overall statistical properties of the turbulence
nevertheless remain unchanged as the resolution increases(Johansen et al., 2012a).
The ability of the streaming instability to concentrateparticles depends on both the particle size and the local
metallicity in the disc. Carrera, Johansen and Davies (inpreparation) show that particles down to St = 0.01 areconcentrated at a metallicity slightly above the solar value,
while smaller particles down to St = 0.001 require significantly increased metallicity, e.g. by photoevaporation ofgas or pile up of particles from the outer disc. Whether the
streaming instability can explain the presence of mm-sizedchondrules in primitive meteorites is still not known. Alternative models based on small-scale particle concentration

While non-stratified simulations are excellent for testing
the robustness of a numerical algorithm and comparing toresults obtained with different codes, stratified simulations
including the component of the stellar gravity towards the
mid-plane are necessary to explore the role of the streaming instability for planetesimal formation. In the stratified
case the mid-plane particle mass-loading is no longer a parameter that can be set by hand. Rather its value is determined self-consistently in a competition between sedimentation and turbulent diffusion. The global metallicity
Z = Σp /Σg is now a free parameter and the mass-loadingin the mid-plane depends on the scale height of the layer
and thus on the degree of turbulent stirring.
Stratified simulations of the streaming instability display
a binary behaviour determined by the heavy element abundance. At Z around the solar value or lower, the mid-planelayer is puffed up by strong turbulence and shows little particle concentration. Bai and Stone (2010a) showed that theparticle mid-plane layer is stable to Kelvin-Helmholtz instabilities thriving in the vertical shear of the gas velocity,so the mid-plane turbulence is likely a manifestation of thestreaming instability which is not associated with particle
concentration. Above solar metallicity very strong particle
concentrations occur in thin and dominantly axisymmetricfilaments (Johansen et al., 2009a; Bai and Stone, 2010a).
The metallicity threshold for triggering strong clumping
may be related to reaching unity particle mass-loading in
the mid-plane, necessary for particle pile-ups by the streaming instability. The particle density can reach several thou-

14

0.10

0.10
250
200

0.05

0.05

0.00

−0.05

R/km

0.00

y/H

y/H

150

100

−0.05

50
−0.10
−0.10

−0.05

0.00
x/H

0.05

0.10

−0.10

−0.05

0.00
x/H

0.05

−0.10
0.10

325

330
t/Ω−1

335

Fig. 9.— Planetesimal formation in a particle filament formed by the streaming instability. The left panel shows the column density ofparticles with St = 0.3, with dots marking the newly formed planetesimals. The middle panel shows the positions and Hill spheres ofthe planetesimals, while the right panel shows their contracted radii when the scale-free simulation is applied to the asteroid belt.
(Cuzzi et al., 2008) and particle sedimentation (Youdin andShu, 2002; Chiang and Youdin, 2010; Lee et al., 2010) may
be necessary.
4.3.

cold Kuiper belt (Noll et al., 2008). The limited mass resolution which is allowed in direct N -body simulations nevertheless highlights the necessity for further studies of thefate of collapsing pebble clumps. Information from laboratory experiments on coagulation, bouncing and fragmentation of particles is directly applicable to self-gravitating
particle clumps as well (Güttler et al., 2010), and including
realistic particle interaction in collapse simulations should
be a high priority for future studies.

Gravitational collapse

Planetesimals forming by gravitational contraction andcollapse of the overdense filaments are generally found to
be massive, corresponding to contracted radii between 100
km and 1000 km, depending on the adopted disc model (Johansen et al., 2007, 2009a, 2011, 2012a; Kato et al., 2012).
Fig. 9 shows the results of a high-resolution simulation ofplanetesimal formation through streaming instabilities, withplanetesimals from 50 to 200 km radius forming when themodel is applied to r = 3 AU. The formation of large planetesimals was predicted already by Youdin and Goodman(2005) based on the available mass in the linear modes. The
planetesimals which form increase in size with increasing
disc mass (Johansen et al., 2012a), with super-Ceres-sizedplanetesimals arising in massive discs approaching Toomre
instability in the gas (Johansen et al., 2011). Increasing thenumerical resolution maintains the size the largest planetesimals, but allows smaller-mass clumps to condense out of
the turbulent flow alongside their massive counterparts (Johansen et al., 2012a). Although most of the particle massenters planetesimals with characteristic sizes above 100 km,
it is an open question, which can only be answered withvery-high-resolution computer simulations, whether there
is a minimum size to the smallest planetesimals which canform in particle concentration models.
Strictly speaking simulations of particle concentrationand self-gravity only find the mass spectrum of gravitationally bound clumps. Whether these clumps will contract to form one or more planetesimals is currently not
known in details. Nesvorný et al. (2010) took bound particle clumps similar to those arising in numerical simulations and evolved in a separate N -body code with collisions
and merging. The rotating clumps typically collapse to abinary planetesimal, with orbital properties in good agreement with the observed high binary fraction in the classical5.

PARTICLE GROWTH

The size distribution of solid particles in protoplanetary
discs evolves due to processes of coagulation, fragmentation, sublimation and condensation. The particle concentration mechanisms discussed in the previous section are
only relevant for particle sizes from millimeters to meters, much larger than the canonical sub-micron-sized dust
and ice grains which enter the protoplanetary disc. Hence
growth to macroscopic sizes must happen in an environment of relatively uniform particle densities.
Primitive meteorites show evidence of thermal processing of protoplanetary disc solids, in the form of CAIs which
likely condensed directly from the cooling gas near the
proto-Sun and chondrules which formed after rapid heating
and melting of pre-existing dust aggregates (section 2.2).
Thermal processing must also have taken place near the water snow line where the continued process of sublimationand condensation can lead to the formation of large haillike ice particles (Ros and Johansen, 2013); or more distantly around the CO snow line identified observationally at
30 AU around the young star TW Hya (Qi et al., 2013).
Direct sticking is nevertheless the most general mechanism for dust growth, as the densities in protoplanetary
discs are high enough for such coagulation to be important
over the life-time of the gaseous disc. Physically, the collisional evolution is determined by the following factors:• The spatial density of particles. This will depend on
both the location in the disc and on the degree of par-15

However, binning in additional dimensions besides mass
may render this approach impractical. To circumvent themulti-dimensional binning Okuzumi et al. (2009) outlined
an extension to Smoluchowski’s equation, which treats themean value of the additional parameters at every mass m.
For example, the dust at a mass bin m is in addition characterized by a filling factor φ, which is allowed to change withtime. If the distribution in φ at a certain mass m is expected
to be narrow, such an approach is advantageous as one may
not have to deviate from Smoluchoski’s 1D framework.
A more radical approach is to use a direct-simulation
Monte Carlo (MC) method, which drops the concept of the
distribution function altogether (Ormel et al., 2007). Instead the MC-method computes the collision probability between each particle pair of the distribution (essentially the
kernel function Kij ), which depends on the properties ofthe particles. Random numbers then determines which particle pair will collide and the time-step that is involved. The
outcome of the collision between these two particles mustbe summarized in a set of physically-motivated recipes (theanalogy with equation 28 is simply m1 + m2 → m). A new
set of collision rates is computed, after which the procedure
repeats itself.
Since there are many more dust grains in a disc than any
computer can handle, the computational particles should bechosen such to accurately represent the physical size distribution. A natural choice is to sample the mass of the distribution (Zsom and Dullemond, 2008). This method also allows fragmentation to be naturally incorporated. The drawback however is that the tails of the size distribution are notwell resolved. The limited dynamic range is one of the key
drawbacks of MC-methods. Ormel and Spaans (2008) havedescribed a solution, but at the expense of a more complex
algorithm.

ticle concentration (see section 4).
• The collision cross section among particles, σcol . An
appropriate assumption is to take spherical particles,
for which σcol,ij = π(Ri + Rj )2 where Ri and Rj
are simply the radii of two particles in the absence ofgravitational focusing.
• The relative velocity among the particles, δvij . The
relative velocity affects both the collision rate and thecollision outcome.
• The collision outcome, which is determined by thecollision energy, porosity, and other collision parameters like the impact parameter.
The product of δvij and σij , which together with the particle density determine the collision rate, is most often referred to as the collision kernel Kij .
5.1.

Numerical approaches the coagulation equation

The evolution of the dust size distribution is most commonly described by the Smoluchowski equation (Smoluchowski, 1916), which reads:
Z
dn(m)
1
dm′ K(m′ , m − m′ )n(m′ )n(m − m′ )=dt2Z−n(m) dm′ K(m, m′ )n(m′ ), (28)where n(m) is the particle number density distribution (i.e.,n(m)dm gives the number density of particles between
mass m and m + dm). The terms on the right-hand sideof equation (28) simply account for the particles of mass
m that are removed due to collisions with any other particle and those gained by a collision between two particles
of mass m′ and m − m′ (the factor of 1/2 prevents double
counting). In numerical applications one most often subdivides the mass grid in bins (often logarithmically spaced)
and then counts the interactions among these bins.
However, equation (28), an already complex integrodifferential equation, ignores many aspects of the coagulation process. For example, collisional fragmentation is
not included and the dust spatial distribution is assumed tobe homogeneous (that is, there is no local concentration).
Equation (28) furthermore assumes that each particle pair(m, m′ ) collides at a single velocity and ignores the effects
of a velocity distribution two particles may have. Indeed,
perhaps the most important limitation of equation (28) isthat the time-evolution of the system is assumed to be only
a function of the particles’ masses. For the early coagulation phases this is certainly inappropriate as the dust internal structure is expected to evolve and influence the collision outcome. The collision outcomes depend very strongly
on the internal structure (porosity, fractal exponent) of dust
aggregates as well on their composition (ices or silicates,
Wada et al., 2009; Güttler et al., 2010).
In principle the Smoluchowski coagulation equation canbe extended to include these effects (Ossenkopf , 1993).

5.2.

Laboratory experiments

Numerical solutions to the coagulation equation rely crucially on input from laboratory experiments on the outcomeof collisions between dust particles. Collision experimentsin the laboratory as well as under microgravity conditions
over the past 20 years have proven invaluable for the modeling of the dust evolution in protoplanetary discs. Here
we briefly review the state-of-the-art of these experiments(more details are given in the chapter by Testi et al.). A detailed physical model for the collision behavior of silicate
aggregates of all masses, mass ratios and porosities can befound in Güttler et al. (2010) and recent modifications are
published by Kothe et al. (2013).
The three basic types of collisional interactions between
dust aggregates are:
1. Direct collisional sticking: when two dust aggregatescollide gently enough, contact forces (e.g., van der
Waals forces) are sufficiently strong to bind the aggregates together. Thus, a more massive dust aggregate has formed. Laboratory experiments summarized by Kothe et al. (2013) suggest that there is16

a mass-velocity threshold v ∝ m−3/4 for spherical medium-porosity dust aggregates such that only
dust aggregates less massive and slower than thisthreshold can stick. However, hierarchical dust aggregates still stick at higher velocities/masses thangiven by the threshold for homogeneous dust aggregates (Kothe et al., 2013).
2. Fragmentation: when the collision energy is sufficiently high, similar-sized dust aggregates break up
so that the maximum remaining mass decreases under the initial masses of the aggregates and a powerlaw type tail of smaller-mass fragments is produced.
The typical fragmentation velocity for dust aggregates consisting of µm-sized silicate monomer grainsis 1 m s−1 .
3. Bouncing: above the sticking threshold and below
the fragmentation limit, dust aggregates bounce off
one another. Although the conditions under whichbouncing occurs are still under debate (see sections 6
and 7), some aspects of bouncing are clear. Bouncing
does not directly lead to further mass gain so that thegrowth of dust aggregates is stopped (“bouncing barrier”). Furthermore, bouncing leads to a steady compression of the dust aggregates so that their porosity
decreases to values of typically 60% (Weidling et al.,2009).

Fig. 10.— Experimental example of mass transfer in fragmenting collisions. All experiments were performed in vacuum. (a)A mm-sized fluffy dust aggregate is ballistically approaching the
cm-sized dusty target at a velocity of 4.2 m/s. Projectile and target consist of monodisperse SiO2 spheres of 1.5 µm diameter. (b)Shortly after impact, most of the projectile’s mass flies off the target in form of small fragments (as indicated by the white arrows);
part of the projectile sticks to the target. (c) - (e) The same target
after 3 (c), 24 (d), 74 (e) and 196 (f) consecutive impacts on thesame spot. Image credit: Stefan Kothe, TU Braunschweig.

Although the details of the three regimes are complex –regime boundaries are not sharp, and other more subtle processes exist – this simplified physical picture shows that inprotoplanetary discs, in which the mean collision velocity
increases with increasing dust-aggregate size (see Fig. 4), a
maximum aggregate mass exists. Detailed numerical simulations using the dust-aggregate collision model by Güttler
et al. (2010) find such a maximum mass (Zsom et al., 2010;
Windmark et al., 2012b). For a minimum mass solar nebulamodel at 1 AU, the maximum dust-aggregate size is in therange of millimeters to centimeters (e.g. Zsom et al., 2010).
Even if the bouncing barrier can be overcome, anotherperhaps even more formidable obstacle will present itself at
the size scale where collision velocities reach ∼1 m/s (Fig.
4). At this size collisions among two silicate, similar-size,
dust particles are seen to fragment, rather than accrete. Consequently, growth will stall and the size distribution will,
like with the bouncing case, settle into a steady-state (Birnstiel et al., 2010, 2011, 2012). Highly porous icy aggregates
may nevertheless still stick at high collision speeds (see section 7).
6.

depending mainly on the collision velocity and the particle
size ratio. Using the most updated laboratory knowledge ofcollision outcomes, Zsom et al. (2010, 2011) solved the coagulation equation with a Monte Carlo method and foundthat in the ice-free, inner disc regions the dust size distribution settles into a state dominated by mm-size particles.
Further growth is impeded because collisions among two
such mm-size particles will mainly lead to bouncing andcompactification.
The real problem behind this “bouncing barrier” may,counter-intuitively, not be the bouncing but rather the absence of erosion and fragmentation events. Based upon laboratory experiments (Wurm et al., 2005; Teiser and Wurm,
2009b,a; Kothe et al., 2010; Teiser et al., 2011; Meisneret al., 2013), it has become clear that above the fragmentation limit, the collision between two dust aggregates canlead to a mass gain of the larger (target) aggregate if thesmaller (projectile) aggregate is below a certain threshold.
For impact velocities in the right range and relatively small
projectile aggregates, up to 50% of the mass of the projectile can be firmly transferred to the target. This process
has been shown to continue to work after multiple collisions
with the target and under a wide range of impact angles (seeFig. 10).
Growth by mass transfer requires a high fraction of thedust mass present in small dust aggregates which can be
swept up by the larger particles. The number of small dust

GROWTH BY MASS TRANSFEROver the last two decades a large effort has been invested
to study experimentally the collisions among dust particles
and to assess the role of collisions in planetesimal formation. As was seen in section 5 a plethora of outcomes –sticking, compaction, bouncing, fragmentation – is possible

17Fig. 3 shows that compact ice particles grow to Stokes
number 10 in the Epstein regime outside of r = 5 AU
and in the Stokes or non-linear regimes inside of this radius. Fluffy ice particles with a very low internal density ofρ• = 10−5 g cm−3 grow to Stokes number 10 in the nonlinear or quadratic regimes in the entire extent of the disc.
In the Epstein regime the time to grow to Stokes number 10,
τ10 = R10 /Ṙ, is independent of both material density andgas density, ǫ −1  r 5/4
40H
(Ep)
τ10 =
. (33)≈ 25000 yr
ǫZv12
0.5
AU

aggregates can be maintained by catastrophic collisions between larger bodies. A simple two-component model forsuch coagulation-fragmentation growth was developed by
Johansen et al. (2008). They divided the size distribution
into small particles and large particles and assumed that (1)a small particle will stick to a large particle, (2) a small particle will bounce off another small particle, and (3) a large
particle will shatter another large particle. Small particles
collide with large particles at speed v12 , while large particles collide with each other at speed v22 . Under these conditions an equilibrium in the surface density ratio of large andsmall particles, Σ1 /Σ2 , can be reached, with equal massflux from large to small particles by fragmentation as fromsmall to large particles by mass transfer in collisions.
In the equilibrium state where the mass flux from large
to small particles balances the mass flux from small to large
particles, the growth rate of the large particles is√Σp /( 2πH1 )
Ṙ =
×
ρ•
v22
. (29)[2/(1 + H12 /H22 )]1/2 + 4v22 /v12

The growth-time is much shorter in the Stokes regime
where the Stokes number is proportional to the squared particle radius, givingp
360Hρ• (m/σ)(St)≈ 500 yr ×τ10=ǫZρg v12
1/2  ǫ −1  ρr 27/8
•
=
,(34)
−3
0.5
1 g cmAU

with a direct dependence on both material density and gas
density. Here m and σ are the mass and collisional crosssection of a hydrogen molecule. Regarding the sticking
coefficient, Wurm et al. (2005) find mass transfer efficiencies of ǫ ≈ 0.5 for collision speeds up to v12 = 25m/s, hence we have used ǫ = 0.5 as a reference value inthe above equations. It is clearly advantageous to have a
low internal density and cross the radial drift barrier in the
Stokes regime (Okuzumi et al., 2012). The non-linear andquadratic regimes are similarly beneficial for the growth,with the time-scale for crossing the meter barrier in thequadratic drag force regime approximately a factor 6cs /∆v
faster than in the Epstein regime (equation 33).
The two-component model presented in Johansen et al.
(2008) is analytically solvable, but very simplified in thecollision physics. Windmark et al. (2012b) presented solutions to the full coagulation equation including a velocitydependent mass transfer rate. They showed that artificiallyinjected seeds of centimeter sizes can grow by sweeping up
smaller mm-sized particles which are stuck at the bouncing barrier (Zsom et al., 2010), achieving at 3 AU growth
to 3 meters in 100,000 years and to 100 meters in 1 million
years. The time-scale to grow across the radial drift barrier
is broadly in agreement with the simplified model presented
above in equation (34). This time is much longer than theradial drift time-scale which may be as short as 100–1,000
years. It is therefore necessary to invoke pressure bumps toprotect the particles from radial drift while they grow slowly
by mass transfer. However, in that case the azimuthal driftof small particles vanishes and a relatively high turbulent
viscosity of α = 10−2 must be used to regain the flux onto
the large particles. The pressure bumps must additionallybe very long-lived, on the order of the life-time of the protoplanetary disc.
While the discovery of a pathway for growth by masstransfer at high-speed collisions is a major experimental

Here H1 and H2 are the scale heights of the small particles
and the large particles, respectively, and Σp = Σ1 + Σ2 is
the total column density of the dust particles. The complexity of the equation arises from the assumption that the small
grains are continuously created in collisions between thelarge particles. Setting H1 = H and H2 ≪ H1 the equations are valid for a very general sweep-up problem wherea few large particles sweep up a static population of small
particles, resulting in growth of the large particles at rateṘ

 ǫ   r −2.75ǫZρg×v12 = 2.7 mm yr−1
4ρ•
0.5
AU

−1 

ρ•
v12 
Z
.
(30)
0.01
1 g cm−3
50 m s−1=

We introduced here a sticking coefficient ǫ which measures either the sticking probability or the fraction of mass
which is transferred from the projectile to the target. The
growth rate is constant and depends only on the particle size
through v12 . If the small particles are smaller than Stokes
number unity and the large particles are larger, then we canassume that the flux of small particles is carried with the
sub-Keplerian wind, v12 ≈ ∆v, and solve equation (30) for
the time-scale to cross the radial drift barrier. The end of theradial drift barrier can be set at reaching roughly St = 10.
This Stokes number corresponds in the Epstein regime tothe particle size
(Ep)

R10

=

10Hρg
.
ρ•

In the Stokes regime the particle size iss
90Hρg λ
(St)
.
R10 =
4ρ•

(31)

(32)

18

breakthrough, it seems that this effect can not in itself leadto widespread planetesimal formation, unless perhaps in the
innermost part of the protoplanetary disc where dynamical time-scales are short (equation 34). The planetesimal
sizes which are reached within a million years are also quite
small, on the order of 100 meters. These sizes are too small
to undergo gravitational focusing even in weakly turbulent
discs (see Fig. 5). An additional concern is the erosion of
the preplanetesimal by tiny dust grains carried with the subKeplerian wind (Schräpler and Blum, 2011; Seizinger et al.,2013).
Particles stuck at the bouncing barrier are at the lower
end of the sizes that can undergo particle concentration andgravitational collapse. However, if a subset of “lucky” particles experience only low-speed, sticking collisions andmanage to grow past the bouncing barrier, this leads eventually to the formation of a population of cm-m-sized particles
that could concentrate in pressure bumps or participate in
the streaming instability (Windmark et al., 2012a; Garaud
et al., 2013).
7.

FLUFFY GROWTH
Fig. 11.— Examples of collision outcomes of icy PCA clusters

The last years have seen major improvements in N -bodymolecular-dynamics simulations of dust aggregate collisions by a number of groups (e.g. Wada et al., 2008, 2009,2011; Paszun and Dominik, 2008, 2009; Seizinger and Kley,
2013). These N -body simulations show that fluffy aggregates have the potential to overcome barriers in dust growth.
In N -body simulations of dust aggregates, all surface interactions between monomers in contact in the aggregates arecalculated, by using a particle interaction model.
N -body simulations of dust aggregates have some merits
compared to laboratory experiments. Precise informationsuch as the channels through which the collisional energy
is dissipated can be readily obtained and subsequent data
analysis is straightforward. Another advantage of N -bodysimulations is that one can study highly fluffy aggregates,
which are otherwise crushed under the Earth’s gravity. On
the other hand, an accurate interaction model of constituent
particles is required in N -body simulations of dust aggregates. In the interaction model used in most N -body simulations, the constituent particles are considered as adhesive
elastic spheres. The adhesion force between them is described by the JKR theory (Johnson et al., 1971). As fortangential resistive forces against sliding, rolling, or twisting motions, Dominik and Tielens’s model is used (Dominik
and Tielens, 1995, 1996, 1997; Wada et al., 2007). In orderto reproduce the results of laboratory experiments better,the interaction model must be calibrated against laboratory
experiments (Paszun and Dominik, 2008; Seizinger et al.,2012; Tanaka et al., 2012).
In protoplanetary discs, dust aggregates are expected tohave fluffy structures with low bulk densities if impact compactification is negligible (e.g. Okuzumi et al., 2009, 2012;
Zsom et al., 2010, 2011). Especially in the early stage ofdust growth, low-speed impacts result in hit-and-stick ofconsisting of 8000 particles for two values of the impact parameterb (lower panels). The upper panels represent initial aggregates.
The collision velocity of 70 m/s in both cases and R is the radiusof the initial aggregates. The head-on collision (left panel) resultsin sticking with minor fragmentation while aggregates tend to passby each other in off-set collisions (right panels). (This figure is
reproduced based on Fig. 2 of Wada et al. 2009 by permission of
the AAS.)

dust aggregates, with little or no compression. This growthmode makes fluffy aggregates with a low fractal dimensionof ∼2. Thus it is necessary to examine the outcome of collisions between fluffy dust aggregates and the resulting compression. The results of N -body simulations can then be
used in numerical models of the coagulation equation including the evolution of the dust porosity.
7.1.

Critical impact velocity for dust growth

Dominik and Tielens (1997) first carried out N -bodysimulation of aggregate collisions, using the particle interaction model they constructed. They also derived a simple recipe for outcomes of aggregate collisions from theirnumerical results, though they only examined head-on collisions of two-dimensional small aggregates containing as
few as 40 particles. The collision outcomes are classified
into hit-and-stick, sticking with compression, and catastrophic disruption, depending on the impact energy. Bouncing of two aggregates was not observed in their simulations.
According to the DT (Dominik and Tielens, 1997) recipe,
growth is possible at collisions with Eimp < Ank Ebreak ,where A∼10 and nk is the total number of contacts in two
colliding aggregates. For relatively fluffy aggregates, nk
is approximately equal to the total number of constituent
19

particles in the two aggregates, N . The energy for breaking
one contact between two identical particles, Ebreak , is givenby (e.g. Chokshi et al., 1993; Wada et al., 2007)
Ebreak = 23[γ 5 r4 (1 − ν 2 )2 /E 2 ]1/3 ,silicate aggregates (e.g. Blum and Wurm, 2000; Güttler
et al., 2010; Paszun and Dominik, 2009; Seizinger and Kley,
2013).
In the case of CCA clusters, Wada et al. (2009) foundthat the constant A is the same as in the DT recipe (≃10).
Considering much smaller volume filling factors of CCA
clusters than PCA, it indicates that the critical velocity is
only weakly dependent on the volume filling factor.

(35)

where r, γ, E, and ν are the radius, the surface energy,
Young’s modulus, and Poisson’s ratio of constituent particles, respectively. For icy particles the breaking energy is
obtained as 5.2 × 10−10 (r/0.1 µm)4/3 erg while the value
for silicate particles is 1.4 × 10−11 (r/0.1µm)4/3 erg, so icy
aggregates are much more sticky than silicate aggregates
Using the impact velocity, vimp , the impact energy is ex2
pressed as (N m/4)vimp/2, where m is the mass of a constituent particle and N m/4 is the reduced mass of two
equal-sized aggregates. Substituting this expression into
the above energy criterion for dust growth and noting thatnk ≃ N , we obtain the velocity criterion asp
vimp < 8AEbreak /m.
(36)

7.2.

Compression of dust aggregates

The DT recipe does not describe the amount of changes
in the porosity (or the volume filling factor) at aggregate
collisions. The first attempt to model porosity changes wasdone by Ormel et al. (2007), using simple prescriptions forthe collision outcome of porosity.
Wada et al. (2008) examined compression at head-on
collisions of two equal sized CCA clusters, using a highnumber of N -body simulations. Compression (or restructuring) of an aggregate occurs through rolling motions between constituent particles. Thus it is governed by therolling energy Eroll (i.e., the energy for rolling of a particle over a quarter of the circumference of another particle in contact). At low-energy impacts with Eimp . Eroll ,aggregates just stick to each other, as indicated by the DTrecipe. For higher-energy impacts, the resultant aggregatesare compressed, depending on the impact energy. The
radius of resulting aggregates, R, is fitted well with thepower-law function

Here A is a dimensionless parameter that needs to be calibrated with experiments. Note that this velocity criterion
is dependent only on properties of constituent particles and,thus, independent of the mass of the aggregates.
To evaluate the critical impact velocity for growth moreaccurately, Wada et al. (2009) performed N -body simulations of large aggregates made of up to ∼104 sub-micron
icy particles, including off-set collisions (see Fig. 11). They
consider two kinds of aggregate structure, the so-called
CCA and PCA clusters. The CCA (Cluster-Cluster Agglomeration) clusters have an open structure with a fractal
dimension of 2 while the PCA (Particle-Cluster Agglomeration) clusters have a fractal dimension of 3 and a volume
filling factor of 0.15. The PCA clusters are rather compact,
compared to the CCA. Since dust aggregates are expected to
be much more compact than CCA clusters due to compression (see section 7.2), the growth and disruption process of
PCA clusters is of particular importance for elucidating theplanetesimal formation.
For PCA clusters (composed of 0.1µm-sized icy particles), the critical impact velocity is obtained as 60 m/s fromtheir simulations, independent of the aggregate mass within
the mass range examined in the simulations of Wada et al.
(2009). This indicates that icy dust aggregates can circumvent the fragmentation barrier and grow towards planetesimal sizes via collisional sticking. Note that the critical velocity actually increases with the aggregate mass when only
head-on collisions are considered, as seen in Fig. 11. For an
accurate evaluation of the critical velocity for growth, offset collisions should be considered as well, as in Wada et al.
(2009).
The above result also fixes the corresponding constant
A in the critical energy at 30. By using Ebreak of silicate particles in equation (36) with A = 30, the criticalvelocity for growth of silicate aggregates is obtained as
vimp = 1.3 (r/0.6µm)−5/6 m/s. It agrees well with thelaboratory experiments and numerical simulations for small

R ≃ 0.8[Eimp /(N Eroll )]−0.1 N 1/2.5 r ,(37)

where r is the radius of a monomer and N the number of
constituent particles. Analysing the structure of the compressed aggregates, it can be shown that the compressed
aggregates have a fractal dimension of 2.5. This fractal dimension is consistent with equation (37) since equation (37)gives the relation of N ∝ R2.5 for maximally compressed
aggregates (with Eimp ∼N Eroll ). The compression observedin the N -body simulation is less much extreme than in thesimple porosity model used by Ormel et al. (2007) becauseof the low fractal dimension of 2.5.
While Wada et al. (2008) examined the compression ina single collision, dust aggregates will be gradually compressed by successive collisions in realistic systems. In order to examine such a gradual compression process during
growth, Suyama et al. (2008) performed N -body simulations of sequential collisions of aggregates. Even after multiple collisions, the compressed aggregates maintain a fractal dimension of 2.5. Suyama et al. (2012) further extend
their porosity model to unequal-mass collisions.
Okuzumi et al. (2012) applied the porosity evolution
model of Suyama et al. (2012) to dust growth outside of
the snow line in a protoplanetary disc. As a first step, collisional fragmentation is neglected because a relatively highimpact velocity (∼60 m/s) is required for significant disruption of icy dust aggregates. They found that dust particles
evolve into highly porous aggregates (with bulk densities
20much less than 0.1 g/cm3 ) even if collisional compression
is taken into account. This is due to the ineffective compression at aggregate collisions in the porosity model of Suyama
et al. (2012). Another important aspect of fluffy aggregatesis that they cross the radial drift barrier in the Stokes or nonlinear drag force regimes (see Fig. 3 and discussion in section 6) where it is easier to grow to sizes where radial drift
is unimportant. This mechanism accelerates dust growth at
the radial drift barrier effectively and enables fluffy aggregates to overcome this barrier inside 10 AU (Okuzumi et al.,2012; Kataoka et al., 2013).
7.3.

the collision speeds reach a maximum of the escape speedof the forming body. In the cometesimal formation zone
any subsequent impacts by high-speed particles, brought inwith the sub-Keplerian flow, will at most add a few metersof compact debris to the pebble pile over the life-time of thesolar nebula (equation 30 applied to r = 10 AU).
This difference in impact velocity leads to substantial differences in the tensile strengths of the planetesimals/cometesimals. Due to the relatively high impact
speeds in models (1) and (2), the growing bodies are compacted and possess tensile strengths on the order of 1–10
kPa (Blum et al., 2006). Skorov and Blum (2012) pointed
out that these tensile-strength values are too high to explain the continuous dust emission of comets approaching the Sun and favor model (3) above, for which the tensile strengths of loosely packed mm-cm pebbles should bemuch smaller.
Skorov and Blum (2012) base their comet-nucleus model
on a few assumptions: (a) dust aggregates in the inner partsof the protoplanetary disc can only grow to sizes of mmcm; (b) turbulent diffusion can transport dust aggregates tothe outer disc; (c) in the outer disc, dust and ice aggregates
become intermixed with a dust-to-ice ratio found in comets;(d) cometesimals form via gravitational instability of overdense regions of mm-cm-sized particles which will form
planetesimals with a wide spectrum of masses. While kmscale planetesimals have not yet been observed to form inhydrodynamical simulations of planetesimal formation, thismay be an artifact of the limited numerical resolution (seediscussion in section 4.3).
If the comet nuclei as we find them today have been dormant in the outer reaches of the Solar System since theirformation and have not been subjected to intense bombardment and aqueous or thermal alteration, then they today represent the cometesimals of the formation era of the Solar
System. Skorov and Blum (2012) derive for their model
tensile strengths of the ice-free dusty surfaces of comet nuclei r −2/3,(38)
T = T01 mm
with T0 = 0.5 Pa and r denoting the radius of the dust andice aggregates composing the cometary surface.
These extremely low tensile-strength values are, according to the thermophysical model of Skorov and Blum(2012), just sufficiently low to explain a continuous outgassing and dust emission of comet nuclei inside a criticaldistance to the Sun. This model provides strong indication that km-sized bodies in the outer Solar System wereformed by the gravitational contraction of smaller sub-units(pebbles) at relatively low (order of m/s) velocities. Suchpebble-pile planetesimals are in broad agreement with thepresence of large quantities of relatively intact mm-sizedchondrules and CAIs in chondrites from the asteroid belt.
The cometesimals could thus represent the smallest bodies which formed by gravitational instability of mm-cmsized icy/rocky pebbles. These pebbles will continue tocollide inside the gravitationally collapsing clump. Low-

Bouncing condition in N -body simulations

Until 2011, bouncing events were not reported in N-bodysimulations of aggregate collisions, while bouncing events
have been frequently observed in laboratory experiments(e.g. Blum and Münch, 1993; Langkowski et al., 2008; Weidling et al., 2009, 2012). Wada et al. (2011) found bouncing events in their N-body simulations of both icy and silicate cases for rather compact aggregates with filling factor
φ ≥ 0.35. Such compact aggregates have a relatively large
coordination number (i.e., the mean number of particles incontact with a particle), which inhibits the energy dissipation through the rolling deformation and helps bouncing.
Seizinger and Kley (2013) further investigated the bouncing
condition and proposed that the more realistic condition is
φ > 0.5.
This critical volume filling factor for bouncing is a few
times as large as that in the laboratory experiments (e.g.
Langkowski et al., 2008). The origin of this discrepancy
between these two approaches is not clear yet. From a qualitative point of view, however, laboratory experiments showthat fluffy silicate aggregates with φ < 0.1 tend to stick toeach other (Blum and Wurm, 2000; Langkowski et al., 2008;
Kothe et al., 2013). This qualitative trend is consistent with
the N-body simulations. In protoplanetary discs, dust aggregates are expected to be highly porous. In the growth ofsuch fluffy aggregates, the bouncing barrier would not be astrong handicap.
8.

TOWARDS A UNIFIED MODEL

Three main scenarios for the formation of planetesimals have emerged in the last years: (1) formation via
coagulation-fragmentation cycles and mass transfer fromsmall to large aggregates (Wurm et al., 2005; Windmark
et al., 2012b), (2) growth of fluffy particles and subsequent compactification by self-gravity (Wada et al., 2008,2009), and (3) concentration of pebbles in the turbulent gas
and gravitational fragmentation of overdense filaments (Johansen et al., 2007, 2009a; Kato et al., 2012).
The collision speeds at which planetesimal-sized bodies
form are typically 50 m/s in models (1) and (2), as small
dust aggregates are carried onto the growing planetesimalwith the sub-Keplerian wind, or even higher if the turbulent density fluctuations are strong. During the gravitational
collapse of pebble clouds, formation mechanism (3) above,21

mass clumps experience low collision speeds and contract
as kinetic energy is dissipated in inelastic collisions, forming small planetesimals which are made primarily of pristine pebbles. Clumps of larger mass have higher collision
speeds between the constituent pebbles, and hence growthby coagulation-fragmentation cycles inside massive clumps
determines the further growth towards one or more solid
bodies which may go on to differentiate by decay of shortlived radionuclides. This way coagulation is not only important for forming pebbles that can participate in particle
concentration, but continues inside of the collapsing clumps
to determine the birth size distribution of planetesimals.
There is thus good evidence that pebbles are the primarybuilding blocks of planetesimals both inside the ice line (today’s asteroids) and outside the ice line (today’s comets andKuiper belt objects). This highlights the need to understandpebble formation and dynamics better. The bouncing barrier is a useful way to maintain a high number of relatively
small pebbles in the protoplanetary disc (Zsom et al., 2010).
As pebbles form throughout the protoplanetary disc, radial
drift starts typically when reaching mm-cm sizes. Suchsmall pebbles drift relatively fast in the outer disc but slow
down significantly (to speeds of order 10 cm/s) in the inner
few AU. Pebbles with a high ice fraction fall apart at theice line around 3 AU, releasing refractory grains, which can
go on to form new (chondrule-like) pebbles inside the iceline (Sirono, 2011), as well as water vapor which boosts
formation of large, icy pebbles outside the ice line (Rosand Johansen, 2013). Similar release of refractories andrapid pebble growth will occur at ice lines of more volatile
species like CO (Qi et al., 2013). In the optically thin very
inner parts of the protoplanetary disc, illumination by thecentral star leads to photophoresis which causes chondrules
to migrate outwards (Loesche et al., 2013). Radial drift ofpebbles will consequently not lead to a widespread depletion of planetesimal building blocks from the protoplanetary disc. Instead pebble formation can be thought of as acontinuous cycle of formation, radial drift, destruction andreformation.
We propose therefore that particle growth to planetesimal sizes starts unaided by self-gravity, but after reaching
Stokes numbers roughly between 0.01 and 1 proceeds inside self-gravitating clumps of pebbles (or extremely fluffy
ice balls with similar aerodynamic stopping times). In thispicture coagulation and self-gravity are not mutually exclusive alternatives but rather two absolutely necessary ingredients in the multifaceted planetesimal formation process.
Acknowledgments. We thank the referee for usefulcomments that helped improve the manuscript. AJ was
supported by the Swedish Research Council (grant 20103710), the European Research Council under ERC Starting
Grant agreement 278675-PEBBLE2PLANET and by the
Knut and Alice Wallenberg Foundation. MB acknowledges
funding from the Danish National Research Foundation(grant number DNRF97) and by the European ResearchCouncil under ERC Consolidator grant agreement 616027STARDUST2ASTEROIDS. CWO acknowledges support

for this work by NASA through Hubble Fellowship grant
No. HST-HF-51294.01-A awarded by the Space Telescope
Science Institute, which is operated by the Association of
Universities for Research in Astronomy, Inc., for NASA,
under contract NAS 5-26555. HR was supported by the
Swedish National Space Board (grant 74/10:2) and the Polish National Science Center (grant 2011/01/B/ST9/05442).
This work was granted access to the HPC resources of JuRoPA/FZJ, Stokes/ICHEC and RZG P6/RZG made available within the Distributed European Computing Initiativeby the PRACE-2IP, receiving funding from the European
Community’s Seventh Framework Programme (FP7/20072013) under grant agreement no. RI-283493.
REFERENCES
A’Hearn M. F. et al. (2012) Astrophys. J., 758, 29.
Andrews S. M. and Williams J. P. (2005) Astrophys. J., 631, 1134.
Asphaug E. and Benz W. (1996) Icarus, 121, 225.
Bai X.-N. and Stone J. M. (2010a) Astrophys. J., 722, 1437.
Bai X.-N. and Stone J. M. (2010b) Astrophys. J. Suppl., 190, 297.
Bai X.-N. and Stone J. M. (2010c) Astrophys. J. Lett., 722, L220.
Bai X.-N. and Stone J. M. (2013) Astrophys. J., 769, 76.
Baker J. et al. (2005) Nature, 436, 1127.
Baker J. A. et al. (2012) Geochim. Cosmochim. Acta, 77, 415.
Balbus S. A. and Hawley J. F. (1991) Astrophys. J., 376, 214.
Barge P. and Sommeria J. (1995) Astron. Astrophys., 295, L1.
Bell K. R. et al. (1997) Astrophys. J., 486, 372.
Benz W. (2000) Space Sci. Rev., 92, 279.
Binzel R. P. and Xu S. (1993) Science, 260, 186.
Birnstiel T. et al. (2010) Astron. Astrophys., 513, A79.
Birnstiel T. et al. (2011) Astron. Astrophys., 525, A11.
Birnstiel T. et al. (2012) Astron. Astrophys., 544, A79.
Bizzarro M. et al. (2005) Astrophys. J. Lett., 632, L41.
Blandford R. D. and Payne D. G. (1982) Mon. Not. R. Astron. Soc.,
199, 883.
Blum J. and Münch M. (1993) Icarus, 106, 151.
Blum J. and Wurm G. (2000) Icarus, 143, 138.
Blum J. et al. (2006) Astrophys. J., 652, 1768.
Bockelée-Morvan D. et al. (2004) Icarus, 167, 113.
Bottke W. F. et al. (2005) Icarus, 175, 111.
Brasser R. and Morbidelli A. (2013) Icarus, 225, 40.
Brauer F. et al. (2007) Astron. Astrophys., 469, 1169.
Brauer F. et al. (2008a) Astron. Astrophys., 480, 859.
Brauer F. et al. (2008b) Astron. Astrophys., 487, L1.
Brown M. E. (2008) The Largest Kuiper Belt Objects, pp. 335–
344.
Brownlee D. et al. (2006) Science, 314, 1711.
Bus S. J. and Binzel R. P. (2002) Icarus, 158, 146.
Carballido A. et al. (2006) Mon. Not. R. Astron. Soc., 373, 1633.
Chiang E. and Youdin A. N. (2010) Annual Review of Earth andPlanetary Sciences, 38, 493.
Chokshi A. et al. (1993) Astrophys. J., 407, 806.
Ciesla F. J. et al. (2004) Meteoritics and Planetary Science, 39,
1809.
Connelly J. N. et al. (2012) Science, 338, 651.
Cuzzi J. N. and Zahnle K. J. (2004) Astrophys. J., 614, 490.
Cuzzi J. N. et al. (2001) Astrophys. J., 546, 496.
Cuzzi J. N. et al. (2008) Astrophys. J., 687, 1432.
D’Alessio P. et al. (2005) in: Chondrites and the Protoplanetary
Disk, vol. 341 of Astronomical Society of the Pacific Con-

22

Lee A. T. et al. (2010) Astrophys. J., 718, 1367.
Lesur G. and Papaloizou J. C. B. (2010) Astron. Astrophys., 513,
A60.
Levison H. F. et al. (2010) Astron. J., 139, 1297.
Levison H. F. et al. (2011) Astron. J., 142, 152.
Li H. et al. (2001) Astrophys. J., 551, 874.
Lissauer J. J. et al. (2011) Astrophys. J. Suppl., 197, 8.
Loesche C. et al. (2013) Astrophys. J., 778, 101.
Luu J. X. and Jewitt D. C. (2002) Annu. Rev. Astron. Astrophys.,
40, 63.
Lyra W. and Klahr H. (2011) Astron. Astrophys., 527, A138.
Lyra W. et al. (2008a) Astron. Astrophys., 491, L41.
Lyra W. et al. (2008b) Astron. Astrophys., 479, 883.
Lyra W. et al. (2009) Astron. Astrophys., 493, 1125.
Matzel J. E. P. et al. (2010) Science, 328, 483.
McKeegan K. D. et al. (2000) Science, 289, 1334.
McKinnon W. B. et al. (2008) Structure and Evolution of Kuiper
Belt Objects and Dwarf Planets, pp. 213–241.
Meisner T. et al. (2013) Astron. Astrophys., in press.
Morbidelli A. et al. (2009) Icarus, 204, 558.
Nakagawa Y. et al. (1986) Icarus, 67, 375.
Nelson R. P. and Gressel O. (2010) Mon. Not. R. Astron. Soc., 409,
639.
Nelson R. P. and Papaloizou J. C. B. (2004) Mon. Not. R. Astron.
Soc., 350, 849.
Nesvorný D. et al. (2010) Astron. J., 140, 785.
Nesvorný D. et al. (2011) Astron. J., 141, 159.
Noll K. S. et al. (2008) Icarus, 194, 758.
Ogihara M. et al. (2007) Icarus, 188, 522.
Oishi J. S. et al. (2007) Astrophys. J., 670, 805.
Okuzumi S. and Ormel C. W. (2013) Astrophys. J., 771, 43.
Okuzumi S. et al. (2009) Astrophys. J., 707, 1247.
Okuzumi S. et al. (2012) Astrophys. J., 752, 106.
Olsen M. B. et al. (2013) Astrophys. J. Lett., 776, L1.
Ormel C. W. and Cuzzi J. N. (2007) Astron. Astrophys., 466, 413.
Ormel C. W. and Okuzumi S. (2013) Astrophys. J., 771, 44.
Ormel C. W. and Spaans M. (2008) Astrophys. J., 684, 1291.
Ormel C. W. et al. (2007) Astron. Astrophys., 461, 215.
Ossenkopf V. (1993) Astron. Astrophys., 280, 617.
Pan L. and Padoan P. (2010) Journal of Fluid Mechanics, 661, 73.
Pan L. et al. (2011) Astrophys. J., 740, 6.
Paszun D. and Dominik C. (2008) Astron. Astrophys., 484, 859.
Paszun D. and Dominik C. (2009) Astron. Astrophys., 507, 1023.
Prialnik D. et al. (2004) Modeling the structure and activity ofcomet nuclei, pp. 359–387.
Qi C. et al. (2013) Science, 341, 630.
Raettig N. et al. (2013) Astrophys. J., 765, 115.
Ragozzine D. and Brown M. E. (2009) Astron. J., 137, 4766.
Rice W. K. M. et al. (2004) Mon. Not. R. Astron. Soc., 355, 543.
Rivkin A. S. et al. (2000) Icarus, 145, 351.
Ros K. and Johansen A. (2013) Astron. Astrophys., 552, A137.
Russell C. T. et al. (2012) Science, 336, 684.
Sanders I. S. and Taylor G. J. (2005) in: Chondrites and the Protoplanetary Disk, vol. 341 of Astronomical Society of the PacificConference Series, (edited by A. N. Krot, E. R. D. Scott, and
B. Reipurth), p. 915.
Schiller M. et al. (2010) Geochim. Cosmochim. Acta, 74, 4844.
Schiller M. et al. (2011) Astrophys. J. Lett., 740, L22.
Schräpler R. and Blum J. (2011) Astrophys. J., 734, 108.
Seizinger A. and Kley W. (2013) Astron. Astrophys., 551, A65.
Seizinger A. et al. (2012) Astron. Astrophys., 541, A59.
Seizinger A. et al. (2013) Astron. Astrophys., in press.

ference Series, (edited by A. N. Krot, E. R. D. Scott, and
B. Reipurth), p. 353.
Davidsson B. J. R. et al. (2007) Icarus, 187, 306.
Desch S. J. and Connolly Jr. H. C. (2002) Meteoritics and Planetary Science, 37, 183.
Dominik C. and Tielens A. G. G. M. (1995) Philosophical Magazine, Part A, 72, 783.
Dominik C. and Tielens A. G. G. M. (1996) Philosophical Magazine, Part A, 73, 1279.
Dominik C. and Tielens A. G. G. M. (1997) Astrophys. J., 480,
647.
Dominik C. et al. (2007) Protostars and Planets V, pp. 783–800.
Dra̧żkowska J. et al. (2013) Astron. Astrophys., 556, A37.
Dubrulle B. et al. (1995) Icarus, 114, 237.
Dzyurkevich N. et al. (2010) Astron. Astrophys., 515, A70.
Dzyurkevich N. et al. (2013) Astrophys. J., 765, 114.
Fang J. and Margot J.-L. (2012) Astrophys. J., 761, 92.
Fessler J. R. et al. (1994) Physics of Fluids, 6, 3742.
Fleming T. and Stone J. M. (2003) Astrophys. J., 585, 908.
Fromang S. and Nelson R. P. (2005) Mon. Not. R. Astron. Soc.,
364, L81.
Fuentes C. I. and Holman M. J. (2008) Astron. J., 136, 83.
Garaud P. et al. (2013) Astrophys. J., 764, 146.
Gladman B. J. et al. (2009) Icarus, 202, 104.
Goldreich P. and Ward W. R. (1973) Astrophys. J., 183, 1051.
Goldstein J. I. et al. (2009) Chemie der Erde / Geochemistry, 69,
293.
Goodman J. and Pindor B. (2000) Icarus, 148, 537.
Gradie J. and Tedesco E. (1982) Science, 216, 1405.
Greenberg J. M. and Hage J. I. (1990) Astrophys. J., 361, 260.
Gressel O. et al. (2012) Mon. Not. R. Astron. Soc., 422, 1140.
Gustavsson K. and Mehlig B. (2011) Phys. Rev. E, 84, 4, 045304.
Güttler C. et al. (2010) Astron. Astrophys., 513, A56.
Haghighipour N. and Boss A. P. (2003) Astrophys. J., 598, 1301.
Hayashi C. (1981) Progress of Theoretical Physics Supplement,
70, 35.
Holst J. C. et al. (2013) Proc. Natl. Acad. Sci. USA., 110, 8819.
Ida S. et al. (2008) Astrophys. J., 686, 1292.
Johansen A. and Klahr H. (2005) Astrophys. J., 634, 1353.
Johansen A. and Youdin A. (2007) Astrophys. J., 662, 627.
Johansen A. et al. (2007) Nature, 448, 1022.
Johansen A. et al. (2008) Astron. Astrophys., 486, 597.
Johansen A. et al. (2009a) Astrophys. J. Lett., 704, L75.
Johansen A. et al. (2009b) Astrophys. J., 697, 1269.
Johansen A. et al. (2011) Astron. Astrophys., 529, A62.
Johansen A. et al. (2012a) Astron. Astrophys., 537, A125.
Johansen A. et al. (2012b) Astrophys. J., 758, 39.
Johnson K. L. et al. (1971) Royal Society of London Proceedings
Series A, 324, 301.
Kataoka A. et al. (2013) Astron. Astrophys., 554, A4.
Kato M. T. et al. (2012) Astrophys. J., 747, 11.
Klahr H. H. and Bodenheimer P. (2003) Astrophys. J., 582, 869.
Kleine T. et al. (2004) Geochim. Cosmochim. Acta, 68, 2935.
Kleine T. et al. (2009) Geochim. Cosmochim. Acta, 73, 5150.
Kothe S. et al. (2010) Astrophys. J., 725, 1242.
Kothe S. et al. (2013) ArXiv:1302.5532K.
Kowalik K. et al. (2013) ArXiv:1306.3937K.
Kretke K. A. and Lin D. N. C. (2007) Astrophys. J. Lett., 664, L55.
Krot A. N. et al. (2003) Treatise on Geochemistry, 1, 83.
Krot A. N. et al. (2005) Nature, 436, 989.
Langkowski D. et al. (2008) Astrophys. J., 675, 764.
Laughlin G. et al. (2004) Astrophys. J., 608, 489.

23

Sekanina Z. (1997) Astron. Astrophys., 318, L5.
Shakura N. I. and Sunyaev R. A. (1973) Astron. Astrophys., 24,
337.
Shankman C. et al. (2013) Astrophys. J. Lett., 764, L2.
Sheppard S. S. and Trujillo C. A. (2010) Astrophys. J. Lett., 723,
L233.
Siegler N. et al. (2007) Astrophys. J., 654, 580.
Simon J. B. et al. (2012) Mon. Not. R. Astron. Soc., 422, 2685.
Sirono S.-i. (2011) Astrophys. J. Lett., 733, L41.
Skorov Y. and Blum J. (2012) Icarus, 221, 1.
Smoluchowski M. V. (1916) Zeitschrift fur Physik, 17, 557.
Spivak-Birndorf L. et al. (2009) Geochim. Cosmochim. Acta, 73,
5202.
Squires K. D. and Eaton J. K. (1991) Physics of Fluids, 3, 1169.
Stevenson D. J. and Lunine J. I. (1988) Icarus, 75, 146.
Suyama T. et al. (2008) Astrophys. J., 684, 1310.
Suyama T. et al. (2012) Astrophys. J., 753, 115.
Tanaka H. and Ward W. R. (2004) Astrophys. J., 602, 388.
Tanaka H. et al. (2012) Progress of Theoretical Physics Supplement, 195, 101.
Teiser J. and Wurm G. (2009a) Astron. Astrophys., 505, 351.
Teiser J. and Wurm G. (2009b) Mon. Not. R. Astron. Soc., 393,
1584.
Teiser J. et al. (2011) Icarus, 215, 596.
Thomas P. C. et al. (2005) Nature, 437, 224.
Tremaine S. and Dong S. (2012) Astron. J., 143, 94.
Turner N. J. et al. (2006) Astrophys. J., 639, 1218.
van der Marel N. et al. (2013) Science, 340, 1199.
Voelk H. J. et al. (1980) Astron. Astrophys., 85, 316.
Wada K. et al. (2007) Astrophys. J., 661, 320.
Wada K. et al. (2008) Astrophys. J., 677, 1296.
Wada K. et al. (2009) Astrophys. J., 702, 1490.
Wada K. et al. (2011) Astrophys. J., 737, 36.
Weidenschilling S. J. (1977a) Mon. Not. R. Astron. Soc., 180, 57.
Weidenschilling S. J. (1977b) Astrophys. Space Sci., 51, 153.
Weidenschilling S. J. (1980) Icarus, 44, 172.
Weidenschilling S. J. (2011) Icarus, 214, 671.
Weidenschilling S. J. and Cuzzi J. N. (1993) in: Protostars andPlanets III, (edited by E. H. Levy and J. I. Lunine), pp. 1031–
1060.
Weidling R. et al. (2009) Astrophys. J., 696, 2036.
Weidling R. et al. (2012) Icarus, 218, 688.
Weissman P. R. et al. (2004) Structure and density of cometary
nuclei, pp. 337–357.
Wetherill G. W. (1992) Icarus, 100, 307.
Whipple F. L. (1972) in: From Plasma to Planet, (edited by
A. Elvius), p. 211.
Windmark F. et al. (2012a) Astron. Astrophys., 544, L16.
Windmark F. et al. (2012b) Astron. Astrophys., 540, A73.
Wurm G. et al. (2005) Icarus, 178, 253.
Wyatt M. C. (2008) Annu. Rev. Astron. Astrophys., 46, 339.
Yang C.-C. et al. (2009) Astrophys. J., 707, 1233.
Yang C.-C. et al. (2012) Astrophys. J., 748, 79.
Youdin A. and Johansen A. (2007) Astrophys. J., 662, 613.
Youdin A. N. and Goodman J. (2005) Astrophys. J., 620, 459.
Youdin A. N. and Lithwick Y. (2007) Icarus, 192, 588.
Youdin A. N. and Shu F. H. (2002) Astrophys. J., 580, 494.
Zsom A. and Dullemond C. P. (2008) Astron. Astrophys., 489, 931.
Zsom A. et al. (2010) Astron. Astrophys., 513, A57.
Zsom A. et al. (2011) Astron. Astrophys., 534, A73.
This 2-column preprint was prepared with the AAS LATEX macros v5.2.

24