A&A 461, 215-232 (2007)
C. W. Ormel1 - M. Spaans1 - A. G. G. M. Tielens1,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
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 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 and on the position in the nebula, aggregates can grow to (porous) sizes of 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
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 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 CH4 (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 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.
|Figure 1: Sketch of the turbulence induced relative velocities, , as function of friction time, , for equal friction times, , (solid line) and different friction times (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 , and .|
|Open with DEXTER|
|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., ) and the enlargement factor 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 between two particles of equal
in general, however,
and results are
derived numerically (Voelk et al. 1980; Markiewicz et al. 1991).
Weidenschilling (1984a,b) has presented analytical fits to these numerical
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,
), less than t0, these read,
then corresponds to the enlargement factor of a compact particle with specific density . 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 we refer to as the enlargement parameter or enlargement factor. is related to the structural parameter, x, of Ossenkopf (1993), as . Note, the close relationship between and porosity; a higher means more open aggregates, i.e., higher porosity. Therefore we will also use this more familiar reference for , but only in a qualitative sense.
Since our enlargement parameter, , is defined in order to match the geometrical cross-section, A, of a particle, we refer to the resulting radii, , as geometrical. The other cross-section of importance in coagulation models is , the collisional cross-section. Ossenkopf (1993) has consequentially defined a "toothing radius'', , such that and provides expressions for for PCA and CCA aggregates. 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 and have therefore simply used a for the calculation of the collisional cross-section. Likewise, we have also ignored augmentations of 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
More subtle are the effects of
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
increases linearly with size (
). 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.,
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 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 , 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 , 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 and .
|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 ( ) and for the line we have assumed equal masses (m1=m2) and . Arrows indicate how the critical energy lines shift with increasing monomer size and mass-ratio, m1/m2.|
|Open with DEXTER|
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,
|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 is the enlargement factor of the most massive aggregate, i.e., m1 > m2, and plots are shown for (grey lines) and (black lines). At low mass ratios the curves converge.|
|Open with DEXTER|
|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 correction term (Eq. (19)).|
|Open with DEXTER|
In the discussion on the
recipes in Sect. 2.3.1 we have used the CCA fractal exponent (
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.,
with m0 the monomer mass, the change in
after addition of a monomer becomes
Although, with Eq. (19) for , Eq. (15) does achieve the right PCA/CCA fractal limits, we do not claim they actually provide a model for . 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 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.
|Figure 6: Illustration of the adopted Monte Carlo technique. (A) N-particles are assumed to be uniformly distributed in a box of volume . The collision rate between particle i and j is . 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 (heavily exaggerated in this figure) to a volume at which contains N-particles again. The hypothesis of this method is that the collisional evolution within is representative for the coagulation of the total space under consideration.|
|Open with DEXTER|
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 -factor in Eq. (20) gives a further complication since there is no such thing as "conservation of porosity'' and the -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 (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, in which Kij is the collision kernel. A random number determines which of the possible collisions will be the next. The collision then creates a new particle, after which the 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 ; instead, - 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 m0 in a fixed volume, these will over time pile up in one agglomerate of . 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 goes proportional to N2 (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, (see Fig. 6B), under the assumption that the collisional evolution outside of 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- 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 Kij 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.
|Figure 7: Flow chart of the MC-coagulation method. One cycle corresponds to one collision. "Func'' like in 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 (gi -1) copies of particle i, it would be a waste of time to calculate the (same) rates gi times. Rather, gi can be absorbed in the calculation of the (combined) coagulation rate. is then the rate at which one of the i-particles collides with one of the j-particles () and between duplicates. The CPU time per step is now proportional to , the total number of distinct particles, i.e., excluding duplicates, with . To think of it in biological terms: gi gives the multiplicity of species i; the total number of species; and N the total number of living creatures.
|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 9: Tests of the Monte Carlo coagulation code. Linear kernel . - (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 . To obtain better correspondence with theory, the number of particles, N, is increased as the simulation progresses, such that more volume is sampled. 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., . 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 has been kept constant, which ensures a constant growth (in exponential terms) per cycle. We have fixed 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 N2 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, , as variables, can now be put into an evolutionary setting.
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 .
The 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 is to relate it to the observed accretion rate, , (Cuzzi et al. 2005) and then values in the range of seem plausible. Yet, despite its uncertainty, appears in key expressions as, e.g., and . Therefore, models are run that cover a large range in : . 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 and surface energy density of . The `outer' models correspond to conditions at 5 AU, where the coagulation is that of ices ( and ) 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 runs. Compact models (denoted by the C-suffix in Table 1) are models where the internal structure does not evolve, i.e., by definition. An overview of all the models is given in Table 1.
Table 1: Overview of all the runs.
|Figure 10: The mass function plotted at various times for the models. The panels compare the coagulation of the compact models (, 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 until . In the first 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 ( ), 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 . 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 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 models the spread in the runs is very large, causing the error bars to overlap. In the 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, , over the mass-weighted compact size, , ( in the text), against . 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
is defined analogously to Eq. (27), but then for the compact particle size.
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
This ratio is plotted vs.
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, (solid line), and the mean size, (dashed line), are calculated for the default models ( ) at 1 AU (top lines) and 5 AU (bottom lines). Lines are averaged over the 50 simulation runs. After 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 , Brownian motion dominates. The size-distribution is therefore rather narrow, because the Brownian collision kernel favours the lighter particles. Quickly ( ), however, turbulent velocities become dominant and relative velocities are now highest for the more massive (high ) particles. Once the first compaction event occurs (dotted, blue line) 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 (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 models friction times are always in the regime and relative velocities between similar-sized particles are suppressed. For these reasons, particles in the 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 2 order of magnitude in size and 4 orders of magnitude in mass. Particles only rain-out at 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 ) 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, also determines the pace of coagulation. In the models (Fig. 11A) coagulation is rapid. Also, a large degree of stochasticity is seen among different models. In the 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 -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, model a significant compaction of aggregates can clearly be observed (descending line). After , this is followed by a slight increase in ; 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 () of ices. Although the differences are subtle, one can see, e.g., in Fig. 12 that the curves of the porous models level-off at higher -to- ratio than the corresponding 1 AU curves, indicating compaction is achieved "easier'' at 1 AU. Also, at the 1 AU particles that rain-out are strongly compacted (an even higher 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 () of ices translates into a higher rolling energy and, subsequently, to decreased compaction.
|Figure 14: The m- relation for several aggregation models. Plotted is 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 () models. The dashed lines indicate points of equal friction times, and . 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|
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 ( ), the "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 , "PCA grains'' at , porosity-evolving particles of our model at , 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
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 . 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 (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 , or in the mid-plane regions with the subsequent upwards diffusion of small grains. While it might be difficult to sustain 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 ( ), 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 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 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, . 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.
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:
Table A.1: List of symbols.
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.