Issue 
A&A
Volume 574, February 2015



Article Number  A83  
Number of page(s)  16  
Section  Planets and planetary systems  
DOI  https://doi.org/10.1051/00046361/201425222  
Published online  29 January 2015 
Erosion and the limits to planetesimal growth
^{1}
Leiden Observatory, Leiden University,
Niels Bohrweg 2,
2333 CA
Leiden,
The Netherlands
email:
krijt@strw.leidenuniv.nl
^{2} Astronomy Department, University of California, Berkeley, CA
94720, California
^{3}
Anton Pannekoek Institute, University of Amsterdam,
Science Park 904, 1098 XH
Amsterdam, The
Netherlands
Received: 27 October 2014
Accepted: 11 December 2014
Context. The coagulation of microscopic dust into planetesimals is the first step towards the formation of planets. The composition, size, and shape of the growing aggregates determine the efficiency of this early growth. In particular, it has been proposed that fluffy ice aggregates can grow very efficiently in protoplanetary disks, suffering less from the bouncing and radial drift barriers.
Aims. While the collision velocity between icy aggregates of similar size is thought to stay below the fragmentation threshold, they may nonetheless lose mass from collisions with much smaller projectiles. As a result, erosive collisions have the potential to terminate the growth of preplanetesimal bodies. We investigate the effect of these erosive collisions on the ability of porous ice aggregates to cross the radial drift barrier.
Methods. We develop a Monte Carlo code that calculates the evolution of the masses and porosities of growing aggregates, while resolving the entire mass distribution at all times. The aggregate’s porosity is treated independently of its mass, and is determined by collisional compaction, gas compaction, and eventually selfgravity compaction. We include erosive collisions and study the effect of the erosion threshold velocity on aggregate growth.
Results. For erosion threshold velocities of 20−40 m s^{1}, highvelocity collisions with small projectiles prevent the largest aggregates from growing when they start to drift. In these cases, our local simulations result in a steadystate distribution, with most of the dust mass in particles with Stokes numbers close to unity. Only for the highest erosion threshold considered (60 m s^{1}) do porous aggregates manage to cross the radial drift barrier in the inner 10 AU of MMSNlike disks.
Conclusions. Erosive collisions are more effective in limiting the growth than fragmentary collisions between similarsize particles. Conceivably, erosion limits the growth before the radial drift barrier, although the robustness of this statement depends on uncertain material properties of icy aggregates. If erosion inhibits planetesimal formation through direct sticking, the sea of ~10^{9} g, highly porous particles appears suitable for triggering streaming instability.
Key words: protoplanetary disks / planets and satellites: formation / circumstellar matter / methods: numerical
© ESO, 2015
1. Introduction
Despite the apparent ease with which nature is forming planets, current models of planet and even planetesimal formation have problems growing large bodies within the typical gas disk lifetime of ~10^{6} years (Haisch et al. 2001). The process of planetesimal formation is a complex one, with many different processes acting on a variety of lengths and timescales (see Testi et al. 2014; and Johansen et al. 2014, for recent reviews).
The first step towards planetesimal formation is the coagulation of small dust aggregates that stick together through surface forces. As aggregates collide and stick to form larger aggregates, these aggregates have to overcome several hurdles on their way to becoming planetesimals. One important obstacle faced by a growing dust aggregate is the radial drift barrier (Whipple 1972; Weidenschilling 1977). When aggregates grow to a certain size (about a meter at 1 AU and a millimeter at 100 AU, assuming compact particles) they will decouple from the pressuresupported gas disk, and start to lose angular momentum to the gas around them. As a result, said particles will drift inward.
Even before radial drift becomes problematic, the coagulation of aggregates can be frustrated by catastrophic fragmentation or bouncing (Blum & Wurm 2008; Güttler et al. 2010; Zsom et al. 2010), which prevents colliding aggregates from gaining mass. These problems are alleviated somewhat by including velocity distributions between pairs of particles (Windmark et al. 2012; Garaud et al. 2013) in combination with mass transfer in highvelocity collisions (Wurm et al. 2005; Kothe et al. 2010), though these solutions require the presence of relatively compact targets.
Recently, it has been proposed that icy aggregates, if they can manage to stay very porous, suffer less from these barriers, and might be able to form planetesimals locally and on relatively short timescales (Okuzumi et al. 2012; Kataoka et al. 2013a). Very porous, or fluffy, aggregates are less likely to bounce (Wada et al. 2011; Seizinger & Kley 2013), and icy particles have much higher fragmentation threshold velocities than refractory ones (Dominik & Tielens 1997; Wada et al. 2013), but perhaps most surprising was the finding that porous aggregates can outgrow the radial drift barrier, by growing very rapidly as a result of their enhanced collisional cross section (Okuzumi et al. 2012). However, Okuzumi et al. (2012) assumed perfect sticking between colliding aggregates, neglecting possible massloss in aggregateaggregate collisions.
We study the effects of the existence of an erosive regime for icy aggregates, where collisions at low mass ratios will produce erosive fragments at velocities below a critical erosion threshold velocity (Schräpler & Blum 2011; Seizinger et al. 2013; Gundlach & Blum 2014). Our goal is to quantify how erosion influences the direct formation of planetesimals through coagulation. To this end, we develop a local Monte Carlo coagulation code, capable of simulating the verticallyintegrated dust population, tracing both the evolution of the mass and the porosity of the entire mass distribution selfconsistently. Section 2 describes the models we use for the protoplanetary disk and the dust aggregates. In Sect. 3, we present the numerical method, which is based on the work of Ormel & Spaans (2008). Then, we test our model against the results of Okuzumi et al. (2012) (Sect. 4.1.1), after which we expand the model to include compaction from gas pressure and selfgravity according to Kataoka et al. (2013a) (Sect. 4.1.2), and erosive collisions (Sect. 4.2). In Sect. 5, we compare the results to a simple semianalytical model, and describe which processes can limit coagulation in different parts of protoplanetary disks. Discussion of the results and implications takes place in Sect. 6, and conclusions are presented in Sect. 7.
2. Disk and dust models
The disk model and collisional compaction prescription are based on Okuzumi et al. (2012), to which we add noncollisional compaction processes (Sect. 2.4.2) and a model for erosive collisions (Sects. 2.3.2 and 2.3.3).
2.1. Disk structure
The disk model used in this work is based on the minimummass solar nebula (MMSN) of Hayashi (1981). The evolution of the gas surface density and temperature as a function of radial distance R from the Sunlike central star are given as The gas sound speed is given by (3)with k_{B} the Boltzmann constant and m_{g} = 3.9 × 10^{24} g the mean molecular weight. The Kepler frequency equals (4)Assuming an isothermal column, the gas density drops with increasing distance from the midplane z according to (5)with the relative vertical scale height of the gas h_{g}/R = 0.05(R/ 5 AU)^{1/4}. The turbulent viscosity is parametrized as following Shakura & Sunyaev (1973), and α is assumed to be constant in both the radial and the vertical direction. The eddie turnover time of the largest eddies equals t_{L} = Ω^{1}.
In our local model, the surface density of the dust is related to the gas surface density through Σ_{d}/ Σ_{g} = 10^{2}, but the vertical distribution of dust depends on its aerodynamic properties. The dust is described by a Gaussian, with the dust scale height h_{d} set by the stopping time t_{s} of the dust particle through (Youdin & Lithwick 2007) (6)Thus, settling becomes important when a dust particle reaches Ωt_{s} ~ α.
2.2. Dust properties
Initially, all dust particles are assumed to be spherical (sub)micronsize monomers. In time, these monomers coagulate through collisions, and aggregates of considerable mass can be formed. Any aggregate is described by two parameters: the mass m, and the filling factor φ. Since aggregates are made up of monomers the mass can be written as m = Nm_{0}, with N the number of monomers and m_{0} the monomer mass. Following Okuzumi et al. (2012), we define the internal density of an aggregate as ρ_{int} = m/V, with V = (4/3)πa^{3} the volume of the aggregate, and a its radius. An aggregate’s radius is defined as , with r_{k} the position of monomer k and r_{CM} the position of the aggregate’s center of mass (Mukai et al. 1992; Suyama et al. 2008; Okuzumi et al. 2009). By definition, monomers have an internal density of ρ_{int} = m_{0}/V_{0} = ρ_{0}, while aggregates can have ρ_{int} ≪ ρ_{0}. Since we are interested in region beyond the snowline, we focus here on monomers composed of mostly ice, and use a density of ρ_{0} = 1.4 g cm^{3}. For the monomer radius we use a_{0} = 0.1 μm. We define the filling factor as (7)as a measure for the internal density.
In the rest of this section, we describe the main ingredients for the simulations presented in Sect. 3. These are: the relative velocities between aggregates, the equations governing the evolution of ρ_{int} through mutual collisions as well as gas ram pressure and selfgravity, and models for the destructive processes of erosion and fragmentation.
2.2.1. Relative velocities
We take into account relative velocities arising from Brownian motion, turbulence, settling, radial drift and azimuthal drift (see Sect. 2.3.2 of Okuzumi et al. 2012). The relative contribution of the velocity components depends strongly on the size and aerodynamic properties of the dust grains in question. More specifically, the relative velocity is a function of the stopping times of the particles. Depending on the size of the particle, the stopping time is set either by Epstein or Stokes drag (8)where is the mean thermal velocity of the gas molecules, and λ_{mfp} = m_{g}/ (σ_{mol}ρ_{g}) is the gas molecule mean free path. Taking σ_{mol} = 2 × 10^{15} cm^{2}, we obtain λ_{mfp} = 120(R/ 5 AU)^{11/4} cm at the disk midplane. In Eq. (8), a = a_{0}(V/V_{0})^{1/3} refers to the dust particle radius, while A is the projected cross section of the particle averaged over all orientations, which can be obtained using the formulation of Okuzumi et al. (2009).
The above equation is accurate when the particle Reynolds number Re_{p} = 4av_{dg}/ (v_{th}λ_{mfp}) < 1, with v_{dg} the relative velocity between the gas and the dust particle. The Reynolds number can become large when aggregates grow very big or their velocity relative to the gas is very large. In general, the stopping time can be written as (9)In the Stokes regime the drag coefficient equals C_{D} = 24/Re_{p}, and the stopping time becomes independent of v_{dg}. However, for larger Reynolds number the stopping time becomes a function of the velocity relative to the gas. This regime is called the Newton drag regime. Since the relative velocity depends in turn on the stopping time, we have to iterate to find the corresponding stopping time. Following Weidenschilling (1977), we use (10)Figure 1 shows Stokes numbers (Ωt_{s}) for different particles in the midplane of a MMSN disk at 5 AU. Different lines show compact particles (red), porous aggregates with constant φ = 10^{4} (yellow), and aggregates with a constant fractal dimension of 2.5 (green). For the solid lines, all drag regimes (Epstein, Stokes and Newton) have been taken into account, while the dashed lines indicate the results using only Epstein and Stokes drag, i.e., assuming that Re_{p}< 1. Focussing on the D_{f} = 2.5 aggregates, we can clearly distinguish the different drag regimes. The smallest particles are in the Epstein regime, and switch to the Stokes regime around Ωt_{s} = 10^{3}. Then, at a mass of m/m_{0} ~ 10^{21}, the Reynolds number exceeds unity and we enter the second regime of Eq. (10). We note that this transition occurs before Ωt_{s} = 1. The most massive particles, m/m_{0}> 10^{26} are in the regime where C_{D} = 0.44. Compact particles on the other hand, reach Ωt_{s} = 1 while still in the Epstein drag regime.
Fig. 1 Particle Stokes numbers as a function of mass, in the midplane of an MMSN disk at 5 AU. Different lines show compact particles (red), porous aggregates with constant φ = 10^{4} (yellow), and aggregates with a constant fractal dimension of 2.5 (green). For the solid lines, all drag regimes (Epstein, Stokes and Newton) have been taken into account, while the dashed lines indicate the results using only Epstein and Stokes drag. Horizontal lines indicate Ωt_{s} = 1 (where drift is fastest) and Ωt_{s} = α = 10^{3} (where particles start to settle to the midplane). 

