A&A 461, 215-232 (2007)
DOI: 10.1051/0004-6361:20065949

Dust coagulation in protoplanetary disks: porosity matters

C. W. Ormel¹ - M. Spaans¹ - A. G. G. M. Tielens^1,2

1 - Kapteyn Astronomical Institute, University of Groningen, PO box 800, 9700 AV Groningen, The Netherlands
2 - NASA Ames Research Center, Mail Stop 245-3, Moffett Field, CA 94035, USA

Received 30 June 2006 / Accepted 27 September 2006

Abstract
Context. Sticking of colliding dust particles through van der Waals forces is the first stage in the grain growth process in protoplanetary disks, eventually leading to the formation of comets, asteroids and planets. A key aspect of the collisional evolution is the coupling between dust and gas motions, which depends on the internal structure (porosity) of aggregates.
Aims. To quantify the importance of the internal structure on the collisional evolution of particles, and to create a new coagulation model to investigate the difference between porous and compact coagulation in the context of a turbulent protoplanetary disk.
Methods. We have developed simple prescriptions for the collisional evolution of porosity of grain-aggregates in grain-grain collisions. Three regimes can then be distinguished: "hit-and-stick'' at low velocities, with an increase in porosity; compaction at intermediate velocities, with a decrease of porosity; and fragmentation at high velocities. This study has been restricted to physical regimes where fragmentation is unimportant. The temporal evolution has been followed using a Monte Carlo coagulation code.
Results. This collision model is applied to the conditions of the (gas dominated) protoplanetary disk, with an $\alpha _{\rm T}$ parameter characterising the turbulent viscosity. We can discern three different stages in the particle growth process. Initially, growth is driven by Brownian motion and the relatively low velocities involved lead to a rapid increase in porosity of the growing aggregate. The subsequent second stage is characterised by much higher, turbulent driven velocities and the particles compact. As they compact, their mass-to-surface area increases and eventually they enter the third stage, the settling out to the mid-plane. We find that when compared to standard, compact models of coagulation, porous growth delays the onset of settling, because the surface area-to-mass ratio is higher, a consequence of the build-up of porosity during the initial stages. As a result, particles grow orders of magnitudes larger in mass before they rain-out to the mid-plane. Depending on the precise value of $\alpha _{\rm T}$ and on the position in the nebula, aggregates can grow to (porous) sizes of $\sim$ ${\rm 10\ cm}$ in a few thousand years. We also find that collisional energies are higher than in the limited PCA/CCA fractal models, thereby allowing aggregates to restructure. It is concluded that the microphysics of collisions plays a key role in the growth process.

Key words: ISM: dust, extinction - planetary systems: formation - planetary systems: protoplanetary disks - accretion, accretion disks

1 Introduction

Understanding the formation of planetary systems is one of the central themes of modern astrophysics. New stars form in molecular cloud cores when these cores contract under the influence of gravity. This contraction leads to the formation of a central object surrounded by a disk (Shu et al. 1987). Planets are generally thought to form in these disks, but neither the precise physical conditions required for, nor the processes involved in planetary body assemblage are well understood. Generally, it is thought that grain growth starts with the sticking of sub-micron-sized grains colliding at low velocities (Weidenschilling & Cuzzi 1993). The sticking is then driven by weak van der Waals interaction forces between the grains. Relative velocities may reflect Brownian motion or differences in coupling to turbulence in the disk. Eventually, when the grains have grown to $\sim$cm-sizes, they will settle to the mid-plane of the disk, forming a thin dust layer where further growth to planetesimal sizes can take place (Weidenschilling & Cuzzi 1993; Safronov 1969).

There is much observational support for the growth of dust grains in protoplanetary disks from sub-micron to millimetre size scales. In particular, the contrast of the $10~ \mu{\rm m}$ silicate emission feature relative to the local continuum shows that the grain size in the disk's photosphere - where these features originate - has increased from sub-micron sizes characteristic for interstellar dust, to the micron-sized range (Kessler-Silacci et al. 2006; Przygodda et al. 2003; Meeus et al. 2003; van Boekel et al. 2003). Furthermore, observations of the continuum (sub)millimetre emission - originating from the deeper layers of protoplanetary disks - show typical grain sizes in the range of millimetres, outside the range of interstellar grain sizes by many orders of magnitudes (Beckwith & Sargent 1991). Additional support for the importance of collisional grain growth follows from analytical studies of solar system dust. In particular, many interplanetary dust particles (IDPs) collected at stratospheric altitudes, consist of a large number of small grains assembled in a very open structure as expected for collisional aggregates (Brownlee 1979). These types of IDPs are thought to derive from comets and, indeed, comets may consist largely of such loosely bound aggregates. In addition, chondrules recovered from meteorites show dust rims which are generally attributed to collisional accretion processes in the solar nebula before the meteorite and its parent body were assembled (Metzler et al. 1992; Cuzzi 2004).

The properties of the dust are of key importance to the evolution of protoplanetary disks. First and foremost, planet formation starts at the dust size and the dust characteristics at the smallest sizes will set the table for further growth. Second, the opacity is dominated by absorption and scattering by dust grains. Hence, the radiative transfer, temperature structure, as well as emission spectrum of protoplanetary disks are controlled by the characteristics of the dust (Beckwith et al. 2000; Bouwman et al. 2000). Third, in turn, the temperature structure will dominate the structure of the disk, including such aspects as flaring. Fourth, dust grains provide convenient surfaces that can promote chemical reactions. Specifically, ice mantles formed by accretion and reactions between simple precursor molecules are widespread in regions of star formation (Boogert et al. 2004; Gibb et al. 2004). In fact, grain surfaces may be catalytically active in the warm gas of the inner disks around protostars, converting CO into CH₄ (Kress & Tielens 2001).

Most early studies of the coagulation process and the characteristics of the resulting aggregates assumed hit-and-stick collisions where randomly colliding partners stick at the point of initial contact (Meakin 1988; Ossenkopf 1993; Meakin & Donn 1988). The structure of the aggregate then depends on whether the collision is between a cluster and a monomer (PCA) or between two clusters (CCA). The latter leads to very open and fluffy structures with fractal dimensions less than 2, while the former leads to more compact structures and a fractal dimension (for large aggregates) near 3. Ossenkopf (1993) also investigated the pre-fractal limit in which aggregates consist of $\lesssim$1000 monomers. He provides simple analytical expressions for, e.g., the geometrical and collisional cross-section in the case of PCA and CCA aggregation. These expressions include a structural parameter, x, which describes the openness (or fluffiness) of the particle. In this study we also introduce a structural parameter and confirm its importance in coagulation studies.

Over the last decade much insight has been gained in the structure of collisional aggregates through extensive, elegant, experimental studies (Blum 2004; Blum et al. 2002) supported by theoretical analysis (Dominik & Tielens 1997; Chokshi et al. 1993). These studies have revealed the importance of rolling of the constituent monomers for the resulting aggregate structure. At low collision velocities, the hit-and-stick assumption is well justified but at high collision velocities, aggregates will absorb much of the collision energy by restructuring to a more compact configuration. At very high velocities, collisions will lead to disruption, fragmentation, and a decrease in particle mass. The critical velocities separating these collisional regimes are related to material properties such as surface free energy and Young's modulus as well as monomer size and the size of the clusters.

While the porous and open structure of collisional aggregates is well recognised, their importance for the evolution of growing aggregates in a protoplanetary setting is largely ignored. Most theoretical studies represent growing aggregates either by an equivalent sphere (e.g., Nomura & Nakagawa 2006; Mizuno et al. 1988; Weidenschilling 1984a; Tanaka et al. 2005) or adopt the fractal dimension linking mass and size characteristic for CCA or PCA growth (e.g., Weidenschilling 1997; Dullemond & Dominik 2005; Suttner & Yorke 2001). Ossenkopf (1993) and Kempf et al. (1999) explicitly account for aggregates' internal structure in their numerical models, although, due to computational reasons, only a limited growth can be simulated. Indeed, the internal structure of collisional aggregates is the key to their subsequent growth. The coupling of aggregates to the turbulent motion of the gas is controlled by their surface area-to-mass ratio, while the relative velocity between the collision partners dictates in turn the restructuring of the resulting aggregate. Moreover, as a result of the growth process from sub-micron-sized monomers to cm-sized aggregates, the coupling to the gas velocity field may well change due to compaction. Indeed, compaction can be an important catalyst for aggregates to settle out in a mid-plane dust layer. Despite its importance for the collisional growth of aggregates in a protoplanetary environment, the evolution of porosity has not yet been theoretically investigated. The present paper focuses precisely on this aspect of grain growth in protoplanetary disks.

This paper is organised as follows. In Sect. 2 a model is presented for the treatment of the porosity as a dynamic variable. This entails defining how porosity, or rather the openness of aggregates, is related to the surface area-to-mass ratio (Sect. 2.2) as well as quantifying how it is affected by collisions (Sect. 2.3). In Sect. 3 a Monte Carlo model is presented to compute the collisional evolution, taking full account of the collisional aspect and all features of the porosity model. Section 4 then applies the model to the upper regions of the protoplanetary disks. Results from the porous model of this paper are compared to traditional, compact models. In Sect. 5 we discuss the differences of the new collision model with respect to PCA and CCA models of aggregation and also discuss our results from an observational perspective, before summarising the main results in Sect. 6.

  
2 Collision model

Dust grains are dynamically coupled to the turbulent motions of the gas and "suspended'' in the (slowly accreting) protoplanetary nebula. Upon collisions, these small dust grains can stick due to van der Waals forces, forming larger aggregates. Eventually, when these agglomerates have grown very large ($\sim$cm-sized), they will decouple from the gas motions and settle in a thin disk at the mid-plane. Further collisional growth in the mid-plane can then lead to the formation of planetesimals ($\sim$km-sized). Up to this point, growth is driven by weak van der Waals forces, but for km-sized planetesimals gravitational forces take over and rapid growth to a planet takes place. Here we focus on this process of small grains suspended in the nebula forming larger conglomerates. Coagulation is driven by the relative grain velocities. Velocities and the kinematics of dust in a turbulent nebula are discussed in Sect. 2.1. The frictional coupling between dust and gas depends largely on the area-to-mass ratio of the grains and hence on the internal structure of the dust. Section 2.2 describes the relation between the area-to-mass ratio and the porosity of the dust agglomerates. In Sect. 2.3, we discuss experimental and theoretical studies on the microphysics of dust coagulation and develop a simple model, given in the form of easily applicable recipes, which describes how the mass and porosity of the dust evolve in collisions between two dust agglomerates. In Sect. 2.4, finally, we compare these recipes in the fractal limit to the frequently used formulations of Particle-Cluster Aggregation (PCA) and Cluster-Cluster Aggregation (CCA).

Figure 1: Sketch of the turbulence induced relative velocities, $\Delta {v}$, as function of friction time, $\tau _1$, for equal friction times, $\tau _1 = \tau _2$, (solid line) and different friction times $\tau _2 \ll \tau _2$ (dashed line) according to the analytic fits of Eqs. ((5), (5)) ( after Weidenschilling 1984a). t0 and v0 are set by the global properties of the disk, while the range over which turbulence is important is determined by the Reynolds number, Re. Values between brackets denote typical values for $t_0 = 1~ {\rm yr}$, $c_{\rm g} = 10^5~{\rm cm~ s}^{-1}$ and $\alpha _{\rm T} = 10^{-4}$.

Open with DEXTER

  
2.1 The turbulent protoplanetary disk

For the characterisation of the gas parameters of the protoplanetary disk we use the minimum-mass solar nebula (MSN) model as described by Hayashi (1981) and Nakagawa et al. (1986). The surface gas density of the disk, $\Sigma_{\rm g}$, is assumed to fall off as a -1.5 power-law with heliocentric radius (R) and the temperature scales as R^-1/2. The vertical structure of the disk is assumed isothermal, resulting in a density distribution which is Gaussian in the height above the mid-plane, z. The scaleheight of the disk, $H_{\rm g}$, is an increasing function of radius, $H_{\rm g} = c_{\rm g}/\Omega \propto R^{5/4}$, with $c_{\rm g}$ the local isothermal sound speed and $\Omega$ the Keplerian rotation frequency. The thermal gas motions will induce relative (Brownian) velocities between dust particles with a mean rms-relative velocity of

