A&A 513, A57 (2010)

The outcome of protoplanetary dust growth: pebbles, boulders, or planetesimals?^,

II. Introducing the bouncing barrier

A. Zsom¹ - C. W. Ormel¹ - C. Güttler² - J. Blum² - C. P. Dullemond¹

1 - Max-Planck-Institute für Astronomie, Königstuhl 17, 69117 Heidelberg, Germany
2 - Institut für Geophysik und extraterrestrische Physik, Technische Universität Braunschweig, Mendelssohnstr. 3, 38106 Braunschweig, Germany

Received 24 July 2009 / Accepted 4 January 2010

Abstract
Context. The sticking of micron-sized dust particles caused by surface forces within circumstellar disks is the first stage in the production of asteroids and planets. The key components describing this process are the relative velocity between the dust particles in this environment and the complex physics of dust aggregate collisions.
Aims. We present the results of a collision model based on laboratory experiments of these aggregates. We investigate the maximum aggregate size and mass that can be reached by coagulation in protoplanetary disks.
Methods. We use the results of laboratory experiments to establish the collision model previously published by Güttler et al. The collision model is based on the assumptions that we model the aggregates as spheres with compact and porous ``phases'' and that there is a continuous transition between these two. We apply this collision model to the Monte Carlo method developed previously by Zsom & Dullemond and include Brownian motion, radial drift, and turbulence as contributors of relative velocity between dust particles.
Results. We model the growth of dust aggregates at 1 AU in the midplane for three different gas densities. We find that the evolution of the dust does not follow the previously assumed growth-fragmentation cycles. Catastrophic fragmentation hardly occurs in the three disk models. Furthermore, we see long-lived, quasi-steady states in the distribution function of the aggregates caused by bouncing. We explore how the mass and the porosity depend on both the turbulence parameter and the critical mass ratio of dust particles. Upon varying the turbulence parameter, the system behaves in a non-linear way, and we find that the critical mass ratio has a strong effect on the particle sizes and masses. Particles reach Stokes numbers of roughly 10^-4 during the simulations.
Conclusions. The particle growth is stopped by bouncing rather than fragmentation in these models. The final Stokes number of the aggregates is rather insensitive to the variations in the gas density and the strength of turbulence. The maximum mass of the particles is limited to 1 g (chondrule-sized particles). Planetesimal formation can proceed by the means of the turbulent concentration of these aerodynamically size-sorted, chondrule-sized particles.

Key words: planets and satellites: formation - accretion, accretion disks - methods: numerical

1 Introduction

In the core accretion paradigm of planet formation (Pollack et al. 1996; Mizuno 1980), planets are the outcome of an accretion process that starts with micron-size dust grains and covers 40 mag in mass. The paradigm can be divided into three stages. The first stage involves the formation of rocky planets and the rocky cores of gas giant planets and begins with the coagulation of dust in the protoplanetary disks surrounding many pre-main-sequence stars (Safronov 1969; Blum & Wurm 2008; Weidenschilling & Cuzzi 1993). The next stage of planet formation is the formation of protoplanetary cores from the planetesimals. The idea is that the kilometer-size planetesimals are so large that gravity begins to dominate leading to the gravitational agglomeration of these bodies to rocky planets. This scenario was studied by Safronov (1969) and modeled using numerical methods by Weidenschilling (1980), Nakagawa et al. (1983), Mizuno et al. (1988), Schmitt et al. (1997), Wetherill (1990), Nomura & Nakagawa (2006), Garaud & Lin (2004), Tanaka et al. (2005), and several additional authors. These models solve for the size distribution of dust aggregates in the disk as a function of time, and investigate if, where, and how larger dusty bodies form, and how long this takes. Finally, in the third stage, gas accretes onto these protoplanets forming giant planets or - in the absence of gas - gravitational encounters occur between these protoplanets result in a chaotic, giant impact phase, until orbital stability is achieved (Kokubo et al. 2006; Thommes et al. 2008; Chambers 2001).