Open with DEXTER 
The turbulenceinduced relative velocity between two particles with stopping times t_{s,1} and t_{s,2} ≤ t_{s,1} has three regimes (Ormel & Cuzzi 2007) (11)where δv_{g} = α^{1/2}c_{s} is the mean random velocity of the largest turbulent eddies, and t_{η} = Re_{t}^{1/2}t_{L} is the turnover time of the smallest eddies. The turbulence Reynolds number is given by , with the molecular viscosity ν_{mol} = v_{th}λ_{mfp}/ 2. We will refer to the first two cases of Eq. (11) as the first and second turbulence regimes. Relative velocities between similar particles (similar in the sense that they have comparable stopping times) are very small^{1} in the first turbulence regime because of the (t_{s,1} − t_{s,2}) term, but considerably larger in the second regime.
Fig. 2 Relative velocities between compact (left) or very porous (right) particles with masses m_{i} and m_{j} at the midplane of an MMSNdisk at 5 AU with α = 10^{3}. The masses range from single monomers to aggregates containing 10^{32} monomers. The contours give relative velocities in m s^{1}, and the colors indicate the dominating source for the relative velocity: Brownian motion (BrM), turbulence (Trb), settling (Set), radial drift (Rad), or azimuthal drift (Azi). Epstein, Stokes, and Newton drag regimes have been taken into account. 