$$\Delta {v_{\rm Brownian}} (m_1, m_2) = \sqrt{ \frac{8k_{\rm B}T(m_1+m_2)}{\pi m_1m_2} },$$ (1)

with m₁,m₂ the masses of the particles and $k_{\rm B}$ Boltzmann's constant. Except for low mass particles of size $\lesssim 10\ \mu{\rm m}$, these velocities are negligibly small when compared to the relative velocities induced by the coupling with the turbulent gas. We assume that the turbulence is characterised by the Shakura & Sunyaev (1973) $\alpha _{\rm T}$-parameter,

$$\nu_{\rm T} = \alpha_{\rm T} c_{\rm g} H_{\rm g} = \alpha_{\rm T} c_{\rm g}^2 / \Omega \approx {v}_0 \ell_0 = {v}_0^2 t_0,$$ (2)

with ${v}_0, \ell_0$ and t₀, respectively, the velocity, the size and the turn-over time of the largest eddies. According to the standard (Kolmogorov) turbulence theory, turbulent energy is injected on the largest scales and transported to and eventually dissipated at the smallest eddies - characterised by turn-over time (or dissipation timescale), $t_{\rm s}$, velocity, ${v}_{\rm s}$, and scale size, $\ell_{\rm s}$, given by $\nu_{\rm m} = {v}_{\rm s} \ell_{\rm s}$, with $\nu_{\rm m}$ the molecular viscosity. We can then define the Reynolds number as, ${\rm Re} = \nu_{\rm T}/\nu_{\rm m}$, giving ${v}_{\rm s} = {\rm Re}^{-1/4} {v}_0$ and $t_{\rm s} = {\rm Re}^{-1/2} t_0$. If t₀ is assumed to be (of the order of) the Kepler time (Dubrulle & Valdettaro 1992), $t_0 = 2\pi/\Omega^{-1}$, the turbulence is fully characterised by the $\alpha_T$-parameter (see Fig. 1):

$$t_0 = 2\pi \Omega^{-1}, \qquad t_{\rm s} = {\rm Re}^{-1/2} t_0$$ (3a)

$${v}_0 = \alpha_{\rm T}^{1/2} c_{\rm g}, \qquad {v}_{\rm s} = {\rm Re}^{-1/4} {v_0}$$ (3b)

$$\ell_0 = \alpha_{\rm T}^{1/2} H_{\rm g}, \qquad \ell_{\rm s} = {\rm Re}^{-3/4} \ell_0.$$ (3c)

This specification of turbulence is of importance, for, together with the friction times of two particles, it determines the (average root-mean-square) velocity between the particles, $\Delta {v}_{ij}$, which in turn plays a crucial role in both the collision model of this section as well as in the time-evolution model of Sect. 2.3. These relative velocities are a function of the friction time ( $\tau_{\rm f}$) of the particles - the time it takes a particle to react to changes in the motion of the surrounding gas - which in the Epstein limit is

$$\tau_{\rm f} = \frac{3}{4c_{\rm g} \rho_{\rm g}} \frac{m}{A},$$ (4)

where $\rho_{\rm g}$ is the local mass density of the gas, m the mass of the particle and A its cross-section. In particular, if the friction time of a particle is much smaller than the turnover time of the smallest eddy ( $\tau_{\rm f} \ll t_{\rm s}$), the particle is coupled to all eddies and flows with the gas. Therefore, grains with $\tau_{\rm f} \ll t_{\rm s}$ have highly correlated velocities. Eventually, however, due to growth and compaction, agglomerates will cross the Kolmogorov "barrier'' (i.e., $\tau_{\rm f} > t_{\rm s}$) and are then insensitive to eddies with turnover times smaller than $\tau_{\rm f}$, since these eddies will have disappeared before they are capable of "aligning'' the particle's motion. At this point, grains can develop large relative motions (Fig. 1).

Figure 2: Projections of fractal aggregates illustrating the relation between A, V and V*. All fractals consist of 1000 monomers and have consequently the same compact volume, V* (inner circle). The black, outer circle gives the area, A, equal to the total projected surface area of the agglomerate. This area subsequently defines the volume, V (e.g., $V \sim A^{3/2}$) and the enlargement factor $\psi$ of the aggregate (Eq. (6)). (left) Compact aggregate, defining the compact volume V*. (centre) Porous aggregate. The geometric cross-section, A (outer circle), is larger than its compact equivalent (inner circle). (right) Even more porous aggregate. Arrows denote the compact and porous radii, a* and a.

Open with DEXTER

To calculate relative velocities accurately, contributions from all eddies have to be included. Voelk et al. (1980) started quantifying these effects by dividing eddies into several classes and subsequently integrated over the full Kolmogorov power spectrum. Cuzzi & Hogan (2003) provide analytical expressions for $\Delta {v}$between two particles of equal $\tau_{\rm f}$; in general, however, $\tau_1\neq\tau_2$ and results are derived numerically (Voelk et al. 1980; Markiewicz et al. 1991). Weidenschilling (1984a,b) has presented analytical fits to these numerical result, which have been frequently applied in subsequent coagulation models (e.g., Tanaka et al. 2005; Dullemond & Dominik 2005; Suttner & Yorke 2001). For particle friction times, $\tau _1$ and $\tau_2$ ( $\tau_1 \ge \tau_2$), less than t₀, these read,

$$\Delta {v}_{12}\ (\tau_1, \tau_2) = \left\{ \begin{array}{ll... ...rm s}}} & \qquad t_{\rm s} < \tau_1 <t_0, \end{array} \right.$$ (5a)

where $\tau _1$ is the larger of the two friction times. If $\tau _1$ exceeds t₀ the relative velocities become,

$$\Delta {v}_{12} = \left\{ \begin{array}{ll} {v}_0 & \qquad \... ...au_1\tau_2} & \qquad \tau_1, \tau_2 > t_0. \end{array}\right.$$ (5b)

For example, in the regime where both friction times are small ( $\tau_1, \tau_2 < t_{\rm s}$) the turbulence induced relative velocity is negligible when $\tau_2 \approx \tau_1$, but scales linearly with $\tau _1$ when particle 2's mass-to-area ratio is much larger than that of particle 1 (Fig. 1). Thus, in the $\tau_1, \tau_2 < t_{\rm s}$ regime particles with very different A/m-ratios will preferentially collide, because differential velocities are then highest. When one of the particles' friction time exceeds $t_{\rm s}$ the dependence on the absolute difference in $\tau_{\rm f}$ vanishes and the relative velocities now scale with the square root of m/A of the largest $\tau_{\rm f}$. As can be seen in Fig. 1 relative velocities are rather insensitive to the precise ratio of the particles' friction times in this regime. (Because of the simplicity of the expressions for $\Delta {v}$ in Eqs. ((5), (5)) the lines do not connect at $\tau_1 = t_{\rm s}$ and $\tau_1 = t_0$.) In the $\tau_{\rm f}>t_0$ regime, (Eq. (5)) relative velocities would stop increasing and eventually become only a minor perturbation to the motion of the particle. These large $\tau_{\rm f}$ regimes are, however, not reached in the early stages of coagulation considered in this paper. Here, it is clear that the relative grain velocities are regulated by the area-to-mass ratio of the colliding grains which sets the friction timescale (Eq. (4)) and hence the velocities (Eqs. (5), (5)).

  
2.2 Porosity of agglomerates

For compact, solid particles, the area-to-mass ratio scales linearly with the size. However, for porous aggregates the A/m ratio depends on the internal structure of the aggregates. In this section, we discuss how the internal structure of the aggregates affects collisions and, conversely, how collisions affect the aggregates' internal structure. The internal structure is modelled using the enlargement parameter. Here, we define an enlargement parameter $\psi$ that is the ratio between its extended volume, V, and its compact volume, V^*, i.e.,

$$\psi = \frac{V}{V^\ast}, \qquad \psi \ge 1.$$ (6)

V^* is the volume the material occupies in its compacted state, i.e., when particles are packed most efficiently, and V the total (extended) volume of the particle. While V^* reflects the mass of the particle, V is related to the spatial extent of the aggregate. Here, V has been defined as the volume corresponding to the geometric cross-section of the aggregate (see Fig. 2). Figure 2 shows three aggregates of 1000 monomers, such that V^* is the same, though with different degrees of compaction. The (from left to right) decreasing compaction translates to an increased geometric cross-section, A (outer circle), of the aggregates and hence to an increased $\psi$. Using the linearity between V^* and m and Eq. (6), two parameters, e.g., m and $\psi$, determine A and, consequently, also determine the coupling with the gas (Eq. (4)).

$\psi =1$ then corresponds to the enlargement factor of a compact particle with specific density $\rho_{\rm s}=m/V^*$. This does not necessarily correspond to a homogeneous particle (of zero porosity); it corresponds, however, to the lowest porosity that can be achieved by re-arranging the constituent particles (monomers) within the aggregate, without destroying this substructure through, e.g., melting. Since $\psi \ge 1$ we refer to $\psi$ as the enlargement parameter or enlargement factor. $\psi$ is related to the structural parameter, x, of Ossenkopf (1993), as $x \propto \psi^2$. Note, the close relationship between $\psi$ and porosity; a higher $\psi$ means more open aggregates, i.e., higher porosity. Therefore we will also use this more familiar reference for $\psi$, but only in a qualitative sense.

Since our enlargement parameter, $\psi$, is defined in order to match the geometrical cross-section, A, of a particle, we refer to the resulting radii, $a(\psi)$, as geometrical. The other cross-section of importance in coagulation models is $\sigma$, the collisional cross-section. Ossenkopf (1993) has consequentially defined a "toothing radius'', $a_{\rm tooth}$, such that $\sigma = \pi (a_{\rm tooth, 1} + a_{\rm tooth, 2})^2$ and provides expressions for $a_{\rm tooth}$ for PCA and CCA aggregates. $\sigma$ is larger than A by a factor of order unity (see also Krause & Blum 2004) and also depends slightly on the structure of the aggregates, i.e., whether PCA or CCA. In this paper, however, we are mainly concerned with obtaining a model for $\psi$ and have therefore simply used a for the calculation of the collisional cross-section. Likewise, we have also ignored augmentations of $\sigma$ due to effects such as, e.g., charges and grain asymmetries (see again Ossenkopf 1993) and rotation of grains (Paszun & Dominik 2006).

The aggregates' internal structure has important consequences for the coagulation rate. The geometric cross-section scales, for example, as $\psi^{2/3}$. More subtle are the effects of $\psi$ on relative velocities, which, as seen above, depend critically on the A/m-ratio of both aggregates. In coagulation models where grains are represented as compact bodies $\tau_{\rm f}$ increases linearly with size ( $A\propto m^{2/3}; \tau_{\rm f} \propto m^{1/3}$). However, in the initial stages of coagulation aggregates stick where they meet - a process characterised by a build-up of porous, fluffy structures, which can be described by fractals. Meakin & Donn (1988) have computed the A/m ratio of an initially monodisperse population and found a profound flattening of this ratio as compared to compact models in which it decreases as m^-1/3. Often, in the "hit-and-stick'' regime the growth shows fractal behaviour and the cross-section can be directly parameterised as a power-law, i.e.,

$$A \propto m^\delta, \qquad \frac{2}{3} \le \delta \le 1.$$ (7)

The lower limit $\delta=\frac{2}{3}$ corresponds to the growth of compact particles, whereas the upper limit of $\delta=1$ denotes the aggregation of chains or planar structures. In the cluster-cluster coagulation (CCA) model of Ossenkopf (1993) $\delta=\delta_{\rm CCA}=0.95$, a result that agrees well with the findings of Kempf et al. (1999), using N-particle simulations in the Brownian motion regime. More recent models, which also include rotation of aggregates (Paszun & Dominik 2006) also confirm this exponent. The scaling of friction time and enlargement factor with mass then become