In this study, we focus on the first phase and study the effectiveness of the dust growth by surface force, that is, how large particles become due to simple sticking processes only. It is known that initially, for micron-size grains, the growth is driven by Brownian motion. This typically leads to slow collisions and forms aggregates of fractal structure (Kempf et al. 1999; Blum et al. 1996). In the current picture of dust growth, as these aggregates grow, at some point the growth will leave the fractal regime, and collisions will start to lead to the compaction and breaking of the aggregates (Blum & Wurm 2000), embedding of small bodies into larger aggregates (leading to ``filling up'' of these larger aggregates and compaction caused by the force of the collision (Ormel et al. 2007). As the size of the dust aggregates increases, differential vertical settling (Safronov 1969), radial drift (Whipple 1972), and turbulence (Mizuno et al. 1988; Ormel & Cuzzi 2007; Völk et al. 1980) will become important new mechanisms driving relative velocities between aggregates. The increasing relative velocities caused by these mechanisms will at least partly compensate the lower collision probability due to lower surface-over-mass ratio of large aggregates. When the aggregates grow to sizes of between millimeter and meter, however, the sticking efficiency drops strongly (e.g., Blum & Münch 1993) and the relative velocities become so large that aggregates can fragment (Blum & Wurm 2008, so-called ``fragmentation barrier''). Another hurdle that the particles have to circumvent is the ``drift barrier'' (Weidenschilling 1977a), namely that millimeter or centimeter-sized particles are lost to the star because of radial drift that occurs on a short timescale. Okuzumi (2009) pointed out the existence of a ``charge barrier'', which possibly halts the particle growth at an early stage of fractal aggregates. Despite many years of efforts, it is unknown whether the coagulation process can overcome these barriers. These barriers have been and remain the main open question about the initial stages of planet formation, i.e., the growth from dust to planetesimals.

Several mechanisms have been proposed to overcome this problem, among which are the trapping of dust in vortices (Barge & Sommeria 1995; Klahr & Henning 1997; Lyra et al. 2009), the trapping of decimeter-sized boulders in turbulent eddies and the subsequent gravitational collapse of swarms of these trapped boulders (Johansen et al. 2007), the trapping of particles in a pressure bump caused by the evaporation front of water (Kretke & Lin 2007; Brauer et al. 2008b) and many other scenarios. However, the correct modeling of any of these requires detailed knowledge of the collisional physics, and these models have so far relied on either simplified input phyisics or simplified initial conditions.

Because of their complexity, collisional evolution models have to make simplifying assumptions about the outcome of dust aggregate collisions, for example that collisions always result in sticking, or otherwise use simple recipes for the collisional outcome. Ideally, one requires detailed knowledge of the outcome of every collision. But modeling this microphysics within an evolution model is simply impractical. There are computer programs that model these individual collisions in detail (e.g., Dominik & Tielens 1997; Suyama et al. 2008; Geretshauser et al., in prep.), but each collision model takes anywhere between hours and weeks to complete on a computer. They are therefore impractical to use at run-time in a model that computes the overall time-dependent evolution of the dust size distribution inside protoplanetary disks. These collision models themselves often depend on poorly known input physics.

Another approach to obtaining the collisional outcome of dust aggregates is to model these collisions in the laboratory. From the many experiments that have been performed, a picture emerges of the outcome of dust aggregate collision under a variety of conditions in the protoplanetary disk (PPD). In a companion paper (Güttler et al. 2009, henceforth Paper I), we collected data from over 19 experiments, and compiled a set of formulae to describe reasonably well the outcomes of these collisions in such a way that they can be used as input to models that address the temporal evolution in the dust size distribution.

In this paper, we directly rely on the outcome of these laboratory experiments to model the dust aggregate size distribution. As described in Paper I, we have produced a mapping of all available collision experiments for silicate-like particles. The velocity range of these experiments is also sufficiently wide to cover various disk models that roughly correspond to the conditions at 1 AU in the PPD. For details of the collisional mapping, we refer to Paper I but summarize the elements of our new collision model in Sect. 3.1.

We build this collision kernel into a Monte Carlo code for modeling the size and porosity distribution of dust in a protoplanetary disk (Zsom & Dullemond 2008, hereafter ZsD08). The outcome of our laboratory-driven dust coagulation model is difficult to predict a priori since the key variables involved depend on a non-trivial interplay between the collision kernel (Paper I) and the velocity field. We can, however, propose two scenarios. In the first, particle growth proceeds beyond the meter-size barrier, all the way to forming planetesimals. In the second scenario, growth terminates at an intermediate size. In this case, additional growth to planetesimal sizes may proceed by the means of the concentration and subsequent gravitational collapse of these particles (Johansen et al. 2007; Cuzzi et al. 2008). Thus, our model providea the starting conditions for these concentration models. We emphasize, however, that in this work we do not in any way ``optimize'' the outcome by laboriously scanning all the parameter space or treating environments that may be more conducive to growth, such as nebula pressure bumps or the trapping of dust in vortices (Kretke & Lin 2007; Lyra et al. 2009). These are obvious extensions of this present work. But by considering the sensitivity of a few key parameters (e.g., gas density, and turbulence strength) to the outcome of the growth process, we obtain a picture of where the arrow of coagulation typically points to in protoplanetary environments: pebbles, boulders, or planetesimals.

In this paper, we describe the three nebulae models used in this work and the sources of relative velocity between the aggregates (Sect. 2), how we developed our coagulation/fragmentation model of Paper I into a Monte Carlo code (Sect. 3), and our first results (Sect. 4). We also test the sensitivity of these results to variations in gas density, the velocity field, and other key model parameters. Section 5 reflects the importance of our result to planetesimal formation and provides suggestions for future experiments. Finally, Sect. 6 lists our main conclusions.

2 The nebulae model

2.1 Disk models

We briefly describe the disk models considered in this paper.

The low density model:

Resolved millimeter emission maps of protoplanetary disks seem to indicate a shallow surface density profile (Andrews & Williams 2007) given by . The systematic effects of some of their assumptions, such as the disk inclinations or the simplified treatment of the temperature distribution, may produce steeper profiles. Therefore, Brauer et al. (2008a) adopted the following profile:

Here we assumed that the central star is of solar mass, the disk extends from 0.03 AU to 150 AU, and the total mass of the disk is 0.01 . Assuming that the pressure scale-height is and the vertical structure is Gaussian, we obtain:

the density at 1 AU in the midplane (z=0) being  g cm^-3, approximately two orders of magnitude lower than the minimum mass solar nebulae (MMSN) value.

MMSN model:

The minimum mass solar nebulae model (MMSN) was introduced by Weidenschilling (1977b) and Hayashi et al. (1985). From the present state of the Solar System today, it is possible to infer a lower limit to the mass in the solar nebulae from which the planets were formed. The model assumes that the planets were formed where they are currently located (no migration included). It also assumes that all the solid material presented in the solar nebula had been incorporated into the planets. The loss of solid material because of radial drift is not taken into account. Despite these uncertainties, the MMSN model is frequently used as a benchmark. The surface density of the MMSN disk is given by

which corresponds to a total disk mass of 0.01  between 0.4 and 30 AU (between the orbits of Mercury and Neptune). Assuming that the vertical structure of the gas follows a Gaussian distribution, we infer a midplane density at 1 AU of  g cm^-3.

The high density model:

Desch (2007) introduced a ``revised MMSN model'' by adopting the starting positions of the planets in the Nice model of planetary dynamics (Tsiganis et al. 2005) thus taking into account planetary migration. The revised MMSN model predicts that the Solar System was initially in a far more compact configuration and its surface density profile is given by

This model is consistent with that of a decretion disk that is being photoevaporated by the central star. Although the model of Desch (2007) was defined for the outer Solar System, we extrapolate the profile to 1 AU to cover a broad range of surface density values in our calculations. Assuming, as in the MMSN model, a Gaussian vertical distribution, the density at 1 AU in the midplane is  g cm^-3.

For simplicity, we adopt a midplane temperature of 200 K (isothermal sound speed of  cm s^-1) in all three models.

2.2 Relative velocities

We consider three contributors to the relative velocities between dust aggregates: Brownian motion, radial drift and turbulence. In the following, we discuss each of these sources.

The average relative velocity of two particles with mass m₁ and m₂ in a region of a disk with temperature T due to Brownian motion is

For micron-sized particles, the relative velocity is of the order of 0.1 cm s^-1, but for cm-sized particles this value drops several orders of magnitude. Therefore, Brownian motion is only effective for collisions between small particles during the initial stages of growth. Coagulation caused by Brownian motion results in fluffy aggregates of fractal dimensions of around 2 and 3 (Kempf et al. 1999; Blum et al. 1996; Krause & Blum 2004; Blum et al. 2000). In practice, no growth is caused by Brownian motion for aggregates larger than 100 micron.

The second contributor to the relative velocity is turbulence. The relative velocity of aggregates produced by the random motion of turbulent eddies was calculated numerically by Völk et al. (1980), Mizuno et al. (1988), and Markiewicz et al. (1991). We use the closed form expressions presented by Ormel & Cuzzi (2007). We assume that turbulence is parameterized by the Shakura & Sunyaev (1973)  parameter

where is the turbulent viscosity, is the isothermal sound speed, and is the pressure scale height of the disk. The value of the parameter reflects the strength of the turbulence in the disk. Typical values of  in this paper range between 10^-3 and 10^-5. The turbulent relative velocity is a function of the stopping times of the two colliding particles. The stopping time (or friction time) is the time the particle needs to react to the changes in the motion of the surrounding gas. As long as the radius of the particle is smaller than the mean free path of the gas ( ), the particle is in the Epstein regime, where the stopping time is (Epstein 1924):

where m and A are the mass and the cross-section of the particle, and and are the gas density and the thermal velocity. At high gas densities where the mean free path is low or in the case of larger particles, the first Stokes regime applies and the stopping time is

In the first Stokes regime, the stopping time is independent of the gas particle relative velocity as well as the gas density. This regime can be used as long as the particle Reynolds number is smaller than unity. The particle Reynolds number is calculated to be (Weidenschilling 1977a)

where is the relative velocity between the particle and the gas, and is the gas viscosity. For particles outside the Epstein regime, we can assume that the systematic velocity (radial drift) dominates over the random velocities (turbulence); therefore, , where is the drift velocity of the particle, defined in the next paragraph. The particle Reynolds number never exceeds unity in our simulations. Therefore, we do not include additional Stokes regimes.

Radial drift also contributes to relative velocities between aggregates. Radial drift () has two sources: the drift of individual particles () and drift caused by accretion processes of the gas ( ), thus the total radial drift velocity is .

The radial drift of individual dust aggregates with mass m is (Weidenschilling 1977a)

where St is the Stokes number of the aggregate ( , where is the orbital frequency) and v_N is the maximum radial drift velocity (Whipple 1972).

The second contribution to the radial velocity is produced by the accretion of the gas. This part of the radial velocity is calculated as follows (Kornet et al. 2001)

where is the accretion velocity of the gas (Takeuchi & Lin 2002).

The relative velocity due to radial drift is then simply the difference between the radial velocity of particles 1 and 2. However, as the Stokes number of the aggregates is always smaller than 10^-3 (see Sect. 4), the second term of the radial velocity ( ) can be safely neglected

This study uses two quantities to describe the porosity of the aggregates. The volume filling factor is

where V^* is the volume occupied by the monomers and is the total volume of the aggregate, including pores, and A and A^* are the surface area equivalents of these quantities. In this way, the filling factors also enters the definition of the friction time (Eqs. (7) and (8)). The density of aggregates then follows as , where  g cm^-3 is the material density of the silicate. In this study, we also use the reciprocal parameter of the filling factor, which is denoted by the enlargement parameter, .

We illustrate the relative velocity between equal-sized and different-sized aggregates with ( ) in Fig. 1 for the disk models considered in this work. Adopting a threshold (fragmentation) velocity of 1 m s^-1, the maximum particle size that can be reached in the models are 0.025 cm in the low density model, 1.4 cm in the MMSN model, and 1.7 cm in the Desch model. The Stokes numbers of these particles are identical in all three models, . The constant fragmentation velocity of 1 m s^-1 is the typical velocity at which silicate particles fragment. In our collision model this is not the case for all combinations of mass ratio and porosity (Paper I), but the m s^-1 threshold remains a useful proxy for the point where fragmentation processes become important.

Figure 2 shows the Stokes number as a function of particle radii in the three models. Initially, particles are in the Epstein regime, where the stopping time, thus the Stokes number, depends on the gas density. When the particles enter the Stokes regime, the stopping time becomes independent of the gas density (see Eq. (8)). One can see that particles in the Desch model are in the Stokes regime at a Stokes number of (when the particles have relative velocities of 1 m s^-1), while the aggregates in the MMSN model are close to it, which explains why the maximum particle size is almost the same in these two models.

As discussed in Ormel & Cuzzi (2007), particles are initially in the ``tightly coupled particle'' regime, where the eddies are all of class I type, meaning that the turnover time of all eddies is longer than the friction time of the particles (Völk et al. 1980). Upon entering a class I eddy, a particle therefore forgets its initial motion and aligns itself with the gas motions of the eddy before the eddy decays or the particle leaves it. This regime is presented in Figs. 1a and b. Different-sized particles are found in this relative velocity regime as long as their masses are less than 10<sup>-8</sup> g in the low density model, 10<sup>-3</sup> g in the MMSN model, and 10<sup>-2</sup> g in the high density model assuming fluffy particles (). If the particles leave this regime and enter the ``intermediate particle'' regime, their relative velocity increases. This transition affects the particle evolution, as discussed in e.g., Sect. 4.3.

3 Collision model and implementation

We use a statistical or ``particle in a box'' method to compute the collisional evolution, that is, we assume that all particles are homogeneously distributed within a certain volume (the simulation volume). In reality, the particles could however leave the simulated volume or new particles could enter from outside because of radial drift or random motions (turbulence and Brownian motion). Since we do not resolve the spatial dependence of the aggregates, we simply assume that local conditions hold during the run. The gas and dust densities are kept constant and particles cannot leave or enter the simulation volume (hereafter ``local approach'').

3.1 Short overview of the collision model

Many laboratory experiments on dust aggregate collisions have been performed (see Blum & Wurm 2008). The growth begins as fractal growth and we use the recipe of Ormel et al. (2007) to describe this initial stage. However, once aggregates have restructured into non-fractal, macroscopic aggregates (e.g., m in size), laboratory experiments show that the collisional outcomes become very diverse. In this regime, many new experiments were performed with dust aggregates consisting of 1.5 m diameter SiO₂ monomers of either high porosity (Blum & Schräpler 2004), or intermediate porosity (). Paper I compiled 19 experiments with different aggregate masses, collision velocities, and aggregate porosities.

From these experiments, we identified nine different collisional outcomes involving sticking, bouncing, or fragmentation (see Fig. 3). The occurrence of these regimes depends mainly on aggregate masses and collision velocities. However, it also depends on both the porosity of the particles and the critical mass ratio. For example, Paper I finds fragmentation in collisions between a porous aggregate and a solid wall, whereas Langkowski et al. (2008) find sticking of a porous projectile by penetrating a target that is also porous. Likewise, Heißelmann et al., in prep. detect the bouncing of two similar-sized, porous dust aggregates, while Langkowski et al. (2008) uncover sticking of the same velocity where one collision partner (target) is significantly bigger. To address the importance of the mass ratio and porosity, we identified eight different collision regimes (look-up tables) based on a binary treatment of porosity and mass ratio i.e., (i) similar-sized or different-sized collision partners; and (ii) porous or compact collision partners. The additional distinction between target, which we always define as the heavier collision partner, and projectile then defines eight different collision regimes. We denote these regimes as ``pP'' (porous projectile, porous target; target significantly bigger than the projectile) and ¸ (porous projectile, compact target; target of similar size than the projectile).

In Paper I, we classified each of these 19 experiments into one or more of these eight regimes (see Fig. 10 of Paper I). Based on extrapolation of experimental findings, we decide in which mass and velocity range collisions result in sticking, bouncing, or fragmentation. These results are presented in Fig. 11 in Paper I.

It should be noted that the critical mass ratio of the equal-size (e.g., ``pp'', ``cc'') to the different-size regimes (e.g., ``pP'', ``cC'') is ill-constrained by experiments. Therefore, we use critical mass ratios of , 100, and 1000 to explore the effect of this parameter.

3.2 Porosity

Paper I defined a binary representation of the porosity, particles are either porous or compact. Following the simple model of Weidling et al. (2009), we include a continuous transition between these two ``phases''. These authors showed that the compaction of particles caused by bouncing can be described by porous and compacted sites on the surface of the aggregate. A site of the aggregate is porous if it has not yet experienced any collisions (e.g., bouncing); a compacted site has encountered at least one collision but any additional collision at that part of the surface cannot change the porosity of the site. We define the probability of hitting a passive site of the aggregate with the experssion

where is the volume-filling factor of the aggregate, is the critical porosity ( , see Paper I), and is the volume filling factor of the porous site, which is chosen to be 0.15. If is between 0.15 and 0.4, a random number decides whether the particle collided with a porous or a compact site. This treatment of the porosity ensures a continuous transition from porous to compact aggregates.

During the initial hit & stick (S1) phase, particles are in the fractal growth regime (Ossenkopf 1993; Blum et al. 2000; Krause & Blum 2004). Particles grow initially because of Brownian motion and later due to turbulence. The structure of the aggregate depends on whether the collision happened between a cluster and a monomer (PCA) or between two clusters (CCA). The latter type of collision produces fluffy aggregates with fractal dimensions of 2, while the former leads to more compact structures of fractal dimension 3. The hit and stick recipe of Ormel et al. (2007) attempts to interpolate between these two fractal models. At one point, collisional energies become high enough to invalidate this assumption. This occurs when the collisional energy is five times higher than the rolling energy of monomers (Dominik & Tielens 1997; Blum & Wurm 2000). The internal structure then becomes homogenous. In our model we assume that once the fractal growth caused by S1 is over, bouncing with compaction (B1) restructures the aggregates producing compact structures of fractal dimension 3. Since our model cannot follow the exact shape of the particles, we assume that the aggregates are spheres and can be described with a single density () thus neglecting the effects of e.g., the ``toothing radius'' of Ossenkopf (1993) or the craters forming during penetration (S3).

3.3 The Monte Carlo method

Using the expressions for the relative velocity, the collisional cross-section between the dust particles, and the collisional outcome, we solve for the temporal evolution of the dust size distribution. Traditionally, the Smoluchowski equation is used to determine the evolution of the mass distribution function (e.g., Dullemond & Dominik 2004; Dullemond & Dominik 2005; Tanaka et al. 2005; Brauer et al. 2008a). The continuous form of the Smoluchowski equation used in these works lacks the stochasticity of the coagulation problem (Safronov 1969). All bodies of mass m will grow in the same way thus the spatial and temporal fluctuations of the particle ensemble are averaged out. In reality, however, particles with similar masses might follow a different evolutionary path depending on which other particles they collide with. The collision model typically used in these works is, by necessity, rather simple because in the Smoluchowski formulation the collision and time evolution steps are linked together. These collision models consist of sticking and fragmentation and only the mass of the particles is followed. The advantage of such a model is that it is not computationally too expensive: the entire disk can be practically modeled.

Ormel et al. (2007) introduced a new Monte Carlo method to solve for the mass and the porosity distribution function simultaneously. Their collision model consists of sticking and compaction; ZsD08 also added a simple fragmentation model. Although these models are more detailed, one can see that they still lack the full complexity observed in ``the zoo'' of laboratory collision experiments.

The MC-approach used in this study was previously presented by ZsD08. It can be characterized by two key properties: (1) the number of MC-particles (also referred to as representative particles) is kept constant; (2) the method follows the mass of the particle distribution.

Property (1) is required to ensure statistics. Because of the noise of MC-methods, a large fluctuation in N would severely affect the accuracy of the method (Ormel & Spaans 2008). The second property states that our primary interest lies in the particles that contain most of the mass of the system. Moreover, it has been shown that following the particle's mass distribution - rather than its number distribution - is also a prerequisite to preserving a close correspondence with systems that experience strong growth (Ormel & Spaans 2008).

Property (2) ensures that the MC method samples only the parameter space where a signi�cant portion of the total dust mass is. However, this is not always desirable. For instance, radiative transfer calculations require as input the surface area distribution of the aggregates, which determines the opacity. If most of the particle mass is contained in large particles (which are not observable), the number of small particles (which could contain most of the surface area and determines the IR appearance of the disk) might be resolved with poor statistics. But if we are interested in following the evolution of the dominant portion of the dust, then MC methods naturally focus on these parts of phase space.

A necessary condition for the ZsD08 method to work is that the number of the representative particles N is much less than the number of actual aggregates present in the system being considered - a condition that is safely met in any of our simulation runs. A representative particle will then collide only with the non-representative particles, whose distribution is assumed to be identical to that of the representative particles. We refer to ZsD08 for details of the precise implementation and accuracy of the method; here we concentrate further on how the method operates in the new collisional setup.

The collision kernel is defined as the product of the cross-section of the colliding particles and their relative velocity:

where the index i corresponds to the representative particle and k is the index of the non-representative particle. The kernel is proportional to the probability of a collision. The value of K_i,k is calculated for every possible particle pair, and random numbers determine which of the collisions occur first and at which time interval.

The above properties and conditions specify the essence of the ZsD08 method: one of the two collision particles is a representative particle and, because of property (1), only one of the collisional products becomes the new representative particle. Because of property (2), the choice for the new representative particle is weighed by the mass of the collision products. A very helpful analogy here is that of the representative ``atom'', which is contained within the representative particle. The choice of new representative particle after the collision is then proportional to the probability of the representative ``atom'' ending up in the collision products. If, for instance, a collision leads to the production of an entire distribution of debris particles, the probability that a particular debris fragment becomes the new representative particle is proportional to the likelihood of this fragment containing the representative ``atom''.

3.4 Implementation of the collision types

We describe the implementation of the collision model using the representative ``atom'' concept. We refer to Paper I for details of the various collision types mentioned below.

Hit & stick (S1), sticking through surface effects (S2), penetration (S3):

All three of these collision types cause sticking and increase the mass of the aggregate by that of the projectile, but the porosity changes in a different manner (see Paper I). The new mass of the representative particle i is then the sum of the original particle masses, , where m_i is the mass of the representative particle and m_k is the mass of the non-representative particle.

Mass transfer (S4):

In the case of mass transfer (S4), a certain percentage of the mass of the projectile sticks to the target, while the leftover mass of the projectile fragments according to a power-law distribution (see Paper I).

There are two situations to consider:

1.

The representative ``atom'' is part of the target. The mass of the new aggregate will be the mass of the original aggregate plus the transferred mass from the non-representative particle ( , where is the transferred mass calculated according to Paper I).

2.

The representative ``atom'' is part of the projectile. We again have two situations.

(a)

The representative ``atom'' will be transferred to the non-representative particle. The mass of the new representative particle will be the mass of the non-representative particle plus the transferred material ( ). The probability of transferring (removing) the representative atom from the projectile is simply , the ratio of the transferred mass to the mass of the projectile.

(b)

The representative ``atom'' remains in one of the fragments. The probability of this event is , the ratio of the fragmented mass to the original mass of the representative particle. As discussed in Paper I, the fragments follow a power-law mass distribution. The distribution is defined by the maximum mass of the fragments, which is a function of the relative velocity and the total mass of the fragments. The total mass of the fragments is . We randomly choose from the fragment distribution to determine the new mass of the representative particle (to find which of the fragments will contain the representative ``atom'').

Bouncing with compaction (B1):

In this process, particles collide and bounce. Bouncing itself does not change the mass of the particles, but it compactifies them according to Paper I. As observed in laboratory experiments (Weidling et al. 2009), there is a small probability ( ) that the bouncing particle will break apart. If this happens, we break the particle into two equal-mass pieces.

Bouncing with mass transfer (B2):

From the implementation point of view, this is similar to mass transfer (S4). The recipe to define the new representative particle is as in mass transfer (S4). The difference is that the projectile does not fragment during the collision, and that the porosity changes differently (see Paper I).

Fragmentation (F1):

Fragmentation only happens between similar-sized aggregates in the ``pp'' and ``cc'' regimes. The fragments follow a power-law mass distribution, where the maximum mass of the fragments is determined by the relative velocity of the particles and the total mass that goes into the fragments (Paper I). We randomly choose from this distribution to determine the new mass of the representative particle.

Erosion (F2):

Erosion (F2) happens between only different-sized particles. During the collision, the projectile ``kicks out'' pieces from the target aggregate. These pieces follow a power-law distribution (see Paper I). We have to consider two cases:

1.

The representative ``atom'' is in the target. Again, we have two possibilities.

(a)

The representative ``atom'' remains in the target after the collision. The mass of the new particle will be , where is the eroded mass. The probability of this event is , which is the ratio of the left-over mass (which does not erode) to the mass of the original particle.

(b)

The representative ``atom'' originates in the eroded particles. Since the eroded particles follow a power-law distribution, we randomly draw from this distribution to determine the new mass of the representative particle. The likelihood of this event is the ratio of the eroded mass to the original mass of the particle ( ).

2.

The representative ``atom'' is part of the small particle that caused the erosion. As the particles do not stick and the small particle does not fragment, the representative particle remains unaffected.

Fragmentation with mass transfer (F3):

The porous particle becomes destroyed by the compact one and transfers a certain amount of mass to the compact particle. fragmentation with mass transfer (F3) only happens in the ``cp'' regime. We again have two possibilities.

1.

The representative ``atom'' is part of the compact particle. In this case, the representative ``atom'' cannot leave the particle. The new mass of the representative particle will be , the sum of the original mass plus the transferred mass.

2.

The representative ``atom'' was part of the porous aggregate.

(a)

The representative ``atom'' is part of the material that is transferred to the compact particle. In this case, the new mass of the particle will be that of the compact (non-representative) particle plus the transferred material ( ). The probability of this event is .

(b)

The representative ``atom'' is among the fragments. As before, the mass distribution follows a power-law and we draw randomly from this distribution to determine the new mass of the representative particle. The probability of this event is .

3.5 Evolving the particle properties in time

We summarize how the particle properties evolve in time using the aforementioned kernel. We begin with the size and porosity distribution of the particles at a given time, t. At t=0, we provide the initial size and porosity distribution (see Sect. 4.1). Knowing these:

• We calculate both the cross-sections of all possible collision partners and their relative velocities using the equations described in Sect. 2.2. Both sets of quantities are used to determine the collision rates between the particle pairs.
• By using random numbers, we identify from the collision rates the representative particle involved in the collision, the non-representative particle it collides with, and the time at which the collision takes place ( ).
• Knowing the masses (mass ratio) and porosities of the collision partners, we identify in which of the eight regimes the collision takes place (e.g., ``pP'', or ``pC'', etc.).
• Next, we identify which of the nine collision types materializes (Fig. 3) using the relative velocity of the particles and the mass of the projectile (see Paper I).
• Based on the collision recipe described in Paper I and Sect. 3.4, the new mass and new porosity of the representative particle is calculated and the new size and porosity distribution of the particles at time is obtained.
• In the final step, we update the collision rates.

3.6 Numerical issues

As mentioned in ZsD08, a sufficiently high number of representative particles is needed to properly reproduce the physics of the collision kernel. We performed simulations with an increasing number of representative particles and found that for more than 200-300 particles the results of the simulations do not change significantly. For all simulations described in the following sections, 500 representative particles are used and we average the results of 20 simulations to decrease the numerical uncertainty of the code.

The required computational time strongly depends on the collision rate of the particles thus determined by the dust density and the relative velocity (the turbulence parameter mostly) and the length of the simulation. On a 2.83 GHz CPU, performing the simulations on a single core, the CPU time varies between nine hours (the low density model with ) and three days (the high density model with ). Both simulations modeled 10⁶ years of particle evolution. The high density simulation with , critical mass ratio of 100, and t=10⁷ years of evolution takes twelve days to complete.

4 Results

4.1 Initial conditions and setup of simulations

All simulations begin with silicate monomers of 1.5 m diameter and 2 g cm^-3 material density (monodisperse size distribution). We simulate the dust evolution at the midplane of our disk models at a distance of 1 AU from the central star. The gas density is obtained from the disk models described in Sect. 2.1. We assume a typical 1:100 dust-to-gas ratio. We follow the history of each collision: the mass and porosity of the colliding particles, their relative velocity, the occurred collision type, and the new mass and porosity of the particles. In this way we are able to reconstruct the history of the dust evolution.

The parameters that we vary in this study are the gas density and the turbulence parameter . We also treat the critical mass ratio as a free parameter to explore its effect on the dust evolution.

We provide a detailed description of the low density model with and critical mass ratio of 100 in Sect. 4.2. We then compare this with the MMSN model and the high density model using the same turbulence parameter and the critical mass ratio (Sects. 4.3 and 4.4). In Sects. 4.5 and 4.6, we discuss the effects of changing the turbulence parameter and critical mass ratio by comparing those results with the two example runs.

4.2 The low density model

In this disk model, the gas density at 1 AU is  g cm^-3, the turbulence parameter is , and the critical mass ratio is . As shown in Fig. 1, the particles reach the fragmentation velocity (1 m s^-1) at sizes smaller than a millimeter because the particles in low gas density environment decouple from the gas at these small radii.

Figure 4 shows the evolution of the mass distribution (a), the porosity distribution (b), and the collision frequency of the various collision types (c). The x-axis shows the time on a logarithmic scale. The y-axis of Figs. 4a,b shows the mass and enlargement distributions, respectively. The intensity of the color reflects the number density of representative particles, which, as explained in ZsD08, indicates the mass density of the distribution. Thus, in Fig. 4a the intensity levels directly reflect the mass density, while in Fig. 4b the colors indicate the mass-weighted enlargement parameter. The black lines show the average of these quantities over the particle distribution. The y-axis in Fig. 4c represents the nine collision types used in this paper. Every stripe shows the total collision rate of the collision types at a given time.

Figure 5 represents the collision history in the eight collision regimes. The x-axis is the velocity, the y-axis shows the mass of the projectile. A mass-velocity grid is created and we calculate how many collisions occurred within each grid cell during which our ``local approach'' assumption is correct. The different collision types and their border lines, as well as the areas covered by the laboratory experiments (indicated by grey colors) are plotted. For more details on the experiments, we refer to Paper I.

In the ``cc'' panel, we indicate two curves with dotted lines. These curves are evolution tracks. The left curve is obtained by calculating the relative velocity between equal-sized particles with an enlargement parameter of 2.5 (volume-filling factor of 40%). The right curve represents the relative velocity between particles with a mass ratio of 100. These two curves serve as a guide to our results, as collisions should occur between these two curves in the ``cc'' panel. The lower part of the left curve, where the relative velocity decreases with increasing mass, is a sign that relative velocities between equal-sized particles are dominated by Brownian motion. For higher masses, the relative velocity is dominated by turbulence. These curves do not precisely match the contours because we assumed a constant enlargement parameter of 40% when calculating the evolution tracks, whereas  is a free parameter in the simulation.

4.2.1 Early evolution

We discuss the evolution of the distribution functions until the ``local approach'' assumption becomes invalid (  yr). The long-term evolution of the dust is discussed in Sect. 4.2.3.

We distinguish between two distinct phases. During the first 300 yr, particles grow by the means of the hit & stick (S1) mechanism. The second phase is dominated by bouncing with compaction (B1); the particles leave the S1 regimes. During this phase, the mass of the particles slowly decreases and the enlargement parameter asymptotically reaches a minimum value of 2.23. As discussed in Paper I, keeping the bouncing velocity of a particle constant, the porosity of the aggregate will asymptotically reach a maximum value, (see Paper I). The relative velocity of a particle is a function of the friction time (Eq. (7)), which depends on the ratio of the mass to surface area, m/A. Since particle growth is halted at this point in the simulation (m remain constant), only a decrease in A caused by compaction can further increase the velocity between particles. The particle radius can decrease until either for the given relative velocity is reached or particles reach the maximum compaction possible. The latter limit, random close packing (RCP), corresponds to an enlargement parameter of 1.6 (volume-filling factor of 60%).

We find that fragmentation does not play a role during the evolution of these particles indicated by Fig. 4c. As can be seen in Fig. 5, their evolution is halted by bouncing before the particles are able to reach the fragmentation barrier. The two dominant collision types are hit & stick (S1) and bouncing with compaction (B1).

4.2.2 Termination of growth

As we can see from Fig. 5, sticking at higher energies than the hit & stick (S1) limit is possible only inside the ``pP'' regime. As soon as we no longer have collisions inside this regime or the S1 regimes, the growth is halted. There can be two reasons why this occurs: 1.) all particles become compact, i.e., there are simply no collisions in the ``pP'' regime; 2.) the width of the particle mass distribution is smaller than the critical mass ratio (), such that all collisions take place in the equal-size regimes (e.g., ``pp'', ¸).