Open with DEXTER 
Figure 2 shows the midplane relative velocity in m s^{1} (contours), and its dominant source (color), for a range of combinations of masses m_{i} and m_{j}. The velocities have been calculated for the disk properties of Sect. 2.1, at 5 AU, and assuming a turbulence α = 10^{3}. The left plot corresponds to two compact particles (ρ_{int} = ρ_{0}), and the right plot to two very porous ones (ρ_{int} = 10^{4}ρ_{0}). The general picture is the same for all porosities: Brownian motion dominates the relative velocity at the smallest sizes, followed by turbulence for larger particles, and systematic drift for bodies that have Ωt_{s} ~ 1. However, the masses at which various transitions occur can vary by several orders of magnitude depending on the particle porosity. For this particular location and turbulence strength, there is no combination of particle masses whose relative velocity is dominated by differential settling.
2.3. Collisional outcomes
A collision between porous aggregates can have a number of outcomes, ranging from perfect sticking to catastrophic fragmentation. For silicates, Blum & Wurm (2008) and Güttler et al. (2010) offer reviews of the various outcomes as observed in laboratory experiments. For porous ices, experimental investigations are scarce, and we have to turn to numerical simulations when predicting the outcome (e.g., Dominik & Tielens 1997; Wada et al. 2007; Suyama et al. 2008; Wada et al. 2009).
In general, a collision can result in sticking, erosion, or fragmentation, depending on the relative velocity and the mass ratio R^{(m)} ≡ m_{i}/m_{j} ≤ 1 of the colliding bodies. Collisions between particles with comparable masses result in catastrophic fragmentation if they collide above the fragmentation velocity (Sect. 2.3.1). When colliding bodies have a mass ratio R^{(m)} ≪ 1, catastrophic fragmentation of the larger body is difficult, but the collision can result in erosion if the velocity is high enough. The transition from erosion to the fragmentation regime occurs at a mass ratio , specified in Sect. 2.3.3. In an erosive event, the larger body will lose mass. From Fig. 2 it is clear that the highest velocities are reached between particles with very different masses, and thus erosion might very well be a common collisional outcome. We discuss erosion in more detail in Sect. 2.3.2. We should note at this point that we do not consider bouncing collisions. For relatively compact silicate particles, bouncing is frequently observed in the laboratory (e.g., Güttler et al. 2010), and indeed can halt growth in protoplanetary disks (Zsom et al. 2010). However, in porous aggregates, the average coordination number (the number of contacts per monomer) is much lower than in compact ones. As a result, collision energy is more easily dissipated, and it is safe to neglect bouncing (Wada et al. 2011; Seizinger & Kley 2013).
2.3.1. Catastrophic fragmentation
For collisions between roughly equal icy aggregates (mass ratio R^{(m)} ≥ 1/64), Wada et al. (2013) find a critical fragmentation velocity of (12)The quantity E_{break} represents the energy needed to break a single mononermonomer contact (Dominik & Tielens 1997). Collisions below this critical velocity result in sticking, while collisions at or above v_{frag} result in fragmentation of the collision partners.
2.3.2. The case for erosion
The relative velocity between similarsized aggregates will generally not reach the fragmentation velocity (Eq. (12)) behind the snow line in a protoplanetary disk, especially not if the turbulence is weak. However, relative velocities between particles with very different masses can be much larger than velocities between similar particles, especially when radial and azimuthal drift are important (Fig. 2). In this paper, we study the effects of an erosive regime, where collisions at low mass ratios will produce erosive fragments at velocities below a critical erosion threshold velocity v_{eros} ≲ v_{frag}. Here, we briefly revisit numerical and experimental studies of erosion, before outlining the erosion model used in this work. The process of erosion can be described by two main quantities: the erosion threshold velocity, v_{eros}, above which erosion takes place, and the (normalized) erosion efficiency, ϵ_{eros}, that indicates how much mass is eroded in units of projectile mass.
For silicate particles, Güttler et al. (2010) summarize a number of experimental investigations and describe a threshold velocity of a few m s^{1}, and an erosion efficiency that increases roughly linearly with collision velocity. Similar trends were observed by Schräpler & Blum (2011), who found an erosion threshold velocity of a few m s^{1} using micronsize silicate projectiles. We note that the threshold velocity is comparable to the monomer sticking velocity of micronsize silicate particles (Poppe et al. 2000). In the experiments of Schräpler & Blum (2011), the erosion efficiency also increased with impact velocity, reaching ~10 for the highest velocity of 60 m s^{1} (their Fig. 5). Seizinger et al. (2013) used molecular dynamics simulations, based on a new viscoelastic model (Krijt et al. 2013), to reproduce the experimental results. In addition, Seizinger et al. (2013) studied the variation on the threshold velocity and erosion efficiency with projectile mass, showing a trend of decreasing erosion threshold with decreasing mass ratio (e.g., their Fig. 11). For monomer projectiles, the threshold velocity equals the monomermonomer sticking velocity , after which it increased linearly with velocity to eventually flatten off around 10 m s^{1}. This flattening off indicates the onset of catastrophic fragmentation, and occurs at a mass ratio of ~10^{2}.
For ice particles, Gundlach & Blum (2014) present recent experimental results on the sticking and erosion threshold of (sub)micronsize particles. For a projectile distribution between 0.2−6 μm (with a mean value of 1.5 μm) impinging an icy target with a filling factor φ ≃ 0.5, an erosion threshold of 15.3 m s^{1} was found. These results confirm the increased stickiness of ice compared to silicate particles, and indicate v_{eros} could indeed be very high for (monodisperse) 0.1μm monomers, possibly even >60 m s^{1}. However, the aggregates acting as targets in the simulations presented here have a much higher porosity (φ ~ 10^{3}), and the lower coordination number is expected to reduce the erosion threshold (Dominik & Tielens 1997). Lastly, while Gundlach & Blum (2014) used a distribution of grain sizes, numerical investigations (e.g., Seizinger et al. 2013; Wada et al. 2013), for computational reasons, often employ a monodisperse monomer distribution, making a direct comparison difficult. For a single grain size, the size significantly influences the strength of the aggregates, with larger monomers leading to weaker aggregates (Eq. (12)). Little is known about the expected grain sizes in the icy regions of protoplanetary disks, let alone their size distribution, or about the effect a monomer size distribution has on the strength and collisional behavior of porous aggregates.
For these reasons, we believe that the existence of an erosive regime for icy aggregates is plausible. However, at present the data are unfortunately ambiguous with other simulations indicating the opposite trend: that the massloss in low mass ratio collisions is relatively small. Using molecular dynamics Nbody simulations Wada et al. (2013) find that the threshold velocity (where fragmentary collisions become more numerous than sticky collisions) increases for smaller mass ratios, suggesting that only similarsize particles colliding at v_{frag} fragment efficiently. This trend of an increased erosion threshold for smaller size ratios is corroborated by recent simulations by Tanaka et al. (in prep.). This would imply that for monodisperse submicron grains, both threshold velocities might not be reached (cf. Eq. (12) and Fig. 2). In this paper, we bury the uncertainty of the erosion threshold velocity in the parameter v_{eros}, which we vary to investigate the implications of effective versus ineffective erosion.
2.3.3. Erosion model
Erosive collisions only occur below a mass ratio , and their outcome is parametrized in terms of a (velocitydependent) erosion efficiency. In accordance with Güttler et al. (2010) and Seizinger et al. (2013) we will use . For smaller mass ratios, we will assume a constant value for v_{eros}, that does not depend on mass ratio or projectile/target porosity. We vary v_{eros} between 20 and 60 m s^{1}, corresponding to (1/4)v_{frag} and (3/4)v_{frag} for 0.1 μm monomers (Eq. (12)). In the erosive regime, the normalized erosion efficiency can be written as (13)with c_{1} ~ 1 (Güttler et al. 2010; Seizinger et al. 2013). While in supersonic cratering collisions γ = 16/9 (Tielens et al. 1994), the velocities encountered in this work are not that high and at most comparable to the sound speed in porous aggregates (Paszun & Dominik 2008). Hence, we will use γ = 1, in agreement with both numerical and experimental work in the appropriate velocity range (Güttler et al. 2010; Schräpler & Blum 2011; Seizinger et al. 2013).
Lastly, we need a prescription for the filling factors after an erosive collision. We assume that i) the filling factor of the target remains unchanged; and ii) the filling factor of the fragments is found by assuming they have the same fractal dimension as the target, where the target’s fractal dimension D_{f} is estimated as (14)The assumptions of the erosion model employed in this work are discussed further in Sect. 6.
2.4. Aggregate compaction
An aggregate’s porosity can be altered through collisions, or through noncollisional mechanisms. In this section, we first describe how porosity can increase and decrease as the result of sticking collisions. Then, we discuss gas and selfgravity compaction.
2.4.1. Collisional compaction
When two particles i and j collide at a relative velocity v_{rel} that is below the thresholds for fragmentation or erosion, the particles stick, and form a new aggregate with mass m_{i} + m_{j}. The internal density of the new particle depends on how the impact energy compares to the energy needed for restructuring. When the impact energy is not enough to cause significant restructuring, particles grow by hitandstick collisions, and very fractal aggregates can be formed (Kempf et al. 1999). When the impact energy is much larger, significant restructuring can take place, reducing the internal density of the dust aggregates. In this work, we will make use of the model presented in Suyama et al. (2012) and Okuzumi et al. (2012). Specifically, we use Eq. (15) of Okuzumi et al. (2012) to calculate the volume of the a newlyformed aggregate, as a function of the masses and volumes of particles i and j, the impact velocity, and the rolling energy E_{roll}; the energy needed to roll two monomers over an angle of 90° (Dominik & Tielens 1997).
Gundlach et al. (2011) measured the rolling force between ice particles with radii of ~1.5 μm to be 1.8 × 10^{3} dyn, implying a rolling energy of 1.8 × 10^{7} erg. Assuming the rolling force is sizeindependent (Dominik & Tielens 1995), the rolling energy is then often extrapolated using E_{roll} ∝ a_{0}. Recently however, Krijt et al. (2014) showed that the rolling force scales with the size of the area of the monomers that is in direct contact, resulting in , and , leading to significantly smaller rolling energies when extrapolating down to monomer radii well below a micrometer. In this work, we use the scaling law of Krijt et al., resulting in a rolling energy of 4 × 10^{9} erg for 0.1μm radius ice particles. Physically, a lower rolling energy means less energy is needed to start restructuring of an aggregate. As a result, a lower rolling energy will lead to compacter aggregates.
2.4.2. Gas and selfgravity compaction
Aggregates can also be compressed by the ram pressure of the gas, or their own gravity, if they become very porous or massive. For low internal densities, Kataoka et al. (2013b) found that the external pressure a dust aggregate can just withstand equals (15)This pressure can then be compared to the pressure arising form the surrounding gas and from selfgravity (16)with G the gravitational constant, in order to see whether an aggregate will be compacted as a result of these noncollisional processes (Kataoka et al. 2013a). In this work, we will take these effects into account in a selfconsistent way, while calculating the collisional evolution of the dust distribution.
3. Monte Carlo approach
Numerical techniques for studying coagulation can be divided in two categories^{2}: integrodifferential methods (e.g., Weidenschilling 1980; Dullemond & Dominik 2005; Birnstiel et al. 2010), and Monte Carlo (MC) methods (Gillespie 1975; Ormel et al. 2007; Zsom & Dullemond 2008). Tracing particle porosity as well as mass becomes computationally expensive in the integrodifferential approach. A solution to this issue was presented by Okuzumi et al. (2012), who assumed the porosity distribution for a given mass bin was narrow, but could vary in time. Since we are interested in including erosive processes, this assumption is not expected to hold, and for this reason we opt for the MC method.
The approach to calculate the collisional evolution is based on the distribution method as described in Ormel & Spaans (2008). In this section we briefly revisit the method, focussing on what is new in this work.
We let f(x) be the (timedependent) particle distribution function, with x_{i} the unique parameters describing dust particle i, in our case mass and filling factor^{3}. For every pair of particles i and j, one can determine the collision rate as (17)with the surface area of the column^{4}, and K_{ij} the collision kernel, which in this case equals (18)where h_{d,i} is given by Eq. (6), and and σ_{ij} = π(a_{i} + a_{j})^{2} equals the collisional cross section (Okuzumi et al. 2012). This rate equation takes into account that particles with different properties inhabit different vertical scale heights, and is correct as long as the coagulation timescale is longer than the vertical settling/diffusion timescale. In this work, we approximate the integral over z by assuming the midplane relative velocity is a good indication for v_{rel} throughout the column. This allows us to solve the integral analytically and write (19)For the purpose of this paper, this approximation is sufficiently accurate, since most of the growth is expected to take place near the midplane.
Then, we can define the total collision rate for particle C_{i} = ∑ _{j>i}C_{ij}, and the total collision rate C_{tot} = ∑ _{i}C_{i}. With all these rates known, 3 random numbers are used to identify which particles collide, and the time Δt after which this collision occurs. The colliding particles are then removed from f, and the collision product is added. As a result, all collision rates C_{i} have to be adjusted, since the particle distribution f has changed. This cycle is then repeated.
The simple method has two main drawbacks. First, the time needed for updating the rates in between collisions scales with N^{2}, where N is the total number of particles. Second, this method describes 1 collision per cycle, which can become a problem whenever the mass distribution is broad.
3.1. Grouping method
Rather than following every particle individually, identical particles can be grouped together. In our approach, the dust distribution is described by N_{f} particle families. Within a single family, all particles have identical properties, in our case mass and internal density. In every family i, there are w_{i} particle groups, each containing 2^{zi} individual particles, where we call z_{i} the zoom factor. The total number of particles in a single family therefor equals g_{i} = w_{i}2^{zi}, and the total number of particles is N = ∑ _{i}g_{i}. Instead of 2 particles colliding per cycle, collisions now happen between groups of particles (see Ormel & Spaans 2008, for details about this method). Letting i refer to the group with the lower zoom factor, we obtain for the group collision rates (20)where the i = j case in Eq. (20) describe socalled ingroup collisions. Like before, we can define the total collision rate per family λ_{i} = ∑ _{j ≥ i}λ_{ij}, and the total collision rate λ_{tot} = ∑ _{i}λ_{i}, which can be used to determine which groups collide and when. This grouped approach has tremendous advantages, but there are also pitfalls, which we discuss in the following section.
3.2. Sequential collisions
Imagine the collision between a group of large bodies i with a group of much smaller bodies j, such that m_{i} ≫ m_{j}. Thus, a total of 2^{zi}iparticles will collide with 2^{zj}jparticles. Assuming z_{j} ≫ z_{i}, every iparticle in the group will collide with 2^{zj − zi}jparticles in a single sequence before the collision rates are updated and the next groups to collide are chosen. We are assuming that the collision rates and the relative velocity between i and j particles are constant during this sequence, but this is only true if the properties of particle i do not change significantly. For this reason we include the group splitting factor N_{ε}, that limits the number of collisions to 2^{zj − zi − Nε}.
We let δm_{i} be the change in the mass of the larger particle i, after a single collision with a jparticle. Assuming the changes are small, we can then extrapolate to find the total change after the full sequence of collisions (21)Now, by imposing that (Δm_{i}/m_{i}) ≤ f_{m}, we obtain (22)where the square brackets indicate that is truncated to integers ≥ 0, which has the effect of particles with mass ratios ≥ f_{m} always colliding 1on1. In the case of perfect sticking, obviously δm_{i} = m_{j}, and Eq. (22) reduces to Eq. (12) of Ormel & Spaans (2008). We write an equivalent expression for the filling factor of the bigger grain (23)where δφ_{i} denotes the change in φ after a single collision. The two limits are combined by writing (24)and ensure that neither the filling factor, nor the mass of the larger particle change by too much during a single MC cycle. We note that N_{ε} is not only a function of the masses and densities of both particles, but also of the relative velocity, since this influences δφ_{i} (and δm_{i}, when erosion is present). In this work, we will typically use f_{m} = f_{φ} = 0.1.
Imposing this limit has two consequences. First, since the group of iparticles can now only collide with part of the group of jparticles, this needs to be taken into account when the group collision rates are calculated, changing Eq. (20) into (25)Second, since it can occur that only part of a group collides, group numbers w_{i} can now become fractional. This is fine as long as w_{i} ≥ 1, ensuring that at least one full group collision can occur in the future (Ormel & Spaans 2008).
3.3. The distribution method
For a given number of family members g_{i}, we have some freedom in choosing z_{i}; either creating many groups with a few members (low z_{i}) or a few groups with many members (high z_{i}). This choice for the zoom factors is crucial because it determines how many groups of a certain mass exist, which is related to the numerical resolution in that part of the mass range. Two approaches for determining the zoom factors have been proposed by Ormel & Spaans (2008).
One approach is the socalled equal mass method, in which one strives to have groups of equal total mass. This method is essentially identical to the method of Zsom & Dullemond (2008). With this approach, the peak of the mass distribution is very well traced, but parts of the particle distribution that carry little mass are described by a few groups, resulting in larger uncertainties. The second option is the distribution method, where one strives to have an equal number of groups per mass decade, independent of the total mass present in that interval. The difference between the two methods is nicely illustrated in Fig. 4 of Ormel & Spaans (2008). Since we are interested in erosion, it is crucial to resolve the particle distribution over the entire mass range. It is for that reason that we adopt the distribution method.
In practice, this means that at certain times during the simulation, we calculate the total number of particles N_{10} in every mass decade. The optimal zoom number for families in that mass range then equals (26)where w^{∗} is the desired number of groups per mass decade. In this way, we construct a function z^{∗}(m), which gives the desired zoom number for a family with particle mass m. We then check every existing family: if a certain zoom number is too big, we magnify the group (z_{i} → z_{i} − 1, w_{i} → 2w_{i}) until z_{i} = z^{∗}(m_{i}). Similarly, if the zoom number is too small, we demagnify (z_{i} → z_{i} + 1, w_{i} → w_{i}/ 2). The (de)magnification process conserves particle number, but does force one to update the various collision rates. A more detailed description of (de)magnification is given by Ormel & Spaans (2008). In the rest of this work, we calculate and update the zoom factors after every 10^{2} collision cycles, whenever the peak or average mass has changed by > 5%, or when the maximum mass has changed by > 50%, which we found to ensure a smooth evolution of the zoom factors. We will use w^{∗} = 60 for the perfect sticking calculations, and w^{∗} = 40 for the ones including erosion.
3.4. Merging
Lastly, we have to address the merging of families. It can occur that demagnification results in a group number w_{i}< 1, which is not allowed. When this occurs, the family does not contain enough individual particles to adopt z_{i} = z^{∗}(m_{i}). At this point, the family is insignificant. As we are simulating a fixed volume and the total mass needs to be conserved, we merge the family with another, healthy (meaning w_{i}> 1) one. First, we find the family j that resembles family i the most. In order to do so, we find the family that gives the largest product (R^{(m)})(R^{(φ)})^{3}, where R^{(φ)} ≤ 1 is the ratio of the filling factors^{5}. Then, we merge the families into a new family k with properties (27)The new zoom and group numbers are chosen such that z_{k} = z^{∗}(m_{k}). Merging is necessary to suppress the total number of groups.
3.5. Noncollisional compaction
Noncollisional compaction is implemented as follows: whenever a new aggregate is created in a collision, we calculate its compressive strength using Eq. (15), and compare this to the external pressures from gas ram pressure and selfgravity, calculated with Eq. (16) (Kataoka et al. 2013a). If either one of the external pressures exceeds P_{c}, we compactify the dust grain (i.e., increase φ) until the aggregate can withstand the external pressures.
Fig. 3 a) Schematic of a collision between unequal particles with a mass ratio R^{(m)} = (m_{proj}/m_{target}) ≪ 1. b) Sticking occurs when v_{rel}<v_{eros}. The mass of the projectile is added to the target. c) Collisions above the erosion threshold velocity lead to erosion. The mass loss of the target is given by the erosion efficiency ϵ_{eros} and the mass of the projectile. 