$$\tau_{\rm f} \propto m^{1-\delta}, \qquad \psi \propto m^{\frac{3}{2}\delta -1},$$ (8)

such that for compact aggregation models ( $\delta=2/3$) $\psi$ stays constant. On the other hand, in models of cluster-cluster aggregation ( $\delta \approx 1)$, area scales roughly as mass and $\tau_{\rm f}$ stays constant or is only weakly dependent on mass, while $\psi$ increases monotonically. Thus, in CCA relative velocities are rather insensitive to growth. As a result, collisions are also less energetic in CCA models, e.g., as compared to compact aggregation models.

The key to the coagulation process in protoplanetary disks is the coupling of the dust to the turbulent motions of the gas and the resulting velocity distribution. In essence, the enlargement parameter $\psi$ provides a relationship between mass and surface area for growing aggregates which controls this gas-dust coupling. Equation (8) provides a relation for the evolution of $\psi$, but this relation is only valid in certain specific aggregation models, e.g., CCA-coagulation, where similar, equally sized aggregates meet. In reality, however, collisions between particles of all kinds of sizes will occur, although, dependent on the parameters that determine $\Delta {v}$, some are just more probable than others. In the end, the growth of grains in protoplanetary disks is controlled by the individual collisions between two aggregates. Therefore, we have to provide prescriptions for the outcome of all possible collisional encounters, i.e., all relevant combinations of $m, \psi$ and  $\Delta {v}$.

  
2.3 The collision model

We consider a collision between two particles with the aim of applying the results to a true coagulation model. Essentially, we have to provide a recipe for the enlargement factor, $\psi$, after the collision. Two relevant parameter sets enter into the new $\psi$: the progenitor masses and enlargement factors (i.e., the $(m_i, \psi_i)$ of the colliding particles). Moreover, the collision energy,

$$E = \frac{1}{2}\frac{m_1m_2}{m_1+m_2} \Delta {v}^2 = \frac{1}{2} \mu \Delta {v}^2,$$ (9)

with $\mu$ the reduced mass, is of key importance. At very low velocities, collisions between two aggregates will lead to sticking without internal restructuring, i.e., the particles will stick where they make first contact. At moderate velocities, the internal structure of the resulting aggregate will adjust - dissipating a major fraction of the kinetic collision energy - resulting in a compaction of the aggregates. Finally, at very high collision velocities, the colliding aggregates will fragment into smaller units and this can lead to substantial erosion. Following the numerical model of Dominik & Tielens (1997) and the experimental studies by Blum & Wurm (2000), these collisional regimes are distinguished by the following critical (collision) energies:

• $E_{\rm restr} = 5 E_{\rm roll}$, the energy needed for the onset of compaction;
• $E_{\rm max-c} \simeq N_{\rm c} \cdot E_{\rm roll}$, the energy at which aggregates are maximum compressed. Here, $N_{\rm c}$ is the total number of contact surfaces (between monomers) of the agglomerate. For a very open, fluffy agglomerate, $N_{\rm c} = N$. With increasing compaction the number of contacts will increase and for a cubic packing $N_{\rm c}=3N$. Here, for simplicity, we will adopt $N_{\rm c} = N$;
• $E_{\rm frag} \simeq 3 N_{\rm c} \cdot E_{\rm break}$, the energy needed for (the onset of) fragmentation of the agglomerate. Here $E_{\rm break}$ is the energy required to break a bond between two monomers. Its magnitude is of similar order as $E_{\rm roll}$.

For monomers of the same size $E_{\rm roll}$ is given by (Blum & Wurm 2000; Dominik & Tielens 1997)

$$E_{\rm roll} = 3\pi^2 \gamma a_0 \xi_{\rm crit} = \frac{1}{2} \pi a_0 F_{\rm roll},$$ (10)

with a₀ the radius of the monomer, $\gamma$ the specific surface adhesion energy and $\xi_{\rm crit}$ the yield displacement at which rolling commences. $F_{\rm roll}$, the rolling force, was measured by Heim et al. (1999) to be $F_{\rm roll} = (8.5\pm1.6)\times10^{-5}\ {\rm dyn}$ for uncoated ${\rm SiO}_2$-spheres of surface energy density $\gamma = 14\pm2\ {\rm erg~ cm}^{-2}$. We adopt this value for $F_{\rm roll}$ and assume proportionality with $\gamma$ when applying it to other materials. The monomer size, a₀, is also an important model parameter directly affecting the outcome of the collision at a given mass; a smaller a₀ provides more contacts (for the same mass) and the agglomerate is more resistant to breakup. With these energy thresholds, three qualitatively different collision regimes can then be discerned (see Fig. 3):

Figure 3: The collision regimes as function of total particle mass and relative velocity. Thick dashed lines indicate the critical energies for the onset of rolling and fragmentation. Parameters are that of quartz particles ( $\rho _{\rm s}=3.0\ {\rm g\ cm^{-3}}, \gamma =25.0\ {\rm ergs\ cm^{-2}}, a_0=0.1\ \mu {\rm m}$) and for the $E_{\rm frag}$ line we have assumed equal masses (m1=m2) and $E_{\rm break} = E_{\rm roll}$. Arrows indicate how the critical energy lines shift with increasing monomer size and mass-ratio, m1/m2.

Open with DEXTER

i)

$E < E_{\rm restr}$: no internal restructuring. The particles stick where they meet ("hit-and-stick'' growth) as in traditional, lattice-based, aggregation models (e.g., Meakin 1988) - a process leading to fractal aggregates.

ii)

$E_{\rm restr} < E < E_{\rm frag}$: restructuring (compaction) of aggregates.

iii)

$E > E_{\rm frag}$. Initially, loss of monomers, but for high energies complete fragmentation (e.g., catastrophic collision). This phase requires a recipe for the mass distribution of the fragments. In this paper, fragmenting collisions are avoided by, e.g., "choosing'' moderate $\alpha _{\rm T}$ and stopping coagulation for particles that reach a critical friction time (see Sect. 4). Within these constraints, our results are therefore not compromised by ignoring the fragmentation issue.

From Fig. 3 it is clear that, when starting with small particles, growth will initially be in the fractal regime. This fractal growth will be followed by a period in which the collisions will promote the restructuring of the growing agglomerate. Fragmentation becomes only important for velocities in excess of $10^3\ {\rm cm\ s^{-1}}$. Since we assume that every contact can absorb a unit energy  $E_{\rm roll}$, the $E_{\rm frag}$ line in Fig. 3 is independent of total mass. The $\Delta {v} \simeq 10^3\ ~{\rm cm\ s^{-1}}$ limit translates to a critical value for the turbulent $\alpha _{\rm T}$ parameter: $\alpha_{\rm T} \lesssim 10^{-2}$. Within Fig. 3, the precise "growth path'' of the agglomerate will be controlled by evolution of the relative velocities and hence A/m or, equivalently, $\psi$. We will now construct recipes for $\psi$ in, respectively, i) the fractal and ii) the compaction regime.

  
2.3.1 ${E}<{E}f{_{\rm restr}}$: hit-and-stick

Besides the usual CCA and PCA formalisms, there have been a few attempts to give prescriptions for the evolving internal structure of aggregates in the hit-and-stick regime. Kostoglou & Konstandopoulos (2001) discuss a formalism for obtaining the new fractal dimension in terms of the sizes and fractal dimensions of the two colliding aggregates. One point is, however, that, apart from the fractal dimension, another parameter - the prefactor - is needed to fully describe the fractal, although it is usually of order unity. Ossenkopf (1993), like this study, introduces only one structural parameter and interpolates between the CCA and PCA limits. We will also follow this idea, but use a different interpolation mechanism.

We recognise that in the pure sticking regime most collisions are between evolved, fluffy aggregates, since the size distribution evolves progressively toward larger sizes. For low velocity-mass combinations (Fig. 3), where restructuring is unimportant, the collisional growth then resembles the CCA growth process the most. We therefore simply rewrite the fractal law in terms of the individual masses of the particles and keep the CCA exponent,

$$A = A_1 \left( \frac{m_1+m_2}{m_1} \right)^{\delta_{\rm CCA}}, \qquad m_1 > m_2.$$ (11)

Although collisions between different particles are included in Eq. (11), we still adopt the CCA-characteristic exponent ( $\delta_{\rm CCA} = 0.95$) to ensure that for the "pure'' CCA case ( $m_1=m_2;\ A_1=A_2$) this prescription is in accordance with detailed numerical studies (Paszun & Dominik 2006; Kempf et al. 1999; Ossenkopf 1993; Meakin & Donn 1988). There is, however, a modification to Eq. (11) that must be made. The term in brackets in Eq. (11) determines the amount of increase in A in the fractal regime. Because fractal growth results from inefficient packing of voluminous objects, it is clear that this term must include parameters describing the spatial extent of the collision partners. These cannot, however, be given by the masses of the particles, since m (alone) does not reflect a spatial size. For example, if we would replace one of the aggregates by one of the same mass, but of lower porosity (i.e., a more compact aggregate), its volume is obviously smaller and packing becomes more efficient. These effects, however, are not reflected in Eq. (11). For these reasons, we replace m by V in Eq. (11),

$$A = A_1 \left( \frac{V_1+V_2}{V_1} \right)^{\delta_{\rm CCA}}, \qquad V_1 > V_2.$$ (12)

Note that for particles of the same internal density (porosity) Eqs. (12) and (11) agree, such that Eq. (12) also can be seen as an extrapolation from the CCA case, but one that takes account of the different internal structures of the collision partners. Using the relation $A\propto V^{2/3}$ and Eq. (6), Eq. (12) can be re-written in terms of m and $\psi$ only

$$\psi = \langle \psi \rangle_m \left( 1 + \frac{m_2\psi_2}{m_1\psi_1} \right)^{\frac{3}{2}\delta_{\rm CCA}-1},$$ (13)

with $\langle \psi \rangle_m$ the mass-averaged enlargement factor of the collision partners,

$$\langle \psi \rangle_m \equiv \frac{m_1 \psi_1 + m_2 \psi_2}{m_1 + m_2}\cdot$$ (14)

In CCA coagulation (m₁=m₂ and $\psi_1 = \psi_2$) we recover Eq. (8), but Eq. (13) now includes all collisions in the hit-and-stick regime. For example, if a large, fluffy aggregate sticks to a compact one, the enlargement factor of the newly formed aggregate is higher than the mass-averaged enlargement factor of the progenitor particles, $\langle \psi \rangle_m$, but smaller than that of the fluffy collision partner. In Sect. 2.4, it will be shown, however, that the $\langle \psi \rangle_m$ term underestimates the porous growth when one of the particles is very small, i.e., in PCA-like collisions. This is solved by adding to Eq. (13) a term that compensates for these cases and our final recipe then becomes

$$\psi = \langle \psi \rangle_m \left( 1 + \frac{m_2\psi_2}{m_1\psi_1} \right)^{\frac{3}{2}\delta_{\rm CCA}-1} + \psi_{\rm add},$$ (15)

where $\psi _{\rm add}$, a term important only for small m₂, is explained in Sect. 2.4.

  
2.3.2 ${E}>{E}_{\rm restr.}$: compaction

In the compaction limit monomers are restructured at the expense of the porous volume of the aggregates. Following the theoretical study of Dominik & Tielens (1997) we will assume that the (relative) amount of compaction, $\Delta V_{\rm p}/V_{\rm p}$, scales linearly with collisional energy, i.e., $\Delta V_{\rm p}/V_{\rm p} = -f_C = E/E_{\rm max-c} = -E/(N \cdot E_{\rm roll})$, where $V_{\rm p} = V - V^*$ denotes the porous volume within the aggregate and N is the total number of monomers present. Essentially, this implies that the collision energy is used to set individual monomers in an agglomerate rolling and that this rolling is only stopped when an additional contact is made, resulting in compaction. Recalling that $V = \psi V^*$ and $V_{\rm p} = V - V^* = (\psi - 1) V^*$ with V^* proportional to mass, we then find that the porous volume after colliding is