In the case of the current simulation, the small particles have been ``consumed''. Once the heavy particles grow into the bouncing with compaction (B1) area of the ``pp'' regime, their growth in the ``pp'' regime stops. The heavy particles collect the small ones via collisions in the ``pP'' regime and by doing so, the width of the distribution is reduced to a value smaller than . Therefore, before particles are able reach the fragmentation barrier, growth is halted. Because of B1, particles become compact and collisions in the ``cc'', ``cp'', and ¸ regimes appear.

4.2.3 Long-term evolution

Before discussing the long-term evolution of the distribution functions, we must consider how long our initial assumptions (``local approach'' and constant gas density) hold true.

Using Eq. (11), we calculate that a particle with Stokes number 10^-4 drifts a distance of 1 AU in roughly  yr. This is the drift timescale beyond which the ``local approach'' assumption (discussed in Sect. 3) is no longer valid and particles become separated from each other on this timescale.

Another process by means of which particles separate is viscous spreading. We assume that the viscous timescale of the disk at 1 AU is given by:

where r is the distance from the central star (1 AU), is defined in Eq. (6). The viscous timescale in our model, using , is of the order of 10⁶ yr.

One has to consider the results of the simulation with caution for longer times than the drift or viscous timescales. The equilibrium or final state of the particles is reached when mass decrease during bouncing with compaction (B1) and bouncing with mass transfer (B2), and mass increase by hit & stick (S1) and sticking through surface effects (S2) are in equilibrium. In other words, the final state is reached when the evolution of the average mass and the enlargement parameter can only be determined by the stochastic fluctuations of the simulation. We find that the equilibrium state of the particles is hardly reached within these timescales. Upon neglecting these warnings, we find that the equilibrium state of the dust is reached at  yr. The equilibrium is reached between the bouncing collisions, resulting in breakage and hit & stick (S1) (see Fig. 4c). The equilibrium average mass and porosity of the particles are , and .