Open with DEXTER 
3.6. Erosion
For every collision, we check if the conditions for erosion are met (i.e., v_{rel}>v_{eros} and ), and if so, we determine the erosion efficiency using Eq. (13). After a single erosive event, the mass that does not end up in the target body equals (1 + ϵ_{eros})m_{proj}, see Fig. 3. To limit the number of new families, we redistribute this mass over fragments with a mass of m_{frag} = m_{proj}/ 10.
4. Results
In this section we show the results of our simulations for different erosion recipes, compaction mechanisms, turbulence strengths, and disk locations. When discussing the particle distribution at a given time, we shall use a number of quantities. These are the average mass and porosity (28)which trace the properties of the average particle, and the peak mass and filling factor (29)which trace the properties of the massdominating particle. We will also use the maximum mass m_{max}, which is simply the mass of most massive particle.
4.1. Perfect sticking
4.1.1. Collisionalcompactiononly scenario
As a test for the MC approach, we attempt first to match the trends observed in Okuzumi et al. (2012), who assumed perfect sticking between the dust grains. We adopt a turbulence strength parameter of α = 10^{3}, and focus on a vertical column at 5 AU in a typical MMSN disk. At this point, we only include collisional compaction and omit erosion. To allow a direct comparison to the work of Okuzumi et al., we do not include the effects of Newton drag for particles with large Reynolds numbers in this simulation. In the rest of this work, Newton drag is always included selfconsistently.
Fig. 4 Evolution of the normalized particle mass distribution at 5 AU with α = 10^{3}, assuming perfect sticking and without compaction through gas and selfgravity. Only Epstein are Stokes drag are considered. Solid lines indicate averages over 4 MC runs with identical starting conditions, and the shaded areas represent a spread of 1σ. 

Open with DEXTER 
Figure 4 shows the evolution of the normalized mass distribution m^{2}f(m) as a function of time. Solid lines mark the average over 4 MC runs. Thanks to the distribution method described in Sect. 3.3, the sampling of the mass distribution is very good over the entire mass range: even at later times, when most of the mass is located in particles with masses of ~10^{15} g, the distribution of particles all the way down to 10^{9} g is resolved remarkably well, despite these particles only making up a very small fraction of the total mass.
When we compare our Fig. 4 to Fig. 7 of Okuzumi et al. (2012), it is clear that our local MC method yields very similar results. We recognize the familiar narrow mass peak when growth is governed by Brownian motion, followed by a broader distribution once turbulence kicks in. Once particles reach Ωt_{s} ~ 1 (m_{j} ≃ 10^{10} g in this case), systematic drift greatly increases their collision rate, and very rapid growth ensues. The slight difference in timescales is attributed to i) the slightly different value for the rolling energy, ii) our approximation of Eq. (18), and iii) our use Eq. (8) to calculate the stopping times, while Okuzumi et al. used to ensure a smooth transition between Epstein and Stokes drag (Okuzumi, priv. comm.).
Fig. 5 Like Fig. 4, but with compaction through gas and selfgravity and Newton drag for particles with Re_{p}> 1. 

Open with DEXTER 
Fig. 6 Evolution of the growth and radial drift timescale of the peak mass for the perfect sticking model at 5 AU with α = 10^{3}. The dotted line indicates (t_{drift}/ 30). Only collisional compaction has been taken into account. 

Open with DEXTER 
Fig. 7 Evolution of the growth and radial drift timescale of the peak mass for the perfect sticking model at 5 AU with α = 10^{3}. The dotted line indicates (t_{drift}/ 30). Compaction from gas and selfgravity, and Newton drag have been taken into account. 

Open with DEXTER 
While we take into account driftinduced relative velocities, the dust particles are bound to our simulated column and cannot move radially through the disk. To test the validity of this assumption, we compare the growth timescale of the peak mass, defined as (30)to the radial drift timescale at that mass (31)The radial drift velocity is given by (Weidenschilling 1977) (32)where v_{K} = RΩ is the Keplerian orbital velocity, and η can be written as (Nakagawa et al. 1986) (33)Figure 6 shows both the growth and radial drift timescales during the complete evolution of the peak mass. Initially, relative velocities are dominated by Brownian motion. Since this velocity drops with increasing particle mass, the growth timescale increases. Around a mass of 10^{9} g, turbulent velocities start to dominate the relative velocity, and the growth timescale stays approximately constant. Particles larger than 10^{3} g enter the second turbulent regime as t_{s}(m_{p}) >t_{η}. In this regime, velocities between similar particles are increased (see Eq. (11)), which leads to a decrease in the growth timescale. Since the growth timescale is always much smaller than the drift timescale, the aggregates in this simulation do indeed outgrow the radial drift barrier.
4.1.2. Including gas and selfgravity compaction
The next step is to include compaction by gas pressure and selfgravity, as described in Sect. 2.4.2. In addition, we now take into account Newton drag for particles with large Reynolds numbers. Figure 5 shows the results for the same disk parameters as before. The general shape of the evolution looks similar to Fig. 4 initially, but from the corresponding times it is clear that the growth is slower for the largest aggregates. The main reason for this is that the largest dust grains are compacted by the gas and selfgravity, resulting in a smaller collisional cross section. In addition, the aerodynamic properties are different, which affects the relative velocities.
The growth and drift timescales are plotted in Fig. 7. When we compare Figs. 6 and 7, we confirm that the growth close to the drift barrier is slower when using the full compaction recipe. For the largest particles, the growth timescale is increased by more than 2 orders of magnitude. In addition, including Newton drag has broadened the drift barrier somewhat. Nonetheless, the growth is still fast enough to prevent particles from drifting significant distances.
4.1.3. Evolution of internal densities
It is interesting to compare the evolution of the internal densities of the particles for the models with and without noncollisional compaction. In Fig. 8, the peak filling factor is plotted versus the peak mass for the simulations described so far. The symbols correspond to important points in the evolution of the aggregates: open circles are related to the stopping time of the aggregates, and closed symbols indicate the onset of various compaction mechanisms^{6}.
Initially, aggregates grow through hitandstick collisions, and evolve along a line of constant fractal dimension close to 2. In the collisionalcompactiononly scenario, particles reach a filling factor of ~10^{5} during hit and stick growth, before collisional compaction kicks in, after which φ stays almost constant. When Ωt_{s}(m_{p}) > 1, the internal density drops even further. The general picture, as well as the location of the various turnover points, is consistent with the top panel of Fig. 10 of Okuzumi et al. (2012). When noncollisional compaction is included, the filling factor, in general, is much higher at later times, and follows the boundaries that have been described by Kataoka et al. (2013a; e.g., their Fig. 3). For this particular combination of turbulence, rolling energy, and monomer size, compacting by gas ram pressure actually occurs before the first collisional compaction event takes place^{7}. Significant settling occurs when Ωt_{s}>α, which corresponds to m ~ 10^{3} g for compact particles (see Fig. 1). From Fig. 8 however, we see that porous particles only begin to settle when their masses reach ~10^{4}−10^{5} g. Lastly, aggregates with masses above ~10^{10} g are compacted by selfgravity, causing the filling factor for the largest bodies to be several orders of magnitude higher. In the remainder of this work, we include both collisional and noncollisional compaction mechanisms, and Epstein, Stokes, and Newton drag selfconsistently.
Fig. 8 Evolution of the internal structure of the massdominating particles, for the perfect sticking models at 5 AU, for the models with and without noncollisional compaction mechanisms. Aggregates start out as monomers in the top left corner, and grow towards larger sizes and porosities. Lines show individual simulations. Open symbols correspond to points where the mass dominating particles reach a = λ_{mfp} (°); t_{s} = t_{η} (◇); Ωt_{s} = α (▿); and Ωt_{s} = 1 (□). Filled symbols show peak mass and filling factor at the times of first: collisional compaction (⋆); gaspressure compaction (◆); and selfgravity compaction (•). 