$$(V_1^* + V_2^*) (\psi - 1) = (1 - f_C) \left[ V_1^* (\psi_1 -1) + V_2^* (\psi_2 -1) \right].$$ (16)

And the new enlargement factor

 

$$\psi-1 (1-f_C) \cdot \frac{1}{m_1+m_2} \left( m_1 (\psi_1-1) + m_2 (\psi_2-1) \right)$$ =

$$(1 - f_C) \left( \langle \psi \rangle_m -1 \right),$$ = (17)

with $\langle \psi \rangle_m$ again the mass-averaged enlargement factor of the (two) collision partners (Eq. (14)). We illustrate Eq. (17) in Fig. 4 for the limiting cases of equal porosity collisions ( $\psi_1 = \psi_2$; black lines) and very porous vs. compact particle collisions ( $\psi_1 \gg \psi_2$; grey lines). Higher mass-ratios give higher collision energies (higher $\mu$) and hence more compaction. If velocities are low, ( ${v} < 100\ {\rm cm\ s^{-1}}$) the net compaction occurs primarily through the $\langle \psi \rangle_m$ factor and the curves converge on the $0\ {\rm cm\ s^{-1}}$ (thick) line. Only when ${v} > 100~ {\rm cm\ s^{-1}}$ does the f_C-factor start to become important. Collisions at velocities higher than $1000\ {\rm cm\ s^{-1}}$ can result in fragmentation if the mass-ratios are similar. For low mass ratios the $\langle \psi \rangle_m$ is determined by $\psi _1$ (the highest mass) and there is little difference between the two limiting cases.

Figure 4: Relative compaction as a function of the mass ratio, m2/m1, at various collisional velocities. The expression on the y-axis measures the compaction relative to particle 1 and is a function of mass ratio, m2/m1, only (see Eq. (17)). Here $\psi _1$ is the enlargement factor of the most massive aggregate, i.e., m1 > m2, and plots are shown for $\psi _2 \ll \psi _1$ (grey lines) and $\psi _2 = \psi _1$ (black lines). At low mass ratios the curves converge.

Open with DEXTER

  
2.4 Porosity increase in the PCA and CCA limits

Figure 5: The enlargement factor in the PCA and CCA limits as a function of the number of monomers (N) of the particle. The thick grey lines show the fits of Ossenkopf (1993) for CCA ( top) and PCA ( bottom) coagulation. The solid line corresponds to the limiting cases of Eq. (13). The dashed lines show the same limits, but now the zero point lies at N=40 (after which the PCA/CCA curves of Ossenkopf (1993) diverge in porosity) and includes the $\psi _{\rm add}$ correction term (Eq. (19)).

Open with DEXTER

In the discussion on the $\psi$ recipes in Sect. 2.3.1 we have used the CCA fractal exponent ( $\delta_{\rm CCA} = 0.95$) as the starting point and extrapolated this empirical relation to Eq. (13): a general recipe for all collisions in the "hit-and-stick'' regime. By definition, CCA-collisions are incorporated into this recipe and we may expect Eq. (13) also to account for collisions between aggregates having about the same size. The PCA model, where one monomer collides with an agglomerate, is the opposite extreme. If we take the PCA-limit of Eq. (13), i.e., $m_1 \gg m_2 \sim m_0$, with m₀ the monomer mass, the change in $\psi$ after addition of a monomer becomes

$$\psi_{N+1} \approx \psi_N + \frac{N_2}{N}\left( \frac{3}{2}\delta_{\rm CCA}\psi_2 - \psi_N \right); \qquad N\gg N_2,$$ (18)

with N=N₁=m/m₀ the number of monomers of the PCA agglomerate and N₂=1 for monomers. Thus, $\psi_N$ goes to $\frac{3}{2} \delta_{\rm CCA} \psi_2 \approx 1.5$ in the PCA asymptotic limit ( $N_2 =1; \psi_2=1$). The fact that an asymptotic limit is reached can be understood intuitively, since there must be a point at which the inward penetration of monomers into the centre of the aggregate, which decreases $\psi$, starts to balance the porous growth due to hit-and-stick collisions at the surface. The asymptotic limit of $\psi \approx 1.5$, however, is much lower than typical PCA models indicate ( $\psi \approx 10$) as is illustrated in Fig. 5, where the PCA/CCA limits of our model (solid lines) are compared to detailed numerical simulations of Ossenkopf (1993) (thick lines). Equation (13) thus underpredicts the porous growth for PCA-like collisions in which one of the particles is small; a result not too surprising since it originates from the CCA fractal law (Eq. (11)), which is constructed to apply only for similar (i.e., equal-sized) particles. For these reasons we add to Eq. (13) a term that compensates the $-N_2 \psi_N / N$ term in Eq. (18),

$$\psi_{\rm add} = B \frac{m_2}{m_1} \psi_1 \exp \left[ -\mu/ m_{\rm F} \right],$$ (19)

where the exponential ensures $\psi _{\rm add}$ is unimportant for collisions between particles well above a certain mass-scale, $m_{\rm F}$. With B=1.0 and $m_{\rm F}=10\ m_0$ we find good correspondence with the results of Ossenkopf (1993). In Fig. 5 the new CCA and PCA limits are shown by the dashed curves, where we have shifted the "zero point'' from N = 1 to N = 40, i.e., the starting point for $\psi$ after which we recursively apply Eq. (15). The transition toward fractal behaviour emerges only after this point and we therefore directly take $\psi$ from the results of Ossenkopf (1993) when $N \lesssim 40$.

Although, with Eq. (19) for $\psi _{\rm add}$, Eq. (15) does achieve the right PCA/CCA fractal limits, we do not claim they actually provide a model for $\psi$. The collision recipes are based on empirical findings and extrapolations from these. However, in contrast to the CCA and PCA limiting collisional growth models, where $\psi$ can be directly parameterised in a single exponent, we recognise that the internal structure of the aggregates is changed during - and caused by - the collisional growth process. This is the main qualitative difference of our collision model captured in Eq. (15). This equation, together with Eq. (17), provides recipes for the collisional evolution of the enlargement factor, which can be easily incorporated into time-dependent coagulation models. We acknowledge this is an active area of research and future model efforts may well improve on the present formulation.

  
3 Monte Carlo coagulation

3.1 Outline

Figure 6: Illustration of the adopted Monte Carlo technique. (A) N-particles are assumed to be uniformly distributed in a box of volume ${\cal V}$. The collision rate between particle i and j is $C_{ij} = \sigma _{ij} \Delta {v}_{ij} / {\cal V}$. A random number determines which two particles will collide. (B) After the collision the number of particles is restored by duplicating one of the remaining N-1 particles. Physically, this can be interpreted as an expansion of ${\cal V}$ (heavily exaggerated in this figure) to a volume at which ${\cal V}$ contains N-particles again. The hypothesis of this method is that the collisional evolution within ${\cal V}$ is representative for the coagulation of the total space under consideration.

Open with DEXTER

To determine the implications the new collision model has for the solar nebula, e.g., as compared to collision models where mass is the only parameter, it must be embedded in a coagulation model that evolves the particle distribution function, $f(\vec{x})$. $f(\vec{x},t)$ gives the number density of particles with a set of properties (parameters) $\left\{ x_i \right\}$ at time t. In compact coagulation models all properties depend only on mass, so $f(\vec{x},t)=f(m,t)$. In the model described in Sect. 2, however, the particle's enlargement factor has been included as an independent parameter such that f becomes a function of three variables, $f(\vec{x},t) = f(m,\psi,t)$. The coagulation equation which describes the evolution of  $f(\vec{x})$ is

 