To be able to compare the distribution functions of different runs, we define some quantities using the mean of the distribution functions indicated with black lines in Figs. 4a and b: , the maximum of the mean mass; , the maximum of the mean enlargement parameter; , the minimum mean enlargement parameter when particles can no longer be compacted anymore; , the time when is reached, which is when the time derivative of becomes zero ( ); and , the maximum average Stokes number reached during the simulation. The values of these quantities are listed in Table 1 (model id ``Lt1d-4m100''). In this table, Col. 1 describes the model names, ``L'' representing the low density model, ``M'' being the MMSN model, ``H'' being the high density model, the letter ``t'' and the following number indicates the value of the turbulence parameter, and the letter ``m'', and the number providing the used critical mass ratio values. Columns 2-4 show the gas density, turbulence parameter, and the critical mass ratio respectively. Columns 5-9 list the parameters defined to characterize the distribution functions. These are in Col. 5, in Col. 6, in Col. 7, in Col. 8, and finally in Col. 9.

4.3 The MMSN model

In the MMSN model, the gas density at 1 AU at the midplane is  g cm^-3, , and the critical mass ratio is 100. As shown in Fig. 1, the particles become larger than in the low density model, because they are more tightly coupled to the gas and the relative velocities are suppressed. As in the previous section, we first discuss the evolution of the distribution functions for as long as the ``local approach'' assumption holds true (  yr in this model).