Open with DEXTER 
4.2. Erosion
With this framework in place, the final step is to include the erosion model of Sect. 2.3.3 in the simulations, and calculate the evolution of the particle distribution selfconsistently. Figure 9 shows the mass distribution at various times for v_{eros} = 20 m s^{1}. Initially, the evolution proceeds just like in 5, but as the largest aggregates approach Ωt_{s} = 1, their velocity relative to smaller particles is high enough for erosion, and their growth stalls. As a direct consequence of the erosion, the amount of small particles increases, and after ~4000 yr a steadystate is reached, with a significant amount of mass residing in particles smaller than a few grams.
Fig. 9 Evolution of the normalized particle mass distribution at 5 AU with α = 10^{3}, assuming v_{eros} = 20 m s^{1}. The full compaction model is used. 

Open with DEXTER 
To investigate how erosion halts the growth of the largest bodies, it is instructive to plot socalled projectile mass distributions (Okuzumi et al. 2009). For a certain particle mass m_{t}, these distributions show the contribution to the growth of that particle as a function of projectile mass m ≤ m_{t}. An example of such a plot is shown in Fig. 9 of Okuzumi et al. (2012), where the distribution function is plotted at various times for m_{t} = m_{p}. For our distribution plots, we make two important changes: first, since we are interested in the growth of the largest bodies, we plot projectile distributions for m_{t} = m_{max}, with m_{max} the largest mass in the simulation at a given time. Second, to illustrate the effect of erosive collisions, we calculate the mass loss for every erosive collision, taking into account the correct erosion efficiency^{8}. As a result, the sign of the distribution function can be both positive and negative. Figure 10 shows the distribution for one of the simulations of Fig. 9 (colors correspond to the same times). When we examine the rightmost projectile mass distribution, corresponding to a time t = 10^{4} yr, it is immediately clear how erosion affects the evolution of particles with Ωt_{s} ~ 1. While these aggregates grow by collisions with similarsized bodies, they lose mass by colliding with particles that have a mass below 10^{2}m_{t}. This could have been predicted by looking at Fig. 2, from which it is clear that the highest velocities are attained between particles with mass ratios well below unity. The importance of this erosion however, depends on the current mass distribution, and can only be tested through dedicated simulations like the ones presented here. Since the area under the negative part of the projectile distribution outweighs the positive part, the erosion is so effective that it stops the growth of the largest bodies, resulting in the behavior seen in Fig. 9.
Fig. 10 Projectile distribution mass functions for simulation E1, constructed for the maximum masses (m_{t} =•) at various times. Colors and times correspond to Fig. 9. For each distribution, the stopping time of the m_{t}particle is given, and the weights of the total positive and negative area are plotted. The distributions have been normalized in such a way, that the absolute sum of the contributions equals 1. 

Open with DEXTER 
We define a parameter ζ using the positive and negative areas under the projectile mass distributions (34)with ∑ C_{+} and ∑ C_{−} the sums of the positive and negative part of the projectile mass distribution respectively. The parameter ζ ranges from 1 (no erosion) to –1 (only erosion), and equals 0 when there is a balance between growth and erosion.
The top panel of Fig. 11 shows the evolution of ζ for the most massive particle during one of the simulations of Fig. 9, plotted as a function of Ωt_{s} of the maximum mass. Early on, there is no erosion present and ζ = 1, but as the largest bodies grow towards Ωt_{s} = 1, erosion increases and ζ drops. When 0 <ζ< 1, the massive particles still grow faster then they are eroded, but the erosion can be significant in that it results in the creation of more small particles, thus increasing its destructive effect. When ζ< 0, erosion dominates over growth and the most massive particles are loosing considerable mass. This causes the curve in Fig. 11 to turn around. As bodies shrink, there is less erosion and ζ increases again. A quasi steadystate is reached with ζ just below unity and Ωt_{s}(m_{max}) ~ 0.6. The reason ζ ≠ 0 during the steady state, is that it is not the same particle that is the most massive at all times. Instead, particles take turn at being the most massive body. Since the largest particles are stuck at a mass and size for which drift is fastest, they will move radially towards the central star. For this combination of parameters, we conclude that growth beyond the drift barrier is impeded by erosion.
Fig. 11 Evolution of ζ(m_{max}) for erosive simulations with v_{eros} = 20,40,60 m s^{1} as a function of Ωt_{s}(m_{max}), showing the impact of erosion on the ability of the largest bodies to grow. In the upper two panels, the steadystate is indicated by the ⊚symbol. 

Open with DEXTER 
The other panels of Fig. 11 show similar plots but for different erosion threshold velocities. For v_{eros} = 40 m s^{1} (middle panel), erosion is less efficient and the largest bodies grow to Ωt_{s} ≃ 10 before they start to lose mass rapidly. The reason particles can grow larger is twofold. First, the threshold velocity itself is somewhat higher, causing erosion to start for higher masses. Second, since the erosion efficiency is proportional to (v_{rel}/v_{eros}), the highvelocity projectile are less efficient in excavating mass from the targets. Both effects together cause the largest mass in the steady state to be about a factor of 10 larger than in the top panel of Fig. 11. Finally, the bottom panel shows the results for v_{eros} = 60 m s^{1}. This is a special case, since now the erosion threshold velocity can only be reached around Ωt_{s} = 1, with radial drift, azimuthal drift, and turbulence contributing (see Fig. 2). Indeed, erosion is strongest around Ωt_{s} = 1, but it is inefficient and ζ never drops below 0. When Ωt_{s}> 20, erosion reappears, as a result of smaller particles drifting into the larger bodies, but since ζ ~ 1, bodies can continue to grow relatively unaffected.
Fig. 12 Masses and filling factors of all unique families at different times, for simulation E1. One dot corresponds to one family, and does not provide information about the total mass or number of members in that family. 

Open with DEXTER 
4.2.1. Variation in porosity
One of the biggest advantages of the MC method is that aggregate mass and porosity are treated truly independently. In other words, aggregates of identical mass can have a very different porosity. However, the collision model used in this work immediately implies that the spread in porosities (for a given particle mass) will be narrow, when sticking collisions dominate the evolution. For example, the collision model, at the moment, does not include an impactparameter dependence in collisions, or a random component in the relative velocity. As a result, collisions between particles with certain properties always occur at the same relative velocity, and always result in the same collision product(s). Moreover, when gas compaction (or selfgravity compaction) limits the porosity of an aggregate, bodies will evolve along P_{c} = P_{gas} (or P_{c} = P_{grav}), according to Eqs. (15) and (16). As a result, massporosity relations as shown in Fig. 8 accurately represent the internal structure of the majority of aggregates.
This picture changes when erosion starts to play a role. Figure 12 shows the evolution of the properties of each family in one of the E1 simulations. We note that each dot corresponds to a single family, and that the total masses and number of family members can vary significantly between families. Nonetheless, Fig. 12 gives a good indication of the spread in porosity. For the reasons described above, the spread in porosity is very small during the first 3000 years of the evolution. After 3400 years, the first erosive collisions have occurred, and created a population of fragments with a fractal dimension set by the parent body. At this point, the porosity distribution becomes bimodal, and the assumption of a single porosity parameter – which only depends on aggregate mass – is untenable. Later, after ~6000 years, the original population of aggregates, whose porosity was set by their growth history, has disappeared. A steadystate is reached in which the internal structure of the fragments is dominated by the porosity of the particles that act as targets for erosion, i.e., the large bodies with Ωt_{s} ~ 1.
5. Semianalytical model
The evolution of the massdominating particles can be captured in a simple semianalytical model. Assuming the entire dust mass is located in particles of identical mass m_{p}, the growth rate can be written as (Okuzumi et al. 2012) (35)The collisional cross section depends directly on the particle porosity, and the relative velocity and dust scale height depend on φ through the particle stopping time. As a simple model for the aggregate’s internal structure, we assume the aggregates initially grow with a constant fractal dimension of ~2, until the kinetic energy in samesized collisions exceeds E_{roll}. After that, the internal structure can be calculated through Eq. (31) of Okuzumi et al. (2012), but in practice is always dominated by the gas/selfgravity compression of Kataoka et al. (2013a), see Sect. 2.4.2.
This approach, similar to Kataoka et al. (2014, Sect. 5.3), is valid when particles grow primarily through collisions with similarsized particles. This is valid in most regimes, but not true in the first turbulence regime. Here, relative velocities between identical particles are suppressed, and aggregates grow by collecting smaller particles. However, it can be shown that in this regime the growth timescale is approximately constant (Okuzumi et al. 2009). Hence, we will assume that t_{grow} is constant in the regime where turbulence dominates v_{rel}, and t_{s}<t_{η}.
At the same time, the radial drift of the particles is governed by (36)with the drift velocity a function of Ωt_{s}. Assuming a fixed dust to gas ratio of 10^{2} throughout the disk, we can solve Eqs. (35) and (36) to obtain the evolution of the vertically integrated peak mass. Catastrophic fragmentation is taken into account by setting (dm_{p}/ dt) = 0 when v_{turb}>v_{frag} for two particles of mass m_{p}. Figure 13 shows lines along which the dust evolves, starting from m = m_{0} at various locations in the disk. The left plot shows the results for compact growth (i.e., φ = 1 at all times), after 10^{6} yr. (For the compact case, we have temporarily set v_{frag} = 10 m s^{1}.) Initially, growing aggregates are not moving radially, resulting in vertical lines in Fig. 13. As the particles’ Stokes numbers increase, collision velocities and drift speeds increase. In the inner regions of the disk, the maximum size is limited by fragmentation through samesized collisions. Particles cannot grow larger than ~cm, and will inevitably drift inwards. In the intermediate region, from 20−100 AU, the fragmentation velocity is not reached. Here, the maximum size is set by radial drift. In the outermost disk (beyond 10^{2} AU), growth is very slow because of the low dust densities, and 10^{6} yr is not enough to reach the size necessary to start drifting. The general behavior is identical to what is observed in full compact coagulation models (see Fig. 3 in Testi et al. 2014).
Fig. 13 Evolution of m_{p}(t) and R(t) for dust coagulation as obtained from the semianalytical model (Eqs. (35) and (36)), for an MMSN disk and α = 10^{3}. Lines indicate different starting conditions R(t = 0), and are evolved for 10^{6} yrs. Left: compact growth: φ = 1 at all times, and v_{frag} = 10 m s^{1}. Right: porous growth: the internal structure of the aggregates is set by hit and stick growth, followed by collisional compaction or gas and selfgravity compaction. Gray lines have no erosion, while black lines show the results for v_{eros} = 40 m s^{1}. Colored lines and ⊚symbols indicate the evolution and steady state peak mass obtained through local MC simulations (Sect. 4.2). 