$$\frac{\partial f(\vec{x},t)}{\partial t} = - f(\vec{x},t) \int {\... ...delta_{\rm k} \left( \vec{x} - \Gamma \left(\vec{x'},\vec{x''} \right) \right),$$ (20)

with $K = \sigma \Delta {v}$ the collision kernel, $\delta_{\rm k}$ the Kronecker $\delta$-function and $\Gamma$ the collision function, which maps in the case of sticking 2n parameters (those of $\vec{x'}$ and $\vec{x''}$) to n, where n is the number of independent parameters with which a particle is characterised. Equation (20) of course is just an extension of the Smoluchowski equation (Smoluchowski 1916) to more than one dimension. Applied to the formalism in Sect. 2 it describes the loss of particles in state $\vec{x} = (m,\psi)$ due to collisions with any other particle (first term) and the gain of "$\vec{x}$-particles'' that happen to be formed out of any suitable collision between two other particles (second term). Applied to the findings in Sect. 2, $\Gamma$ symbolises the collision recipes with $\Gamma(m_1,\psi_1;m_2,\psi_2) = (m_1+m_2, \psi)$ and $\psi$ is given by Eq. (15) or Eq. (17), dependent on the collisional energy.

One approach to implement coagulation is to numerically integrate the Smoluchowski equation. However, it is immediately clear that numerically integrating Eq. (20) becomes a daunting exercise. Integrating the ordinary (1-dimensional) Smoluchowski equation is already a non-trivial matter. Problems of near cancellation (the gain terms often equal the loss terms), mass conservation (systematic propagation of errors) and the problem concerning the determination of a time-step must all be tackled. (Dullemond & Dominik (2005) in their Appendix B give an overview of the subtleties involved.) The $\Gamma$-factor in Eq. (20) gives a further complication since there is no such thing as "conservation of porosity'' and the $\delta$-factor cannot be easily integrated away. Although there is no fundamental reason against the binning method - see, e.g., Ossenkopf (1993) who solves the Smoluchowski equation in two dimensions - these issues make the whole procedure quite elaborate. We felt that much of the simplicity of the collision model of Sect. 2 would be "buried'' by numerical integration of a 2d-Smoluchowski equation and therefore have found it suitable to opt for an approach using direct Monte-Carlo simulation (DSMC) techniques.

The simplicity of using Monte-Carlo methods for coagulation problems is appealing. Basically, N-particles are distributed over a volume ${\cal V}$ (see Fig. 6A). The evolution then boils down to the determination of the two particles which are involved in the next collision. We hereafter assume that the particles are well mixed, i.e., no potential is present, such that the determination of the next collision is governed by basic stochastic principles. Then the probability of a collision between particles i and j is given by the collision rate, $C_{ij} = K_{ij} / {\cal V}$ in which K_ij is the collision kernel. A random number determines which of the $\frac{1}{2}N(N-1)$ possible collisions will be the next. The collision then creates a new particle, after which the $\{C_{ij}\}$ must be updated and the procedure repeats itself.

The advantages of such an approach are obvious. Most striking perhaps is the "physical character'' of Monte-Carlo simulations. The growth-evolution of individual particles is directly traced and can be studied. The algorithm does not use the distribution function, f, in a direct way; it is obtained indirectly by binning the particles. Secondly, the above described method is exact, i.e., no "time vs. accuracy'' considerations have to be made in choosing the time-step $\Delta t$; instead, $\Delta t$ - the inter-collision time - is an outcome of the stochastic coagulation process as it is in nature. Furthermore, due to its stochastic nature, no Monte-Carlo simulation is the same. A series of (independent) runs gives at once a measure of the statistical spread in the distribution. Note that the fluctuations around the average are a combination of real stochastic noise and random noise, but it is qualitatively different from the Smoluchowski approach, which describes the evolution of the mean of all possible realisations and is therefore completely deterministic (Gillespie 1977). From a practical point of view the straightforwardness of the DSMC-method makes that there is no need for resorting to "control parameters'' like those required in the numerical-integration method.

The DSMC-method, however, has its limitations. It can be immediately seen, for instance, that when having started with N (say identical) particles of mass m₀ in a fixed volume, these will over time pile up in one agglomerate of $m_{\rm final}=Nm_0$. The accuracy then steadily decreases during the simulation (in MC-simulations the statistical error scales proportional to N^-1/2) and most of the computing time is spent in the first few (quite uninteresting) collisions. Consequently, to achieve orders of magnitude growth, the initial number of particles must be extremely large. And because the calculation of the $\{C_{ij}\}$ goes proportional to N² (every particle can collide with each other), it becomes clear that this method becomes impracticable. To counter the dependence on large initial particle numbers, one can also try to preserve the number of particles during the simulation. This can be done, for instance, by "tossing-up'' the particle-distribution when the number of particles reaches N/2 as described by Liffman (1991). A more natural approach perhaps, given by Smith & Matsoukas (1998), and adopted here, keeps the number of particles constant at each step; every time a collision takes place one of the remaining N-1 particles is randomly duplicated such that the number of particles throughout the simulation stays the same. This procedure can graphically be represented as an increase in the simulated volume, ${\cal V}$ (see Fig. 6B), under the assumption that the collisional evolution outside of ${\cal V}$ progresses identically. Smith & Matsoukas (1998) have shown that the error in the coagulation-scheme now scales logarithmically with the extent of growth, - or growth factor (GF), defined here as the mean mass over the initial mean mass of the population - much improved over the constant-${\cal V}$ case, where the error has a square-root dependence on GF. We might worry though about the consequences of the duplication mechanism. It causes a certain degree of "inbreeding'', which effects we cannot quantify directly. Smith & Matsoukas (1998) show it is unimportant for the constant or Brownian kernels used in their studies. However, these kernels are known to behave gently, i.e., they are rather insensitive to irregular changes in the population since the variations in the K_ij are small. Perhaps, more erratic kernels are more sensitive to the "duplication mechanism'', but we might equally well attribute this sensitivity to the stochastic nature of the coagulation process. At any rate, these consequences are best quantified by running the code multiple times.

  
3.2 Implementation

In implementing the DSMC approach we follow the "full conditioning method'' of Gillespie (1977). This involves the calculation and updating of partial sums, C_i,

$$C_i \equiv \sum_{j=i+1}^{N} C_{ij}, \qquad i=1,\dots,N-1,$$ (22)

with $C_{ij} = K_{ij}/{\cal V} = \sigma_{ij}(a_i, a_j) \Delta {v}_{ij} (\tau_{i}, \tau_{j}) / {\cal V}(\kappa)$ (here $\kappa$ is the total number of collisions since the start of the simulation). These N-1 quantities are stored in the memory of the computer. $C_{\rm tot} = \sum_i C_i$ is the total coagulation rate. The probability density function, P(t,i,j), i.e., the probability that the next collision will occur in time-interval $(t,t+{\rm d}t)$ and involves particles i and j (i<j) can then be written as (Gillespie 1977)

$$C_{ij} \exp \left[ -C_{\rm tot} t \right]$$                      P(t,i,j) =

$$\Big( C_{\rm tot} \exp \left[ -C_{\rm tot} t \right] \Big) \times \Big( C_i/C_{\rm tot} \Big) \times \Big( C_{ij}/C_i \Big).$$ = (23)

Three random deviates, $r_i \in [0,1]$, then determine successively: I) the time it takes until the next collision takes place, $t = C_{\rm tot}^{-1} \ln (1/r_1)$; II) the first particle (i) to collide, by summing over the C_i-s (starting with i=1) until $r_2 C_{\rm tot}$ is exceeded (this fixes i); and III) its collision partner (j) by summing the C_ij-s over the j-index (starting with j=i+1) until the value r₃ C_i is exceeded (Gillespie 1977, Eq. (19)). The outcome of the collision is evaluated using the relevant equations in Sect. 2 and the new particle is stored in the i-slot. Another random number then determines which of the N-1 particles (excluding j) will be duplicated and this one is stored in the j-slot. Having created, removed and duplicated particles, all of the C_i-s need to be updated. This implies only the subtraction/addition of the C_ij-s that have changed, not the re-computation of Eq. (22). Moreover, the `duplication procedure' entails a rescaling of the simulated volume, ${\cal V}$, such that the spatial density of solids, $\rho_{\rm d} = \sum_i m_i /{\cal V(\kappa)}$ remains constant. The algorithm can then be repeated. All these steps are order-N calculations at worst; most time-consuming are the determination of j (for which C_ij has to be calculated) and the update of the $\{ C_i \}$. To achieve a given GF another factor N in computation time is needed, bringing the total CPU-time proportional to N². These procedures are graphically summarised in Fig. 7.

	Figure 7: Flow chart of the MC-coagulation method. One cycle corresponds to one collision. "Func'' like in $i={\rm Func}(\ \{C_{i}\}, r\ )$ indicates that i is a function of the Ci values and a random deviate, r. The procedure is further explained in the text.
Open with DEXTER

It is possible, however, to save some CPU time by taking the duplicates together in the calculation of the collision rates. If there are (g_i -1) copies of particle i, it would be a waste of time to calculate the (same) rates g_i times. Rather, g_i can be absorbed in the calculation of the (combined) coagulation rate. $\tilde{C}_{ij} = g_i g_j C_{ij}$ is then the rate at which one of the i-particles collides with one of the j-particles ($i\neq j$) and $\tilde{C}_{ii} = \frac{1}{2}g_i(g_i-1) C_{ii}$ between duplicates. The CPU time per step is now proportional to $N_{\rm s}$, the total number of distinct particles, i.e., excluding duplicates, with $\sum_{i=1}^{N_{\rm s}} g_i = N$. To think of it in biological terms: g_i gives the multiplicity of species i; $N_{\rm s}$ the total number of species; and N the total number of living creatures.

3.3 Tests

The Monte-Carlo coagulation model described above is tested with kernels that have an analytic solution. These are I) the constant kernel, K_ij=1, and II) the linear kernel, $K_{ij}= \frac{1}{2}(m_i+m_j)$. The evolution of the mean mass of the distribution, $\langle m \rangle$ for these kernels is (Tanaka & Nakazawa 1994; Silk & Takahashi 1979)

$$\langle m \rangle = \left\{ \begin{array}{ll} m_0(1+\frac{1... ...\right] & \qquad {\rm linear\ kernel}, \\ \end{array}\right.$$ (25)

where the distribution at t=0 is monodisperse of mass m₀. Well-known coagulation models have either $K \propto m^{1/3}$ (Brownian coagulation) or $K\propto m^{2/3}$ (geometric area), but here, due to changes in $\psi$ and $\Delta {v}$, we should be prepared for K to vary with time. Thus, it is important that the Monte-Carlo code passes both these tests. Initial conditions for these test cases are monodisperse with all relevant parameters put at 1 at t=0 (i.e., m₀=1 and f(t=0)=1, $\rho_{\rm d}=1$, $\rho_{\rm s}=1$) and we do not take porosity into account ($\psi =1$ always). At various times the particles are binned by mass and the distribution function f(m) is determined by summing over the masses in the bin and dividing by the width of the bin (to get the spectrum) and the volume of the simulation (to get the density). Multiple runs of the simulation then provide the spread in f. Figures 8 and 9 present the results. On the y-axis f(m) is multiplied by m² to show the mass-density per logarithmic bin. Analytical solutions (Tanaka & Nakazawa 1994) are overplotted by solid curves, while the dotted line shows the (hypothetical) distribution function if all the bins would be occupied by 1 particle. Thus, the dotted line corresponds to $m^2f(m) = m^2/ {\cal V} w_{\rm b} \approx m/{\cal V}$, because the widths of the bins, $w_{\rm b}$, are also exponentially distributed. In a single simulation, the distribution function should lie above this (auxiliary) line and the vertical distance to this line is a measure for the number of particles in a bin.

Figure 8: Test of the Monte Carlo coagulation code - constant kernel (Kij=1 ), 20 000 particles. At (dimensionless) times t=1,102,104,106,108 particles are binned and the distribution function f is computed (symbols). The analytical solutions at these times are overplotted by the solid lines. The error bars (hardly visible) show the spread averaged over 10 runs. The dotted lines show the distribution function if all the bins would be occupied by only 1 particle - it is an auxiliary line with slope 1 showing the lower limit of the (individual) distribution function.

Open with DEXTER

Figure 8 shows that the code passes the constant kernel test with flying colours. The spread in the data is limited, and it does not noticeably increase with time. The linear kernel (Fig. 9A), on the other hand, shows a different story. Here the mean mass, $\langle m \rangle$, as well as the peak mass, $m_{\rm p}$ - defined as the peak of the m² f(m) size distribution - evolve exponentially with time. Note that the position of $m_{\rm p}$ only depends on the particles that contain most of the mass, while $\langle m \rangle$ is also sensitive to the total number of particles. Therefore, $\langle m \rangle$ lags $m_{\rm p}$ at any time and one can show that the gap between the two also increases exponentially with time. Inevitably, at some point in time, the theoretical value of the m²f(m) mass-peak becomes larger than the total mass present inside ${\cal V}$. This corresponds to the crossing of the dotted line at, e.g., $t\simeq 20$ in Fig. 9A. In other words, the duplication mechanism, needed to conserve N but which has the additional effect of enlarging ${\cal V}$, is incapable of keeping up with the exponential growth: "runaway particles'' could have formed, but the simulated volume ${\cal V}$ is just not large enough to take them into account. The postulate of the "duplication mechanism'' - the particle distribution evolves similarly in and outside ${\cal V}$ - then breaks down. The only way to avoid this effect is to enlarge ${\cal V}$ by having more particles in the simulation, i.e., to improve the "numerical resolution''. In Fig. 9B, we show the results, in which $N_{\rm s}$, instead of N, is held constant (Sect. 3.2). In these simulations N increases with time, starting with N=20 000 and ends with more than a million particles. The distribution now represents more closely the theoretical curve. Growth factors of 14 orders of magnitude in mass ($\simeq$5 in size) can then be accurately simulated.

Figure 9: Tests of the Monte Carlo coagulation code. Linear kernel $\left( K_{ij}= \frac{1}{2} (m_i+m_j) \right)$. - (A) fixed N (N=200 000). The distribution function f is computed at times t=(1,10,20,30). At t=30 insufficient mass is present to provide a good fit to the analytical distribution. - (B) fixed $N_{\rm s}$. To obtain better correspondence with theory, the number of particles, N, is increased as the simulation progresses, such that more volume is sampled. $N_{\rm s}$ is fixed at 20 000 and N ends with over a million particles. The average corresponds well to the theory, yet the amount of CPU time is disproportionally larger at the later times.

Open with DEXTER

The drawback, however, is that the computation slows down as N increases, since the relative increase in the average mass is inversely proportional to N, i.e., $\Delta \langle m \rangle/\langle m \rangle \propto N^{-1}$. These simulations therefore take much more CPU time. It is clear that a fundamental limit is reached, in which, given a certain CPU power, the calculation of the mass-distributions can only be achieved for a limited range. A way to overcome this problem is to collide multiple particles per event. In such an algorithm collisions are no longer between two particles but between two groups of particles. Although this approximates the collisional process, the coagulation can be speeded up by grouping especially small particles into a single unit. We will discuss the grouping algorithm and its implication in a future paper. For the simulations in Sect. 4, the product of N and $N_{\rm s}$ has been kept constant, which ensures a constant growth (in exponential terms) per cycle. We have fixed $\sqrt{N\times N_{\rm s}}$ at 20 000 but made sure that the numerical resolution issues as discussed here for the linear kernel, did not occur. Fortunately, realistic mass-distributions are not that broad as in the analytic, linear kernel; e.g., the smaller particles are quickly removed due to Brownian coagulation.

In summary, we have built an efficient Monte Carlo code to follow coagulation. The advantage of this code (above other numerical methods) is that it is intuitive, simple to implement and expand, and that it takes full account of the stochastic nature of coagulation. Minor disadvantages are the N² dependence on CPU time, the somewhat artificial duplication procedure, and the resolution problems shown for the linear kernel at high m as the simulation progresses. We have addressed these concerns and developed methods suitable for the scope of this work. DSMC-coagulation methods are very appropriate to work in conjunction with multi-parameter models. The collision model of Sect. 2, with mass, m, and enlargement parameter, $\psi$, as variables, can now be put into an evolutionary setting.

  
4 Results

4.1 Application to a protoplanetary disk

The collision model of Sect. 2 is embedded in the Monte Carlo code of the previous section and applied to the protoplanetary disk (Sect. 2.1). The coagulation is restricted to a turbulent environment in which particles are well mixed with the gas, yet do not fragment upon collision. For these reasons, we start out with a monodisperse population of sub-micron-sized grains (we choose them to be $a_0=0.1\ \mu{\rm m}$) present at $z=H_{\rm g}$, i.e., at the scaleheight of the disk, where the density is a factor $\exp [-0.5]$ lower than at the mid-plane. Its collisional evolution within the gaseous nebula is followed until the point that systematic motions, i.e., settling to the mid-plane, start to dominate. Thus, particles are either present and well mixed by turbulence, or have started to settle and are no longer in the region of interest. Settling is then modelled as a sudden phenomenon. The reality is, of course, a more gradual transition, but the discrete picture here is not a bad approximation, since the vertical structure is quickly established (Youdin & Chiang 2004). Settling occurs at the point when the friction time of a particle has exceeded a critical friction time, $\tau_{\rm rain}$, such that the scaleheight of the particle, $h(\tau_{\rm f})$, becomes smaller than that of the gas, i.e., $h(\tau_{\rm f}) < H_{\rm g}$. If self-gravity is neglected (valid in the gaseous nebula) and Schmidt numbers, Sc, are close to unity, then h can be obtained by equating the particle diffusion timescale, $t_{\rm diff} = h^2/\nu_{\rm T}$, with the particle settling timescale, $t_{\rm settl}=(\Omega^2 \tau_{\rm f})^{-1}$, i.e., $h(\tau_{\rm f}) = H_{\rm g} \sqrt{\alpha_{\rm T} / \Omega \tau_{\rm f} }$ ( cf. Schräpler & Henning 2004; Youdin 2005). The critical friction time is then

$$\tau_{\rm f} = \frac{\rho_{\rm s}}{c_{\rm g} \rho_{\rm g}} \f... ...\psi^{2/3}} = \tau_{\rm rain} = \frac{\alpha_{\rm T}}{\Omega},$$ (26)

with a^* the compact size of particles (see Fig. 2). As expected, higher values of $\alpha _{\rm T}$ and higher gas densities (lower $\tau_{\rm f}$) delay the onset of settling in the sense that the particle has to grow further in size before it arrives at the critical friction time. Alternatively, increasing $\psi$ also delays settling.

In the code, settling is implemented as a "rain-out'': the particle is removed from the simulation and the spatial dust density decreases. The evolution of these "rain-out particles'' during settling is not traced anymore. The focus stays on the particles that remain in the layers above the mid-plane. Their evolution is followed for $10^7~ {\rm yr}$.

The $\alpha _{\rm T}$ parameter is one of the major uncertainties concerning the characterisation of protoplanetary disks. One of the prime candidates for turbulence is the magneto-rotational instability (Balbus & Hawley 1991; Hawley & Balbus 1991), which seems to be most robust in well-ionised regions, i.e., in the upper layers of the disk. Another way to characterise $\alpha _{\rm T}$ is to relate it to the observed accretion rate, ${\rm d}M/{\rm d}t$, (Cuzzi et al. 2005) and then values in the range of $10^{-4} \lesssim \alpha_{\rm T} \lesssim 10^{-2}$ seem plausible. Yet, despite its uncertainty, $\alpha _{\rm T}$ appears in key expressions as, e.g., $\Delta {v}_{ij}$ and $\tau_{\rm rain}$. Therefore, models are run that cover a large range in $\alpha _{\rm T}$: $\alpha_{\rm T} = 10^{-6}{-}10^{-2}$. Furthermore, we divide the runs in two categories, which reflect the spatial position in the solar nebula. The `inner' models correspond to conditions at 1 AU where the monomers are quartz with internal density $\rho_{\rm s} = 3.0\ {\rm g\ cm}^{-3}$ and surface energy density of $\gamma = 25\ {\rm ergs\ cm}^{-2}$. The `outer' models correspond to conditions at 5 AU, where the coagulation is that of ices ( $\rho_{\rm s} = 1\ {\rm g\ cm}^{-3}$ and $\gamma = 300~ {\rm ergs\ cm}^{-2}$) and with an enhanced surface density of a factor 4.2 over the minimum solar nebula (Nakagawa et al. 1986). For comparison, we also run compact models for the $\alpha _{\rm T} = 10^{-4}$ runs. Compact models (denoted by the C-suffix in Table 1) are models where the internal structure does not evolve, i.e., $\psi =1$ by definition. An overview of all the models is given in Table 1.

Table 1: Overview of all the runs.

4.2 Particle growth and compaction

In Fig. 10 mass distributions are shown at various times during their collisional evolution. On the y-axis the mass function is plotted in terms of $m\cdot a^\ast \cdot f(a^*)$, which shows the mass-density, i.e., the mass of grains of compact size a^* in logarithmic bins. The panels compare the results of compact coagulation (panels A and B) with those of porous coagulation (panels C and D) for a turbulent strength parameter of $\alpha _{\rm T} = 10^{-4}$. The coagulation is calculated at 1 AU (quartz; panels A and C) and at 5 AU (ices, panels B and D). In Fig. 11 the $\alpha _{\rm T} = 10^{-4}$ model at 1 AU is compared to other $\alpha _{\rm T}$ models at 1 AU (see Table 1). In Fig. 13 averages of the size distributions of Figs. 10A, C are shown explicitly with time. Here $\langle a \rangle_m$ is the mass-weighted size,

$$\langle a \rangle_m = \frac{1}{\sum_i m_i} \sum_i m_i a_i,$$ (27)

of the population. Thus, while $\langle a \rangle$ gives the average particle size, $\langle a \rangle_m$ corresponds to the (average) size to which a unit of mass is confined. Because of this weighting, $\langle a \rangle_m \ge \langle a \rangle$, with the equality valid only for monodisperse distributions.

Figure 10: The mass function plotted at various times for the $\alpha _{\rm T} = 10^{-4}$ models. The panels compare the coagulation of the compact models ($\psi =1$, top panels) with those where porosity effects are included ( bottom panels). Left ( right) panels show the coagulation of quartz (ice) particles at 1 AU (5 AU). Each plot shows the mass function at every logarithmic interval in time from $t=10~ {\rm yr}$ until $t = 10^7~ {\rm yr}$. In the first $\sim$ $10^2\ ~{\rm yr}$ Brownian motion dominates the coagulation. Subsequent evolution is driven by turbulence-induced velocity differences and includes the moment of first compaction (blue, dotted line) and first rain-out (red, dashed line). After rain-out ( $t \gtrsim 10^4~ {\rm yr}$), the mass density in the gaseous nebula decreases and the mass function collapses. In the compact models the blue, dotted curve also corresponds to the first time that $E > E_{\rm roll}$. Greyscales indicate the spread in the 50 realisations of the simulation.

Open with DEXTER

Figure 11: Effects of turbulence on the coagulation. Panels show the collisional evolution at 1 AU of the porous models, yet with $\alpha _{\rm T}$ values of 10-2  A), 10-4 B) and 10-6  C). The scaling of the axis is the same throughout the panels. In the $\alpha _{\rm T}=10^{-2}$ models the spread in the runs is very large, causing the error bars to overlap. In the $\alpha _{\rm T}=10^{-6}$ model the particles rain-out without compacting.