Figure 6 again shows the time evolution of the mass (a), the enlargement parameter (b), and the collision frequency (c). Figure 7 shows the collision history. These figures depict a rather different evolution than that of the previous model.

4.3.1 Early evolution

During the fractal growth regime, we find that the collision rate of hit & stick (S1) is much higher than in the low density model (Fig. 6c). This is because of the higher dust densities. We can see from Fig. 7, for the ``cc'' regime, that growth begins with Brownian motion because the relative velocity decreases with increasing particle mass for particle masses less than 10<sup>-9</sup> g. As a result of these low velocity collisions, some particles reach enlargement parameter values of higher than 30 (volume-filling factors less than of 3.3%). At 200 yr, some particles grow above the border line of hit & stick (S1) and enter the area of sticking through surface effects (S2) in the ``pP'' plot, and bouncing with compaction (B1) in the ``pp'' plot. Growth caused by both S1 and S2 continues until different-sized particles enter the transition regime in the ``pP'' plot. One can see in Fig. 6a, that some particles reach 1 g in mass. However, when particle collisions enter the transition regime between bouncing with mass transfer (B2), and sticking through surface effects (S2) in the ``pP'' plot, their masses are equalized because of the mass transfer that occurs during the B2 collisions and the collisions shift to the similar-sized regime (B1). After roughly 10<sup>4</sup> yr, we find that particles mostly bounce and become compact. The enlargement parameter reaches a minimum value of 1.85 (54% volume filling factor), the mass distribution function slowly decreases because of the low probability of breakage. Collisions at this point occur mainly in the ``cc'' regime.

