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. BR  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. BR  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 aggregates 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. BR  Results For erosion threshold velocities of 2040 m sSUP1SUP 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 sSUP1SUP do porous aggregates manage to cross the radial drift barrier in the inner 10 AU of MMSNlike disks. BR  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 10SUP9SUP g highly porous particles appears suitable for triggering streaming instability.
['protoplanetary disks', 'planets and satellites: formation', 'circumstellar matter', 'methods: numerical', 'Astrophysics - Earth and Planetary Astrophysics']
Erosion and the limits to planetesimal growth
{'Erosion and the limits to planetesimal growth': 'S. Krijt 1 , C.W. Ormel 2 ; 3 , C. Dominik 3 , and A.G.G.M. Tielens 1  \n- 1 Leiden Observatory, Leiden University, Niels Bohrweg 2, 2333 CA Leiden, The Netherlands e-mail: krijt@strw.leidenuniv.nl\n- 2 Astronomy Department, University of California, Berkeley, CA 94720\n- 3 Anton Pannekoek Institute, University of Amsterdam, Science Park 904, 1098 XH Amsterdam, The Netherlands', 'ABSTRACT': "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 e GLYPH<14> ciency of this early growth. In particular, it has been proposed that flu GLYPH<11> y ice aggregates can grow very e GLYPH<14> ciently in protoplanetary disks, su GLYPH<11> ering less from the bouncing and radial drift barriers.  \nAims. 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 pre-planetesimal bodies. We investigate the e GLYPH<11> ect of these erosive collisions on the ability of porous ice aggregates to cross the radial drift barrier.  \nMethods. 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 self-gravity compaction. We include erosive collisions and study the e GLYPH<11> ect of the erosion threshold velocity on aggregate growth.  \nResults. For erosion threshold velocities of 20 GLYPH<0> 40 m s GLYPH<0> 1 , high-velocity collisions with small projectiles prevent the largest aggregates from growing when they start to drift. In these cases, our local simulations result in a steady-state distribution, with the majority of the dust mass in particles with Stokes numbers close to unity. Only for the highest erosion threshold considered (60 m s GLYPH<0> 1 ), do porous aggregates manage to cross the radial drift barrier in the inner 10 AU of MMSN-like disks.  \nConclusions. Erosive collisions are more e GLYPH<11> ective in limiting the growth than fragmentary collisions between similar-size 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 GLYPH<24> 10 9 g, highly porous particles appears well-suited for triggering streaming instability.  \nKey words. protoplanetary disks - planet and satellites: formation - stars: circumstellar matter - methods: numerical", '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 GLYPH<24> 10 6 years (Haisch et al. 2001). The process of planetesimal formation is a complex one, with many di GLYPH<11> erent processes acting on a variety of length and timescales (see Testi et al. (2014) and Johansen et al. (2014) for recent reviews).  \nThe 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.  \nBut even before radial drift becomes problematic, coagulation of aggregates can be frustrated by catastrophic fragmentation or bouncing (Blum & Wurm 2008; Güttler et al. 2010; Zsom et al. 2010), which prevent colliding aggregates from gaining mass. These issues are alleviated somewhat by including velocity distributions between pairs of particles (Windmark et al.  \n2012; Garaud et al. 2013) in combination with mass transfer in high-velocity collisions (Wurm et al. 2005; Kothe et al. 2010), though these solutions require the presence of relatively compact targets.  \nRecently, it was proposed that icy aggregates, if they can manage to stay very porous, su GLYPH<11> er 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 flu GLYPH<11> y, 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 due to their enhanced collisional cross section (Okuzumi et al. 2012). However, Okuzumi et al. (2012) assumed perfect sticking between colliding aggregates, neglecting possible mass-loss in aggregateaggregate collisions.  \nIn this paper, we study the e GLYPH<11> ects 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 co-  \ngulation code, capable of simulating the vertically-integrated dust population, tracing both the evolution of the mass and the porosity of the entire mass distribution self-consistently. Section 2 describes the models we use for the protoplanetary disk and the dust aggregates. In section 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) (Section 4.1.1), after which we expand the model to include compaction from gas pressure and self-gravity according to Kataoka et al. (2013a) (Section 4.1.2), and erosive collisions (Section 4.2). In Section 5, we compare the results to a simple semi-analytical model, and describe which processes can limit coagulation in di GLYPH<11> erent parts of protoplanetary disks. Discussion of the results and implications takes place in Section 6, and conclusions are presented in Section 7.', '2. Disk and dust models': 'The disk model and collisional compaction prescription are based on Okuzumi et al. (2012), to which we add non-collisional compaction processes (Section 2.4.2) and a model for erosive collisions (Sections 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 Sun-like central star are given as  \nGLYPH<6> g = 152 GLYPH<18> R 5 AU GLYPH<19> GLYPH<0> 3 = 2 g cm GLYPH<0> 2 ;  \n(1) GLYPH<30> GLYPH<17> GLYPH<26> int GLYPH<26> 0 ; (7)  \nas a measure for the internal density.  \nT = 125 GLYPH<18> R 5 AU GLYPH<19> GLYPH<0> 1 = 2 K : (2)  \nThe gas sound speed is given by  \ncs = q k B T = m g = 6 : 7 GLYPH<2> 10 2 GLYPH<18> R 5 AU GLYPH<19> GLYPH<0> 1 = 4 ms GLYPH<0> 1 ; (3)  \nwith k B the Boltzmann constant and m g = 3 : 9 GLYPH<2> 10 GLYPH<0> 24 g the mean molecular weight. The Kepler frequency equals  \nGLYPH<10> = p GM GLYPH<12> = R 3 = 1 : 8 GLYPH<2> 10 GLYPH<0> 8 GLYPH<18> R 5 AU GLYPH<19> GLYPH<0> 3 = 2 s GLYPH<0> 1 : (4)  \nAssuming an isothermal column, the gas density drops with increasing distance from the mid plane z according to  \nGLYPH<26>g = GLYPH<6> g p 2 GLYPH<25> h g exp 0 B B B B @ GLYPH<0> z 2 2 h 2 g 1 C C C C A ; (5)  \nwith the relative vertical scale height of the gas h g= R = 0 : 05( R = 5 AU) 1 = 4 . The turbulent viscosity is parametrized as GLYPH<23> turb = GLYPH<11> c 2 s = GLYPH<10> following Shakura & Sunyaev (1973), and GLYPH<11> is assumed to be constant in both the radial and the vertical direction. The eddie turn-over time of the largest eddies equals tL = GLYPH<10> GLYPH<0> 1 .  \nIn our local model, the surface density of the dust is related to the gas surface density through GLYPH<6> d = GLYPH<6> g = 10 GLYPH<0> 2 , but the vertical distribution of dust depends on its aerodynamic properties. The dust is described by a Gaussian, with the dust scale height hd set  \nIn the rest of this section, we describe the main ingredients for the simulations presented in Section 3. These are: the relative velocities between aggregates, the equations governing the evolution of GLYPH<26> int through mutual collisions as well as gas ram pressure and self-gravity, 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 Section 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  \nts = 8 > > > > > > < > > > > > > : t (Ep) s = 3 m 4 GLYPH<26>g v th A for a < 9 4 GLYPH<21> mfp ; t (St) s = 4 a 9 GLYPH<21> mfp t (Ep) s for a > 9 4 GLYPH<21> mfp ; (8)  \nwhere v th = p 8 =GLYPH<25> cs is the mean thermal velocity of the gas molecules, and GLYPH<21> mfp = m g= ( GLYPH<27> mol GLYPH<26>g ) is the gas molecule mean free path. Taking GLYPH<27> mol = 2 GLYPH<2> 10 GLYPH<0> 15 cm 2 , we obtain GLYPH<21> mfp = 120( R = 5 AU) 11 = 4 cm at the disk mid plane. In Equation 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,  \nby the stopping time ts of the dust particle through (Youdin & Lithwick 2007)  \nhd h g = 1 + GLYPH<10> ts GLYPH<11> 1 + 2 GLYPH<10> ts 1 + GLYPH<10> ts ! GLYPH<0> 1 = 2 : (6)  \nThus, settling becomes important when a dust particle reaches GLYPH<10> ts GLYPH<24> GLYPH<11> .', '2.2. Dust properties': "Initially, all dust particles are assumed to be spherical (sub)micron-size 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 GLYPH<30> . 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 GLYPH<26> int = m = V , with V = (4 = 3) GLYPH<25> a 3 the volume of the aggregate, and a its radius. An aggregate's radius is defined as a = [5 = (3 N ) P N k = 1 ( r k GLYPH<0> r CM) 2 ] 1 = 2 , 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 GLYPH<26> int = m 0 = V 0 = GLYPH<26> 0, while aggregates can have GLYPH<26> int GLYPH<28> GLYPH<26> 0. Since we are interested in region beyond the snow-line, we focus here on monomers composed of mostly ice, and use a density of GLYPH<26> 0 = 1 : 4 g cm GLYPH<0> 3 . For the monomer radius we use a 0 = 0 : 1 GLYPH<22> m. We define the filling factor as  \nwhich can be obtained using the formulation of Okuzumi et al. (2009).  \nThe above equation is accurate when the particle Reynolds number Rep = 4 a v dg = ( v th GLYPH<21> 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  \nts = 2 m CD GLYPH<26>g v dg A : (9)  \nIn the Stokes regime the drag coe GLYPH<14> cient equals CD = 24 = Rep, 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  \nCD = 8 > > > > > > < > > > > > > : 24(Rep) GLYPH<0> 1 for Rep < 1 ; 24(Rep) GLYPH<0> 3 = 5 for 1 < Rep < 800 ; 0 : 44 for 800 < Rep : (10)  \nFigure 1 shows Stokes numbers ( GLYPH<10> ts ) for di GLYPH<11> erent particles in the mid plane of a MMSN disk at 5 AU. Di GLYPH<11> erent lines show compact particles (red), porous aggregates with constant GLYPH<30> = 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 Epstein and Stokes drag only, i.e. assuming that Rep < 1. Focussing on the Df = 2 : 5 aggregates, we can clearly distinguish the di GLYPH<11> erent drag regimes. The smallest particles are in the Epstein regime, and switch to the Stokes regime around GLYPH<10> ts = 10 GLYPH<0> 3 . Then, at a mass of m = m 0 GLYPH<24> 10 21 , the Reynolds number exceeds unity and we enter the second regime of Equation 10. Note that this transition occurs before GLYPH<10> ts = 1. The most massive particles, m = m 0 > 10 26 are in the regime where CD = 0 : 44. Compact particles on the other hand, reach GLYPH<10> ts = 1 while still in the Epstein drag regime.  \nThe turbulence-induced relative velocity between two particles with stopping times ts ; 1 and ts ; 2 GLYPH<20> ts ; 1 has three regimes (Ormel & Cuzzi 2007)  \nv turb ' GLYPH<14>vg GLYPH<2> 8 > > > > > > > > < > > > > > > > > : Ret 1 = 4 GLYPH<10> ( ts ; 1 GLYPH<0> ts ; 2) for ts ; 1 GLYPH<28> t GLYPH<17> ; 1 : 4 : : : 1 : 7 GLYPH<0> GLYPH<10> ts ; 1 GLYPH<1> 1 = 2 for t GLYPH<17> GLYPH<28> ts ; 1 GLYPH<28> GLYPH<10> GLYPH<0> 1 ; 1 1 + GLYPH<10> ts ; 1 + 1 1 + GLYPH<10> ts ; 2 ! 1 = 2 for ts ; 1 GLYPH<29> GLYPH<10> GLYPH<0> 1 ; (11)  \nwhere GLYPH<14>vg = GLYPH<11> 1 = 2 cs is the mean random velocity of the largest turbulent eddies, and t GLYPH<17> = Ret 1 = 2 tL is the turn-over time of the smallest eddies. The turbulence Reynolds number is given by Ret = GLYPH<11> c 2 s = ( GLYPH<10> GLYPH<23> mol), with the molecular viscosity GLYPH<23> mol = v th GLYPH<21> mfp = 2. We will refer to the first two cases of Equation 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  \nFigure 2 shows the mid plane relative velocity in m s GLYPH<0> 1 (contours), and its dominant source (color), for a range of combinations of masses mi and mj . The velocities have been calculated for the disk properties of Section 2.1, at 5 AU, and assuming a turbulence GLYPH<11> = 10 GLYPH<0> 3 . The left plot corresponds to two compact particles ( GLYPH<26> int = GLYPH<26> 0), and the right plot to two very porous ones ( GLYPH<26> int = 10 GLYPH<0> 4 GLYPH<26> 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 GLYPH<10> ts GLYPH<24> 1. However, the masses at which various transitions occur can vary 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 di GLYPH<11> erential settling.  \n<!-- image -->  \n0  \nFig. 1. Particle Stokes numbers as a function of mass, in the mid plane of an MMSN disk at 5 AU. Di GLYPH<11> erent lines show compact particles (red), porous aggregates with constant GLYPH<30> = 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 Epstein and Stokes drag only. Horizontal lines indicate GLYPH<10> ts = 1 (where drift is fastest) and GLYPH<10> ts = GLYPH<11> = 10 GLYPH<0> 3 (where particles start to settle to the mid plane).  \nregime because of the ( ts ; 1 GLYPH<0> ts ; 2) term, but considerably larger in the second regime.", '2.3. Collisional outcomes': 'Acollision 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) o GLYPH<11> er 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).  \nIn general, a collision can result in sticking, erosion, or fragmentation, depending on the relative velocity and the mass ratio R ( m ) GLYPH<17> mi = mj GLYPH<20> 1 of the colliding bodies. Collisions between particles with comparable masses result in catastrophic fragmen-  \nlocity. We treat this dispersion in the same way as Okuzumi et al. (2012), by taking into account the standard deviation in the mass-to-area ratio of a porous aggregate (Okuzumi et al. 2011). The size of this standard deviation, normalized by the mean mass-to-area ratio, is parametrized as " , which we take to equal 0.1, following Okuzumi et al. (2011).  \nFig. 2. Relative velocities between compact (left) or very porous (right) particles with masses mi and mj at the mid plane of an MMSN-disk at 5 AU with GLYPH<11> = 10 GLYPH<0> 3 . The masses range from single monomers to aggregates containing 10 32 monomers. The contours give relative velocities in ms GLYPH<0> 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.  \n<!-- image -->  \n<!-- image -->  \ni  \ntation if they collide above the fragmentation velocity (Section 2.3.1). When colliding bodies have a mass ratio R ( m ) GLYPH<28> 1, catastrophic fragmentation of the larger body is di GLYPH<14> cult, 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 R ( m ) crit , specified in Section 2.3.3. In an erosive event, the larger body will lose mass. From Figure 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 Section 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 ) GLYPH<21> 1 = 64), Wada et al. (2013) find a critical fragmentation velocity of  \nv frag ' 20 E break m 0 ! 1 = 2 ' 80 a 0 0 : 1 GLYPH<22> m ! GLYPH<0> 5 = 6 ms GLYPH<0> 1 : (12)  \nThe quantity E break represents the energy needed to break a single mononer-monomer 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 similar-sized aggregates will generally not reach the fragmentation velocity (Equation 12) behind  \nthe snow line in a protoplanetary disk, especially not if the turbulence is weak. However, relative velocities between particles with very di GLYPH<11> erent masses can be much larger than velocities between similar particles, especially when radial and azimuthal drift are important (Figure 2). In this paper, we study the e GLYPH<11> ects 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 e GLYPH<14> ciency, GLYPH<15> eros, that indicates how much mass is eroded in units of projectile mass.  \nFor silicate particles, Güttler et al. (2010) summarize a number of experimental investigations and describe a threshold velocity of a few m s GLYPH<0> 1 , and an erosion e GLYPH<14> ciency 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 GLYPH<0> 1 using micron-size silicate projectiles. Note that the threshold velocity is comparable to the monomer sticking velocity of micron-size silicate particles (Poppe et al. 2000). In the experiments of Schräpler & Blum (2011), the erosion e GLYPH<14> ciency also increased with impact velocity, reaching GLYPH<24> 10 for the highest velocity of 60 m s GLYPH<0> 1 (their Figure 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 e GLYPH<14> ciency with projectile mass, showing a trend of decreasing erosion threshold with decreasing mass ratio (e.g. their Figure 11). For monomer projectiles, the threshold velocity equals the monomer-monomer sticking velocity v s ' p E break = m 0, after which it increased linearly with velocity to eventually flatten o GLYPH<11> around 10 m s GLYPH<0> 1 . This flattening o GLYPH<11> indicates the onset of catastrophic fragmentation, and occurs at a mass ratio of GLYPH<24> 10 GLYPH<0> 2 .  \nFor ice particles, Gundlach & Blum (2014) present recent experimental results on the sticking and erosion threshold of  \nFig. 3. (a): Schematic of a collision between unequal particles with a mass ratio R ( m ) = ( m proj = m target) GLYPH<28> 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 e GLYPH<14> ciency GLYPH<15> eros and the mass of the projectile.  \n<!-- image -->  \n(b) After sticking  \n(  \nv  \nrel <  \nv  \neros)  \n(c) After erosion  \n(  \nv  \nrel >  \nv  \neros)  \nm target -ϵ eros m proj  \n<!-- image -->  \n(sub)micron-size particles. For a projectile distribution between 0 : 2 GLYPH<0> 6 GLYPH<22> m(with a mean value of 1 : 5 GLYPH<22> m) impinging an icy target with a filling factor GLYPH<30> ' 0 : 5, an erosion threshold of 15 : 3 m s GLYPH<0> 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.1GLYPH<22> m monomers, possibly even > 60 m s GLYPH<0> 1 . However, the aggregates acting as targets in the simulations presented here have a much higher porosity ( GLYPH<30> GLYPH<24> 10 GLYPH<0> 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 di GLYPH<14> cult. For a single  \nm target + m proj  \n<!-- image -->  \ngrain size, the size significantly influences the strength of the aggregates, with larger monomers leading to weaker aggregates (Equation 12). Little is known about the expected grain sizes in the icy regions of protoplanetary disks, let alone their size distribution, or about the e GLYPH<11> ect a monomer size distribution has on the strength and collisional behavior of porous aggregates.  \nFor 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 mass-loss in low mass ratio collisions is relatively small. Using molecular dynamics N-body 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 similar-size particles colliding at v frag fragment e GLYPH<14> -ciently. 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. Equation 12 and Figure 2). In this paper, we take the agnostic view by 'burying' the uncertainty of the erosion threshold velocity in the parameter v eros, which we vary to investigate the implications of e GLYPH<11> ective versus ine GLYPH<11> ective erosion.", '2.3.3. Erosion model': 'Erosive collisions occur only below a mass ratio R ( m ) crit , and their outcome is parametrized in terms of a (velocity-dependent) erosion e GLYPH<14> ciency. In accordance with Güttler et al. (2010) and Seizinger et al. (2013) we will use R ( m ) crit = 10 GLYPH<0> 2 . 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 GLYPH<0> 1 , corresponding to (1 = 4) v frag and (3 = 4) v frag for 0 : 1 GLYPH<22> m monomers (Equation 12). In the erosive regime, the normalized erosion e GLYPH<14> ciency can be written as  \nGLYPH<15> eros = c 1 v rel v eros ! GLYPH<13> ; (13)  \nwith c 1 GLYPH<24> 1 (Güttler et al. 2010; Seizinger et al. 2013). While in supersonic cratering collisions GLYPH<13> = 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 GLYPH<13> = 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).  \nLastly, 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 Df is estimated as  \nDf \' 3 " 1 GLYPH<0> log( GLYPH<30> ) log( m = m 0) # GLYPH<0> 1 : (14)  \nThe assumptions of the erosion model employed in this work are discussed further in Section 6.', '2.4. Aggregate compaction': "An aggregate's porosity can be altered through collisions, or through non-collisional mechanisms. In this Section, we first describe how porosity can increase and decrease as the result of  \nsticking collisions. Then, we discuss gas- and self-gravity 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 mi + mj . 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 hit-and-stick 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 Equation (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 GLYPH<14> (Dominik & Tielens 1997).  \nGundlach et al. (2011) measured the rolling force between ice particles with radii of GLYPH<24> 1 : 5 GLYPH<22> mto be 1 : 8 GLYPH<2> 10 GLYPH<0> 3 dyn, implying a rolling energy of 1 : 8 GLYPH<2> 10 GLYPH<0> 7 erg. Assuming the rolling force is size-independent (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 F roll / a 2 = 3 0 , and E roll / a 5 = 3 0 , 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 GLYPH<2> 10 GLYPH<0> 9 erg for 0.1GLYPH<22> 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 self-gravity 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  \nPc = E roll a 3 0 GLYPH<30> 3 : (15)  \nThis pressure can then be compared to the pressure arising form the surrounding gas and from self-gravity  \nP gas = v dg m GLYPH<25> a 2 ts ; P grav = Gm 2 GLYPH<25> a 4 ; (16)  \nwith G the gravitational constant, in order to see whether an aggregate will be compacted as a result of these non-collisional processes (Kataoka et al. 2013a). In this work, we will take these e GLYPH<11> ects into account in a self-consistent 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 : integro-di GLYPH<11> erential 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 integro-di GLYPH<11> erential 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 Monte Carlo method.  \nThe 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.  \nLet f ( x ) be the (time-dependent) 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  \nCij = Kij = S ; (17)  \nwith S the surface area of the column 4 , and Kij the collision kernel, which in this case equals  \nKij = GLYPH<27> i j 2 GLYPH<25> hd ; ihd ; j Z 1 GLYPH<0>1 v rel( z ) exp 0 B B B B B @ GLYPH<0> z 2 2 h 2 d ; i j 1 C C C C C A dz ; (18)  \nwhere hd ; i is given by Equation 6, and hd ; i j = ( h GLYPH<0> 2 d ; i + h GLYPH<0> 2 d ; j ) GLYPH<0> 1 = 2 and GLYPH<27> i j = GLYPH<25> ( ai + aj ) 2 equals the collisional cross section (Okuzumi et al. 2012). This rate equation takes into account that particles with di GLYPH<11> erent properties inhabit di GLYPH<11> erent vertical scale heights, and is correct as long as the coagulation timescale is longer than the vertical settling / di GLYPH<11> usion timescale. In this work, we approximate the integral over z by assuming the mid plane relative velocity is a good indication for v rel throughout the column. This allows us to solve the integral analytically and write  \nKij \' GLYPH<27> i jhd ; i j p 2 GLYPH<25> hihj v rel( z = 0) : (19)  \nFor the purpose of this paper, this approximation is su GLYPH<14> ciently accurate, since most of the growth is expected to take place near the mid plane.  \nThen, we can define the total collision rate for particle Ci = P j > i Cij , and the total collision rate C tot = P i Ci . With all these rates known, 3 random numbers are used to identify which particles collide, and the time GLYPH<1> 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 Ci have to be adjusted, since the particle distribution f has changed. This cycle is then repeated.  \nThe 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 Nf 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 zi 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 = P 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  \nGLYPH<21> i j = ( w i w j 2 zi Cij for i , j ; w i ( w i 2 zi GLYPH<0> 1) Cii for i = j ; (20)  \nwhere the i = j case in Equation 20 describe so-called in-group collisions. Like before, we can define the total collision rate per family GLYPH<21> i = P j GLYPH<21> i GLYPH<21> i j , and the total collision rate GLYPH<21> tot = P i GLYPH<21> 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 mi GLYPH<29> mj . Thus, a total of 2 zi i -particles will collide with 2 zj j -particles. Assuming zj GLYPH<29> zi , every i -particle in the group will collide with 2 zj GLYPH<0> zi j -particles 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 GLYPH<0> zi GLYPH<0> N " .  \nLet GLYPH<14> mi be the change in the mass of the larger particle i , after a single collision with a j -particle. Assuming the changes are small, we can then extrapolate to find the total change after the full sequence of collisions  \nGLYPH<1> mi mi = 2 zj GLYPH<0> zi GLYPH<0> N " GLYPH<14> mi mi : (21)  \nNow, by imposing that ( GLYPH<1> mi = mi ) GLYPH<20> fm , we obtain  \nN ( m ) " = " GLYPH<0> log 2 fmmi 2 zi GLYPH<14> mi 2 zj !# ; (22)  \nwhere the square brackets indicate that N ( m ) " is truncated to integers GLYPH<21> 0, which has the e GLYPH<11> ect of particles with mass ratios GLYPH<21> fm always colliding 1-on-1. In the case of perfect sticking, obviously GLYPH<14> mi = mj , and Equation 22 reduces to Equation 12 of Ormel & Spaans (2008). We write an equivalent expression for the filling factor of the bigger grain  \nN ( GLYPH<30> ) " = " GLYPH<0> log 2 f GLYPH<30> GLYPH<30> i 2 zi GLYPH<14>GLYPH<30> i 2 zj !# ; (23)  \nwhere GLYPH<14>GLYPH<30> i denotes the change in GLYPH<30> after a single collision. The two limits are combined by writing  \nN " = max GLYPH<16> N ( m ) " ; N ( GLYPH<30> ) " GLYPH<17> ; (24)  \nand ensure that neither the filling factor, nor the mass of the larger particle change by too much during a single Monte Carlo 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 GLYPH<14>GLYPH<30> i (and GLYPH<14> mi , when erosion is present). In this work, we will typically use fm = f GLYPH<30> = 0 : 1.  \nImposing this limit has two consequences. First, since the group of i -particles can now only collide with part of the group of j -particles, this needs to be taken into account when the group collision rates are calculated, changing Equation 20 into  \nGLYPH<21> i j = ( w i w j 2 zi + N " Cij for i , j ; w i ( w i 2 zi GLYPH<0> 1) Cii for i = j : (25)  \nSecond, 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 GLYPH<21> 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 zi ; either creating many groups with few members (low zi ) or a few groups with many members (high zi ). 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).  \nOne approach is the so-called "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 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 di GLYPH<11> erence between the two methods is nicely illustrated in Figure 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.  \nIn 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  \nz GLYPH<3> = GLYPH<20> log 2 GLYPH<18> N 10 w GLYPH<3> GLYPH<19>GLYPH<21> ; (26)  \nwhere w GLYPH<3> is the desired number of groups per mass decade. In this way, we construct a function z GLYPH<3> ( 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 ( zi ! zi GLYPH<0> 1, w i ! 2 w i ) until zi = z GLYPH<3> ( mi ). Similarly, if the zoom number is too small, we "demagnify" ( zi ! zi + 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 GLYPH<3> = 60 for the perfect sticking calculations, and w GLYPH<3> = 40 for the ones including erosion.', '4.1. Perfect sticking': '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 zi = z GLYPH<3> ( mi ). 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 ( GLYPH<30> ) ) 3 , where R ( GLYPH<30> ) GLYPH<20> 1 is the ratio of the filling factors 5 . Then, we merge the families into a new family k with properties  \ng k = g i + g j ; mk = mi g i + mj g j g i + g j ; GLYPH<30> k = GLYPH<30> i g i + GLYPH<30> j g j g i + g j : (27)  \nThe new zoom- and group numbers are chosen such that zk = z GLYPH<3> ( mk ). Merging is necessary to suppress the total number of groups.', '3.5. Non-collisional compaction': 'Non-collisional compaction is implemented as follows: whenever a new aggregate is created in a collision, we calculate its compressive strength using Equation 15, and compare this to the external pressures from gas ram pressure and self-gravity, calculated with Equation 16 (Kataoka et al. 2013a). If either one of the external pressures exceeds Pc , we compactify the dust grain (i.e. increase GLYPH<30> ) until the aggregate can withstand the external pressures.', '3.6. Erosion': 'For every collision, we check if the conditions for erosion are met (i.e. v rel > v eros and R ( m ) < R ( m ) crit ), and if so, we determine the erosion e GLYPH<14> ciency using Equation 13. After a single erosive event, the mass that does not end up in the target body equals (1 + GLYPH<15> eros) m proj, see Figure 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 di GLYPH<11> erent 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  \nma = h mi i ; GLYPH<30> a = h GLYPH<30> i i ; (28)  \nwhich trace the properties of the average particle, and the peak mass and filling factor  \nmp = h m 2 i i h mi i ; GLYPH<30> p = h mi GLYPH<30> i i h mi i ; (29)  \nwhich trace the properties of the mass-dominating particle. We will also use the maximum mass m max, which is simply the mass of most massive particle.', '4.1.1. Collisional compaction only': "As a test for the Monte Carlo 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 GLYPH<11> = 10 GLYPH<0> 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 for a direct comparison to the work of Okuzumi et al., we do not include the e GLYPH<11> ects of Newton drag for particles with large Reynolds numbers in this simulation. In the rest of this work, Newton drag is always included self-consistently.  \nFigure 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 Monte Carlo runs. Thanks to the distribution method described in Section 3.3, the sampling of the mass distribution is  \n<!-- image -->  \n0  \nFig. 4. Evolution of the normalized particle mass distribution at 5 AU with GLYPH<11> = 10 GLYPH<0> 3 , assuming perfect sticking and without compaction through gas and self-gravity. Only Epstein are Stokes drag are considered. Solid lines indicate averages over 4 Monte Carlo runs with identical starting conditions, and the shaded areas represent a spread of 1 GLYPH<27> .Fig. 5. Like Figure 4, but with compaction through gas and self-gravity and Newton drag for particles with Rep > 1.  \n<!-- image -->  \nFig. 6. Evolution of the growth- and radial drift timescale of the peak mass for the perfect sticking model at 5 AU with GLYPH<11> = 10 GLYPH<0> 3 . The dotted line indicates ( t drift = 30). Only collisional compaction has been taken into account.  \n<!-- image -->  \nFig. 7. Evolution of the growth- and radial drift timescale of the peak mass for the perfect sticking model at 5 AU with GLYPH<11> = 10 GLYPH<0> 3 . The dotted line indicates ( t drift = 30). Compaction from gas and self-gravity, and Newton drag have been taken into account.  \n<!-- image -->  \nvery good over the entire mass range: even at later times, when most of the mass is located in particles with masses of GLYPH<24> 10 15 g, the distribution of particles all the way down to 10 GLYPH<0> 9 g is resolved remarkably well, despite these particles only making up a very small fraction of the total mass.  \nWhen we compare Figure 4 to Figure 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 GLYPH<10> ts GLYPH<24> 1 ( mj ' 10 10 g in this case), systematic drift greatly increases their collision rate, and very rapid growth ensues. The slight di GLYPH<11> erence in timescales is attributed to i ) the slightly di GLYPH<11> erent value for the rolling energy, ii ) our approximation of Equation 18, and iii ) our use Equation 8 to calculate the stopping times, while Okuzumi et al. used ts = t (Ep) s + t (St) s to ensure a smooth transition between Epstein and Stokes drag (S. Okuzumi, private communication).  \nWhile we take into account drift-induced 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,  \ndefined as  \nt grow GLYPH<17> mp ( dmp = dt ) ; (30)  \nto the radial drift timescale at that mass  \nt drift GLYPH<17> R v drift( mp ) : (31)  \nThe radial drift velocity is given by (Weidenschilling 1977)  \nv drift = GLYPH<0> 2 GLYPH<10> ts 1 + ( GLYPH<10> ts ) 2 GLYPH<17>v K ; (32)  \nwhere v K = R GLYPH<10> is the Keplerian orbital velocity, and GLYPH<17> can be written as (Nakagawa et al. 1986)  \nGLYPH<17> GLYPH<17> GLYPH<0> 1 2 cs v K ! 2 @ ln( GLYPH<26>g c 2 s ) @ ln R = 4 GLYPH<2> 10 GLYPH<0> 3 GLYPH<18> R 5 AU GLYPH<19> 1 = 2 : (33)  \nFigure 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 GLYPH<0> 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 ts ( mp ) > t GLYPH<17> . In this regime, velocities between similar particles are increased (see Equation 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 out-grow the radial drift barrier.", '4.1.2. Including gas and self-gravity compaction': 'The next step is to include compaction by gas pressure and selfgravity, as described in Section 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 Figure 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 self-gravity, resulting in a smaller collisional cross section. In addition, the aerodynamic properties are di GLYPH<11> erent, which affects the relative velocities.  \nThe growth- and drift timescales are plotted in Figure 7. When we compare Figures 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 non-collisional compaction. In Figure 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,  \nFig. 8. Evolution of the internal structure of the mass-dominating particles, for the perfect sticking models at 5 AU, for the models with and without non-collisional 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 = GLYPH<21> mfp ( GLYPH<14> ); ts = t GLYPH<17> ( ^ ); GLYPH<10> ts = GLYPH<11> ( O ); and GLYPH<10> ts = 1 ( GLYPH<3> ). Filled symbols show peak mass and filling factor at the times of first: collisional compaction ( ? ); gas-pressure compaction ( GLYPH<7> ); and self-gravity compaction ( GLYPH<15> ).  \n<!-- image -->  \nand closed symbols indicate the onset of various compaction mechanisms 6 .  \nInitially, aggregates grow through hit-and-stick collisions, and evolve along a line of constant fractal dimension close to 2. In the collisional-compaction-only scenario, particles reach a filling factor of GLYPH<24> 10 GLYPH<0> 5 during hit and stick growth, before collisional compaction kicks in, after which GLYPH<30> stays almost constant. When GLYPH<10> ts ( mp ) > 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 Figure 10 of Okuzumi et al. (2012). When non-collisional 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 Figure 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 GLYPH<10> ts > GLYPH<11> , which corresponds to m GLYPH<24> 10 GLYPH<0> 3 g for compact particles (see Figure 1). From Figure 8 however, we see that porous particles only begin to settle when their masses reach GLYPH<24> 10 4 GLYPH<0> 10 5 g. Lastly, aggregates with masses above GLYPH<24> 10 10 g are compacted by self-gravity, 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 non-collisional compaction mechanisms, and Epstein, Stokes, and Newton drag selfconsistently.  \nFig. 9. Evolution of the normalized particle mass distribution at 5 AU with GLYPH<11> = 10 GLYPH<0> 3 , assuming v eros = 20 m s GLYPH<0> 1 . The full compaction model is used.  \n<!-- image -->  \n0', '4.2. Erosion': "With this framework in place, the final step is to include the erosion model of Section 2.3.3 in the simulations, and calculate the evolution of the particle distribution self-consistently. Figure 9 shows the mass distribution at various times for v eros = 20 m s GLYPH<0> 1 . Initially, the evolution proceeds just like in 5, but as the largest aggregates approach GLYPH<10> ts = 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 GLYPH<24> 4000 yr a steady-state is reached, with a significant amount of mass residing in particles smaller than a few grams.  \nTo investigate how erosion halts the growth of the largest bodies, it is instructive to plot so-called projectile mass distributions (Okuzumi et al. 2009). For a certain particle mass mt , these distributions show the contribution to the growth of that particle as a function of projectile mass m GLYPH<20> mt . An example of such a plot is shown in Figure 9 of Okuzumi et al. (2012), where the distribution function is plotted at various times for mt = mp . 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 mt = m max, with m max the largest mass in the simulation at a given time. Second, to illustrate the e GLYPH<11> ect of erosive collisions, we calculate the mass loss for every erosive collision, taking into account the correct erosion e GLYPH<14> ciency 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 Figure 9 (colors correspond to the same times). When we examine the right-most projectile mass distribution, corresponding to a time t = 10 4 yr, it is immediately clear how erosion a GLYPH<11> ects the evolution of particles with GLYPH<10> ts GLYPH<24> 1. While these aggregates grow by collisions with similar-sized bodies, they lose mass by colliding with particles that have a mass below 10 GLYPH<0> 2 mt . This could have been predicted by looking at Figure 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 dis-  \nFig. 10. Projectile distribution mass functions for simulation E1, constructed for the maximum masses ( mt = GLYPH<15> ) at various times. Colors and times correspond to Figure 9. For each distribution, the stopping time of the mt -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.  \n<!-- image -->  \ntribution, 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 e GLYPH<11> ective that it stops the growth of the largest bodies, resulting in the behavior seen in Figure 9.  \nWedefine a parameter GLYPH<16> using the positive and negative areas under the projectile mass distributions  \nGLYPH<16> = P C + GLYPH<0> P C GLYPH<0> P C + + P C GLYPH<0> ; (34)  \nwith P C + and P C GLYPH<0> the sums of the positive and negative part of the projectile mass distribution respectively. The parameter GLYPH<16> ranges from 1 (no erosion) to -1 (only erosion), and equals 0 when there is a balance between growth and erosion.  \nThe top panel of Figure 11 shows the evolution of GLYPH<16> for the most massive particle during one of the simulations of Figure 9, plotted as a function of GLYPH<10> ts of the maximum mass. Early on, there is no erosion present and GLYPH<16> = 1, but as the largest bodies grow towards GLYPH<10> ts = 1, erosion increases and GLYPH<16> drops. When 0 < GLYPH<16> < 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 e GLYPH<11> ect. When GLYPH<16> < 0, erosion dominates over growth and the most massive particles are loosing considerable mass. This causes the curve in Figure 11 to turn around. As bodies shrink, there is less erosion and GLYPH<16> increases again. A quasi steady-state is reached with GLYPH<16> just below unity and GLYPH<10> ts ( m max) GLYPH<24> 0 : 6. The reason GLYPH<16> , 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.  \nThe other panels of Figure 11 show similar plots but for different erosion threshold velocities. For v eros = 40 m s GLYPH<0> 1 (middle panel), erosion is less e GLYPH<14> cient and the largest bodies grow to  \nFig. 11. Evolution of GLYPH<16> ( m max) for erosive simulations with v eros = 20 ; 40 ; 60 m s GLYPH<0> 1 as a function of GLYPH<10> ts ( m max), showing the impact of erosion on the ability of the largest bodies to grow. In the upper two panels, the steady-state is indicated by the } -symbol.  \n<!-- image -->  \nGLYPH<10> ts ' 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 e GLYPH<14> ciency is proportional to ( v rel =v eros), the high-velocity projectile are less e GLYPH<14> cient in excavating mass from the targets. Both e GLYPH<11> ects together cause the largest mass in the steady state to be about a factor of 10 larger than in the top panel of Figure 11. Finally, the bottom panel shows the results for v eros = 60 m s GLYPH<0> 1 . This is a special case, since now the erosion threshold velocity can only be reached around GLYPH<10> ts = 1, with radial drift, azimuthal drift, and turbulence contributing (see Figure 2). Indeed, erosion is strongest around GLYPH<10> ts = 1, but is ine GLYPH<14> cient and GLYPH<16> never drops below 0. When GLYPH<10> ts > 20, erosion reappears, as a result of smaller particles drifting into the larger bodies. But, since GLYPH<16> GLYPH<24> 1, bodies can continue to grow relatively una GLYPH<11> ected.", '4.2.1. Variation in porosity': 'One of the biggest advantages of the Monte Carlo method is that aggregate mass and porosity are treated truly independently. In other words, aggregates of identical mass can have a very dif-  \nFig. 12. Masses and filling factors of all unique families at di GLYPH<11> erent 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.  \n<!-- image -->  \nferent 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 impact-parameter 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 Pc = P gas (or Pc = P grav), according to Equations 15 and 16. As a result, mass-porosity relations as shown in Figure 8 accurately represent the internal structure of the majority of aggregates.  \nThis 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. 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, Figure 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 GLYPH<24> 6000 years, the original population of aggregates, whose porosity was set by their growth history, has disappeared. A steady-state 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 GLYPH<10> ts GLYPH<24> 1.', '5. Semi-analytical model': "The evolution of the mass-dominating particles can be captured in a simple semi-analytical model. Assuming the entire dust mass is located in particles of identical mass mp , the growth rate can be written as (Okuzumi et al. 2012)  \ndmp dt = GLYPH<6> d p 2 GLYPH<25> hd GLYPH<27> col v rel : (35)  \nThe collisional cross section depends directly on the particle porosity, and the relative velocity and dust scale height depend on GLYPH<30> 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 GLYPH<24> 2, until the kinetic energy in same-sized collisions exceeds E roll. After that, the internal structure can be calculated through Equation 31 of Okuzumi et al. (2012), but in practice is always dominated by the gas / self-gravity compression of Kataoka et al. (2013a), see Section 2.4.2.  \nThis approach, similar to Kataoka et al. (2014, Section 5.3), is valid when particles grow primarily through collisions with similar-sized 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 ts < t GLYPH<17> .  \nAt the same time, the radial drift of the particles is governed by  \ndR dt = GLYPH<0> v drift ; (36)  \nwith the drift velocity a function of GLYPH<10> ts . Assuming a fixed dust to gas ratio of 10 GLYPH<0> 2 throughout the disk, we can solve Equations 35 and 36 to obtain the evolution of the vertically integrated peak  \nmass. Catastrophic fragmentation is taken into account by setting ( dmp = dt ) = 0 when v turb > v frag for two particles of mass mp . 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. GLYPH<30> = 1 at all times), after 10 6 yr. (For the compact case, we have temporarily set v frag = 10 m s GLYPH<0> 1 .) Initially, growing aggregates are not moving radially, resulting in vertical lines in Figure 13. As 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 same-sized collisions. Particles cannot grow larger than GLYPH<24> cm, and will inevitably drift inwards. In the intermediate region, from 20 GLYPH<0> 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 (cf. Figure 3 of Testi et al. (2014)).  \nThe grey lines in the right-hand panel of Figure 13 show the results of the semi-analytical model for porous growth, where GLYPH<30> is set by collisional, gas pressure, and self-gravity compaction, assuming perfect sticking. Since we are assuming the mp particles carry the total dust mass, we do not have any information about the mass-distribution of smaller particles. Nonetheless, we can mimic the e GLYPH<11> ect of e GLYPH<11> ective erosion, by setting ( dmp = 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 Figure 13 show the results for v eros = 40 m s GLYPH<0> 1 , while the red lines indicate results for the peak mass of the full Monte Carlo models for the same erosion threshold velocity (note that the maximum mass reached in these models can be a factor of GLYPH<24> 10 larger). We have also included a full model run at 200 AU, which we evolved for 10 6 yrs. The results of the semi-analytical model agree with the simulations of the previous section remarkably well.", '6. Discussion': "From the maximum sizes flu GLYPH<11> y aggregates can reach at a given location, we identify three zones in the protoplanetary disk:  \n- -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 out-grow the radial drift barrier, and reach planetesimal sizes without experiencing significant drift. However, when erosion is e GLYPH<14> -cient, mass loss in erosive collisions stalls the growth around GLYPH<10> ts GLYPH<24> 1, preventing the porous aggregates from crossing the radial drift barrier (Figure 11).\n- -10-100 AU: At intermediate radii growth timescales increase and radial drift takes over, even before aggregates reach sizes and stopping times that allow for erosive collisions to take place.\n- -> 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 hit-and-stick regime, and their surface-to-mass 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 surface-to-mass ratio, and the onset of radial drift is already visible for very low particle masses.  \nFor 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 GLYPH<10> ts GLYPH<29> 1, it moves on a Keplerian orbit, and v dg ' GLYPH<17>v K , while bodies with GLYPH<10> ts = 1 have a slightly larger velocity with respect to the gas (Weidenschilling 1977). For the disk model employed in this work (Section 2.1), the quantity GLYPH<17>v K does not depend on R , and thus the maximum drift speed is constant though out the disk. It is clear then from Equation 33 that growing aggregates in colder disks (lower cs ), or disks with a (locally) shallower gas density profile might su GLYPH<11> er less from erosion.  \nIt is clear from Figure 11 that the size of v eros is a key parameter: its value, together with GLYPH<17>v K , determines whether growth beyond GLYPH<10> ts = 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 highly-porous icy bodies in question (Section 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 (Equation 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 . mmand porosities & 10 GLYPH<0> 1 , and cover a soberingly small portion of the parameter space encountered in this work (e.g., Figure 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 re-accreted (Wurm et al. 2004). If e GLYPH<14> cient, 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).  \nSo far we have assumed that while erosion can play an important role, catastrophic fragmentation does not occur. The maximum velocity between same-sized bodies is reached for GLYPH<10> ts = 1, and equals GLYPH<24> (3 = 2) GLYPH<11> 1 = 2 cs (Equation 11). Since the sound speed diminishes for increasing radii, this velocity is highest in the inner disk. For typical turbulence strengths ( GLYPH<11> . 10 GLYPH<0> 3 ) and small icy monomers, this velocity will not exceed the fragmentation threshold velocity (Equation 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 GLYPH<22> m) are removed from the protoplanetary nebula very rapidly, contradicting observational constraints (Dullemond & Dominik 2005; Dominik & Dullemond 2008). Drift-induced erosion can alleviate these issues, since the maximum drift velocity is high throughout the entire disk.  \nIn 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 Monte Carlo model approach like the one presented  \nFig. 13. Evolution of mp ( t ) and R ( t ) for dust coagulation as obtained from the semi-analytical model (Equations 35 and 36), for an MMSN disk and GLYPH<11> = 10 GLYPH<0> 3 . Lines indicate di GLYPH<11> erent starting conditions R ( t = 0), and are evolved for 10 6 yrs. Left: Compact growth: GLYPH<30> = 1 at all times, and v frag = 10 m s GLYPH<0> 1 . Right: Porous growth: the internal structure of the aggregates is set by hit and stick growth, followed by collisional compaction or gas and self-gravity compaction. Grey lines have no erosion, while black lines show the results for v eros = 40 m s GLYPH<0> 1 . Colored lines and } -symbols indicate the evolution and steady state peak mass obtained through local Monte Carlo simulations (Section 4.2).  \n<!-- image -->  \nR  \n[AU]  \nhere, is that it is relatively straightforward to include an additional random number to determine, for example, the impact parameter. The di GLYPH<14> culty lies in obtaining a collision model that describes the collisional outcome as a function of said parameter. A good start would be the work of Wada et al. (2013), who show the growth e GLYPH<14> ciency as a function of impact parameter (Figure 4). Basically, head-on 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.  \nAt the heart of the model of Section 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 for 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 self-gravity 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 Monte Carlo method developed in this paper would be well suited for such studies, since adding parameters describing the aggregates is relatively straightforward.", '6.1.1. Pebble accretion': 'A novel idea in the field of planet formation is the process of pebble accretion, where protoplanets grow very e GLYPH<14> ciently by accreting small pebbles (Ormel & Klahr 2010; Lambrechts & Johansen 2012, 2014; Kretke & Levison 2014). These models  \nFig. 14. Same as the right-hand plot of Figure 13, but with aggregate size on the vertical axis.  \n<!-- image -->  \nrely 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 (Figure 13). While the Stokes numbers of these particles are similar to the drifting pebbles in the compact case, their masses, sizes, and porosities can di GLYPH<11> er by many orders of magnitude (see also Figure 14). In addition, the drag regime that the drifting bodies experience di GLYPH<11> ers from the compact case (Figure 1). Future studies are needed to address the e GLYPH<11> ect of these factors on the e GLYPH<14> ciency of pebble accretion.', '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  \n10  \n10  \n15  \n12  \n9  \n10  \n6  \n10  \n3  \n10  \n0  \n10  \n10  \n10  \n10  \n10  \n-3  \n-6  \n-9  \n-12  \n]  \ng  \n[  \nm  \n1  \n10  \nPorous growth  \nInefficient erosion  \nEfficient erosion  \nMC Models  \n2  \n10  \nmajority of mass needs to reside in particles with high Stokes numbers; the mid plane dust to gas ratio has to be close to unity; and the local vertically integrated dust-to-gas ratio needs to exceed GLYPH<24> 0 : 03 (Drazkowska & Dullemond 2014). The first two conditions can be studied with simulations like the ones presented in this work. For example, for the steady-state distribution reached for v eros = 40 m s GLYPH<0> 1 at 5 AU for GLYPH<11> = 10 GLYPH<0> 3 , approximately 50% of the dust mass resides in particles with GLYPH<10> ts > 10 GLYPH<0> 2 , and the mid plane dust-to-gas ratio is GLYPH<24> 10 GLYPH<0> 1 . For weaker turbulence, the mid plane dust-to-gas ratio will be increased further, since hd GLYPH<24> GLYPH<11> 1 = 2 (Equation 6). Because our simulations are local, the vertically integrated dust-to-gas ratio stays constant at 10 GLYPH<0> 2 . To fulfill the third condition, the dust-to-gas 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, drift-induced erosion appears to be a robust way of concentrating mass around GLYPH<10> ts GLYPH<24> 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 e GLYPH<11> ective in stopping growth at mm-cm 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 di GLYPH<11> erential and Monte Carlo methods (Dra¸˙zkowska et al. 2014). The distribution method used in this work, as outlined in Section 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 well-suited for studying the breakthrough case.', '6.1.4. Opacities of porous grains': 'The optical properties of dust distributions resulting from porous growth are very di GLYPH<11> erent from populations containing exclusively solid particles. Not only are the mass distributions themselves di GLYPH<11> erent (e.g. Figure 13), but the scattering and absorption opacities of the individual grains are a GLYPH<11> ected significantly by porosity (Kataoka et al. 2014; Cuzzi et al. 2014). For simple dust mass distributions, the e GLYPH<11> ect of grain porosity on the appearance of protoplanetary disks has been investigated by Kirchschlager &Wolf (2014). Combining self-consistent coagulation models including erosion and fragmentation - with porosity-dependent dust opacities will reveal the full impact porous growth has on the appearance of protoplanetary disks.', '7. Conclusions': 'Porous growth is very di GLYPH<11> erent from compact growth (Figure 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 di GLYPH<11> er greatly from those of compact particles (Figure 1), causing di GLYPH<11> erences in relative velocities (Figure 2), vertical settling, and radial drift.  \nWe 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.  \nWe consider compaction through collisions, gas pressure, and self-gravity (Figure 8), and include a physical model for erosive collisions (Sections 2.3.2 and 2.3.3). The main findings of this work are:  \n- 1. Porous icy aggregates can outgrow the radial drift barrier in the inner GLYPH<24> 10 AU, despite increased growth timescales resulting from gas- and self-gravity compaction, if the perfect sticking assumption holds (Figures 7 and 13). This is in agreement with Okuzumi et al. (2012) and Kataoka et al. (2013b).\n- 2. While the maximum collision velocity between similar particles ( GLYPH<24> GLYPH<11> 1 = 2 cs ) typically does not exceed the critical fragmentation threshold velocity for icy bodies, the velocity between drifting aggregates (with GLYPH<10> ts GLYPH<21> 1) and smaller bodies is much larger ( GLYPH<24> GLYPH<17>v K), and can exceed the critical threshold velocity for erosion (Figure 2).\n- 3. In these cases, we find that the mass loss through erosive collisions can balance the growth through same-size collisions, halting the growth of the largest bodies (Figures 10 and 11). In our local simulations, this results in a steadystate where the largest bodies have GLYPH<10> ts GLYPH<24> 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 (Figures 9 and 12). Only for the highest erosion threshold velocity we considered ( v eros = 60 m s GLYPH<0> 1 ), do the aggregates with GLYPH<10> ts GLYPH<24> 1 manage to gain mass and grow through the drift barrier.\n- 4. A simple semi-analytical model (Section 5) accurately describes the growth and drift behavior of the mass-dominating 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, or dust-togas ratio; and aggregate properties such as monomer size and erosion / fragmentation threshold velocities.  \nAcknowledgements. 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¸˙zkowska and T. Birnstiel for useful discussions, and S. Okuzumi for comments and helping with the comparison to his work. The authors thank S. 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. 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.', 'References': 'Bai, X.-N. & Stone, J. M. 2010a, ApJ, 722, 1437 Bai, X.-N. & Stone, J. M. 2010b, ApJ, 722, L220 Birnstiel, T., Dullemond, C. P., & Brauer, F. 2010, A&A, 513, A79 Blum, J. & Wurm, G. 2008, ARA&A, 46, 21 Cuzzi, J. N., Estrada, P. R., & Davis, S. S. 2014, ApJS, 210, 21 Dominik, C. & Dullemond, C. P. 2008, A&A, 491, 663 Dominik, C. & Tielens, A. G. G. M. 1995, Philosophical Magazine, Part A, 72, 783 Dominik, C. & Tielens, A. G. G. M. 1997, ApJ, 480, 647 Drazkowska, J. & Dullemond, C. P. 2014, ArXiv e-prints Dra¸˙zkowska, J., Windmark, F., & Dullemond, C. P. 2014, A&A, 567, A38 Dullemond, C. P. & Dominik, C. 2005, A&A, 434, 971 Garaud, P., Meru, F., Galvagni, M., & Olczak, C. 2013, ApJ, 764, 146 Gillespie, D. T. 1975, Journal of Atmospheric Sciences, 32, 1977 Gundlach, B., Kilias, S., Beitz, E., & Blum, J. 2011, Icarus, 214, 717 Gundlach, S. & Blum, J. 2014, Accepted for publication in ApJ Güttler, C., Blum, J., Zsom, A., Ormel, C. W., & Dullemond, C. P. 2010, A&A, 513, A56 Haisch, Jr., K. E., Lada, E. A., & Lada, C. J. 2001, ApJ, 553, L153  \nHayashi, C. 1981, Progress of Theoretical Physics Supplement, 70, 35'}