Open with DEXTER

Figure 12: Evolution curves of various models. In these panels the enlargement factor is characterised by the ratio of the mass-weighted porous size, $\langle a \rangle_m$, over the mass-weighted compact size, $\langle a_{\rm comp} \rangle_m$, ( $\langle a^* \rangle_m$ in the text), against $\langle a_{\rm comp} \rangle_m$. Rising curves correspond to an increase of porosity due to fractal growth, and horizontal or declining lines indicate compaction. Numbers give the temporal stage of the coagulation (i.e., t=10i). The detached, coloured crosses indicate the average porous sizes (x-axis) and the size enhancement of the rain-out particles (y-axis).

Open with DEXTER

Table 2: Properties of the particle distribution at rain-out.

The porosity evolution is also displayed in Fig. 12, where we plot the ratio of $\langle a \rangle_m$ to $\langle a^* \rangle_m$. $\langle a^* \rangle_m$ is defined analogously to Eq. (27), but then for the compact particle size. $\langle a \rangle_m / \langle a^* \rangle_m$ then gives the mass-weighted averaged enhancement of the geometrical radius of the particles (see, Fig. 2). In the case of monodisperse populations or runaway growth one size dominates the average and the ratio is directly related to $\langle \psi \rangle_m$ as $\langle a \rangle_m / \langle a \rangle_m = \langle \psi \rangle_m^{1/3}$ This ratio is plotted vs. $\langle a^* \rangle_m$ in Figs. 12 for various models. It shows the build-up of porosity during the fractal stage, the stabilisation during the compaction stage, up to the stage where the particles rain-out. In Fig. 12 the average sizes of some of these rain-out particles are indicated by the detached crosses. Note that (for illustrative reasons) the x-axis for these particles corresponds to porous size, a, rather than compact size, a^*. Properties of the "rain-out'' population are given in Table 2. In Fig. 12 the temporal stage is indicated by numbered ticks.

Figure 13: Evolution of the size distribution with time. The mass-weighted size, $\langle a \rangle_m$ (solid line), and the mean size, $\langle a \rangle$ (dashed line), are calculated for the default models ( $\alpha _{\rm T} = 10^{-4}$) at 1 AU (top lines) and 5 AU (bottom lines). Lines are averaged over the 50 simulation runs. After $t \gtrsim 10^3~ {\rm yr}$ the coagulation drives most of the mass into the largest particles.

Open with DEXTER

The evolution of the mass-distribution (Figs. 10, 11, 13) can be divided into three stages. Initially, since particles start out as grains with sizes of $0.1\ \mu{\rm m}$, Brownian motion dominates. The size-distribution is therefore rather narrow, because the Brownian collision kernel favours the lighter particles. Quickly ($\sim$ $10^2\ {\rm yr}$), however, turbulent velocities become dominant and relative velocities are now highest for the more massive (high $\tau_{\rm f}$) particles. Once the first compaction event occurs (dotted, blue line) $\tau_{\rm f}$ enters the regime in which it becomes (at least) proportional to size and the pace of coagulation strongly accelerates toward larger sizes. These findings correspond well with the simple analytical model of Blum (2004), where in his Fig. 5 the Brownian motion driven growth is also followed by a stage in which the growth is exponential. This evolution could be called "run-away'', were it not for the fact that (in our case) the particle distribution is truncated at $\tau_{\rm rain}$ (Eq. (26)). The distribution at the first rain-out event is shown by the dashed, red line. Thereafter, the mass function collapses and evolves to a monodisperse population, close to the rain-out size (Figs. 10, 13). Two causes conspire to make these particles favoured: first, large particles can only be (efficiently) removed by a collision with a similar-sized particle (and no longer by a larger particle since these have disappeared); second, in the $\alpha _{\rm T} = 10^{-4}$ models friction times are always in the $\tau_{\rm f} < t_{\rm s}$ regime and relative velocities between similar-sized particles are suppressed. For these reasons, particles in the $\alpha _{\rm T} = 10^{-4}$ models near rain-out deplete the smaller particles faster than they deplete themselves, and the size distribution evolves again to monodispersity. Note, however, that this behaviour is essentially caused by the imposed presence of a sharp cut-off size due to the rain-out criterion. In reality, a more smooth transition can be expected.

Although the qualitative trends between the porous and compact models are essentially the same - fractal growth, compaction and run-away growth, rain-out and depletion - it is unambiguously clear that the porous models evolve to larger particles, as is also seen in Table 2 in which the values for the rain-out particles are tabulated. The size difference at rain-out is $\sim$2 order of magnitude in size and $\sim$4 orders of magnitude in mass. Particles only rain-out at $\tau_{\rm rain}$ and in the porous models particles have to grow much further before the critical friction time is reached (Eq. (26)). The inclusion of porosity in the coagulation models thus allows particles (when $\alpha_{\rm T} > 10^{-4}$) to grow to cm/dm sizes in the gaseous nebula, i.e., before settling to the mid-plane.