A peculiar feature of Fig. 6a is a peak at  yr, which is accompanied by a fast decrease in the enlargement parameter in Fig. 6b and an increased collision rate of sticking through surface effects (S2) and bouncing with mass transfer (B2) in Fig. 6c. At this point, the relative velocity produced by turbulence increases. As discussed in Sect. 2.2, particles leave the ``tightly coupled particle'' regime and enter the ``intermediate particle'' regime (see the relative velocity bump in Fig. 1b). We calculate the growth timescale of the heaviest particle with mass M in the simulation to be:

which is illustrated with a dotted line in Fig. 8. As a comparison, we also calculate the minimum growth timescale that a particle can have (solid line), which

where is the cross-section of the largest particle. Here, we assume that the ``swept up'' particles have masses of M/100, therefore we use the relative velocity curve presented in Fig. 1b, dotted line. The effects of both the relative velocity bump and the increased growth rate can be seen at 0.1 g.

The relative velocity ``boost'' happens shortly after the particles enter the transition regime of S2 and B2 in the ``pP'' plot. The heaviest particle, which experiences the velocity transition the earliest, acquires higher relative velocities, leading to a higher rate of collision with the other particles. As the particles are initially in the lower part of the S2-B2 transition regime (with masses of 10^-3 g, see Fig. 6a), the heaviest particle experiences rapid growth reaching masses of 30 g. The simulated timescale, however, does not reach the minimum growth timescale because of the bouncing with mass transfer (B2) collisions, which reduce the mass of the heaviest particle. The remainder of the particle population increase in mass because of B2, and the growth rate of the heaviest particle decreases. Eventually, the rapid growth of the heaviest particle is halted, and the growth timescale at m=30 g becomes infinity. From this point on, the heaviest particle reduced in mass, and B2 equalizes the masses of the particles.

4.3.2 Long-term evolution

We calculate the drift and viscous timescales to determine how long our assumptions of ``local approach'' and constant gas density remain valid. Assuming Stokes number 10^-4 particles, we find that the drift timescale is of the order of  yr, and the viscous timescale is 10⁶ yr. These timescales are indicated by solid and dotted white lines in Fig. 6.

We find that the final equilibrium is reached at  yr, which is longer than the drift and viscous timescales. The equilibrium is reached between the growth mechanisms of hit & stick (S1), sticking through surface effects (S2) and the destruction mechanisms of bouncing resulting in breakage and bouncing with mass transfer (B2). The final average mass and porosity of the particles are , .

We conclude that the dust evolution is more complex in the MMSN model than in the low density model because the complex interaction of the velocity field and the collision kernel is apparent in this model. As in the previous model, bouncing with compaction (B1) is the most frequent collision type and hit & stick (S1) determines the initial particle growth, but both sticking through surface effects (S2) and bouncing with mass transfer (B2) are of importance in this model. The final equilibrium is not reached within the drift and viscous timescales.

4.4 The high density model

The gas density in this model is  g cm^-3 at the midplane of the disk at a distance of 1 AU from the central star. The values of , , and the dust-to-gas ratio are the same as in the previous models.

Figure 1, dashed line, shows the relative velocity field of fluffy aggregates in this model. As already discussed in Sect. 2.2, the aggregates reach 1 m s^-1 relative velocities at similar masses as the MMSN model, because of the Stokes drag. Therefore, we expect that the final aggregate sizes and masses will be similar to the particles produced in the MMSN model.

Figure 9 shows the time evolution of the mass (a), enlargement parameter (b), and the collision frequency (c). Figure 10 illustrates the collision history.

4.4.1 Early evolution

As seen in Fig. 10, Brownian motion is the dominant source of relative velocity, as long as particles have masses of lower than 10^-8 g (that is an order of magnitude higher than in the MMSN model). Therefore, the enlargement parameter of the aggregates is also higher than in the MMSN model. As the hit & stick (S1) collisions are more frequent than in the MMSN model because of the higher dust densities, the particles reach the transition regime between sticking through surface effects (S2) - bouncing with mass transfer (B2) earlier, at t=200 yr. The peak in the mass distribution is not as pronounced as in the MMSN model. The relative velocity boost occurs at heavier aggregate masses (10^-2 g, see Fig. 1b) because of the higher gas density of the model. When the rapid growth of the heaviest particle begins, most of the projectiles are already in the transition regime. Here, the B2 collisions soon reduce the mass of the heaviest particle and narrow the mass distribution.

In contrast to the MMSN model, the mass of the particles is not reduced because of the low probability of breakage in bouncing with compaction (B1), but is kept nearly constant in time. This is the result of the increased collision rate of sticking through surface effects (S2). The S2 collision rate is increased because of the low velocity collisions, which occur when particles are in the tightly coupled regime and have similar stopping times. These S2 collisions occur in the ``pp'' regime as seen in Fig. 10. These collisions cancel out the effect of breakage in B1.

The maximum Stokes number reached in this model is (see Table 1, model id ``Ht1d-4m100''), which is lower than in the MMSN model. The growth in this model is halted by the bouncing with mass transfer (B2) collisions in the transition regime of the ``pP'' panel. This shows us that particles cannot reach masses much higher than 1 g independently from the exact value of the gas density, because at this point, particles enter the S2-B2 transition regime and the growth is halted. Increasing the gas density yet further would lead to even lower Stokes numbers.

4.4.2 Long-term evolution

The drift and the viscous timescales in the high density model are both 10⁶ yr. As seen in Fig. 9a, the particle masses do not change significantly after t=10³ yr. The porosity is reduced by bouncing with compaction (B1) and reaches a final value of 5.41 at  yr.

4.5 Varying the turbulence parameter

Table 1:   Overview and results of all the simulations.

To explore the effects of turbulence, we perform two more simulations in each of the disk models. We keep the critical mass ratio fixed (100) and vary only the turbulence parameter () to obtain values of 10^-3, 10^-4, and 10^-5. The results are shown in Table 1, for the first nine models, and in Fig. 11a.

The work of Brauer et al. (2008a) suggests that in situations where fragmentation limits the growth, a lower turbulence strength produces larger aggregates. This, of course, directly reflects the shift in the fragmentation threshold (1 m/s) to larger sizes when is lower (Fig. 1). In this study, it is fragmentation that balances the growth, producing a (quasi) steady-state. For the low density models, we do see a decrease in the final particle mass, but it is balanced by bouncing and not fragmentation. In the low density model, particles grow only in the hit & stick (S1) regimes. When particles leave these regimes, the growth stops due to bouncing. The border of the S1 regime corresponds to the collision energy that is lower than , where is the rolling energy of monomers (see Paper I). As the collision energy is , particles in strong turbulence leave the S1 regimes at lower particle masses.