Open with DEXTER 
The gray lines in the righthand panel of Fig. 13 show the results of the semianalytical model for porous growth, where φ is set by collisional, gas pressure, and selfgravity compaction, assuming perfect sticking. Since we are assuming the m_{p} particles carry the total dust mass, we do not have any information about the massdistribution of smaller particles. Nonetheless, we can mimic the effect of effective erosion, by setting (dm_{p}/ dt) = 0 when the relative velocity between the mass dominating particle and small projectiles (taken to be monomers) exceeds v_{eros}. The black solid lines in the right panel of Fig. 13 show the results for v_{eros} = 40 m s^{1}, while the red lines indicate results for the peak mass of the full MC models for the same erosion threshold velocity (although that the maximum mass reached in these models can be a factor of ~10 larger). We have also included a full model run at 200 AU, which we evolved for 10^{6} yrs. The results of the semianalytical model agree with the simulations of the previous section remarkably well.
6. Discussion
From the maximum sizes fluffy aggregates can reach at a given location, we identify three zones in the protoplanetary disk:

3−10 AU: assuming perfect sticking, the combination of Stokes drag and enhanced collisional cross sections allows the porous aggregates in the inner disk to outgrow the radial drift barrier, and reach planetesimal sizes without experiencing significant drift. However, when erosion is efficient, mass loss in erosive collisions stalls the growth around Ωt_{s} ~ 1, preventing the porous aggregates from crossing the radial drift barrier (Fig. 11).

10−100 AU: at intermediate radii growth timescales increase and radial drift takes over, even before aggregates reach sizes and stopping times that allow erosive collisions to take place.

>100 AU: in the outer disk, the disk lifetime is not long enough for particles to grow to sizes where significant drift occurs. In the porous growth scenario, aggregates this far out are in the hitandstick regime, and their surfacetomass ratio does not change when they gain mass. As a result, hardly any drift is visible. In the compact case, an increase in mass automatically results in a decrease in the surfacetomass ratio, and the onset of radial drift is already visible for very low particle masses.
For erosion to start, the collision velocity between target and projectile needs to exceed v_{eros}. In the limit where the projectiles are monomers that couple to the gas extremely well, this collision velocity equals the relative velocity of the large bodies with respect to the gas. When the largest particle has Ωt_{s} ≫ 1, it moves on a Keplerian orbit, and v_{dg} ≃ ηv_{K}, while bodies with Ωt_{s} = 1 have a slightly larger velocity with respect to the gas (Weidenschilling 1977). For the disk model employed in this work (Sect. 2.1), the quantity ηv_{K} does not depend on R, and thus the maximum drift speed is constant though out the disk. It is clear then from Eq. (33) that growing aggregates in colder disks (lower c_{s}), or disks with a (locally) shallower gas density profile might suffer less from erosion.
It is clear from Fig. 11 that the size of v_{eros} is an essential parameter: its value, together with ηv_{K}, determines whether growth beyond Ωt_{s} = 1 is possible or not. Unfortunately, the value of v_{eros}, or even its relation to v_{frag}, is not accurately known for the large and highlyporous icy bodies in question (Sect. 2.3.3). Numerical investigations, showing conflicting trends for erosion efficiency with mass ratio, often employ monodisperse grain sizes (e.g., Seizinger et al. 2013; Wada et al. 2013), and the threshold velocities depend almost linearly on the grain radius (Eq. (12)), a parameter which itself is not well constrained. At the same time, the only available experimental work on erosion for ices used a distribution of grain sizes (Gundlach & Blum 2014). In addition, both numerical and experimental studies are restricted to sizes ≲ mm and porosities ≳ 10^{1}, and cover a soberingly small portion of the parameter space encountered in this work (e.g., Fig. 8). Future studies, numerical as well as experimental, are encouraged to elucidate these matters, and constrain the threshold for erosion and its dependence on target/projectile sizes and porosity. Finally, we assume that material that is eroded locally is removed from the target. In reality, the fate of the fragments will be determined by the local gas flow and the velocity with which they are ejected. For very porous targets, the gas flow through and around the surface of the target might result in these fragments being reaccreted (Wurm et al. 2004). If efficient, this reaccretion might be a way to alleviate the destructive influence of erosive collisions. On the other hand, the flow through a body is likely to be insignificant, unless it is extremely porous (Sekiya & Takeda 2005).
So far we have assumed that while erosion can play an important role, catastrophic fragmentation does not occur. The maximum velocity between samesized bodies is reached for Ωt_{s} = 1, and equals ~(3/2)α^{1/2}c_{s} (Eq. (11)). Since the sound speed diminishes for increasing radii, this velocity is highest in the inner disk. For typical turbulence strengths (α ≲ 10^{3}) and small icy monomers, this velocity will not exceed the fragmentation threshold velocity (Eq. (12)), and, especially in the outer disk, fragmentation of icy bodies through catastrophic fragmentation is very unlikely. However, if all collisions result in sticking, small particles (≲ 100 μm) are removed from the protoplanetary nebula very rapidly, contradicting observational constraints (Dullemond & Dominik 2005; Dominik & Dullemond 2008). Driftinduced erosion can alleviate these issues, since the maximum drift velocity is high throughout the entire disk.
In this work, we have assumed collisions below the fragmentation threshold to result in perfect sticking, i.e., the mass of the resulting aggregate equals the sum of both colliding masses. However, even for collisions below the fragmentation threshold velocity, a significant amount of mass may be ejected during a collision, especially if the collision occurs at a large impact parameter (Paszun & Dominik 2009; Wada et al. 2013). An advantage of a MC model approach like the one presented here, is that it is relatively straightforward to include an additional random number to determine, for example, the impact parameter. The difficulty lies in obtaining a collision model that describes the collisional outcome as a function of this parameter. A good start would be the work of Wada et al. (2013), who show the growth efficiency as a function of impact parameter (Fig. 4). Basically, headon collisions promote growth, while collisions with a large impact parameter result in little mass gain. Unfortunately, much less is known about the porosities of the resulting aggregates.
At the heart of the model of Sect. 3 lies the assumption that an aggregate is adequately described by two quantities: its mass and (average) porosity. While this represents a considerable improvement on the compact coagulation assumption, a single average porosity does not allow a complex internal structure of the aggregates. For small grains, the accuracy of this assumption will depend on their collisional history. For example, one can imagine a porous aggregate with a denser outer shell being formed if the aggregate is compacted through many collisions with small mass ratios (Meisner et al. 2012). Such a compact rim will hardly alter the aggregate’s average porosity, but can influence its sticking and erosion behavior (Schräpler & Blum 2011). Likewise, gas and selfgravity compaction need not result in a homogenous internal structure. With instruments such as CONCERT on board ESA’s Rosetta and Philae capable of probing the internal structure of large solar system objects, studies focussing on the internal structure of the larger bodies, as determined by its growth and compaction history would be very interesting. The MC method developed in this paper would be suitable for such studies, since adding parameters describing the aggregates is relatively straightforward.
6.1. Future work and implications
6.1.1. Pebble accretion
A novel idea in the field of planet formation is the process of pebble accretion, where protoplanets grow very efficiently by accreting small pebbles (Ormel & Klahr 2010; Lambrechts & Johansen 2012, 2014; Kretke & Levison 2014). These models rely on the radial influx of particles drifting in from the outer disk. As in the compact case, porous growth leads to the creation of rapidly drifting bodies in the region between 10 and 10^{2} AU (Fig. 13). While the Stokes numbers of these particles are similar to the drifting pebbles in the compact case, their masses, sizes, and porosities can differ by many orders of magnitude (see also Fig. 14). In addition, the drag regime that the drifting bodies experience differs from the compact case (Fig. 1). Future studies are needed to address the effect of these factors on the efficiency of pebble accretion.
Fig. 14 Same as the righthand plot of Fig. 13, but with aggregate size on the vertical axis. 