Apart from determining the size at which particles rain-out, $\alpha _{\rm T}$ also determines the pace of coagulation. In the $\alpha _{\rm T}=10^{-2}$ models (Fig. 11A) coagulation is rapid. Also, a large degree of stochasticity is seen among different models. In the $\alpha _{\rm T}=10^{-6}$ models, on the other hand, the turbulent velocities are very small and the support against gravity is weak, such that that rain-out happens before any compaction takes place. These models are most reminiscent of the "laminar nebula'', where systematic (i.e., settling) velocities dominate and are therefore most prone to gravitational instability effects (Hubbard & Blackman 2005). The comparison between the various $\alpha _{\rm T}$-models is perhaps best seen in the "evolution tracks'' of Figs. 12. They show that initially all porous models follow the same (fractal) curve, until the moment compaction occurs. In the 1 AU, $\alpha _{\rm T}=10^{-2}$ model a significant compaction of aggregates can clearly be observed (descending line). After $t > 10^3\ {\rm yr}$, this is followed by a slight increase in $\langle a \rangle_m / \langle a_{\rm comp} \rangle_m$; apparently due to the heavy rain-out, most collisions are again in the fractal regime.

Comparing the coagulation of the two materials studied here, i.e., quartz for the 1 AU models and ice for the 5 AU models, one sees similarities in their collisional evolution (see also Fig. 13). It seems, however, that the coagulation at 5 AU is somewhat slower. This can be explained naturally because of the lower density, but also perhaps because compaction is more difficult to achieve due to the higher surface density ($\gamma$) of ices. Although the differences are subtle, one can see, e.g., in Fig. 12 that the curves of the $\alpha_{\rm T} = 10^{-4}, 10^{-2}$ porous models level-off at higher $\langle a \rangle_m$-to- $\langle a \rangle$ ratio than the corresponding 1 AU curves, indicating compaction is achieved "easier'' at 1 AU. Also, at $\alpha _{\rm T}=10^{-2}$ the 1 AU particles that rain-out are strongly compacted (an even higher $\alpha _{\rm T}$ would have led to fragmentation), while the rain-out particles at 5 AU do not compact considerably before rain-out (see also Table 2). Thus, the larger surface energy ($\gamma$) of ices translates into a higher rolling energy and, subsequently, to decreased compaction.

  
5 Discussion

Figure 14: The m-$\psi$ relation for several aggregation models. Plotted is $\psi (m)$ for: the most massive particle in one of the R1Ta4-P models (black curve); the PCA, CCA aggregation models, discussed in Sect. 2.4, (grey curves); and compact ($\psi =1$) models. The dashed lines indicate points of equal friction times, $\tau _{\rm rain} \approx 500\ {\rm s}$ and $t_{\rm s}\approx 1600\ {\rm s}$. The cross shows the point at which the model experiences the first compaction event.

Open with DEXTER

	Figure 15: Geometrical optical depth to the mid-plane as function of time for the 1 AU models ( left) and the 5 AU models ( right). The optical depth is computed at every logarithmic interval in time and also at the point where particles rain-out (red symbol). For illustrative purposes the points are connected by lines. Error bars denote the spread throughout the 50 runs of each model.
Open with DEXTER

It now becomes clear that the internal structure of particles, represented here by the enlargement parameter $\psi$, is a variable of key importance in models of dust-aggregation. To illustrate this point further, Fig. 14 compares the "evolution curve'' of our (default $\alpha _{\rm T} = 10^{-4}$, 1 AU) model with those of compact, PCA and CCA aggregation (grey curves). The superimposed black curve connects the $(m,\psi)$ values of the most massive particle resulting from the collision model. The small, erratic structure in this curve corresponds to the fact that the particles fluctuate stochastically during the simulation. Furthermore, since we compare single particles here, instead of a (weighted) mean of the distribution as, e.g., in Fig. 12, a fixed point in the figure corresponds to one particular friction time and lines of equal friction times lie parallel to the dashed lines indicating $\tau=\tau_{\rm rain}$ and $\tau=t_{\rm s}$. For the specific choice of $\alpha _{\rm T} = 10^{-4}$ we see that $\tau_{\rm rain}$ comes prior to $t_{\rm s}$ and velocities therefore remain modest. Still, in our collision model the coagulation will reach a point, indicated by a cross in Fig. 14, at which some collisions, having the right mass ratio and relative velocity, are energetic enough to cause compaction ( $E > E_{\rm roll}$), which effectively halts any further increase in porosity. However, due to the earlier extensive build-up of porosity in the fractal regime, the particle distribution now evolves to larger sizes as compared to the compact models (Fig. 10), causing the rain-out masses to be orders of magnitudes higher (Table 2). Thus, a major part of the growth takes place in the nebula phase. Large, porous particles are quickly produced, stay in the nebula mixed with the gas and only settle when they are sufficiently compacted, e.g., by energetic collisions (this paper) or shocks (see below).

The curve of our model, with its characteristic bending point due to compaction, is a direct result of treating porosity as a dynamic variable that is altered by the collisional process. In all other models in Fig. 14, on the other hand, compaction is not incorporated, resulting in straight lines when the growth is in the fractal regime. In the PCA/CCA fractal models compaction is of course a priori ruled-out. However, the pre-assumed absence of compaction in the PCA/CCA models is consistent with the low collision energies in these limiting cases. In CCA, friction times barely grow ( $\tau_{\rm f} \propto m^{0.05}$), the $t_{\rm s}$ "threshold'' is not reached and relative velocities vanish for similar particles. In PCA, on the other hand, relative velocities are higher, but the collision energy is now suppressed by the reduced mass. Thus, if the coagulation process would (for some reason) be restricted to these limiting growth modes, aggregates will not restructure. Figure 14 shows how this affects the overall coagulation process. It shows, e.g., that compact grains rain-out at $\sim$ $10^{-4}~ {\rm g}$, "PCA grains'' at $\sim$ $10^{-2}\ {\rm g}$, porosity-evolving particles of our model at $\sim$ $1\ {\rm g}$, and "CCA particles'' will grow forever! It is clear that the porosity evolution of collisional agglomerates is of decisive influence on the coagulation process. Modelling the porosity evolution in combination with a microphysical collision model is therefore a key requirement for a full understanding of the first stages of planet formation.

To quantify the effects of coagulation on the appearance of the disk, we have calculated optical depths to the mid-plane in the MSN. As a first order approximation, the vertical structure is to be taken of constant density and extends over one scaleheight. We assume that this layer is represented by the particles of our simulation and that the rain-out particles (which we do not follow) are below it, i.e., at the mid-plane regions of the disk. The geometrical optical depth (i.e., at visible/UV-wavelengths) to the mid-plane is then calculated as

$$\tau_{\rm geom} = H_{\rm g} \int {\rm d}m\ f(m) \pi a^2.$$ (28)

Results are given in Fig. 15 for the 1 AU and 5 AU models. In Fig. 15A it is seen that the $\alpha _{\rm T}=10^{-6}$ models stay optically thick for most of the disk's evolution. Note also that the $\alpha _{\rm T} = 10^{-4}$ porous models (solid lines) and the $\alpha _{\rm T} = 10^{-4}$ compact models (dashed-dotted line) do not deviate much in $\tau_{\rm geom}$. This shows again the dual effects of porosity on the population: it increases the geometrical cross-section, yet it also speeds up the coagulation, causing more mass to be "locked'' inside large particles. At 5 AU (Fig. 15B), the timescales are longer and the disks only becomes optically thin after $\sim$ $10^6\ {\rm yr}$ when $\alpha_{\rm T} \lesssim 10^{-4}$. In the $\alpha _{\rm T}=10^{-2}$ models, evolution to optical thinness is very fast at both radii. It is clear that, within the frame-work of these models, the inner regions of protoplanetary disks are rapidly depleted of small grains unless $\alpha_{\rm T} \sim 10^{-6}$ or less.

There is a further serious issue hidden here. All our models show a rapid decline of the (sub-)micron size population on a timescale of $\sim$ $10^3\ {\rm yr}$. This is a well-known problem in models for grain coagulation in protoplanetary environments: the densities are high enough for coagulation to proceed rapidly (Dullemond & Dominik 2005). Furthermore, the turbulence induced relative velocities promote collisions between particles with different friction times, i.e., between small and larger particles. In contrast, observations reveal the presence of copious amounts of small grains in the disk-photospheres of isolated Herbig AeBe stars and T-Tauri stars, pointing toward the presence of small grains on timescales of some $10^6~ {\rm yr}$ (Natta et al. 2006; Meeus et al. 2003; van Boekel et al. 2005,2004). This discrepancy implies a continuous replenishment of (sub)micron-sized grains. In particular, it may reflect the importance of vaporisation and condensation processes continuously forming fresh, small grains. Likely, this vaporisation and condensation would be localised in the hot and dense, inner regions of the disk and these grains would then have to be transported upwards and outwards to the disk photosphere through diffusion processes. The high degree of crystallinity of silicates in the inner few AU of protoplanetary disks also points toward the importance of condensation processes in these environments (van Boekel et al. 2004). Likewise, the presence of crystalline silicates in the cold outer regions of protoplanetary disks has been attributed to large scale mixing of materials in these environments (Gail 2004; Bockelée-Morvan et al. 2002). Alternatively, the replenishment of small grains is through collisional fragmentation. These energetic collisions could take place either in the gaseous nebula due to high relative velocities driven by a high $\alpha _{\rm T}$, or in the mid-plane regions with the subsequent upwards diffusion of small grains. While it might be difficult to sustain $\alpha_{\rm T} \gtrsim 10^{-2}$ over a prolonged period of time, fragmentation in the mid-plane regions seems viable since more massive particles will reside here. Furthermore, if the mid-plane becomes dust-dominated ( $\rho_{\rm d} > \rho_{\rm g}$), shear-turbulence will develop, further augmenting the collisional energies (Cuzzi et al. 1993).

Further constraints on the collisional growth of grains in protoplanetary disks are provided by the solar system record; specifically, the chondrules and Ca-Al-rich Inclusions (CAI), which dominate the composition of primitive meteorites. These millimetre-sized igneous spherules are high-temperature components that formed during transient heating events in the early solar system. We realise that the cm-sized fluff-balls formed in our porous coagulation models are in the right mass-range of these meteorite components. It is tempting to identify these fluff-balls as the precursors to the chondrules and CAIs. We might then envision a model where the collisional evolution is terminated by the flash-heating event, for example a shock or lightning, which leads to melting and the formation of a chondrule and subsequent immediate settling. During the settling phase the chondrule may acquire a dust rim by sweeping up small dust grains or other unprocessed fluff-balls still suspended in the nebula (Cuzzi 2004). One key point to recognise here is that chondrules show a spread in age of a few million years (Wood 2005), which indicates that the collisional grain growth process takes place over a much longer timescale than our models would predict (see above).

In this study we have focused on the agglomeration driven by random motions - either Brownian or turbulent - high up in the nebula. Growth is then presumed to be terminated once the aggregate has compacted enough to settle. At that point, the aggregate will drop-down about one scaleheight after which further growth must occur for further settling to continue. In reality, instead of the simple two-component picture of nebula and mid-plane presented in this work, the nebula will acquire a stratified appearance (Dubrulle et al. 1995), where larger particles with higher friction times have smaller scaleheights. Collisional evolution models should be able to include this stratified nature of the disk. This stratification also extends in the radial direction due to radial drift motions and turbulent diffusion. Incorporation of these motions into the Monte Carlo code will be challenging. We expect that the study presented here can serve as the basis for incorporating realistic grain growth in hydrodynamical models, likely in the form of well-selected "collision-recipes''. An obvious step would be to include collisions that exceed the $E_{\rm frag}$ limit in the collision model. Note, for example, that in the Dominik & Tielens (1997) terminology what we have called "fragmentation'' should in fact be sub-divided in a continuous range of aggregate disruptions. At first monomers will be lost and only if $E \gg E_{\rm frag}$ are aggregates completely shattered. Other extensions to the model are to allow for a distribution of monomer sizes and to use monomers of different chemical composition. Both will affect the critical energy for restructuring, $E_{\rm roll}$. However, with an increasing number of parameters characterising the collision, it is worthwhile to verify experimentally - either through laboratory experiments or through detailed numerical calculations - which are of prime importance.

  
6 Conclusions