On the other hand, the MMSN and high density models show that the maximum mass of the particles can even increase with . The precise value of the is determined by the intensity of the peak in the mass-density plots (Sect. 4.3.1) and this may vary between simulations. In the ``Mt1d-4m100'' model, we have argued that the spike is exceptionally pronounced because of the high probability of sticking through surface effects (S2) collisions at the initial part of the rapid growth. However, the main point is that in the MMSN/high-density simulations the maximum particle masses all end up around 1 g, independent of the turbulent strength because of the nature of the S2-B2 transition, which occurs at projectile masses of 10<sup>-4</sup> g in the ``pP'' plot. As explained before, collisions in the ``pP'' plot are the only way by which particles can grow after the hit & stick (S1) phase is finished. Thus, we require a broad distribution with a high growth rate. However, B2 collisions operates in the opposite way: they transfer mass from the target to the projectile, narrowing the distribution and decreasing the overall probability of the ``pP'' process occuring. Thus, once B2 becomes effective, there is a shift from the ``pP'' panel to the ``pp'' panel. For the MMSN/high-density models this behavior is always present and the important quantities involved (i.e., relative probability of B2 over S2) scale with mass but not with velocity. The result is that the maximum masses that particles reach are 1 g, which is rather insensitive to the strength of the turbulence.

4.6 Varying the critical mass ratio

We perform simulations in the disk models with but a varying critical mass ratio. We explore how the dust distributions change upon using , 100, and 1000. Table 1, lines 10 to 18, shows the parameters that describe the distribution functions, and Fig. 11b illustrates the maximum particle mass as a function of the critical mass ratio.

By examining Table 1, we see that using in the low density model (``Lt1d-4m10'') results in heavier and more compact particles. The low critical mass ratio means that the largest particles in the different-sized regimes can sweep up the projectiles and grow to larger sizes, eventually reaching the fragmentation line, where growth stops. As discussed in Sect. 2.2, assuming a fragmentation velocity of 1 m s^-1, the maximum Stokes number of the aggregates is , a value almost reached in this model.

We find that there is no significant difference between the and 1000 simulations in the low density model. The explanation for this can be found by examining the width of the mass distribution in the hit & stick (S1) phase. This initial phase occurs in the same way independently of the critical mass ratio. If the critical mass ratio is equal to or larger than the width of the distribution function, collisions between different-size particles in the ``pP'' regime are inhibited. After the S1 phase, the width of the distribution in the low density regime is approximately 100. Therefore, we do not see any difference when the mass threshold is shifted from to ; in both cases, collisions occur between equal-size particles only, and these are either S1 or (when this stage is over) B1.

For the high density models (MMSN/Desch), we find that the outcome is again similar: growth halts at 0.1 g (within a factor of 10) and no clear dependence on is seen. For the high mass ratios, growth is always in the similar-size regime. Here, it is the gas density that determines the velocity, i.e., whether we have a sticking (S1) or a bouncing (B1) collision. Therefore, if , the high density model produces heavier particles than the MMSN model (see Fig. 11b). For lower , it is again the nature of the S2-B2 transition regime that limits the maximum mass.

Thus, the critical mass ratio is an important parameter because it determines the relative likelihood of collisions occurring in the different-size regime, which are in general more conducive to growth. In contrast, for simulations where B2 collisions are important - which have the effect of narrowing the distribution - the width of the distribution will correspond to the value of the parameter, although we have also seen that the absolute size/mass is rather insensitive to the parameter. Overall, these arguments indicate that a good knowledge of this parameter is important.

5 Discussion

We have performed simulations of varying turbulence parameter and critical mass ratio values in three disk models with low, intermediate and high gas densities. We have found that hit & stick (S1) and bouncing with compaction (B1) are the most dominant collision types. All simulations show evidence of long-lived, quasi-steady states. Fragmentation is rarely present, but even then, for only a limited time period. The absence of fragmentation is caused by the bouncing collisions.

5.1 The sensitivity of the results

As presented in Sect. 4, the outcome of our simulations is determined by the collision kernel, and the relative velocity field. A significant change in one, or both can alter the evolution of the aggregates.

We present the results of a test simulation, where the sticking through surface effects (S2) - bouncing with mass transfer (B2) transition regime in the ``pP'' plot is neglected and replaced by S2 collisions. This alternative transition regime provides a good opportunity to examine the rapid growth presented in Sect. 4.3.1, as the kernel is now simplified. The new kernel also provides us with the possibility to see how much the outcome of our simulations can be altered by changing critical areas of the parameter space. As the transition regime is only constrained by one experiment in a rather small area (see e.g. Fig. 5 or Fig. 11 in Paper I), this part of the parameter space may require changing in the future. We use the same initial conditions as in the ``Mt1d-4m100'' model described in Sect. 4.3.

In this case, the heaviest particle experiences increased relative velocities, as soon as it reaches m=0.1 g, and the particle undergoes a rapid growth period (as in the original MMSN simulation, Sect. 4.3). Figure 13 indicates the growth timescale of the heaviest particle (dotted line) and the minimum growth timescale possible (solid line). Since none of the B2 collisions can reduce the mass of the heaviest particle, the growth timescale reaches its possible maximum. The heaviest particle increases in mass until the rest of the particle population enters the B2 regimes above 0.1 g in the ``pP'' plot. In this simulation, the maximum average mass is 27 g, whereas in the original simulation for the transition regime, the value was 4.18 g.

Together with Paper I, this work is the first attempt to calculate dust growth in protoplanetary disks on an empirical, thus more realistic basis. However, a few additional cycles of the feedback loop between the laboratory experiments, the models of the kind described by Paper I, and the models described in this paper have to be conducted before we can obtain a truly reliable model of dust growth in protoplanetary disks.

5.2 Retention of small grains

Dullemond & Dominik (2005) showed that without a mechanism that reduces the sticking probability of particles in the upper layers of the disk or without a continuous source of small particles, the observed SEDs of T Tauri stars should exhibit a very weak infrared excess. The SEDs of T Tauri stars have strong IR excess (e.g., Furlan et al. 2005; Kessler-Silacci et al. 2006); therefore, some kind of grain-retention mechanism is needed to explain these SEDs. Previous models of grain growth assumed a continuous cycle of growth and fragmentation, which provides the necessary amount of small particles (see e.g., Brauer et al. 2008a; Dullemond & Dominik 2005; Birnstiel et al. 2009). Our simulations, however, have shown that the mass distribution function is narrow. Small, monomer-sized particles are not present and fragmentation is ineffective in providing small particles, which could be transported to disk atmospheres. The question naturally arises: how can small grains be produced in our collision model?

One possible solution might be provided by bouncing. Weidling et al. (2009) performed bouncing experiments by placing an aggregate onto an oscillating metal plate and measuring the porosity of particles due to collisions with the plate. They observed that approximately 10% of the projectile mass was eroded during the experiment (see Table 1 of their paper). This mass loss may be caused by the initial collisions, the eroded mass sticking to the baseplate. It is also possible that small pieces of fragments grind off when the aggregates bounce, which cannot be observed in the experiment. These ground-off particles can then diffuse out of the midplane and provide the small particles that we observe to enter to the upper layers of the disk. Future laboratory experiments are needed to quantify the level of ground-off particles created by bouncing collisions.

The second possible explanation involves dust growth in the upper layers of the disk. We performed two simulations at four pressure scale-heights in both the low density model and the MMSN model for . We find that the relative velocity of two monomers in the Brauer model is 2 m s^-1, thus monomers at these heights do not coagulate, only bounce. The particles in the MMSN model can form aggregates of a maximum of 10 m in size. Using a higher value of  (as usually assumed in the upper layers of the disk), we can completely halt even this limited growth. Therefore, bouncing could be the key mechanism reducing the sticking probability of the particles. However, if substantial vertical turbulent mixing takes place, bouncing may not be able to help, because these monomers would then be ``vacuum-cleaned'' away by the larger particles at the interior of the disk. Additional studies of 1D vertical slices of disk models are needed to investigate this scenario.

5.3 Implications for planetesimal formation models

It is also evident that coagulation alone is unable to produce planetesimals under the conditions presented in this work. Even if the turbulence parameter is assumed to be zero, relative velocity caused by radial drift prevents particles crossing the so called ``meter-size barrier''. An ideal environment for particle growth is a pressure bump in the dead zone, where both the turbulent and radial relative velocities are reduced. This environment is located around the snow line (Kretke & Lin 2007). Brauer et al. (2008b) showed that in these pressure bumps relative velocities stayed below a fragmentation threshold of 10 m s^-1, providing a window through which particles could overcome the m-size barrier, although they assumed perfect sticking (no bouncing) below the fragmentation barrier. Future studies should verify whether planetesimals can be formed with the collision model presented in this study.