Open with DEXTER 
6.1.2. Streaming instability
While – depending on the critical erosion velocity – rapid coagulation into masses as large as planetesimals might be prevented by erosive collisions, the conditions created by this process might be favorable for triggering planetesimal formation by streaming instability (Youdin & Goodman 2005; Johansen et al. 2007; Bai & Stone 2010a,b). To trigger streaming instability, the majority of mass needs to reside in particles with high Stokes numbers; the midplane dust to gas ratio has to be close to unity; and the local vertically integrated dusttogas ratio needs to exceed ~0.03 (Dra¸żkowska & Dullemond 2014). The first two conditions can be studied with simulations like the ones presented in this work. For example, for the steadystate distribution reached for v_{eros} = 40 m s^{1} at 5 AU for α = 10^{3}, approximately 50% of the dust mass resides in particles with Ωt_{s}> 10^{2}, and the midplane dusttogas ratio is ~10^{1}. For weaker turbulence, the midplane dusttogas ratio will be increased further, since h_{d} ~ α^{1/2} (Eq. (6)). Because our simulations are local, the vertically integrated dusttogas ratio stays constant at 10^{2}. To fulfill the third condition, the dusttogas ratio either has to be larger from the beginning, or must increase by material drifting in from the outer disk. To study this, a global model is required, that calculates the evolution of the dust surface density in the presence of radial drift and erosion. In conclusion, driftinduced erosion appears to be a robust way of concentrating mass around Ωt_{s} ~ 1, and is expected to create conditions favorable for streaming instability.
6.1.3. The breakthrough case
For compact silicate bodies in the inner disk, bouncing and fragmentation are very effective in stopping growth at mmcm sizes. The breakthrough scenario, in which a small number of lucky particles still manages to gain mass, might render further growth possible (Windmark et al. 2012; Garaud et al. 2013). The total mass fraction of these lucky particles can be extremely small, making this a challenging process to model for both differential and MC methods (Dra¸żkowska et al. 2014). The distribution method used in this work, as outlined in Sect. 3.3, is capable of resolving the entire mass distribution, including parts that contribute very little to the total dust mass, and appears to be suitable for studying the breakthrough case.
6.1.4. Opacities of porous grains
The optical properties of dust distributions resulting from porous growth are very different from populations containing exclusively solid particles. Not only are the mass distributions themselves different (e.g., Fig. 13), but also the scattering and absorption opacities of the individual grains are affected significantly by porosity (Kataoka et al. 2014; Cuzzi et al. 2014). For simple dust mass distributions, the effect of grain porosity on the appearance of protoplanetary disks has been investigated by Kirchschlager & Wolf (2014). Combining selfconsistent coagulation models  including erosion and fragmentation – with porositydependent dust opacities will reveal the full impact porous growth has on the appearance of protoplanetary disks.
7. Conclusions
Porous growth is very different from compact growth (Fig. 13). For example, porous particles have larger collisional cross sections than compact particles of the same mass. More importantly, the aerodynamical properties of porous aggregates can differ greatly from those of compact particles (Fig. 1), causing differences in relative velocities (Fig. 2), vertical settling, and radial drift.
We have modeled the coagulation of porous icy particles in the outer parts of protoplanetary disks, tracing the evolution of the mass and filling factor of the individual aggregates in time. We consider compaction through collisions, gas pressure, and selfgravity (Fig. 8), and include a physical model for erosive collisions (Sects. 2.3.2 and 2.3.3). The main findings of this work are:

1.
Porous icy aggregates can outgrow the radial drift barrier in the inner ~10 AU, despite increased growth timescales resulting from gas and selfgravity compaction, if the perfect sticking assumption holds (Figs. 7 and 13). This is in agreement with Okuzumi et al. (2012) and Kataoka et al. (2013b).

2.
While the maximum collision velocity between similar particles (~ α^{1/2}c_{s}) typically does not exceed the critical fragmentation threshold velocity for icy bodies, the velocity between drifting aggregates (with Ωt_{s} ≥ 1) and smaller bodies is much larger (~ ηv_{K}), and can exceed the critical threshold velocity for erosion (Fig. 2).

3.
In these cases, we find that the mass loss through erosive collisions can balance the growth through samesize collisions, halting the growth of the largest bodies (Figs. 10 and 11). In our local simulations, this results in a steadystate where the largest bodies have Ωt_{s} ~ 1, and the porosity of the small fragment distribution is dominated by the fact that all fragments have at some point been part of these large eroded particles (Figs. 9 and 12). Only for the highest erosion threshold velocity we considered (v_{eros} = 60 m s^{1}) do the aggregates with Ωt_{s} ~ 1 manage to gain mass and grow through the drift barrier.

4.
A simple semianalytical model (Sect. 5) accurately describes the growth and drift behavior of the massdominating bodies. While no information is obtained about the dust mass distribution, such an approach is very useful for investigating how the size of the largest bodies depends on disk parameters such as the total disk mass, turbulence strength, and dusttogas ratio, or aggregate properties such as monomer size and erosion/fragmentation threshold velocities.
According to Eq. (11), v_{turb} = 0 for aggregates with identical stopping times in the first turbulence regime. In reality, the dispersion in the aggregate’s masstoarea ratio will give rise to a small relative velocity. We treat this dispersion in the same way as Okuzumi et al. (2012), by taking into account the standard deviation in the masstoarea ratio of a porous aggregate (Okuzumi et al. 2011). The size of this standard deviation, normalized by the mean masstoarea ratio, is parametrized as ε, which we take to equal 0.1, following Okuzumi et al. (2011).
See Dra¸żkowska et al. (2014) for a comparison between the two methods in the breakthrough growth case.
In fact, the gas compaction starts when the aggregates are still in the Epstein drag regime. Equation (15) is determined by static compression of porous aggregates, and Eq. (16) assumes the external pressure can be treated as continuous. However, if the collision frequency of gas molecules with individual monomers of the aggregate is low compared to the frequency at which monomermonomer contacts oscillate and dissipate energy, this approach might not be accurate. Future work is encouraged to investigate the effect of collisions between the aggregate and gas molecules in this regime.
Acknowledgments
Dust studies at Leiden Observatory are supported through the Spinoza Premie of the Dutch science agency, NWO. S.K. would like to thank C.P. Dullemond, J. Dra¸żkowska and T. Birnstiel for useful discussions, and S. Okuzumi for comments and helping with the comparison to his work. The authors thank B. Gundlach and J. Blum for sharing a version of their manuscript. C.W.O. acknowledges support for this work by NASA through Hubble Fellowship grant No. HSTHF51294.01A 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 526555.
References
 Bai, X.N., & Stone, J. M. 2010a, ApJ, 722, 1437 [NASA ADS] [CrossRef] [Google Scholar]
 Bai, X.N., & Stone, J. M. 2010b, ApJ, 722, L220 [NASA ADS] [CrossRef] [Google Scholar]
 Birnstiel, T., Dullemond, C. P., & Brauer, F. 2010, A&A, 513, A79 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Blum, J., & Wurm, G. 2008, ARA&A, 46, 21 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Cuzzi, J. N., Estrada, P. R., & Davis, S. S. 2014, ApJS, 210, 21 [NASA ADS] [CrossRef] [Google Scholar]
 Dominik, C., & Dullemond, C. P. 2008, A&A, 491, 663 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Dominik, C., & Tielens, A. G. G. M. 1995, Philosophical Magazine, Part A, 72, 783 [Google Scholar]
 Dominik, C., & Tielens, A. G. G. M. 1997, ApJ, 480, 647 [NASA ADS] [CrossRef] [Google Scholar]
 Drążkowska, J., & Dullemond, C. P. 2014, A&A, 572, A78 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Drążkowska, J., Windmark, F., & Dullemond, C. P. 2014, A&A, 567, A38 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Dullemond, C. P., & Dominik, C. 2005, A&A, 434, 971 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Garaud, P., Meru, F., Galvagni, M., & Olczak, C. 2013, ApJ, 764, 146 [NASA ADS] [CrossRef] [Google Scholar]
 Gillespie, D. T. 1975, J. Atmos. Sci., 32, 1977 [NASA ADS] [CrossRef] [Google Scholar]
 Gundlach, B., & Blum, J. 2014, ApJ, accepted [Google Scholar]
 Gundlach, B., Kilias, S., Beitz, E., & Blum, J. 2011, Icarus, 214, 717 [NASA ADS] [CrossRef] [Google Scholar]
 Güttler, C., Blum, J., Zsom, A., Ormel, C. W., & Dullemond, C. P. 2010, A&A, 513, A56 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Haisch, Jr., K. E., Lada, E. A., & Lada, C. J. 2001, ApJ, 553, L153 [NASA ADS] [CrossRef] [Google Scholar]
 Hayashi, C. 1981, Prog. Theor. Phys. Suppl., 70, 35 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Johansen, A., Oishi, J. S., Mac Low, M.M., et al. 2007, Nature, 448, 1022 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Johansen, A., Blum, J., Tanaka, H., et al. 2014, in Protostars and Planets VI, University of Arizona Press, eds. H. Beuther, et al., in press [arXiv:1402.1344] [Google Scholar]
 Kataoka, A., Tanaka, H., Okuzumi, S., & Wada, K. 2013a, A&A, 557, L4 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kataoka, A., Tanaka, H., Okuzumi, S., & Wada, K. 2013b, A&A, 554, A4 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kataoka, A., Okuzumi, S., Tanaka, H., & Nomura, H. 2014, A&A, 568, A42 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kempf, S., Pfalzner, S., & Henning, T. K. 1999, Icarus, 141, 388 [NASA ADS] [CrossRef] [Google Scholar]
 Kirchschlager, F., & Wolf, S. 2014, A&A, 568, A103 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kothe, S., Güttler, C., & Blum, J. 2010, ApJ, 725, 1242 [NASA ADS] [CrossRef] [Google Scholar]
 Kretke, K. A., & Levison, H. F. 2014, AJ, 148, 109 [NASA ADS] [CrossRef] [Google Scholar]
 Krijt, S., Güttler, C., Heißelmann, D., Dominik, C., & Tielens, A. G. G. M. 2013, J. Phys. D Appl. Phys., 46, 5303 [NASA ADS] [CrossRef] [Google Scholar]
 Krijt, S., Dominik, C., & Tielens, A. G. G. M. 2014, J. Phys. D Appl. Phys., 47, 175302 [NASA ADS] [CrossRef] [Google Scholar]
 Lambrechts, M., & Johansen, A. 2012, A&A, 544, A32 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lambrechts, M., & Johansen, A. 2014, A&A, 572, A107 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Meisner, T., Wurm, G., & Teiser, J. 2012, A&A, 544, A138 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mukai, T., Ishimoto, H., Kozasa, T., Blum, J., & Greenberg, J. M. 1992, A&A, 262, 315 [NASA ADS] [Google Scholar]
 Nakagawa, Y., Sekiya, M., & Hayashi, C. 1986, Icarus, 67, 375 [NASA ADS] [CrossRef] [Google Scholar]
 Okuzumi, S., Tanaka, H., & Sakagami, M.A. 2009, ApJ, 707, 1247 [NASA ADS] [CrossRef] [Google Scholar]
 Okuzumi, S., Tanaka, H., Takeuchi, T., & Sakagami, M.A. 2011, ApJ, 731, 95 [NASA ADS] [CrossRef] [Google Scholar]
 Okuzumi, S., Tanaka, H., Kobayashi, H., & Wada, K. 2012, ApJ, 752, 106 [NASA ADS] [CrossRef] [Google Scholar]
 Ormel, C. W., & Cuzzi, J. N. 2007, A&A, 466, 413 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ormel, C. W., & Klahr, H. H. 2010, A&A, 520, A43 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ormel, C. W., & Spaans, M. 2008, ApJ, 684, 1291 [NASA ADS] [CrossRef] [Google Scholar]
 Ormel, C. W., Spaans, M., & Tielens, A. G. G. M. 2007, A&A, 461, 215 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Paszun, D., & Dominik, C. 2008, A&A, 484, 859 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Paszun, D., & Dominik, C. 2009, A&A, 507, 1023 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Poppe, T., Blum, J., & Henning, T. 2000, ApJ, 533, 454 [NASA ADS] [CrossRef] [Google Scholar]
 Schräpler, R., & Blum, J. 2011, ApJ, 734, 108 [NASA ADS] [CrossRef] [Google Scholar]
 Seizinger, A., & Kley, W. 2013, A&A, 551, A65 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Seizinger, A., Krijt, S., & Kley, W. 2013, A&A, 560, A45 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Sekiya, M., & Takeda, H. 2005, Icarus, 176, 220 [NASA ADS] [CrossRef] [Google Scholar]
 Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337 [NASA ADS] [Google Scholar]
 Suyama, T., Wada, K., & Tanaka, H. 2008, ApJ, 684, 1310 [NASA ADS] [CrossRef] [Google Scholar]
 Suyama, T., Wada, K., Tanaka, H., & Okuzumi, S. 2012, ApJ, 753, 115 [NASA ADS] [CrossRef] [Google Scholar]
 Testi, L., Birnstiel, T., Ricci, L., et al. 2014, in Protostars and Planets VI, University of Arizona Press, eds. H. Beuther, et al., in press [arXiv:1402.1354] [Google Scholar]
 Tielens, A. G. G. M., McKee, C. F., Seab, C. G., & Hollenbach, D. J. 1994, ApJ, 431, 321 [NASA ADS] [CrossRef] [Google Scholar]
 Wada, K., Tanaka, H., Suyama, T., Kimura, H., & Yamamoto, T. 2007, ApJ, 661, 320 [NASA ADS] [CrossRef] [Google Scholar]
 Wada, K., Tanaka, H., Suyama, T., Kimura, H., & Yamamoto, T. 2009, ApJ, 702, 1490 [NASA ADS] [CrossRef] [Google Scholar]
 Wada, K., Tanaka, H., Suyama, T., Kimura, H., & Yamamoto, T. 2011, ApJ, 737, 36 [NASA ADS] [CrossRef] [Google Scholar]
 Wada, K., Tanaka, H., Okuzumi, S., et al. 2013, A&A, 559, A62 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Weidenschilling, S. J. 1977, MNRAS, 180, 57 [NASA ADS] [CrossRef] [Google Scholar]
 Weidenschilling, S. J. 1980, Icarus, 44, 172 [NASA ADS] [CrossRef] [Google Scholar]
 Whipple, F. L. 1972, in From Plasma to Planet, ed. A. Elvius, 211 [Google Scholar]
 Windmark, F., Birnstiel, T., Ormel, C. W., & Dullemond, C. P. 2012, A&A, 544, L16 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Wurm, G., Paraskov, G., & Krauss, O. 2004, ApJ, 606, 983 [NASA ADS] [CrossRef] [Google Scholar]
 Wurm, G., Paraskov, G., & Krauss, O. 2005, Icarus, 178, 253 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Youdin, A. N., & Goodman, J. 2005, ApJ, 620, 459 [NASA ADS] [CrossRef] [Google Scholar]
 Youdin, A. N., & Lithwick, Y. 2007, Icarus, 192, 588 [NASA ADS] [CrossRef] [Google Scholar]
 Zsom, A., & Dullemond, C. P. 2008, A&A, 489, 931 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Zsom, A., Ormel, C. W., Güttler, C., Blum, J., & Dullemond, C. P. 2010, A&A, 513, A57 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