We have presented a model that incorporates the internal structure of aggregates as a variable in coagulation models. We used the enlargement parameter $\psi$ to represent the internal structure. It is seen that the internal structure is key to the collisional evolution, since it crucially affects the dust-gas coupling. However, in the model presented in Sect. 2 $\psi$ is not a static variable. It is altered by the collisional process and we have constructed simple recipes to include this aspect in coagulation models.

Next, we have applied the new collision model to the collisional evolution of the turbulent protoplanetary disk, until particles rain-out to the mid-plane. Our main conclusions are:

• With the treatment of porosity as a variable, three different regimes can be distinguished: fractal growth, compaction and fragmentation (Blum 2004). These regimes are also reflected in the collisional evolution of the size distribution: fractal growth (mostly Brownian motion), compaction (growth accelerates) and rain-out.
• The collisional evolution of our porosity-evolving model is quantitatively different from PCA/CCA aggregation models in which porosity can be parameterised by a fixed exponent. Therefore, a microphysical collision model is a key requirement for coagulation models.
• Due to the porous evolution, particles up to dm-sizes can be suspended in the gaseous nebula, orders of magnitude larger in mass than in models of compact coagulation. Therefore, chondrule precursors could have had their origin in regions above the mid-plane.
• If $\alpha_{\rm T} < 10^{-2}$, no fragmentation occurs in the gaseous nebula. Therefore, if $10^{-6} < \alpha_{\rm T} < 10^{-2}$, the inner nebula will become optically thin on timescales of $\sim$ $10^4\ {\rm yr}$, unless an influx of small grains takes place.

Appendix A: List of symbols

Table A.1: List of symbols.

Acknowledgements

CWO likes to show his appreciation to Carsten Dominik and Dominik Paszun for extensive discussion on the manuscript. These greatly helped to improve the quality of the paper, especially regarding Sect. 2. Carsten Dominik is also thanked for suggesting Figs. 2 and 13 and Dominik Paszun for useful assistance with creating Fig. 2. Finally, the referee, Jürgen Blum, is thanked for providing a thorough and helpful review.

References

• Balbus, S. A., & Hawley, J. F. 1991, ApJ, 376, 214 [NASA ADS] [CrossRef]
• Beckwith, S. V. W., & Sargent, A. I. 1991, ApJ, 381, 250 [NASA ADS] [CrossRef] (In the text)
• Beckwith, S. V. W., Henning, T., & Nakagawa, Y. 2000, Protostars and Planets IV, 533
• Blum, J. 2004, in ASP Conf. Ser. 309: Astrophysics of Dust, ed. A. N. Witt, G. C. Clayton, & B. T. Draine, 369
• Blum, J., & Wurm, G. 2000, Icarus, 143, 138 [NASA ADS] [CrossRef] (In the text)
• Blum, J., Wurm, G., Poppe, T., Kempf, S., & Kozasa, T. 2002, Adv. Space Res., 29, 497 [NASA ADS]
• Bockelée-Morvan, D., Gautier, D., Hersant, F., Huré, J.-M., & Robert, F. 2002, A&A, 384, 1107 [EDP Sciences] [NASA ADS] [CrossRef] [EDP Sciences]
• Boogert, A. C. A., Pontoppidan, K. M., Lahuis, F., et al. 2004, ApJS, 154, 359 [NASA ADS] [CrossRef]
• Bouwman, J., de Koter, A., van den Ancker, M. E., & Waters, L. B. F. M. 2000, A&A, 360, 213 [NASA ADS]
• Brownlee, D. E. 1979, Meteoritics, 14, 358 [NASA ADS] (In the text)
• Chokshi, A., Tielens, A. G. G. M., & Hollenbach, D. 1993, ApJ, 407, 806 [NASA ADS] [CrossRef]
• Cuzzi, J. N. 2004, Icarus, 168, 484 [NASA ADS] [CrossRef]
• Cuzzi, J. N., & Hogan, R. C. 2003, Icarus, 164, 127 [NASA ADS] [CrossRef] (In the text)
• Cuzzi, J. N., Dobrovolskis, A. R., & Champney, J. M. 1993, Icarus, 106, 102 [NASA ADS] [CrossRef] (In the text)
• Cuzzi, J. N., Ciesla, F. J., Petaev, M. I., et al. 2005, in Chondrites and the Protoplanetary Disk, ed. A. N. Krot, E. R. D. Scott, & B. Reipurth, ASP Conf. Ser., 341, 732 (In the text)
• Dominik, C., & Tielens, A. G. G. M. 1997, ApJ, 480, 647 [NASA ADS] [CrossRef]
• Dubrulle, B., & Valdettaro, L. 1992, A&A, 263, 387 [NASA ADS] (In the text)
• Dubrulle, B., Morfill, G., & Sterzik, M. 1995, Icarus, 114, 237 [NASA ADS] [CrossRef] (In the text)
• Dullemond, C. P., & Dominik, C. 2005, ApJ, 434, 971 [NASA ADS]
• Gail, H.-P. 2004, A&A, 413, 571 [EDP Sciences] [NASA ADS] [CrossRef] [EDP Sciences]
• Gibb, E. L., Whittet, D. C. B., Boogert, A. C. A., & Tielens, A. G. G. M. 2004, ApJs, 151, 35 [NASA ADS] [CrossRef]
• Gillespie, D. T. 1977, J. Atmos. Sci., 32, 1977 [NASA ADS] [CrossRef] (In the text)
• Hawley, J. F., & Balbus, S. A. 1991, ApJ, 376, 223 [NASA ADS] [CrossRef]
• Hayashi, C. 1981, Progr. Theor. Phys. Suppl., 70, 35 [NASA ADS] (In the text)
• Heim, L.-O., Blum, J., Preuss, M., & Butt, H.-J. 1999, Phys. Rev. Lett., 83, 3328 [NASA ADS] [CrossRef] (In the text)
• Hubbard, A., & Blackman, E. G. 2005, ArXiv Astrophysics e-prints (In the text)
• Kempf, S., Pfalzner, S., & Henning, T. K. 1999, Icarus, 141, 388 [NASA ADS] [CrossRef] (In the text)
• Kessler-Silacci, J., Augereau, J.-C., Dullemond, C. P., et al. 2006, ApJ, 639, 275 [NASA ADS] [CrossRef]
• Kostoglou, M., & Konstandopoulos, A. G. 2001, J. Aerosol Sci., 32, 1399 (In the text)
• Krause, M., & Blum, J. 2004, Phys. Rev. Lett., 93, 021103 [NASA ADS] [CrossRef] (In the text)
• Kress, M. E., & Tielens, A. G. G. M. 2001, Meteoritics and Planetary Science, 36, 75 [NASA ADS] (In the text)
• Liffman, K. 1991, J. Comput. Phys., 100, 116 [NASA ADS] [CrossRef] (In the text)
• Markiewicz, W. J., Mizuno, H., & Voelk, H. J. 1991, A&A, 242, 286 [NASA ADS]
• Meakin, P. 1988, Ann. Rev. Phys. Chem., 39, 237 [CrossRef]
• Meakin, P., & Donn, B. 1988, ApJ, 329, L39 [NASA ADS] [CrossRef]
• Meeus, G., Sterzik, M., Bouwman, J., & Natta, A. 2003, A&A, 409, L25 [EDP Sciences] [NASA ADS] [CrossRef] [EDP Sciences]
• Metzler, K., Bischoff, A., & Stoeffler, D. 1992, Geochim. Cosmochim. Acta, 56, 2873 [NASA ADS] [CrossRef]
• Mizuno, H., Markiewicz, W. J., & Voelk, H. J. 1988, A&A, 195, 183 [NASA ADS]
• Nakagawa, Y., Sekiya, M., & Hayashi, C. 1986, Icarus, 67, 375 [NASA ADS] [CrossRef] (In the text)
• Natta, A., Testi, L., Calvet, N., et al. 2006, ArXiv Astrophysics e-prints
• Nomura, H., & Nakagawa, Y. 2006, ApJ, 640, 1099 [NASA ADS] [CrossRef]
• Ossenkopf, V. 1993, A&A, 280, 617 [NASA ADS]
• Paszun, D., & Dominik, C. 2006, Icarus, 182, 274 [NASA ADS] [CrossRef] (In the text)
• Przygodda, F., van Boekel, R., Àbrahàm, P., et al. 2003, A&A, 412, L43 [EDP Sciences] [NASA ADS] [CrossRef] [EDP Sciences]
• Safronov, V. S. 1969, Evolution of the protoplanetary cloud and the formation of the Earth and planets; NASA TTF-667
• Schräpler, R., & Henning, T. 2004, ApJ, 614, 960 [NASA ADS] [CrossRef]
• Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337 [NASA ADS] (In the text)
• Shu, F. H., Adams, F. C., & Lizano, S. 1987, ARA&A, 25, 23 [NASA ADS] (In the text)
• Silk, J., & Takahashi, T. 1979, ApJ, 229, 242 [NASA ADS] [CrossRef]
• Smith, M., & Matsoukas, T. 1998, Chem. Eng. Sci., 53, 1777 (In the text)
• Smoluchowski, M. V. 1916, Zeitschrift fur Physik, 17, 557 [NASA ADS] (In the text)
• Suttner, G., & Yorke, H. W. 2001, ApJ, 551, 461 [NASA ADS]
• Tanaka, H., & Nakazawa, K. 1994, Icarus, 107, 404 [NASA ADS] [CrossRef]
• Tanaka, H., Himeno, Y., & Ida, S. 2005, ApJ, 625, 414 [NASA ADS] [CrossRef]
• van Boekel, R., Waters, L. B. F. M., Dominik, C., et al. 2003, A&A, 400, L21 [EDP Sciences] [NASA ADS] [CrossRef] [EDP Sciences]
• van Boekel, R., Min, M., Leinert, C., et al. 2004, Nature, 432, 479 [NASA ADS] [CrossRef]
• van Boekel, R., Min, M., Waters, L. B. F. M., et al. 2005, A&A, 437, 189 [EDP Sciences] [NASA ADS] [CrossRef] [EDP Sciences]
• Voelk, H. J., Morfill, G. E., Roeser, S., & Jones, F. C. 1980, A&A, 85, 316 [NASA ADS] (In the text)
• Weidenschilling, S. J. 1977, MNRAS, 180, 57 [NASA ADS]
• Weidenschilling, S. J. 1984a, Icarus, 60, 553 [NASA ADS] [CrossRef]
• Weidenschilling, S. J. 1984b, in Lunar and Planetary Institute Conference Abstracts, 900
• Weidenschilling, S. J. 1997, Icarus, 127, 290 [NASA ADS] [CrossRef]
• Weidenschilling, S. J., & Cuzzi, J. N. 1993, in Protostars and Planets III, ed. E. H. Levy, & J. I. Lunine, 1031 (In the text)
• Whipple, F. L. 1972, in From Plasma to Planet, ed. A. Elvius, 211
• Wood, J. A. 2005, in Chondrites and the Protoplanetary Disk, ed. A. N. Krot, E. R. D. Scott, & B. Reipurth, ASP Conf. Ser., 341, 953 (In the text)
• Youdin, A. N. 2005 [arXiv:astro-ph/0508659]
• Youdin, A. N., & Chiang, E. I. 2004, ApJ, 601, 1109 [NASA ADS] [CrossRef] (In the text)