Another planetesimal-forming mechanism is the gravitational collapse of swarms of boulders (Johansen et al. 2007). This scenario assumes that a large amount of the solid material is presented in dm-sized boulders ( ) at the midplane of the disk. These boulders then concentrate in long-lived high pressure regions in the turbulent gas and initial over-densities are amplified further by the streaming instability. This mechanism forms 100 km sized objects on a very short timescale (some orbits). However, our simulations produce particles with , which is due to bouncing with compaction (B1) and the low (1 m s^-1) fragmentation velocity of silicates. Using a ``stickier'' material such as ices or particles with organic mantels may produce larger particles. Molecular dynamical simulations (e.g., Dominik & Tielens 1997; Wada et al. 2007; Wada et al. 2008) have shown that icy aggregates could have fragmentation velocities of about 10 m s^-1, although these findings have yet to be confirmed by laboratory experiments. Similarly, it is conceivable that the enhanced sticking capabilities of ices can prevent bouncing, which is so omnipresent for small particles in our simulations, or shifts its proficiency to larger sizes.

Cuzzi et al. (2008) outlined an alternative concentration mechanism to obtain gravitationally unstable clumps of particles, which can then undergo sedimentation and form a ``sandpile'' planetesimal. In this model, turbulence causes dense concentrations of aerodynamically size-sorted, chondrule-size particles (Cuzzi et al. 2001) or more precisely, particles of Stokes numbers in our simulations. Since growth in our models is typically halted at these Stokes numbers, this concentration mechanism is an obvious successor to coagulation - at least where it concerns the conditions adopted in this paper (1 AU, silicates).

However, we emphasize that the formation of a gravitationally unstable clump does not imply that planetesimals will form without impediment. An important question is how collisions will affect the collapse. In the Cuzzi et al. (2008) scenario, the collapse occurs on a sedimentation timescale and at these high densities, collisions between particles are frequent. Likewise, in the Johansen scenario - where the collapse occurs on an orbital timescale and involves particles - collisions can be rather violent. Collisional fragmentation or erosion may change the appearance of the collapse, because the small fragments are carried away by the gas. The role of collisions in these situations is certainly an important question, and our new collision model provides a tool to quantitatively address this issue in future studies.

5.4 Consequences for laboratory experiments

From Fig. 11, in Paper I, one can see that only a small part of the relevant parameter space has been covered by experiments. Although laboratory experiments cannot be performed for every point of the parameter space, we suggest future ones based on Figs. 5, 7, and 10 to understand dust growth in the early stages of planet formation.

• More experiments in the ``cc'' and ``cC'' regimes are needed as particles become compactified toward the end of their evolution. Thus, most of the collisions occur in this regime, at velocities between 0.1 and 100 cm s^-1, and masses between 10^-7 and 10 g.
• As seen in Figs. 5, 7, and 10, the ``hot spots'', where most of the collisions occur, are located in equal-sized regimes, at the left side of the fragmentation line. Therefore, it is important to map these areas of the parameter space in detail.
• We define a sharp border line between the hit & stick (S1) and bouncing with compaction (B1) collisions. If there is a continuous transition between S1 and B1, the growth of particles should not be halted by bouncing at these low particle sizes. Since many collisions occur in the ``pp'' and ``cc'' regimes, even a small probability of growth could increase the particle sizes.
• As seen in Fig. 6b, particles in high gas density environments can have enlargement parameters that are much higher than 6.6 ( ). An interesting question is whether the collision types and regimes are also valid for particles with such low volume-filling factors, or whether these particles have different collision behaviors?
• The sticking through surface effects (S2) - bouncing with mass transfer (B2) transition regime greatly affects the outcome of the simulations (see Sect. 5.1). However, the transition regime is only mapped at the high velocity and low mass regions. Therefore, it is important to constrain more tightly this part of the parameter space.
• The critical-mass ratio affects both particle masses and porosities. Experiments are needed to constrain its value.
• The bouncing model, described in Paper I, has important implications for the evolution of dust aggregates in protoplanetary disks but is unfortunately still constrained by too few experiments. Additional experiments are needed to refine the model, as bouncing with compaction (B1) is the most frequent collision type in all of the simulations.

6 Summary

We have performed simulations of dust growth using the Monte Carlo code of ZsD08 and a dust collision model based on laboratory experiments (Paper I). We have performed simulations in the midplane of three disk models at low (  g cm^-3), intermediate (  g cm^-3), and high (  g cm^-3) gas densities at 1 AU distances from the central star. We have varied the turbulence parameter () and the critical-mass ratio ( ) to explore their effects on the mass and porosity distribution functions. Our main results are:

• Upon using , the low-density/MMSN/high-density model produces particles with maximum mean mass of  g/4.18 g/0.23 g, respectively, the maximum average enlargement parameter of these particles are 7.12/21.9/38.0, respectively. The maximum average Stokes numbers are / / , respectively.
• We find that particle evolution does not follow the previously assumed growth-fragmentation cycles. Although catastrophic fragmentation is present for a short period of time in some models (typically when ), it has a fringe effect. Particles in most of the simulations do not reach the fragmentation barrier because their growth is halted by bouncing.
• We see long-lived, quasi-steady states in the distribution function of the aggregates that are caused by bouncing. The final equilibrium state is not reached within the drift or the viscous timescales.
• We have performed simulations of varying turbulence strength. We find that the system is ``non-linear'': the maximum mass of particles is not a decreasing function of the turbulence parameter and is not an increasing function of the gas density.
• We have explored the effects of the critical mass ratio. We find that different critical mass ratios can affect the particle evolution. Low critical mass ratios can produce heavier particles, while high values of can halt the growth earlier.
• The maximum Stokes number is almost independent of either the gas density or the strength of the turbulence.
• The maximum mass of the aggregates is limited to 1 g because of the S2-B2 transition regime.
• The Stokes number 10^-4 particles can be concentrated by aerodynamical size-sorting, thus planetesimals can form from these particles.

Acknowledgements

A. Zs. would like to thank to Frithjof Brauer, Tilman Birnstiel and Felipe Gerhard for useful discussions. C.G. was funded by the Deutsche Forschungsgemeinschaft within the Forschergruppe 759 ``The Formation of Planets: The Critical First Growth Phase'' under grant Bl 298/7-1. A. Zs. acknowledges the support of the IMPRS for Astronomy & Cosmic Physics at the University of Heidelberg and C.W.O. acknowledges the financial support from the Alexander von Humboldt Foundation. We also thank our referee, Sascha Kempf for providing useful comments which helped to improve the paper.

Appendix A: Time evolution and animations

We present six animations to illustrate the time evolution of the aggregates. These animations are available in the electronic edition of the journal at http://www.aanda.org. The mass_vs_psi_lowd.mpg file shows the time evolution of the masses and enlargement parameters of our representative particles. The x-axis is the particle mass in gram units, the y-axis is the enlargement parameter. The regimes_lowd.mpg file illustrates the time evolution of Fig. 5. The contour levels are collision frequencies/grid cell. The same movies are provided for the MMSN model, and are called mass_vs_psi_mmsn.mpg and regimes_mmsn.mpg respectively. For the high density model the relevant files are mass_vs_psi_highd.mpg and regimes_highd.mpg.

References

Online Material

mass_vs_phi_highd.mpg
mass_vs_phi_MMSN.mpg
regimes_lowd.mpg
mass_vs_phi_lowd.mpg
regimes_highd.mpg
regimes_mmsn.mpg

Footnotes

... planetesimals?

This paper is dedicated to the memory of our dear friend and colleague Frithjof Brauer (14th March 1980-19th September 2009) who developed powerful models of dust coagulation and fragmentation, and thereby studied the formation of planetesimals beyond the meter-size barrier in his Ph.D. Thesis. Rest in peace, Frithjof.

...

Movies are only available in electronique form at http://www.aanda.org

All Tables

Table 1:   Overview and results of all the simulations.

All Figures

---

Copyright ESO 2010