All Figures
Fig. 1 Particle Stokes numbers as a function of mass, in the midplane of an MMSN disk at 5 AU. Different lines show compact particles (red), porous aggregates with constant φ = 10^{4} (yellow), and aggregates with a constant fractal dimension of 2.5 (green). For the solid lines, all drag regimes (Epstein, Stokes and Newton) have been taken into account, while the dashed lines indicate the results using only Epstein and Stokes drag. Horizontal lines indicate Ωt_{s} = 1 (where drift is fastest) and Ωt_{s} = α = 10^{3} (where particles start to settle to the midplane). 

Open with DEXTER  
In the text 
Fig. 2 Relative velocities between compact (left) or very porous (right) particles with masses m_{i} and m_{j} at the midplane of an MMSNdisk at 5 AU with α = 10^{3}. The masses range from single monomers to aggregates containing 10^{32} monomers. The contours give relative velocities in m s^{1}, and the colors indicate the dominating source for the relative velocity: Brownian motion (BrM), turbulence (Trb), settling (Set), radial drift (Rad), or azimuthal drift (Azi). Epstein, Stokes, and Newton drag regimes have been taken into account. 

Open with DEXTER  
In the text 
Fig. 3 a) Schematic of a collision between unequal particles with a mass ratio R^{(m)} = (m_{proj}/m_{target}) ≪ 1. b) Sticking occurs when v_{rel}<v_{eros}. The mass of the projectile is added to the target. c) Collisions above the erosion threshold velocity lead to erosion. The mass loss of the target is given by the erosion efficiency ϵ_{eros} and the mass of the projectile. 

Open with DEXTER  
In the text 
Fig. 4 Evolution of the normalized particle mass distribution at 5 AU with α = 10^{3}, assuming perfect sticking and without compaction through gas and selfgravity. Only Epstein are Stokes drag are considered. Solid lines indicate averages over 4 MC runs with identical starting conditions, and the shaded areas represent a spread of 1σ. 

Open with DEXTER  
In the text 
Fig. 5 Like Fig. 4, but with compaction through gas and selfgravity and Newton drag for particles with Re_{p}> 1. 

Open with DEXTER  
In the text 
Fig. 6 Evolution of the growth and radial drift timescale of the peak mass for the perfect sticking model at 5 AU with α = 10^{3}. The dotted line indicates (t_{drift}/ 30). Only collisional compaction has been taken into account. 

Open with DEXTER  
In the text 
Fig. 7 Evolution of the growth and radial drift timescale of the peak mass for the perfect sticking model at 5 AU with α = 10^{3}. The dotted line indicates (t_{drift}/ 30). Compaction from gas and selfgravity, and Newton drag have been taken into account. 

Open with DEXTER  
In the text 
Fig. 8 Evolution of the internal structure of the massdominating particles, for the perfect sticking models at 5 AU, for the models with and without noncollisional compaction mechanisms. Aggregates start out as monomers in the top left corner, and grow towards larger sizes and porosities. Lines show individual simulations. Open symbols correspond to points where the mass dominating particles reach a = λ_{mfp} (°); t_{s} = t_{η} (◇); Ωt_{s} = α (▿); and Ωt_{s} = 1 (□). Filled symbols show peak mass and filling factor at the times of first: collisional compaction (⋆); gaspressure compaction (◆); and selfgravity compaction (•). 

Open with DEXTER  
In the text 
Fig. 9 Evolution of the normalized particle mass distribution at 5 AU with α = 10^{3}, assuming v_{eros} = 20 m s^{1}. The full compaction model is used. 

Open with DEXTER  
In the text 
Fig. 10 Projectile distribution mass functions for simulation E1, constructed for the maximum masses (m_{t} =•) at various times. Colors and times correspond to Fig. 9. For each distribution, the stopping time of the m_{t}particle is given, and the weights of the total positive and negative area are plotted. The distributions have been normalized in such a way, that the absolute sum of the contributions equals 1. 

Open with DEXTER  
In the text 
Fig. 11 Evolution of ζ(m_{max}) for erosive simulations with v_{eros} = 20,40,60 m s^{1} as a function of Ωt_{s}(m_{max}), showing the impact of erosion on the ability of the largest bodies to grow. In the upper two panels, the steadystate is indicated by the ⊚symbol. 

Open with DEXTER  
In the text 
Fig. 12 Masses and filling factors of all unique families at different times, for simulation E1. One dot corresponds to one family, and does not provide information about the total mass or number of members in that family. 

Open with DEXTER  
In the text 
Fig. 13 Evolution of m_{p}(t) and R(t) for dust coagulation as obtained from the semianalytical model (Eqs. (35) and (36)), for an MMSN disk and α = 10^{3}. Lines indicate different starting conditions R(t = 0), and are evolved for 10^{6} yrs. Left: compact growth: φ = 1 at all times, and v_{frag} = 10 m s^{1}. Right: porous growth: the internal structure of the aggregates is set by hit and stick growth, followed by collisional compaction or gas and selfgravity compaction. Gray lines have no erosion, while black lines show the results for v_{eros} = 40 m s^{1}. Colored lines and ⊚symbols indicate the evolution and steady state peak mass obtained through local MC simulations (Sect. 4.2). 

Open with DEXTER  
In the text 
Fig. 14 Same as the righthand plot of Fig. 13, but with aggregate size on the vertical axis. 

Open with DEXTER  
In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.