Dust coagulation and fragmentation in molecular clouds. I. How collisions between dust aggregates alter the dust size distribution
aa r X i v : . [ a s t r o - ph . S R ] J un Astronomy&Astrophysicsmanuscript no. 11158Orm c (cid:13)
ESO 2018October 24, 2018
Dust coagulation and fragmentation in molecular clouds
I. How collisions between dust aggregates alter the dust size distribution
C.W. Ormel , , D. Paszun , C. Dominik , , and A.G.G.M. Tielens , Kapteyn Astronomical Institute, University of Groningen, PO box 800, 9700 AV Groningen, The Netherlands Max-Planck-Institut f¨ur Astronomie, K¨onigstuhl 17, 69117, Heidelberg, Germany; e-mail: [email protected] Sterrenkundig Instituut Anton Pannekoek, Kruislaan 403, 1098 SJ Amsterdam, The Netherlands; e-mail:
[email protected] Afdeling Sterrenkunde, Radboud Universiteit Nijmegen, Postbus 9010, 6500 GL Nijmegen, The Netherlands Ames Research Center, NASA, Mail Stop 245-3, Mo ff ett Field, CA 94035, USA Leiden Observatory, Leiden University, P.O. Box 9513, 2300 RA Leiden, The Netherlandse-mail: [email protected]
Preprint online version: October 24, 2018
ABSTRACT
The cores in molecular clouds are the densest and coldest regions of the interstellar medium (ISM). In these regions ISM-dust grainshave the potential to coagulate. This study investigates the collisional evolution of the dust population by combining two models: abinary model that simulates the collision between two aggregates and a coagulation model that computes the dust size distributionwith time. In the first, results from a parameter study quantify the outcome of the collision – sticking, fragmentation (shattering,breakage, and erosion) – and the e ff ects on the internal structure of the particles in tabular format. These tables are then used as inputfor the dust evolution model, which is applied to an homogeneous and static cloud of temperature 10 K and gas densities between 10 and 10 cm − . The coagulation is followed locally on timescales of ∼ yr. We find that the growth can be divided into two stages: agrowth dominated phase and a fragmentation dominated phase. Initially, the mass distribution is relatively narrow and shifts to largersizes with time. At a certain point, dependent on the material properties of the grains as well as on the gas density, collision velocitieswill become su ffi ciently energetic to fragment particles, halting the growth and replenishing particles of lower mass. Eventually, asteady state is reached, where the mass distribution is characterized by a mass spectrum of approximately equal amount of mass perlogarithmic size bin. The amount of growth that is achieved depends on the cloud’s lifetime. If clouds exist on free-fall timescales thee ff ects of coagulation on the dust size distribution are very minor. On the other hand, if clouds have long-term support mechanisms,the impact of coagulation is important, resulting in a significant decrease of the opacity on timescales longer than the initial collisiontimescale between big grains. Key words.
ISM: dust, extinction – ISM: clouds – Turbulence – Methods: numerical
1. Introduction
Dust plays a key role in molecular clouds. Extinction of pen-etrating FUV photons by small dust grains allows moleculesto survive. At the same time, gas will accrete on dustgrains forming ice mantles consisting of simple molecules(Tielens & Hagen 1982; Hasegawa et al. 1992). Moreover, sur-face chemistry provides a driving force towards molecularcomplexity (Charnley et al. 1992; Aikawa et al. 2008). Finally,dust is often used as a proxy for the total gas (H ) col-umn density, either through near-IR extinction measurementsor through sub-millimeter emission studies (Johnstone & Bally2006; Alves et al. 2007; Jørgensen et al. 2008). Dust is often pre-ferred as a tracer in these types of studies because it is nowwell established that – except for pure hydrides – all speciescondense out in the form of ice mantles at the high densitiesof prestellar cores (Flower et al. 2006; Bergin & Tafalla 2007;Akyilmaz et al. 2007). Thus, it is clear that our assessment ofthe molecular contents of clouds, as well as the overall state ofthe star and planet formation process, are tied to the propertiesof the dust grains – in particular, its size distribution.The properties of interstellar dust are, however, expected tochange during its sojourn in the molecular cloud phase. First,condensation from the gas phase causes grain sizes to increase, forming ice mantles. This growth is limited, however, becausethere are many small grains – which dominate the total grainsurface area – over which the ice should be distributed; if allthe condensible gas freezes out, the thickness of the ice mantlesis still only 175 Å (Draine 1985). Therefore, in dense clouds,coagulation is potentially a much more important promoter ofdust growth. On a long timescale ( > yr), the interstellar grainsize distribution is thought to reflect a balance between coagula-tion in dense clouds and shattering in interstellar shocks as ma-terial constantly cycles between dense and di ff use ISM phases(Jones et al. 1996; Dominik & Tielens 1997).Infrared and visual studies of the wavelength dependenceof linear polarization and the ratio of total-to-selective extinc-tion were among the first observational indications of the im-portance of grain growth in molecular clouds (Carrasco et al.1973; Wilking et al. 1980; Whittet 2005). Early support for graingrowth by coagulation in molecular clouds was also providedby a Copernicus study that revealed a decreased amount of vi-sual extinction per H-nucleus in the ρ -Oph cloud relative to thevalue in the di ff use interstellar medium (Jura 1980). These typeof visual and UV studies are by necessity limited to the outskirtsof molecular clouds. Subsequent IR missions provided uniquehandles on the properties of dust deep inside dense clouds. In C.W. Ormel et al.: Dust Coagulation and Fragmentation in Molecular Clouds I. particular, comparison of far-IR emission maps taken by IRASand Spitzer and near-IR extinction maps derived from 2MASSindicate grain growth in the higher density regions (Schnee et al.2008). Likewise, evidence for grain coagulation is also providedby a comparison of visual absorption studies (e.g., star counts)and sub-millimeter emission studies which imply that the small-est grains have been removed e ffi ciently from the interstellargrain size distribution (Stepnik et al. 2003). Similarly, a recentcomparison of Spitzer-based, mid-IR extinction and submillime-ter emission studies of the dust characteristics in cloud cores re-veals systematic variations in the characteristics as a function ofdensity consistent with models of grain growth by coagulation(Butler & Tan 2009). Dust-to-gas ratios derived from a compar-ison of line and continuum sub-millimeter data is also consistentwith grain growth in dense cloud cores (Goldsmith et al. 1997).In recent years, X-ray absorption studies with Chandra have pro-vided a new handle on the total H column along a line of sight –that can potentially probe much deeper inside molecular cloudsthan UV studies – and in combination with Spitzer data, the de-creased dust extinction per H-nucleus reveals grain growth inmolecular clouds (Winston et al. 2007). Finally, Spitzer / IRS al-lows studies of the silicate extinction inside dense clouds anda comparison of near-IR color excess with 10 µ m optical depthreveals systematic variations which is likely caused by coagula-tion (Chiar et al. 2007). This is supported by an analysis of thedetailed absorption profile of the 10 µ m silicate absorption bandin these environments (Bowey et al. 1998; van Breemen et al.2009).Because it is the site of planet formation, theoretical coag-ulation studies have largely focused on grain growth in pro-toplanetary disks (Weidenschilling & Cuzzi 1993). In molecu-lar clouds, dust coagulation has been theoretically modeled byOssenkopf (1993) and Weidenschilling & Ruzmaikina (1994).In these studies, coagulation is driven by processes that pro-vide grains with a relative motion. For larger grains ( &
100 Å)turbulence in particularly is important in providing relative ve-locities. These motions – and the outcomes of the collisions –are very sensitive to the coupling of the particles to the turbu-lent eddies, which is determined by the surface area-to-mass ra-tio of the dust particles. At low velocities, grain collisions willlead to the growth of very open and flu ff y structures, while at in-termediate velocities compaction of aggregates occurs. At veryhigh velocities, cratering and catastrophic destruction will haltthe growth (Dominik & Tielens 1997; Paszun & Dominik 2009;Blum & Wurm 2008). Thus, to study grain growth requires us tounderstand the relationships between the macroscopic velocityfield of the molecular cloud, the internal structure of aggregates(which follows from its collision history), and the microphysicsof dust aggregates collisions. In view of the complexity of thecoagulation process and the then existing, limited understand-ing of the coagulation process, previous studies of coagulationin molecular cloud settings have been forced to make a num-ber of simplifying assumptions concerning the characteristics ofgrowing aggregates.Theoretically, our understanding of the coagulation processhas been much improved by the development of the atomic forcemicroscope, which has provided much insight in the binding ofindividual monomers. This has been translated into simple re-lationships between velocities and material parameters, whichprescribe under which conditions sticking, compaction, andfragmentation occur (Chokshi et al. 1993; Dominik & Tielens1997). Over the last decade, a number of elegant experimen-tal studies by Blum and coworkers (see Blum & Wurm 2008)have provided direct support for these concepts and in many ways expanded our understanding of the coagulation process.Numerical simulations have translated these concepts into sim-ple recipes, which link the collisional parameters and the ag-gregate properties to the structures of the evolving aggre-gates (Paszun & Dominik 2009). Together with the developmentof Monte Carlo methods, in which particles are individuallyfollowed (Ormel et al. 2007; Zsom & Dullemond 2008), thesestudies provide a much better framework for modeling the coag-ulation process than hitherto possible.In this paper, we reexamine the coagulation of dust grains inmolecular cloud cores in the light of this improved understand-ing of the basic physics of coagulation with a two-fold goal.First, we will investigate the interrelationship between the de-tailed prescriptions of the coagulation recipe and the structure,size, and mass of aggregates that result from the collisional evo-lution. Therefore, a main goal of this work is to explore the fullpotential of the collision recipes, e.g., by running simulationsthat last long enough for fragmentation phenomena to becomeimportant. Second, we aim to give a simple prescriptions forthe temporal evolution of the total grain surface area in molecu-lar clouds, thereby capturing its observational characteristics, interms of the physical conditions in the core.This paper is organized as follows. Section 2 presents a sim-plified and static model of molecular clouds we adopted in ourcalculations. Section 3 describes the model to treat collisions be-tween dust grains and, more generally, aggregates of dust grains.In Sect. 4 the results are presented: we discuss the imprints ofthe collision recipe on the coagulation and also present a param-eter study, varying the cloud gas densities and the dust materialproperties. In Sect. 5 we review the implications of our result tomolecular clouds and Sect. 6 summarizes the main conclusions.
2. Density and velocity structure of molecularclouds
The physical structure of molecular clouds – the gas density andtemperature profiles – is determined by its support mechanisms:thermal, rotation, magnetic fields, or turbulence. If there is onlythermal support to balance the cloud’s self-gravity and the tem-perature is constant, the density, assuming spherical symmetry,is given by ρ g ∝ r − , where ρ g is the gas density and r the ra-dius. However, the isothermal sphere is unstable as it heraldsthe collapse phase (Shu 1977). The cloud then collapses on afree-fall timescale t ff = s π G ρ g = . × yr (cid:18) n cm − (cid:19) − / , (1)where G is Newton’s gravitational constant, n = ρ g / m H µ thenumber density of the molecular gas, m H the hydrogen mass,and µ = .
34 the mean molecular mass. Thermally supportedcores are stable if the thermal pressure wins over gravity, a sit-uation described by the Bonnor-Ebert sphere, where an externalpressure confines the cloud (still assuming a constant tempera-ture).Magnetic fields in particular are important to support thecloud against the opposing influence of gravity. Ion-moleculecollisions provide friction between the ions and neutrals and in A list of symbols is provided in Appendix D. Note that our definition of density – n , the number of gas moleculesper cm – di ff ers from the density of hydrogen nuclei ( n H ), which ismore commonly referred to as density in the dust / ISM communities.For cosmic abundances n H is related to n as n ≃ . n H ..W. Ormel et al.: Dust Coagulation and Fragmentation in Molecular Clouds I. 3 that way couple the magnetic field to the neutral cloud. Overtime the magnetic field will slowly leak out on an ambipolar dif-fusion timescale t ad ∼ K in πµ m H G (cid:18) n i n (cid:19) ≃ . × yr (cid:18) n cm − (cid:19) − / , (2)where K in is the ion-molecular collision rate and n i the densityin ions. In Eq. (2) we have used an ion-neutral collision rate of K in = × − cm s − and a degree of ionization due to cosmicrays of n i / n = × − / √ n (Tielens 2005).Turbulence is another possible support mechanism of molec-ular cores, but its nature is dynamic – rather than (quasi)static.At large scales it provides global support to molecular clouds,whereas at small scales it locally compresses the gas. If theseoverdensities exist on timescales of Eq. (1), collapse will follow.This is the gravo-turbulent fragmentation picture of turbulence-dominated molecular clouds, where the (supersonic) turbulenceis driven at large scales, but also reaches the scales of quies-cent (subsonic) cores (Mac Low & Klessen 2004; Klessen et al.2005). In this dynamical, turbulent-driven picture both molecu-lar clouds and cores are transient objects.Thus, cloud cores will dynamically evolve due to either am-bipolar di ff usion and loss of supporting magnetic fields or due toturbulent dissipation, or simply because the core is only a tran-sient phase in a turbulent velocity field. In this work, for reasonsof simplicity, we consider only a static cloud model – the work-ing model – in which turbulence is unimportant for the supportof the core, but where (subsonic) turbulence is included in theformalism for the calculation of relative motions between thedust particles. In this exploratory study we will for simplicity adopt an homo-geneous core of mass given by the critical Jeans mass. Moreover,we assume the cloud is turbulent, but neglect the influenceof the turbulence on the support of the cloud. Thus, our ap-proximation is probably applicable for high density, low masscores as velocity dispersions increase towards high mass cores(Kawamura et al. 1998). The homogeneous structure ensuresthat collision timescales are the same throughout the cloud, i.e.,the coagulation and fragmentation can be treated locally. In ourcalculations, the sensitivity of the coagulation on the gas density n will be investigated and the relevant coagulation and fragmen-tation timescales will be compared to the timescales in Eqs. (1)and (2).Starting from the isodense sphere, the cloud outer radius isgiven by the Jeans length (Binney & Tremaine 1987) L J = s π c G ρ g = .
033 pc (cid:18) n cm − (cid:19) − / (cid:18) T
10 K (cid:19) / , (3)where c g = kT /µ is the isothermal sound speed and T the tem-perature of the cloud. A temperature of 10 K is adopted. Thesound crossing time L J / c g is comparable to the free-fall time ofthe cloud.Regarding the driving scales of the turbulence we assume(i) that the largest eddies decay on the sound crossing time, i.e., t L = L J / c g , and (ii) that the fluctuating velocity at the largestscale is given by the sound speed, v L = c g . Thus, the turbulentviscosity is ν t = L v L = v t L = c g L J with L = L J the size of thelargest eddies. Although our parametrization of the large eddy quantities seems rather ad-hoc, we can build some trust in thisrelation by considering the energy dissipation rate v / L per unitmass, which translates into a heating rate of n Γ = v L L ρ g = . × − erg cm − s − (cid:18) T
10 K (cid:19) (cid:18) n cm − (cid:19) / . (4)Based upon observational studies of turbulence in cores, Tielens(2005) gives a heating rate of n Γ = × − n erg s − , with whichEq. (4) reasonably agrees for the range of densities we consider.Additionally, the adoption of the crossing time and sound speedfor the large eddy properties are natural upper limits.The turbulent properties further follow from the Reynoldsnumber,Re = ν t ν m = v L Lc g ℓ mfp / = . × (cid:18) n cm − (cid:19) / (cid:18) T
10 K (cid:19) / , (5)where ν m is the molecular viscosity and ℓ mfp the mean free pathof a gas particle. Assuming a Kolmogorov cascade, the turn-overtime and velocity at the inner scale follow from the Reynoldsnumber t s = Re − / t L = . × yr (cid:18) n cm s − (cid:19) − / (cid:18) T
10 K (cid:19) − / ; (6a) v s = Re − / v L = . × cm s − (cid:18) n cm s − (cid:19) − / (cid:18) T
10 K (cid:19) / . (6b)A collisional evolution model requires a prescription for the relative velocities ∆ v between two solid particles. Apart fromturbulence, other mechanisms, reflecting di ff erences in the ther-mal, electrostatic, and aerodynamic properties of particles, willalso provide particles with a relative motion. However, undermost molecular cloud conditions turbulence will dominate thevelocity field (Ossenkopf 1993) and in this work we only con-sider turbulence. Then, the surface area-to-mass ratio of the dustparticles is of critical importance since this quantity determinesthe amount of coupling between the dust particles and the gas.We use the analytic approximations of Ormel & Cuzzi (2007)for the relative velocity between two particles. These expres-sions only include contributions that arise as a result of the parti-cle’s inertia in a turbulent velocity field and do not contain con-tributions that arise from gyroresonance acceleration (see, e.g.,Yan et al. 2004). See Appendix A for more details.
3. Collision model
Dust grains that collide can stick together, forming aggregates(see Fig. 1). In this work we consider the collisional evolution ofthe distribution of aggregates, f ( N , t ): the number of aggregatesof mass N with time. Many works have studied aggregate growthunder the conditions of perfect sticking upon contact, neglect-ing the e ff ects of the impact energy on the structure of aggre-gates (Meakin 1988; Meakin & Donn 1988; Ossenkopf 1993).Of special interest are the particle-cluster aggregation (PCA) andcluster-cluster aggregation (CCA) modes. In PCA, the aggre-gates collide only with single grains, while in CCA the collisionpartners are of similar size and structure. In CCA, the emergingstructures become true fractals, with a fractal dimension ∼
2. For
C.W. Ormel et al.: Dust Coagulation and Fragmentation in Molecular Clouds I. a out a s Fig. 1.
2D projection of a flu ff y aggregate with indication of thegeometrical radii, a σ , and the outer radii, a out .PCA, however, an homogeneous structure is eventually reachedat a filling factor of ∼ ff erent size. Furthermore, at higher energiesthe assumption of ‘hit-and-stick’ breaks down: aggregate bounc-ing, compaction (in which the constituent grains rearrange them-selves), and fragmentation lead to a rearrangement of the inter-nal structure. These collisional processes, except bouncing, areincluded in our collision model.We quantify the internal structure of aggregates in terms ofthe geometrical filling factor, φ σ , defined as φ σ = N (cid:18) a a σ (cid:19) , (7)where we have assumed that the aggregate contains N equal sizegrains of radius a with a σ = √ σ/π the projected surface equiv-alent radius of the aggregate. For very flu ff y aggregates a σ can bemuch less than the outer radius of the aggregate, a out , see Fig. 1.The definition of the filling factor in terms of the projected areadetermines its aerodynamic behavior, and thereby the relativevelocities ( ∆ v ) between the dust aggregates. Each collisions is classified into one of three groups:1. Hit-and-stick. At low collision energies, the internal struc-ture of the particles is preserved.2. Local. Only a small part of the aggregate is a ff ected by thecollision, as in, e.g., erosion. The mass ratio between the twoparticles is large.3. Global. The collision outcome results in a major change tothe structure or size of the target aggregate. Relevant forequal-size particles or at large energies.Figure 2 provides the adopted decision chain between the threeregimes. Parameters that enter the chain are the collision energy, E = m µ ( ∆ v ) = m m m + m ( ∆ v ) , (8) The compactness parameter φ σ is the inverse of the enlargementfactor ψ , previously used in Ormel et al. (2007). Ossenkopf (1993) uses x = φ − σ as its structural parameter. where m µ the reduced mass, the particle masses ( m and m , or,in number of grains, N and N ), and the critical energies E roll and E br . Here, E br is the energy to break a bond between twograins (the contact) and E roll the energy required to roll the con-tact area over a visible fraction of the grain. These critical ener-gies are defined as (Dominik & Tielens 1997) E br = A br γ / a / µ E ⋆ / ; (9a) E roll = π ξ crit γ a µ , (9b)where a µ is the reduced radius of the grains ( a µ = a / γ the surface energy density of the material,and E ⋆ the reduced elastic modulus. Here, following laboratoryexperiments (Blum & Wurm 2000) we adopt the values ξ crit = × − cm and A br = . × , which are larger than thetheoretically derived values of Dominik & Tielens (1997).Thus, when collisional energies are low enough for aggre-gate restructuring to be unimportant (experimentally determinedto be 5 E roll ; Blum & Wurm 2000) particles are in the hit-and-stick regime. Similarly, when the collision energy is su ffi cient tobreak all contacts the collision falls – obviously – in the globalregime. In the remainder the mass ratio of the colliding particlesdetermines whether the collision is global or local. For mass ra-tios smaller than 0.1 the changes become more localized and itis seen from the simulations that at this point the energy dis-tribution during a collision becomes inhomogeneous. Thus, wetake N / N = . N and φ σ ), the collision outcomeinvolves many other parameters (discussed in Sect. 3.3). Theseparameters are provided in tabulated form as a function ofthree input parameters – dimensionless energy parameter, fill-ing factor, and impact parameter b – for both the local and theglobal recipe. To obtain these collision parameters, direct nu-merical simulations between two colliding aggregates were per-formed, in which the equations of motions for all grains withinthe two colliding aggregates are computed (Paszun & Dominik2009). An example of these quantities is the fraction of miss-ing collisions, which is a result of the fact that we have de-fined the collision cross section σ C in terms of the outer radii, σ C = π ( a out , + a out , ) . Appendix B presents a description of thenumerical simulations with their results, discusses a few auxil-iary relations that are required for a consistent treatment of thecollision model, and treats the format of the collision tables. Two key limitation of these binary aggregate simulations fol-low from computational constraints: (i) the number of grains thatcan be included is limited to N ∼ ; and (ii) the simulationscannot take into account a grain size distribution that spans overorders of magnitude in mass. These limitations are reflected inour collision model and constitute a potential bottleneck for thelevel of realism for the application of our results to molecularclouds. We therefore first motivate our choice for the monomersize and present scaling relations that provide a way to extrapo-late the results beyond the parameter space sampled in the sim-ulation. The tables are provided online..W. Ormel et al.: Dust Coagulation and Fragmentation in Molecular Clouds I. 5 γ, E ⋆ N / N > . E > N + N ) E br E > E roll E br = A br γ / ( a / / ( E ⋆ ) / E roll = π ξ crit γ a E = m N µ ( ∆ u ) mass N ; filling factor f s,2 Particle mass N ; filling factor f s,1 Particle Velocity fieldDust properties material propertiesgrain size a YES
LOCALRECIPEGLOBALRECIPE OR NOYES NO
HIT & STICK
Critical energies breaking energyrolling energy
Collision properties collision energymass ratio N /N
Fig. 2.
Schematic decision chain employed to distinguish between the hit-and-stick, global, and local recipes.
Our recipe is based on simulations of aggregates that are built ofmonomers of a single size. Therefore, we treat a monodispersedistribution of grains. In reality, interstellar dust exhibits a sizedistribution, or a series of size distributions based on the vari-ous grain types (e.g., Desert et al. 1990; Weingartner & Draine2001). For simplicity, we compare our monodisperse approachwith the MRN grain distribution, in which the number of grainsdecreases as a − / n ( a )d a ∝ a − / , be-tween a lower ( a i ) and an upper ( a f ) size (Mathis et al. 1977).Thus, in the MRN-distribution the smallest grains dominate bynumber, whereas the larger grains dominate the mass. For anMRN distribution we adopt, a i =
50 Å and a f = . µ m. Toanswer the question which grain size best represents the MRNdistribution, we consider both its mechanical and aerodynamicproperties.In the monodisperse situation the mechanical properties ofa particle (its strength) can be estimated from the total bindingenergy per unit mass, E br / m , if we assume each grain has oneunique contact. In the literature the strength of a material is usu-ally denoted by the quantity Q . Thus, for a monodisperse modelwe have Q = E br ( a / m = k ( a / / a = k − / a − / , (10)where k = A br γ / / πρ s E ⋆ / is a material constant with ρ s thebulk density of the material. A smaller grain size thus lead tosignificantly stronger aggregates. For the MRN distribution weassume that a typical contact always involves a small grain, i.e., a µ ≃ a i enters in the E br expression. Assuming again that thenumber of contacts is of the order of the number of grains, theiraverage strength is given by Q MRN ≃ E br ( a i ) R a f a i n ( a )d a πρ s / R a f a i n ( a ) a d a ≃ ka − / a − / (11)where we have used that a f ≫ a i . Equating Eqs. (10) and (11)it follows that the grain size at which the monodisperse modelgives aggregates that have the same strength as the MRN is a ≃ / a / i a / f =
560 Å. Apart from the mechanical properties, the aerodynamic properties of aggregates are of crucial importance to the col-lisional evolution. This mainly concerns the initial (fractal)growth stage. For a single grain σ/ m = (3 / ρ s ) a − . For theMRN distribution an upper limit on σ/ m is provided by assum-ing that all of its surface is exposed, i.e., as in a 2D arrangementof grains; then, σ m = R a f a i π n ( a ) a d a R a f a i (4 πρ s / n ( a ) a d a = ρ s ( a i a f ) − / , (12)and the equivalent aerodynamic grain size becomes √ a i a f = σ becomeslower at the same mass; (ii) due to their low rolling energies, thesmallest grains of size a i will already initiate restructure at verylow velocities; (iii) in the case of ice-coating, the lower grainsize a i will be larger by a factor of ∼ σ/ m is nota constant. To nevertheless get a feeling of the trend, we havecalculated the aerodynamic properties of MRN aggregates thatconsists out of a few big grains, such that their total compactvolume is equivalent to a sphere of 0 . − . µ m. These MRNaggregates were constructed through a PCA process, i.e., addingone grain at a time. Because the aggregates contains the largegrains, they fully sample the MRN distribution and can thereforebe regarded as the smallest building blocks for the subsequentcollisional evolution.We observed that, due to the above mentioned self-shielding,the aerodynamic size increases to ∼ . − . µ m, significantlyhigher than the 2D limit of Eq. (12) (see also Ossenkopf 1993).Thus, the initial clustering phase of MRN-grains produces struc-tures that behave aerodynamically as compact grains of ∼ . µ m.We remark that this estimate is approximate – a CCA-like clus-tering will decrease it, whereas the above mentioned preferentialcompaction of the very small grains will increase σ/ m – but thetrend indicates that in 3D the aerodynamic size becomes skewedto the larger grains in the distribution. Therefore, we take 0 . µ mas the equivalent monomer grain size of the MRN-distribution,but to assess the sensitivity of the adopted grain size to the resultswe also consider models with a di ff erent grain sizes. C.W. Ormel et al.: Dust Coagulation and Fragmentation in Molecular Clouds I.
A key limitation of the aggregate-aggregate collision simulationsis the number of grains that can be used; typically, N . isrequired in order to complete a thorough parameter study withina reasonable timeframe. As a consequence, the mass ratio of thecolliding aggregate is also restricted. Furthermore, in the numer-ical experiments all simulations were performed using materialproperties applicable to silicates, whereas in molecular cloudswe expect the grains to be coated with ice mantles. Clearly, scal-ing of the results is required such that the findings of the numer-ical experiments can be applied to aggregates of di ff erent sizeand composition.Therefore, we scale the collisional energy E to the criticalenergies, E br and E roll , since these quantities involve the mate-rial properties. For a collision between silicate aggregates andice-coated aggregates a similar fragmentation behavior may beexpected if the collisional energy in the latter case is a factor E icebr / E silbr higher. Similarly, restructuring is determined by therolling energy, E roll . Thus, the collision energy is normalized to E roll where it concerns the change in filling factor and to E br forall other parameters that quantify the collision outcome.The division between the global and local recipes is alsoclosely linked to scaling arguments. In the global recipe ener-gies are normalized to the total number of monomers, N tot . Thus,a collision taking place at twice the energy and twice the massleads to the same fragmentation behavior. However, in the localrecipe the amount of damage will be independent of the size ofthe bigger particle. In this case we then scale by N µ , essentiallythe mass of the projectile. This information is captured in a sin-gle dimensionless energy parameter ε , ε = EN e ff E crit , (13)where E crit and N e ff depend on the context: the energy E crit canbe either one of E br or E roll , whereas N e ff is one of N tot or N µ (seeTable 3.3). In discussing the collision outcomes, we focus on the local andglobal recipes, which are a direct result of the numerical simu-lations. The hit-and-stick recipe is discussed in Appendix B.2.3.To streamline the recipe for a Monte Carlo approach, the speci-fication of the collision outcomes slightly deviates from our pre-vious study (Paszun & Dominik 2009).In the general case of a collision including fragmentation theemergent mass distribution, f ( m ), consists of two components:(i) a power-law component that describes the small fragments;and (ii) a large fragment component that consist out of one ortwo fragments (see Fig. 3). The separation between the two com-ponents is set, somewhat arbitrarily, at a quarter of the total massof the aggregates, N tot = N + N . (It turns out that for our simu-lations the precise place of the cut is unimportant, because of thelack of severe fragmentation events). Then, the power-law dis-tribution spans the range from monomer mass up to N = N tot / N tot /
4. To obtain the num-ber of large fragments, the recipes provide the mean number oflarge fragments, N f , together with its spread S f .Table 3.3 lists all quantities describing a collision outcome.Apart from N f and S f these include: – The fraction of missing collisions, f miss . This number givesthe fraction of collisions in which no interaction between the N f ( m ) − q power-lawcomponent large fragmentcomponent ← relative contribution of the → components determinedby f pwl N F , S F Fig. 3.
Sketch of the adopted formalism for the size distributionwith which the results of the aggregate collision simulation arequantified. See text and Table 3.3 for the description of the sym-bols. If f pwl = N f and S f parameters essentially indicate whether we have zero, one, ortwo large fragments.aggregates took place. Missing collision are a result from thechoice of normalizing the impact parameter b to the outerradius a out , ˜ b = b / ( a out , + a out , ) (see Appendix B.2.2). Forlarge values of ˜ b and very flu ff y structures f miss becomes sig-nificant. – The mass fraction in the power-law component, f pwl . It givesthe fraction of the total mass ( N tot ) that is contained in thepower-law component. In the local recipe f pwl is defined rel-ative to N µ , because here the amount of erosion is measuredwith respect to the smaller projectile. – The exponent of the power-law distribution, q . It determinesthe distribution of the small fragments, i.e., f ( m ) ∝ m − q . – The relative change in filling factor, C φ . It gives the change infilling factor of the large fragment component, φ new σ = C φ φ ini σ . C φ < C φ > C φ refers to the chance in the filling factorof the large aggregate (for both the local and global recipe),its dimensionless energy parameter ε is always normalizedto the total number of grains, N e ff = N tot . Thus, the com-paction may be local and moderate, but the a ff ected quantity– the filling factor – describes a global property. Moreover, toprevent possible spuriously high values of φ σ , we artificiallyassign an upper limit of 33% to the collisional compactionof aggregates (Blum & Schr¨apler 2004). Figure 4a shows how much mass is ejected during collisionsat di ff erent energies and for di ff erent filling factors (symbols).Recall that in the local recipe the f pwl quantity involves a nor-malization to N µ , rather than N tot . At high energies, then, theexcavated mass may exceed the mass of the small projectile byeven two orders of magnitude. The distribution of the small frag-ments created by the erosion is relatively flat with slopes oscil-lating between q = − . q = − .
3. The number of largefragments N F rarely increases above unity. The exception are the‘lucky projectiles’ that destroy the central contacts of very flu ff yaggregates, causing the two sides of the aggregate to become dis-connected. If energies are su ffi ciently high, fragments producedin a cratering event can result in secondary impacts, enhancingthe erosion e ffi ciency. .W. Ormel et al.: Dust Coagulation and Fragmentation in Molecular Clouds I. 7 Table 1.
Quantities provided by the adjusted collision recipe.Symbol Description Energy scaling parameters in Eq. (13)Global Local(1) (2) (3) (4) f miss Fraction of collisions that resulted in a miss a — — N f Mean number of large fragments N tot E br N µ E br S f Standard deviation of the N f N tot E br N µ E br f pwl The fraction of the mass in the small fragments component.Normalized to N tot (global recipe) or N µ (local recipe). N tot E br N µ E br q Exponent of the power-law distribution of small fragments N tot E br N µ E br C φ = φ σ /φ ini σ Relative change of the geometrical filling factor. N tot E roll N tot E roll N ote . Columns (3)–(4) denote the energy scaling expressions used to obtain the dimensionless energy parameter, ε , see Eq. (13). a Given asfunction of a out / a σ instead of φ σ , see Appendix B.3. Since the influence of the impact is local, the change in fillingfactor is relatively minor (see Fig. 4b). However, increasing thecollision energy results in an increasing degree of compression.Only very flu ff y and elongated aggregates may break in half,causing an artificial increase of the filling factor. This can beobserved in Fig. 4b for aggregates with φ ini σ = .
07 (diamonds),where the change in filling factor shows a strong variation forenergies above E = − N tot E roll . In Fig. 5 a few results from the global recipe are presented, inwhich results of collisions at central impact parameter (˜ b = left panels) and o ff -center collisions (˜ b = right panels) arecontrasted. In Figs. 5a and 5b the number of large particles thatremain after a collision, N f , is given. At low energies the numberof fragments is the same ( N f =
1) in both cases, reflecting stick-ing. At very high energies ( E > N tot E br ), the central collisionresults in catastrophic disruption ( N f = ff -center collisions,on the other hand, tend to produce two large fragments at higherenergies; because they interact only with their outer parts, theamount of interaction is insu ffi cient to let the colliding aggre-gates either stick or fragment.Figures 5c,d show the mass in the power-law component(small fragments). Central collisions result in an equal distri-bution of energy among the monomers. A collision energy of3 N tot E br is su ffi cient to shatter an aggregate. O ff -center colli-sions are more di ffi cult to fully destroy, though, and show, more-over, a strong e ff ect on porosity. In the most compact aggregate(crosses) over 70% of the mass ends up in the power-law compo-nent, whereas the remainder is in one large fragment. However,these are average quantities, and in some experiments all themass ended up in the power-law component as can be seen fromFig. 5b where N f drops below unity. For more flu ff y aggregatesthe fragmentation is much less pronounced, because the redistri-bution of the kinetic energy over the aggregate is less e ff ective.For example, very flu ff y aggregates of filling factor φ σ = . b = .
75 pro-duce small fragments which add up to only 6% of the total mass.The rest of the mass is locked into two large fragments.The degree of damage can also be assessed through the slopeof the power-law distribution of small fragments ( q , not plot-ted in Fig. 5). The steeper the slope, the stronger the damage.Heavy fragmentation produces many small fragments and re-sults in a steepening of the power-law. Although destruction isvery strong in the case of a central impact (the slope reaches val-ues of q = − . E > N tot E br ), it weakens considerably for o ff -center collisions ( q > − . q . However, for erosion thefragments are small in any case, independent of q .At low energies, the amount of aggregate restructuring,quantified in the C φ parameter, is independent of impact param-eter (Fig. 5e,f). This is simply because the collision energy istoo low for restructuring to be significant. The aggregates’ vol-ume then increases in a hit-and-stick fashion, resulting in a de-crease of the filling factor ( C φ < ff ect thefilling factor φ σ . Figure 5e shows that the compression is max-imal at an impact energy of about E = N tot E roll . Aggregatesthat are initially compact are di ffi cult to further compress, be-cause for filling factors above a critical value of 33% the requiredpressures increase steeply (Blum et al. 2006; Paszun & Dominik2008). Any further pressure will preferentially move monomerssideways, causing a flattening of the aggregate and a decreas-ing packing density. O ff -center collisions, however, lead to amuch weaker compression (Fig. 5f). Here, the forces acting onmonomers in the impacting aggregates are more tensile, and tendto produce two large fragments with the filling factor una ff ected, C φ =
4. Results
The formulation of the collision recipe in terms of the six outputquantities enables us to calculate the collisional evolution by aMonte Carlo method (see Appendix C for its implementation).The sensitivity of the collisional evolution to the environment(e.g., gas density, grain size, grain type; see Table 2) is assessed.The coagulation process is generally followed for 10 yr. Whilewe realize that bare silicates and the long timescales may notbe fully relevant for molecular clouds, we have elected here toextend our calculations to fully probe the characteristics of thecoagulation process. In particular, since fragmentation is explic-itly included in the collision model we evolve our runs until asteady-state situation materializes.In Sect. 4.1 the results from the standard model ( n = cm − , a = . µ m, ice-coated silicates) are analyzed.Section 4.2 presents the results of our parameter study. Figure 6 shows the progression of the collisional evolution ofice-coated silicates at a density of n = cm − (the stan-dard model) starting from a monodisperse distribution of 0 . µ m C.W. Ormel et al.: Dust Coagulation and Fragmentation in Molecular Clouds I.
Fig. 4.
Quantities provided by the local recipe. The left panel shows the mass in small fragments of the power-law component,normalized to the reduced mass of the colliding aggregates f pwl = M pwl / m µ . The right panel shows the relative change in thegeometrical filling factor C φ = φ new σ /φ ini σ . Symbols refer to the initial filling factor of the larger aggregate. Fig. 5.
Quantities provided by the global recipe. Left panels correspond to central collisions, while the right panels correspond too ff -center collisions at normalized impact parameter ˜ b = .
75. From top to bottom: number of large fragments N f ( A, B ); mass of thesmall fragments component, M pwl , normalized to the total mass of the two aggregates M tot ( C, D ); relative change in the geometricalfilling factor C φ = φ new σ /φ ini σ (E, F) . .W. Ormel et al.: Dust Coagulation and Fragmentation in Molecular Clouds I. 9 Table 2.
List of the model runs.id Density Type Grain size Figure ref. n [cm − ] a [ µ m](1) (2) (3) (4) (5)1 10 ice 0 .
12 10 silicates 0 . ice 0 . silicates 0 . silicates 16 a ice 0 . ice 1 Fig. 118 10 ice 0 .
03 Fig. 119 10 ice, compact b . ice, head-on c . silicates 0 . ice 0 . ice 0 . N ote . (1) Model number. (2) Number density of the gas. (3) Collisiontype, describing material parameters and collision setup. Here, ‘ice’refers to ice-coated silicates of bulk density identical to bare silicates, ρ s = .
65 g cm − , but di ff erent material properties: γ =
370 erg cm − and E ⋆ = . × dyn cm − . For bare silicates, γ =
25 erg cm − and E ⋆ = . × dyn cm − . (4) Monomer radius. (5) Figure reference. a The standard model; b filling factor of particles restricted to a minimumof 33%; c central impact collisions only ( b = mass, N − − − − − m a ss d e n s i t y , N f ( N ) yr yr yr n = , Ice × yr × yr × yr × yr Fig. 6.
Mass distribution of the standard model ( n = cm − , a = − cm, ice-coated silicates) at several times during itscollisional evolution, until t = × yr. The distribution isplotted at times of 10 i yr (solid lines, except for the 10 yr curve,which is plotted with a dashed line) and 3 × i yr (all dottedlines), starting at t = × yr. The gray shading denotes thespread in 10 runs. Mass is given in units of monomers. The finalcurve (thick dashed line) corresponds to 5 × yr and overlapsthe 3 × yr curve almost everywhere, indicating that steady-state has been reached. grains. Each curve shows the average of 10 simulations, wherethe gray shading denotes the 1 σ spread in the simulations. At t = N =
1. Thedistribution function f ( N ) gives the number of aggregates perunit volume such that f ( N )d N is the number density of particlesin a mass interval [ N , N + d N ]. Thus, at t = f ( N = , t = = n µ m H / R gd m = . × − cm − for n = cm − and a = . µ m. On the y -axis N f ( N ) is plotted, which shows the mass of the distribution perlogarithmic interval, at several times during the collisional evo-lution. The mass where N f ( N ) peaks is denoted the mass peak:it corresponds to the particles in which most of the mass is con-tained. The peak of the distribution curves stays on roughly thesame level during its evolution, reflecting conservation of massdensity.After 10 yr (first solid line) a second mass peak has ap-peared at N ≃
10. The peak at N = φ σ =
1) size and smaller collisional cross-section of monomerscompared with dimers, trimers. Furthermore, the high collisionalcross section of, e.g., dimers is somewhat overestimated, be-ing the result of the adopted power-law fit between the geomet-rical and collisional cross section (Fig. B.3). These e ff ects aremodest, however, and do not a ff ect the result of the subsequentevolution. Meanwhile, the porosity of the aggregates steadilyincreases, initially by hit-and-stick collisions but after ∼ yrmostly through low-energy collisions between equal size parti-cles (global recipe) that do not visibly compress the aggregates.In Fig. 7 the porosity distribution is shown at several times dur-ing the collisional evolution. Initially, due to low-energy colli-sions the filling factor decreases as a power-law with exponent ≃ . φ σ ≃ N − . . This trend ends after N ∼ , at which timecollisions have become su ffi ciently energetic for compaction tohalt the fractal growth. The filling factor then stabilizes and in-creases only slowly. At t = × yr the N ∼ particles arestill quite porous.After 3 × yr collisions have become su ffi ciently ener-getic for particles to start fragmenting, significantly changingthe appearance of the distribution (Fig. 6). Slowly, particles atlow mass are replenished and growth decelerates. When inquir-ing the statistics underlying the fragmenting collisions, we findthat collisions that result in fragmentation show only modest ero-sion: only a tiny amount of the mass of the large aggregate isremoved. Therefore, at the onset of erosion, growth is not imme-diately halted, but it is e ff ective in replenishing the particles atlow mass. Eventually, at N ∼ ( a ∼ µ m) the erosive colli-sions reach a point at which there is no net growth. With increas-ing time and replenishment, the small particles start to reaccreteto produce a nearly flat distribution in terms of N f ( N ). Sincethe final curve ( t = × yr) mostly overlaps the 3 × curve(both in Figs. 6 and 7) it follows that steady state is reached on ∼ yr timescales – much longer than the timescales on whichmolecular clouds are thought to exist.At 10 yr the largest particles have reached the upper limit of33% for the filling factor (see Fig. 7). Compaction increases thecollision velocities between the particles and therefore enhancesthe fragmentation. The presence of a large population of smallparticles in the steady state distribution also hints that they areresponsible for the higher filling factors particles of intermediatemass (i.e., N ∼ − ) have at steady-state compared withthe filling factor of these particles at earlier times. Indeed, theturnover point at N ∼ corresponds to an energy of ∼ E roll these particles have with small fragments. Compaction by smallparticles is thus much more e ffi cient than collisions with larger(but very flu ff y) particles. mass − fi lli n g f a c t o r , φ σ n = , Ice × & × yr yr × yr yr × yr Fig. 7.
The distribution of the filling factor, φ σ , in the standardmodel, plotted at various times. Initially, the porosity decreasesas a power-law, φ ≃ N − . , the fractal regime. Compaction ismost severe for the more massive particles where the fillingfactor reaches the maximum of 33%. Only mean quantities areshown, not the spread in φ σ . To further understand the influence of the porosity on the colli-sional evolution, the progression of a few key quantities as func-tion of time are plotted in Fig. 8: the mean size h a i , the mass-average size h a i m , and the mass-average filling factor h φ σ i m ofthe distribution. Here, mass-average quantities are obtained byweighing the particles of the Monte Carlo program by mass; e.g., h a i m = P i m i a i P i m i , (14)is the mass-weighted size. The weighing by mass has the ef-fect that the massive particles contribute most, because it is usu-ally these particles in which most of the mass resides. On theother hand, in a regular average all particles contribute equally,meaning that this quantity is particularly a ff ected by the particlesthat dominate the number distribution. Thus, initially h a i m = h a i since at t = h a i m > h a i . This picture is consistentwith the distribution plots in Fig. 6.How sensitive is the emergent size distribution to the adoptedcollision recipe? To address this question we ran simulations inwhich the collision recipe is varied with respect to the standardmodel. The distribution functions of these runs are presented inFig. 9, while Fig. 8 also shows the computed statistical quan-tities (until t = yr). In the case of compact coagulation thefilling factor of the particles was restricted to a minimum of 33%(but small particles like monomers still have a higher filling fac-tor). Clearly, Fig. 8 shows that porous aggregates grow duringthe initial stages (cf. also Fig. 9a and Fig. 9b). These results arein line with a simple analytic model for the first stages of the time [yr] − h φ σ i m h a i / a , h a i m / a h φ σ i m h a ih a i m standardcompacthead-on only Fig. 8. ( solid curves) The mean size h a i ( dashed curves), themass-weighted size h a i m ( dotted curves) and the mass-weightedfilling factor, h φ σ i m ( solid curves) of the size distribution as func-tion of time in the standard model ( black curves), the compactmodel ( dark gray curves) and the head-on only model ( light gray curves).growth, presented in Appendix A.2: the collision timescales be-tween similar size aggregates is shorter when they are porous.Figure 9c presents the results of the standard model in whichcollisions are restricted to take place head-on, an assumptionthat is frequently employed in collision studies (e.g., Wada et al.2008; Suyama et al. 2008). That is, except for the missing col-lision probability ( f miss ), the collision parameters are obtainedexclusively from the b = E / N tot E roll ∼
1) central collisions aremuch more e ff ective in compacting than o ff -center collisions.For the same reason growth in the standard model is also some-what faster during the early stages. However, at later times thedi ff erences between Fig. 9b and Fig. 9c become negligible, in-dicating that head-on and o ff -center collisions do not result in adi ff erent fragmentation behavior.Thus, we conclude that inclusion of porosity significantlyboosts the growth rates on molecular cloud relevant timescales( t = − yr). Studies that model the growth by compact par-ticles of the same internal density will therefore underestimatethe aggregation. O ff -center collisions are important to provide a(net) increase in porosity during the restructuring phase but donot play a critical role at later times. Figure 10a-c give the collisional evolution of silicates at sev-eral densities. In most of the models fragmentation is impor-tant from the earliest timescales on. Due to the much lowerbreaking energy of silicates compared with ice, silicates alreadystart fragmenting at relative velocities of ∼
10 m s − . As a re- .W. Ormel et al.: Dust Coagulation and Fragmentation in Molecular Clouds I. 11 mass, N − − − − − m a ss d e n s i t y , N f ( N ) n = , Ice , compact A × yr mass, N n = , Ice B × yr mass, N − − − − − C n = , Ice , head-on × yr Fig. 9.
The e ff ects of the collision recipe on the evolution of the size distribution. The standard model ( b , shown for comparison) isvaried and features: ( a ) compact coagulation, in which the filling factor is restricted to a lower limit of 33%; ( c ) head-on collisionsonly, where the impact parameter is fixed at b = yr.sult, the growth is very modest: only a factor of 10 in size forthe n = cm − model, whereas at lower densities most of themass stays in monomers. For the same reason, silicates reachsteady state already on a timescale of 10 yr, much faster thanice-coated particles.In the case of ice-coated silicates (Fig. 10d-f) much higherenergies are required to restructure and break aggregates andparticles grow large indeed. In all cases the qualitative pic-ture reflects that of our standard model (Fig. 10e), discussed inSect. 4.1: porous growth in the initial stages, followed by com-paction and fragmentation in the form of erosion. The evolutiontowards steady-state is a rather prolonged process and is onlycomplete within 10 yr in Fig. 10f. In the low density modelof Fig. 10d fragmentation does not occur within 10 yr. Steadystate is characterized by a rather flat mass spectrum.In Fig. 11 the collisional evolution is contrasted for threedi ff erent monomer sizes: a =
300 Å (Fig. 11a), 0 . µ m (thestandard model, Fig. 11b), and 1 µ m (Fig. 11c). To obtain aproper comparison, Fig. 11 uses physical units (grams) for themass of the aggregates, rather than the dimensionless number ofmonomers, N . From Fig. 11 it can be seen that the models areextremely sensitive to the grain size. In Fig. 11c, for example,the weaker aggregates result in the dominance of fragmentingcollisions already from the start. These curves, therefore, resem-ble the silicate models of Fig. 10b.Figure 11a, on the other hand, shows that a reduction of thegrain size by about a factor three ( a = . µ m) enhances thegrowth significantly. Despite starting from a lower mass, the 300Å model overtakes the standard model at t ∼ yr. An under-standing of this behavior is provided in Appendix A.2, the keyelement being the persistence of the hit-and-stick regime fromwhich it is very di ffi cult to break out if a is small. Until 4 × yrvisible compaction fails to take place and aggregates becomevery porous indeed ( φ ≃ × − ). The consequence is that frag-mentation is also delayed, and has only tentatively started nearthe end of the simulations. We caution, however, against the rel- evance of the 300 Å model; as explained in Sect. 3.1, the choiceof a =
300 Å is too low to model aerodynamic and mechanicalproperties of MRN aggregates. But Fig. 11 serves the purpose ofshowing the sensitivity of the collisional result on the underlyinggrain properties.
Tables 3 and 4 present the results of the collisional evolutionin tabular format. In Table 3 the mass-weighted size of the dis-tribution ( h a i m , reflecting the largest particles) is given, and inTable 4 the opacity of the distribution is provided, which reflectsthe behavior of the small particles. Here, opacity denotes geo-metrical opacity – the amount of surface area per unit mass –which would be applicable for visible or UV radiation, but notto the IR. Its definition is, accordingly, h κ i = P i π a σ, i P i m i , (15)where the summation is over all particles in the simulation.These tables show, for example, that in order to grow chondrule-size particles ( ∼ − g), dust grains need to be ice-coated and,except for the n = cm − model, coagulation times of ∼ yrare needed.To further assess the impact of grain coagulation we mustcompare the coagulation timescales to the lifetimes of molecularclouds. In a study of molecular clouds in the solar neighborhoodHartmann et al. (2001) hint that the lifetime of molecular cloudis short, because of two key observational constraints: (i) mostcores do contain young stars, rather than being starless; and (ii)the age of the young stars that are still embedded in a cloud is1 − − mass, N − − − − − m a ss d e n s i t y , N f ( N ) n = , Sil. A × yr × yr × yr × yr × yr × yr mass, N n = , Ice D × yr mass, N − − − − − m a ss d e n s i t y , N f ( N ) n = , Sil. B × yr mass, N n = , Ice E × yr mass, N − − − − − m a ss d e n s i t y , N f ( N ) n = , Sil. C × yr mass, N n = , Ice F × yr Fig. 10.
Distribution plots corresponding to the collisional evolution of silicates ( left panels) and ice-coated silicates ( right panels)at densities of n = , and 10 cm − until t = yr. For the silicates a steady-state between coagulation and fragmentationis quickly established on timescales of ∼ yr, whereas ice-coated silicates grow much larger before fragmentation kicks in. Theinitial distribution is monodisperse at a = − cm. Note the di ff erent x -scaling.large particles produced, or (ii) the removal of small particles.This can be seen from Tables 3 and 4 where h a i m and h κ i are alsogiven at the free-fall time of the simulation (col. 6). From Table 3it is seen that the sizes of the largest particles all stay below 1 µ m(except the models that started already with a monomer sizes of a = µ m). Likewise, Table 4 shows that the opacities from the t ff entry are similar to those of the ‘initial’ 10 yr column, i.e.,growth is negligible on free-fall timescales.This information is also displayed in Fig. 12, where the opac-ity with respect to the initial opacity, κ/κ , is plotted againsttime for all densities from the a = − cm ice-coated sili- cate models. In Fig. 12 time is normalized to the initial colli-sion timescale between two grains, t coll , , which is a function ofdensity (see Eq. (A.4)). The similarity of the curves for the first ∼ t coll , is in good agreement with a simple analytic modelpresented in Appendix A.2. In models where small particles arereplenished by fragmentation, κ first obtains a minimum andlater levels-o ff at κ/κ ∼ .
05. Also in Fig. 12, the free-fall andambipolar di ff usion timescales are indicated with diamond andsquare symbols, respectively. Due to the normalization by t coll , these occur within a relatively narrow region, despite the largerange in densities considered. It is then clear that at a free-fall .W. Ormel et al.: Dust Coagulation and Fragmentation in Molecular Clouds I. 13 − − − − mass, m [g] − − − − − m a ss d e n s i t y , m f ( m ) n = , Ice , a = ˚A A × yr − − − − − mass, m [g] B n = , Ice × yr − − − mass, m [g] − − − − − n = , Ice , a = µ m C × yr Fig. 11.
The e ff ects of a di ff erent grain size a to the collisional evolution: ( a ) a = µ m, ( b ) a = . µ m (the default, shown forreasons of comparison), and ( c ) a = µ m. To facilitate the comparison, physical units are used (grams) for the mass of aggregates,rather than the number of monomers ( N ). Table 3.
Mass-weighted size of the distribution, h a i m , at several distinct events during the simulation run. h a i m [cm]model 10 yr 10 yr 10 yr 10 yr t ff ( n ) t ad ( n )(1) (2) (3) (4) (5) (6) (7) n = , ice 1 . −
5) 1 . −
5) 1 . −
5) 8 . −
5) 1 . −
5) 8 . − n = , silicates 1 . −
5) 1 . −
5) 1 . −
5) 1 . −
5) 1 . −
5) 1 . − n = , ice 1 . −
5) 1 . −
5) 4 . −
5) 8 . −
4) 1 . −
5) 8 . − n = , silicates 1 . −
5) 1 . −
5) 4 . −
5) 4 . −
5) 2 . −
5) 4 . − n = , silicates, a = − . −
4) 1 . −
4) 1 . −
4) 1 . −
4) 1 . −
4) 1 . − n = , ice 1 . −
5) 2 . −
5) 6 . −
4) 7 . −
3) 2 . −
5) 3 . − n = , ice, a = − . −
4) 1 . −
4) 2 . −
4) 2 . −
4) 1 . −
4) 2 . − n = , ice, a = × − . −
6) 1 . −
5) 1 . −
3) 4 . . −
5) 2 . − n = , ice, compact 1 . −
5) 1 . −
5) 1 . −
4) 5 . −
3) 1 . −
5) 1 . − n = , ice, head-on 1 . −
5) 2 . −
5) 3 . −
4) 7 . −
3) 2 . −
5) 3 . − n = , silicates 1 . −
5) 1 . −
4) 1 . −
4) 1 . −
4) 4 . −
5) 1 . − n = , ice 1 . −
5) 2 . −
4) 3 . −
2) 2 . −
2) 4 . −
5) 2 . − n = , ice 7 . −
5) 3 . −
2) 5 . −
2) 6 . −
2) 8 . −
5) 7 . − N ote . Column (1) lists the models in terms of the density ( n ) and material properties. The monomer size ( a ) is 0 . µ m, unless otherwise indicated.Cols. (2)–(5) give the mass-weighted size of the distribution at fixed coagulation times. Likewise, cols. (6)–(7) provide h a i m at the free-fall andthe ambipolar di ff usion timescale of the cloud that corresponds to the gas density n . These are a function of density and are given in Eq. (1) andEq. (2), respectively. Values a × b are denoted a ( b ). timescale no significant reduction of the opactiy takes place,since t ff / t coll , . ff usion (AD),and the relevant timescales are much longer than the free-falltimescale (see Eq. (2)), t AD / t coll , ≫ ∼ µ m in the densest models on an AD-timescale. For the highest density models timescales are evensu ffi ciently long for fragmentation to replenish the small grains.(Note that, although t AD decreases with density, the evolutionof the core is determined by the quantity t AD / t coll , , which in-creases with n .) Thus, if cores evolve on AD-timescales, the ob-servational appearance of the core will be significantly a ff ected.Table 4 and Fig. 12 show that the UV-opacity, which is directlyproportional to κ , will be reduced by a factor of ∼
10. Studiesthat relate the A V extinction measurements to column densities Table 4.
Like Table 3 but for the geometrical opacity κ of the particles. h κ i [cm g − ]model 10 yr 10 yr 10 yr 10 yr t ff ( n ) t ad ( n )(1) (2) (3) (4) (5) (6) (7) n = , ice 2 . . . . . . n = , silicates 2 . . . . . . n = , ice 2 . . . . . . n = , silicates 2 . . . . . . n = , silicates, a = − . . . . . . n = , ice 2 . . . . . . n = , ice, a = − . . . . . . n = , ice, a = × − . . . . . . n = , ice, compact 2 . . . . . . n = , ice, head-on 2 . . . . . . n = , silicates 2 . . . . . . n = , ice 2 . . . . . . n = , ice 1 . . . . . . − t / t coll , − − k / k / e / e / e n = n = n = n = n = t ff t ad yr Fig. 12.
The opacity κ normalized to its initial value vs. time inunits of the initial collision time t coll , (Eq. (A.4)) for the ice-coated, a = . µ m silicates models at five di ff erent gas densi-ties n . The decrease in opacity occurs on timescales of ∼ t coll , .In simulations where small grains reemerge due to fragmentation κ starts to increase again. The free-fall ( diamonds ) and ambipo-lar di ff usion timescales ( squares ) are indicated as far as thesefall within 10 yr ( circles ). Points of low density appear at lower t / t coll , .through the standard dust-to-gas ratio therefore could underesti-mate the amount of gas that is actually present.
5. Discussion
In his pioneering work to the study of dust coagulation in molec-ular clouds, Ossenkopf (1993), like our study, follows the inter-nal structure of particles and presents a model for the changein particle properties for collisions in the hit-and-stick regime.Furthermore, the grains are characterized by an MRN size dis- tribution. The model of Ossenkopf (1993) only treats the hit-and-stick collision regime but at the high densities ( n H ≥ cm − )and short timescales ( ∼ yr) he considers any compactionor fragmentation between ice(-coated) particles is indeed of noconcern. The coagulation then proceeds to produce particles ofcompact size ∼ . µ m at n H = cm − . It can be seen fromTable 3 that the growth in the corresponding model of our study(ice, n = cm − ) is higher: 2 . µ m. This large di ff erence(especially in terms of mass) can be attributed to the fact thatOssenkopf (1993) ignores turbulent relative velocities betweenparticles of friction times τ f < t s . As a result, growth is pre-dominantly PCA, because the small grains can only be swept upby bigger aggregates, rendering his coagulation more compactin comparison to our model and therefore slower. Additionally,due to the di ff erent definitions we use for ‘density’ ( n vs. n H , seefootnote 2) our ‘10 cm − model’ is denser by a factor of 1.7,resulting in a lower collision timescale and faster growth.However, at timescales t > t coll , (where t coll , for a dis-tribution would be the collision time between big grains) hit-and-stick growth will turn into CCA. Consequently, fast growthis expected on timescales larger than a collision timescale (seeAppendix A). By 10 yr this condition has clearly been fulfilledin our n = cm − model, but it is likely that, due to theabove mentioned di ff erences, it has not been met, or perhapsonly marginally, in Ossenkopf (1993). Thus, rather than fixingon one point in time, a more useful comparison would be to com-pare the growth curves, a ( t ).On the other hand, Weidenschilling & Ruzmaikina (1994),adopt a Bonnor-Ebert sphere to model the molecular cloud, andcalculate the size distribution for much longer timescales ( t = yr). Like our study, Weidenschilling & Ruzmaikina (1994)include fragmentation in the form of erosion and, at high ener-gies, shattering. Their particles are characterized by a strengthof Q ∼ erg g − , which are, therefore, somewhat weaker thanthe particles of our standard model. Although their work lacks adynamic model for the porosity evolution, it is assumed that theinitial growth follows a fractal law until 30 µ m. At these sizesthe minimum filling factor becomes less than 1%, lower thanour results. On timescales of ∼ yr particles grow to & µ m,comparable to that of our standard model.A major di ff erence between Weidenschilling & Ruzmaikina(1994) and our works concerns the shape of the size dis-tribution. Whereas in our calculations the mass-peak al- .W. Ormel et al.: Dust Coagulation and Fragmentation in Molecular Clouds I. 15 ways occurs at the high-mass end of the spectrum, in theWeidenschilling & Ruzmaikina (1994) models most of the massstays in the smallest particles. Perhaps, the lack of massive parti-cles in the Weidenschilling & Ruzmaikina (1994) models is theresult of the spatial di ff usion processes this work includes; mas-sive particles, produced at high density, mix with less massiveparticles from the outer regions. In contrast, our findings regard-ing steady-state distributions agree qualitatively with the find-ings of Brauer et al. (2008) for protoplanetary disks. Despite thedi ff erent environments, and therefore di ff erent velocity field, wefind that the steady state coagulation-fragmentation mass spec-trum is characterized by a rather flat m f ( m ) mass function.It is also worthwhile to compare the aggregation results fromour study with the constituent particles of meteorites, chon-drules ( a ∼ µ m) and calcium-aluminium inclusions (CAIs, a ∼ cm). Although most meteoriticists accept a nebular originfor these species (e.g., Huss et al. 2001), Wood (1998) suggestedthat, in order to explain Al-26 free inclusions, aggregates thesizes of CAIs (and therefore also chondrules), formed in the pro-tostellar cloud. These large aggregates then were self-shieldedfrom the e ff ects of the Al-26 injection event. However, our re-sults indicate that growth to cm-sizes seems unlikely. Only thedense ( n ≥ ) models can produce chondrule-size progenitorsand only at a (long) ambipolar di ff usion timescale. In our models we observe that the shape of the initially monodis-perse dust size distribution evolves first to a Gaussian-like dis-tribution and eventually to a flat steady-state distribution. Fortimescales longer than the coagulation timescale (Eq. (A.4)) wecan expect that this result is independent of the initial condi-tions, even if the coagulation starts from a power-law distri-bution. Essentially, these distributions are a direct result of thephysics of the coagulation: the Gaussian-like distribution reflectsthe hit-and-stick nature of the agglomeration process at low ve-locities while the ‘flat’ N f ( N ) distribution at later times re-sults from a balance between fragmentation – erosion but notcatastrophic destruction – and growth. In contrast, in interstel-lar shocks grains acquire much larger relative velocities andgrain-grain collisions will then quickly shatter aggregates intotheir constituent monomers (Jones et al. 1996; Hirashita & Yan2009). Hence, the interstellar grain size distribution will be verydi ff erent in the dense phases of the interstellar medium than inthe di ff use ISM and studies of the e ff ects of grains on the opac-ity, ionization state and chemical inventory of molecular cloudswill have to take this into account.As Fig. 12 illustrates, in our calculations, the average surfacearea of the grain mixture – a proxy for the visual and near-IR ex-tinction – decreases by orders of magnitude during the initialcoagulation process. In a general sense this finding is in agree-ment with observational evidence for the importance of graingrowth in molecular clouds as obtained from studies of dust ex-tinction per unit column density of gas, where the latter is mea-sured either through HI / H UV absorption lines, sub-millimeterCO emission lines, or X-ray absorption (cf. Whittet 2005; Jura1980; Winston et al. 2007; Goldsmith et al. 1997). Obviously,this process is faster and therefore can proceed further in denseenvironments (Fig. 12). As a corollary to this, the decrease in to-tal surface area only occurs for timescales well in excess of thefree-fall timescale. Hence, very short lived density fluctuationsdriven by turbulences will not show this e ff ects of coagulationon the total grain surface area of dust extinction. The study by Chiar et al. (2007) is – at first sight – somewhatat odds with this interpretation. They find that the total near-IRextinction keeps rising when probing deeper into dense coreswhile the strength of the 10 µ m feature abruptly levels o ff at anear-IR extinction value which depends somewhat on the cloudsurveyed. The recent study by McClure (2009) also concludesthat the strength of the 10 µ m feature relative to the local contin-uum extinction decreases dramatically when the K-band extinc-tion exceeds 1 magnitude. Clearly, the two features are carriedby di ff erent grain populations (Chiar et al. 2007). Indeed, mod-els for interstellar extinction attribute the near-IR extinction tocarbonaceous dust grains while the 10 µ m feature is a measure ofthe silicate population (Draine & Lee 1984). Hence, these datasuggest that silicates coagulate very rapidly when a certain den-sity (i.e., depth into the cloud) is reached – essentially hiding sil-icates grains in the densest parts of the cloud from view – whilethe carbonaceous grain population is not (as much) a ff ected. Inhis study, McClure (2009) concludes that icy grains are involvedin this change in extinction behavior with A K . Likely, rather thanthe presence of the 13 µ m ice libration band a ff ecting the silicateprofile, this behavior reflects grain growth. Our study shows thatcoagulation in molecular clouds is greatly assisted by the pres-ence of ice mantles. Once grains are covered by ice mantles, theincreased ‘stickiness’ of ice takes over and the precise character-istics of the core become immaterial. Perhaps, therefore, the datasuggest that silicates rapidly acquire ice mantles while carbona-ceous grains do not. However, there is no obvious physical basisfor this suggestion. Further experimental studies on ice forma-tion on di ff erent materials will have to settle this issue.In this study we discussed observational implications of ourmodel in a very coarse way, i.e., by considering the reduction ofthe total geometrical surface area ( κ ) due to aggregation. We thenfind that its behavior can be largely expressed as function of theinitial collision timescale, t coll , . However, for a direct compar-ison with observations, e.g., the 10 µ m silicate absorption fea-ture, it is relevant to calculate the extinction properties of thedust distribution as function of wavelength, and to assess, forexample, the significance of porous aggregates (Min et al. 2008;Shen et al. 2008). This is the subject of a follow up study.
6. Summary and conclusions
We have studied the collisional growth and fragmentation pro-cess of dust in the environments of the molecular cloud (cores).In particular, we have focused on the collision model and theanalysis of the several collisional evolution stages. Except forbouncing, the collision model features all relevant collisionaloutcomes (sticking, breakage, erosion, shattering). Furthermore,we have included o ff -center collisions in the recipe format andalso prescribe the change to the internal structure in terms of thefilling factor. We have treated a general approach, and outcomesof future experiments – either numerical or laboratory – can beeasily included. The collision model features scaling of the re-sults to the relevant masses and critical energies, which allowsthe calculation to proceeds beyond the sizes covered by the origi-nal numerical collision experiments. Our method is, in principle,also applicable to the dust coagulation and fragmentation stagesin protoplanetary disks.We list below the key results that follow from this study:1. The collisional evolution can be divided into three phases: (i) t < t coll , in which the imprints of growth are relatively mi-nor; (ii) a porosity-assisted growth stage, where the N f ( N )mass spectrum peaks at a well-defined size; and (iii) a frag- mentation stage, where the N f ( N ) mass spectrum is rela-tively flat due to the replenishment of small particles by frag-mentation. Fragmentation is primarily caused by erosive col-lisions.2. A large porosity speeds up the coagulation of aggregates inthe early phases. This e ff ect is self-enhancing, because veryporous particles couple very well to the gas, preventing ener-getic collisions capable of compaction. Growth in the secondregime can therefore become very fast. Grazing collisionsare largely responsible for obtaining flu ff y aggregates in thefirst phases, further increasing the porosity.3. Silicate dust grains or, in general, grains without ice-coatingare always in the fragmentation regime. This is a resultof their relatively low breaking energy. Freeze-out of ices,on the other hand, will significantly shift the fragmentationthreshold upwards, fulfilling a prerequisite for significant ag-gregation in molecular clouds.4. Likewise, the (monodisperse) grain size that enters the col-lision model is also critical for the strength of the resultingdust aggregates. Smaller grains will increase the strength sig-nificantly, due to increased surface contacts. Besides, a coag-ulation process that starts with small grains also results in thecreation of very porous aggregates, which further enhancestheir growth. Although a single grain size cannot fully sub-stitute for both the mechanical and aerodynamic propertiesof a grain size distribution, we have argued that for the MRNdistribution a size of 0 . µ m reflects these properties best.5. If cloud lifetimes are restricted to free-fall times, little co-agulation can be expected since the coagulation timescaleis generally longer than t ff . However, if additional supportmechanism are present, like ambipolar di ff usion, and freeze-out of ice has commenced, dust aggregates of ∼ µ m areproduced, which will significantly alter the UV-opacity ofthe cloud. Conversely, our results reveal that the total dustsurface area (and hence the extinction per H-nuclei) pro-vides a convenient clock that measures the lifetime of a densecore in terms of the initial coagulation timescale. As obser-vations typically reveal that the dust extinction per H-nucleiin dense cores has decreased substantially over that in thedi ff use ISM, this implies that such cores are long-lived phe-nomena rather than transient density fluctuations.6. Despite the complexity of the collision model, we find thatthe decrease in (total) dust opacity can be expressed in termsof the initial collision time t coll , only, providing a relationfor the density and lifetime of the cloud to its observationalstate (Fig. 12). Acknowledgements.
The authors thank V. Ossenkopf for discussion on the re-sults of his 1993 paper. C.W.O. appreciates useful discussions with MarcoSpaans, which helped to clarify certain points of this manuscript. The authorsalso acknowledge the significantly contributions the referee, Vincent Guillet, hasmade to the paper by suggesting, for example, Sect. 3.1, Fig. 1, and Fig. A.2.These, together with many other valuable comments, have resulted in a signifi-cant improvement of both the structure and contents of the manuscript.
References
Aikawa, Y., Wakelam, V., Garrod, R. T., & Herbst, E. 2008, ApJ, 674, 984Akyilmaz, M., Flower, D. R., Hily-Blant, P., Pineau Des Forˆets, G., & Walmsley,C. M. 2007, A&A, 462, 221Alves, J., Lombardi, M., & Lada, C. J. 2007, A&A, 462, L17Bergin, E. A. & Tafalla, M. 2007, ARA&A, 45, 339Binney, J. & Tremaine, S. 1987, Galactic dynamics (Princeton, NJ, PrincetonUniversity Press, 1987, 747 p.)Blum, J. 2004, in ASP Conf. Ser. 309: Astrophysics of Dust, ed. A. N. Witt,G. C. Clayton, & B. T. Draine, 369 Blum, J. 2006, Advances in Physics, 55, 881Blum, J. & M¨unch, M. 1993, Icarus, 106, 151Blum, J. & Schr¨apler, R. 2004, Physical Review Letters, 93, 115503Blum, J., Schr¨apler, R., Davidsson, B. J. R., & Trigo-Rodr´ıguez, J. M. 2006,ApJ, 652, 1768Blum, J. & Wurm, G. 2000, Icarus, 143, 138Blum, J. & Wurm, G. 2008, ARA&A, 46, 21Bowey, J. E., Adamson, A. J., & Whittet, D. C. B. 1998, MNRAS, 298, 131Brauer, F., Dullemond, C. P., & Henning, T. 2008, A&A, 480, 859Butler, M. J. & Tan, J. C. 2009, ApJ, 696, 484Carrasco, L., Strom, S. E., & Strom, K. M. 1973, ApJ, 182, 95Charnley, S. B., Tielens, A. G. G. M., & Millar, T. J. 1992, ApJ, 399, L71Chiar, J. E., Ennico, K., Pendleton, Y. J., et al. 2007, ApJ, 666, L73Chokshi, A., Tielens, A. G. G. M., & Hollenbach, D. 1993, ApJ, 407, 806Derjaguin, B. V., Muller, V. M., & Toporov, Y. P. 1975, Journal of Colloid andInterface Science, 53, 314Desert, F.-X., Boulanger, F., & Puget, J. L. 1990, A&A, 237, 215Dominik, C. & N¨ubold, H. 2002, Icarus, 157, 173Dominik, C. & Tielens, A. G. G. M. 1995, Philosophical Magazine A, 72, 783Dominik, C. & Tielens, A. G. G. M. 1996, Philosophical Magazine A, 73, 1279Dominik, C. & Tielens, A. G. G. M. 1997, ApJ, 480, 647Draine, B. T. 1985, in Protostars and Planets II, ed. D. C. Black & M. S.Matthews, 621–640Draine, B. T. & Lee, H. M. 1984, ApJ, 285, 89Filippov, A., Zurita, M., & Rosner, D. 2000, Journal of Colloid and InterfaceScience, 229, 261Flower, D. R., Pineau Des Forˆets, G., & Walmsley, C. M. 2006, A&A, 456, 215Goldsmith, P. F., Bergin, E. A., & Lis, D. C. 1997, ApJ, 491, 615Hartmann, L., Ballesteros-Paredes, J., & Bergin, E. A. 2001, ApJ, 562, 852Hasegawa, T. I., Herbst, E., & Leung, C. M. 1992, ApJS, 82, 167Hirashita, H. & Yan, H. 2009, MNRAS, 394, 1061Huss, G. R., MacPherson, G. J., Wasserburg, G. J., Russell, S. S., & Srinivasan,G. 2001, Meteoritics and Planetary Science, 36, 975Johnson, K. L., Kendall, K., & Roberts, A. D. 1971, Proceeding of the RoyalSociety A, 324, 301Johnstone, D. & Bally, J. 2006, ApJ, 653, 383Jones, A. P., Tielens, A. G. G. M., & Hollenbach, D. J. 1996, ApJ, 469, 740Jørgensen, J. K., Johnstone, D., Kirk, H., et al. 2008, ApJ, 683, 822Jura, M. 1980, ApJ, 235, 63Kawamura, A., Onishi, T., Yonekura, Y., et al. 1998, ApJS, 117, 387Klessen, R. S., Ballesteros-Paredes, J., V´azquez-Semadeni, E., & Dur´an-Rojas,C. 2005, ApJ, 620, 786Kozasa, T., Blum, J., & Mukai, T. 1992, A&A, 263, 423Langkowski, D., Teiser, J., & Blum, J. 2008, ApJ, 675, 764Mac Low, M.-M. & Klessen, R. S. 2004, Reviews of Modern Physics, 76, 125Markiewicz, W. J., Mizuno, H., & V¨olk, H. J. 1991, A&A, 242, 286Mathis, J. S., Rumpl, W., & Nordsieck, K. H. 1977, ApJ, 217, 425McClure, M. 2009, ApJ, 693, L81Meakin, P. 1988, Ann. Rev. Phys. Chem., 39, 237Meakin, P. & Donn, B. 1988, ApJ, 329, L39Min, M., Hovenier, J. W., Waters, L. B. F. M., & de Koter, A. 2008, A&A, 489,135Ormel, C. W. & Cuzzi, J. N. 2007, A&A, 466, 413Ormel, C. W. & Spaans, M. 2008, ApJ, 684, 1291Ormel, C. W., Spaans, M., & Tielens, A. G. G. M. 2007, A&A, 461, 215Ossenkopf, V. 1993, A&A, 280, 617Paszun, D. & Dominik, C. 2006, Icarus, 182, 274Paszun, D. & Dominik, C. 2008, A&A, 484, 859Paszun, D. & Dominik, C. 2009, A&A submittedSchnee, S., Li, J., Goodman, A. A., & Sargent, A. I. 2008, ApJ, 684, 1228Shen, Y., Draine, B. T., & Johnson, E. T. 2008, ApJ, 689, 260Shu, F. H. 1977, ApJ, 214, 488Stepnik, B., Abergel, A., Bernard, J.-P., et al. 2003, A&A, 398, 551Suyama, T., Wada, K., & Tanaka, H. 2008, ApJ, 684, 1310Tassis, K. & Mouschovias, T. C. 2004, ApJ, 616, 283Tielens, A. G. G. M. 2005, The Physics and Chemistry of the InterstellarMedium (The Physics and Chemistry of the Interstellar Medium, byA. G. G. M. Tielens, pp. . ISBN 0521826349. Cambridge, UK: CambridgeUniversity Press, 2005.)Tielens, A. G. G. M. & Hagen, W. 1982, A&A, 114, 245van Breemen, J. M., Min, M., Chiar, J. E., et al. 2009, A&A submittedV¨olk, H. J., Morfill, G. E., Roeser, S., & Jones, F. C. 1980, A&A, 85, 316Wada, K., Tanaka, H., Suyama, T., Kimura, H., & Yamamoto, T. 2008, ApJ, 677,1296Weidenschilling, S. J. & Cuzzi, J. N. 1993, in Protostars and Planets III, ed. E. H.Levy & J. I. Lunine, 1031–1060Weidenschilling, S. J. & Ruzmaikina, T. V. 1994, ApJ, 430, 713Weidling, R., G¨uttler, C., Blum, J., & Brauer, F. 2009, ArXiv e-prints .W. Ormel et al.: Dust Coagulation and Fragmentation in Molecular Clouds I. 17
Weingartner, J. C. & Draine, B. T. 2001, ApJ, 548, 296Whittet, D. C. B. 2005, in Astronomical Society of the Pacific Conference Series,Vol. 343, Astronomical Polarimetry: Current Status and Future Directions, ed.A. Adamson, C. Aspin, C. Davis, & T. Fujiyoshi, 321– + Wilking, B. A., Lebofsky, M. J., Kemp, J. C., Martin, P. G., & Rieke, G. H. 1980,ApJ, 235, 905Winston, E., Megeath, S. T., Wolk, S. J., et al. 2007, ApJ, 669, 493Wood, J. A. 1998, ApJ, 503, L101Yan, H., Lazarian, A., & Draine, B. T. 2004, ApJ, 616, 895Zsom, A. & Dullemond, C. P. 2008, A&A, 489, 931
Appendix A: Analytical background
A.1. Relativevelocitiesandcollisiontimescalesof dustparticles
The friction time, τ f , sets the amount of coupling between parti-cles and gas. In molecular clouds the Epstein regime is applica-ble for which τ f = π c g ρ g m σ , (A.1)where m is the mass of the particle and σ the average projectedsurface area. For compact spheres Eq. (A.1) scales linearly withradius, but for porous aggregates σ can have a much steeper de-pendence on mass (in the case of flat structures, σ ∝ m ) and τ f a much weaker dependence. For spherical grains of size a anddensity ρ s Eq. (A.1) becomes τ = τ f ( a ) = ρ s a c g ρ g (A.2) = . × yr (cid:18) n cm − (cid:19) − (cid:18) T
10 K (cid:19) − / a . µ m ! , where ρ s = .
65 g cm − is used, applicable for silicates.Any coagulation models requires the relative velocities ∆ v between two arbitrary particles. In turbulence, the motions ofparticles can become very correlated, though; e.g., particles re-act in similar ways to the eddy in which they are entrained. Themean relative motion with respect to the gas, therefore, does nottranslate to ∆ v . V¨olk et al. (1980) have pioneered a study to sta-tistically account for the collective e ff ects of all eddies by di-viding the eddies into two classes – ‘strong’ and ‘weak’ – de-pending on the turn-over time of the eddy with respect to theparticle friction time. Ormel & Cuzzi (2007) approximated theframework of V¨olk et al. (1980) and Markiewicz et al. (1991)to provide closed-form expressions for the relative motion be-tween two solid particles. In general three regimes can be distin-guished:1. The low velocity regime, τ ≤ τ ≪ t s . (Here, τ ≥ τ arethe friction times of the particles). Relative velocities scalewith the absolute di ff erence in friction time, ∆ v ∝ τ − τ .2. The intermediate velocity regime, for which t s ≪ τ ≪ t L .Particle velocities scale with the square root of the largestparticle friction time. The particle motion will not align witheddies of shorter turn-over time. These ‘class II’ eddies pro-vide random kicks to the particle motion – an importantsource for sustaining relative velocities of at least ∆ v ∼ v s .3. The heavy particle regime, τ ≫ t L , in which it is τ thatdetermines the relative velocity.Comparing the friction time of the monomer grains(Eq. (A.2)) with the smallest eddy turnover time, t s (Eq. (6a)),it follows that τ > t s under most molecular cloud conditions. We therefore focus on the intermediate velocity regime. In par-ticular, the relative velocity between two equal solid spheres of1 < τ / t s < Re / is (Ormel & Cuzzi 2007) ∆ v ≈ √ v s τ f t s ! / (A.3) = . × cm s − a . µ m ! / (cid:18) n cm − (cid:19) − / (cid:18) T
10 K (cid:19) / . Thus, velocities between silicate dust particles are ∼
10 m / s, andhave a very shallow dependence on density. The same expressionholds when the silicates are coated with ice mantles that are nottoo thick, i.e., ρ s is still the silicate bulk density. Dust monomersthen collide on a collision timescale of t coll , = (cid:16) n d ∆ v π a (cid:17) − = ρ s a R gd ρ g ∆ v (A.4) = . × yr a . µ m ! / (cid:18) n cm − (cid:19) − / (cid:18) T
10 K (cid:19) − / , where n d is the dust number density and R gd =
100 is the stan-dard gas-to-dust density ratio by mass.Equations (A.3) and (A.4) are only valid for solid parti-cles. However, we can scale these relations to two arbitrary butequal aggregates of filling factor φ σ and (dimensionless) mass N . Because m ∝ N and σ ∝ ( N /φ σ ) / it follows that τ f = N / φ / σ τ ; (A.5a) ∆ v ≃ τ f τ ! / ∆ v = N / φ / σ ∆ v ; (A.5b) t coll = (cid:16) n d ∆ vσ C (cid:17) − ≃ N / φ / σ t coll , , (A.5c)where in the latter equation we substituted for simplicity the ge-ometrical cross section σ for the collisional cross-section σ C andused the monodisperse assumption n d ∝ N − for the dust num-ber density n d (this is of course only applicable for narrow dis-tributions). Thus, Eq. (A.5c) shows that the collision timescaledecreases for very porous particles with low filling factors. A.2. Asimpleanalyticalmodelforthe initialstages ofthegrowth
Despite the complexity of the recipe, it is instructive to approx-imate the initial collisional evolution of the dust size distribu-tion with a simple analytical model (cf. Blum 2004). Figure 7suggests that the initial evolution of φ σ can be divided in tworegimes, where the transition point occurs at a mass N . Initially( N < N ), the filling factor is in the fractal regime, which canbe well approximated by a power-law, φ σ ≃ N − / . We referto the fractal regime as including hit-and-stick collisions (no re-structuring) as well as collisions for which E > E roll but whichdo not lead to visible restructuring, i.e., only a small fraction ofthe grains take part in the restructuring. For N > N the fillingfactor starts to flatten-out. It is di ffi cult to assign a trend for φ σ in the subsequent evolution. Following Fig. 7 we may assumethat initially φ σ stays approximately constant for several ordersof magnitude in N , although at some point it will quickly as-sume its compact value of 33%. A sketch of the adopted porositystructure and the resulting scaling of velocities and timescales ispresented in Fig. A.1, in which it is assumed that the collapse N = N N t r e nd r e l a t i v e t o N fractal growth compaction fragment. −
0? 33% ( ∆ v ) ∼ E roll / m ( ∆ v ) ∼ E br / m φ σ ∆ v , t coll ( ∆ v ) Fig. A.1. ( gray solid line ) A simplified model for the behavior ofthe filling factor with growth. Initially, for ∆ v < ∆ v , the poros-ity decreases (fractal growth regime). This phase is followed bya ‘status quo’ phase where filling factors will be approximatelyconstant. The first compaction event is reached when velocitiesreach ∆ v and fragmentation sets in when relative velocities ex-ceeds ∆ v . ( black solid line ) Trend of the collision velocity andcollision timescale. ( dashed line ) Trend of ( ∆ v ) . The numbersdenote the power-law exponents.of the porous structure takes place after the point where the firsterosive collisions occurs, at N = N .From the collision recipe (Sect. 3.3) we identify the criticalenergies at which visible compaction and fragmentation occur.Compaction requires collisions between similar size particles(i.e., the global recipe) and Fig. 5 shows that the transition ( C φ >
1) occurs at an energy of E / N tot E roll ≃ .
2. Regarding fragmen-tation, the simulations clearly show that small particles are re-plenished in the form of erosive collisions (local recipe). FromFig. 4 we assign an energy threshold of E / N µ E br ≃ .
0. Workingout these expressions and using a typical mass ratio of 3 for theglobal recipe ( N µ = N / ∆ v ) ≃ . E roll / m and ( ∆ v ) ≃ . E br / m , respec-tively. These energy thresholds are also indicated in Fig. A.1.From these expressions and the initial expressions for the rel-ative velocity and the collision timescale (Eqs. (A.3) and (A.4)),the turn-over points N and N can be calculated. We assumethat ∆ v < ∆ v such that a fractal growth regime exist. Then, thefirst transition point is reached at a mass N ∼ ∆ v ∆ v ! = . E roll m ( ∆ v ) ! . = (A.6) = × (cid:18) n cm − (cid:19) . γ
370 erg cm − ! . a . µ m ! − . . Unfortunately, the high powers make the numeric evaluationrather unstable. In our simulations we find that N ∼ .Subsequently, we can write for the second transition point, the − time / t coll,0 m a ss , N N N h m i m h m i Eq. (A.9)
Fig. A.2.
Results of the simple analytic model ( solid line ) andcomparison to the h m i and h m i m statistics of the numerical resultsof our standard model.onset of fragmentation, N , N N ∼ ∆ v ∆ v ! = . E br E roll ! (A.7) = × γ
370 erg cm − ! a . µ m ! E ⋆ . × dyn cm − ! − , which corresponds also well to the results from the simulationfor which N ∼ . In our simulations the first fragmentationinvolves particles that are still relatively porous, such that theassumption in Fig. A.1 about the porosity of the N -particles isjustified. However, once steady-state has been reached, particlesof N ∼ will have a 33% filling factor (see Fig. 7).Using again the monodisperse assumption we also obtainthe timescales t , t at which these transition points are reached.Writing,d N d t = Nt coll = N / φ − / σ t coll , , (A.8)where Eq. (A.4) is inserted for t coll , we insert the power-law de-pendence of φ σ on N to obtain t . Straightforward integrationgives tt coll , = Z N N ′− / d N ′ + N − / Z N N N ′− / d N ′ , (A.9)see Fig. A.2, where we plotted the results of Eq. (A.9) togetherwith the averaged mass and the mass-averaged mass of the dis-tribution of the standard model. It follows that the fractal growthstages takes ∼ t coll , , or ∼ × yr (cf. ∼ × yr in thesimulation), while N is reached at ∼ t coll , (cf. ∼ t coll , in thesimulation). At larger time our fit may overestimate the growthrates somewhat because it assumes the filling factor stays fixed atits low value. Overall, the model nicely catches the trends of thegrowth and is also in good agreement with previous approaches(Blum 2004), although, being a monodisperse model, it cannotfit both h m i and h m i m . .W. Ormel et al.: Dust Coagulation and Fragmentation in Molecular Clouds I. 19 Fig. B.1.
Sketch of the initial setup of our simulations. The keyinput parameters ∆ v , b , and φ σ are illustrated. Appendix B: The numerical collision experiments
The skeleton of our collision model consists of the out-comes of many aggregate-aggregate collision simulations. Inthis appendix we briefly review the setup of the simulations(Appendix B.1), discuss some auxiliary relations required tocomplete the collision model (Appendix B.2), discuss the out-put format of the binary collision model (Appendix B.3), andpresent a few limitations that arise due to our reliance on theoutcomes of the numerical simulations (Appendix B.4).
B.1. Collisionsetup and inputparameters
Collisions between aggregates are modeled using the soft ag-gregates numerical dynamics (SAND) code (Dominik & N¨ubold2002; Paszun & Dominik 2008). This code treats interactionsbetween individual grains held together by surface forces in acontact area (Johnson et al. 1971; Derjaguin et al. 1975). TheSAND code calculates the equation of motion for each grain in-dividually and simulates vibration, rolling, twisting, and slidingof the grains that are in contact. These interactions lead to en-ergy dissipation via di ff erent channels. When two grains in con-tact are pulled away, the connection may break, causing loss ofthe energy. Monomers may also roll or slide over each other,through which energy is also dissipated (Dominik & Tielens1995, 1996, 1997). For further details regarding this model andtesting it against laboratory experiments we refer the reader toPaszun & Dominik (2008).To provide both a qualitative and a quantitative descriptionof a collision, Paszun & Dominik (2009) have performed a largenumber of simulations, exploring an extensive parameter space.They formulate a simple collision recipe that quantifies how ki-netic energy, compactness, and mass ratio a ff ect the outcomeof aggregate-aggregate collisions. These results are presented intabular format (Appendix B.3). Figure B.1 sketches the setup ofthese numerical experiments, illustrating three parameters thatshape the outcome of a collision: the initial compactness as rep-resented by the geometrical filling factor φ σ (see below, Eq. (7)),the collision energy E , and the impact parameter b . A collisionfor each parameter set is repeated six times at di ff erent orienta-tions, providing information on the range of outcomes. Becauseof the occasionally flu ff y structure of the aggregates not all ori-entations result in a collision hit, especially not those at largeimpact parameter. An overview of the parameter ranges is givenin Table B.1. The radius of the monomer grains in the simula- Table B.1.
Parameters used in the numerical simulations. ∆ v b / b max φ σ N / N [m / s](1) (2) (3) (4)0.05 0.0 0.070 1.00.30 0.25 0.090 10 − N ote . (1) relative velocity; (2) impact parameter normalized to the sumof the outer radii b max ; (3) geometrical filling factor; (4) mass ratio. tion is a = . µ m and the adopted material properties reflectsilicates ( γ =
25 erg cm − , E = . × dyn cm − ).Relative velocities ∆ v are chosen such that all relevant colli-sion regimes are sampled: from the pure hit-and-stick collisions,where particles grow without changing the internal structure ofthe colliding aggregates, up to catastrophic destruction, wherethe aggregate is shattered into small fragments. In the interme-diate energy regime, restructuring without fragmentation takesplace. For aggregate collisions at large size ratios, high velocityimpacts result in erosion of the large aggregates.The impact parameter b is also well sampled. We probe cen-tral collision ( b = b ≈ b max ), where particles can be stretched due toinertia, and several intermediate cases. In Table B.1 the impactparameter is defined relative to the outer radius of the particles, b max = a out , + a out , . Here the outer radius a out is the radius ofthe smallest sphere centered at the center-of-mass of the particlethat fully encloses it. In the Paszun & Dominik (2009) study theoutcomes of a collision are averaged over the impact parame-ter b . However, in a Monte Carlo treatment, the averaging overimpact is not necessary, and the normalized impact parameter˜ b = b / b max can be determined using a random deviate r , i.e.,˜ b = r / . We have chosen to use the raw data sampled by thePaszun & Dominik (2009) parameter study, explicitly including˜ b as an input parameter for the Monte Carlo model. In this waythe e ff ects of o ff -center impacts can be assessed, i.e., by compar-ing it to models that contain only head-on collisions.The internal structure of the aggregates, or initial com-pactness, is characterized by the geometrical filling factor, φ σ (Eq. (7)). To obtain φ σ , the projected surface area, σ , is aver-aged over a large number of di ff erent orientations of the particle.The parameter space of the filling factor φ σ is chosen such thatwe sample very porous, fractal aggregates that grow due to theBrownian motion (Blum & Schr¨apler 2004; Paszun & Dominik2006), through intermediate compactness aggregates that formby particle-cluster aggregation (PCA), up to compact aggregatesthat may result from collisional compaction.The final parameter that determines a collision outcome isthe mass ratio, N / N ( N being the larger aggregate). ThePaszun & Dominik (2009) experiments sample a mass ratio be-tween 1 and 10 − . In this study, however, we will only use thecollisions corresponding to the two extreme values (i.e., massratios of 1 and 10 − ) as representatives of two distinct classes: global and local . Fig. B.2.
The geometrical filling factor as a function of fragmentmass. Many simulations with di ff erent sets of parameters areoverplotted. The dashed line indicates the least square power-law fit, φ σ ≃ ( m / m ) − . . B.2. Auxiliaryrelationsfor thecollisionrecipe
There are a few more relations required for a consistent approachto the collision recipe. These are presented here for complete-ness.
B.2.1. The filling factor of small fragments
A common filling factor can be assigned to the small fragmentsproduced by erosive or fragmenting collisions that constitute thepower-law component. The compactness of these particles de-pends only on mass and is presented in Fig. B.2, where frag-ments produced in many simulations, reflecting a variety of col-lision properties, are plotted. Almost all particles are locatedalong the power-law with the slope of − .
33. This provides aneasy prescription for the filling factor of small fragments. Thedependence indicates a fractal structure (with the fractal dimen-sion of about D f ≈ .
0) of aggregates formed in a fragmentationevent, since non-fractal aggregates would have a filling factorindependent of mass.As shown by Paszun & Dominik (2009), after reaching themaximum compaction, further increase of the impact energyproduces more restructuring and results in a flattening of theproduced aggregate. Therefore, very flu ff y particles can be pro-duced in collisions of massive aggregates, where the power-lawcomponent extends to larger N . This behavior is also observedin Fig. B.2, where flu ff y, small fragments follow the power-lawrelation, while some large, still compact particles are above thedashed line. B.2.2. Relation between a out and a σ In this study we characterize aggregates by two di ff erent radii:the outer radius a out and the projected surface equivalent radius a σ . The first is used as a reference to the impact parameter b , Fig. B.3.
The geometrical filling factor dependence on the ratioof outer to geometrical radii. In this figure we plot φ σ N . toscale the data for aggregates of di ff erent mass.i.e., the collision o ff set is determined relative to the largest im-pact parameter, b max = a out , + a out , . The cross-section equivalentradius a σ defines our structural parameter φ σ (see Eq. (7)). Wedetermine the relation between the two radii ( a out and a σ ) em-pirically. Both a out and a σ are determined for many aggregatesof various shape and mass. We sample particles with the frac-tal dimension in the range of D f = . . . ff erentaggregates versus the ratio of the outer radius over the cross-section equivalent radius. Diamonds correspond to the producedaggregates of di ff erent fractal dimension and mass. The mass de-pendence in Fig. B.3 is taken into account by plotting φ σ N . .In this way the data for all aggregates are well confined along asingle curve. At small a out / a σ the curve decreases very steeplywith increasing a out / a σ . This corresponds to compact particlesfor which a out / a σ is insensitive to filling factor. The line, how-ever, breaks at about a out / a σ ≈ . − .
3. This shallow relation represents flu ff yaggregates that show a large discrepancy between the projectedsurface equivalent radius and the outer radius.In order to provide a simple relation between the two radii,two power-law functions are fitted to the two regimes: compactparticles below a out / a σ = . ff y aggregates above thatlimit. These two functions are given by φ compact σ = a out a σ ! . − .
21 log N (B.1a) φ flu ff y σ = . a out a σ ! − . N − . . (B.1b)To further verify these relations we use particles produced in sev-eral simulations performed by Paszun & Dominik (2009). These .W. Ormel et al.: Dust Coagulation and Fragmentation in Molecular Clouds I. 21 aggregates are indicated in Fig. B.3 by black squares. They showa similar relation to the one obtained in Eq. (B.1). Points thatare slightly shifted above the fitted lines correspond to aggre-gates that are partly compressed (they did not reach the maxi-mum compaction). Their compact cores are still surrounded bya flu ff y exterior that causes a small increase of the ratio of theouter radius over the projected surface equivalent radius a out / a σ .This behavior, however, occurs at a relatively small value of a out / a σ <
2. At a larger size ratio the filling factor falls backonto the power-law given in Eq. (B.1b).We remark here that, although the fits present a general pic-ture, situations where a out ≫ a σ are not likely to material-ize when N ≫
1. This would corresponds to very open frac-tals of fractal dimension less than two. Instead, in our mod-els φ σ N . & . a out ∼ a σ .Consequently, the fraction of missing collisions f miss is close tozero in most of the cases. B.2.3. Hit and stick
At very low energies ( E ≤ E roll ) two aggregates will stickwhere they meet, without a ff ecting the internal structure of theparticles. This is the ‘hit-and-stick’ regime in which the colli-sional growth can often be described by fractal laws. Two im-portant limits are cluster-cluster coagulation (CCA) and particle-cluster coagulation (PCA). In the former, two particles of equalsize meet, which often results in very flu ff y structures, whereasPCA describes the process in which the projectile particles aresmall with respect to the target. The filling factor then saturatesto a constant value. In the case of monomers, the filling factorwill reach 15% (Kozasa et al. 1992).In general particles do not merely collide with either similar-size particles or monomers. Every size-ratio is possible and leadsto a di ff erent change in filling factor. Ormel et al. (2007) pro-vide an analytical expression, based upon fits to collision exper-iments of Ossenkopf (1993), that give the increase in void spaceas function of the geometrical volume of the collision partners.Here, the geometrical volume V is the volume that correspondsto the geometrical radius, a σ . In this study additional numeri-cal collision experiments were used to further constrain theseanalytical fits. These experiments involved several ‘monomer-bombardments’ of aggregates of di ff erent initial filling factor.Using these data, we fit the volume increase as V void V = max ( V + V ) + V V ! δ/ − − , N . φ exp − V V ! . , (B.2)where V void is the increase in void space (leading to a lower φ σ ), V > V the geometrical volumes of the collision partner, N thenumber of grains in the smaller particle, and V the monomervolume. The first term converges to CCA in the limit of V = V ,and is the same as in Ormel et al. (2007). Here, δ = .
95 is anexponent that reflects the fractal growth in this limit (Ossenkopf1993). The second expression converges to PCA in the limit of V ≪ V . The rationale of providing a second expression is thatin the case of V ≪ V (PCA) the first expression goes to zerovery quickly (no voids are added), which is inconsistent with thePCA limit of 15%. From the results of our new collision exper-iments we have introduced an exponent of 0 .
25 to the PCA-partof Eq. (B.2), which softens the decrease of V void with increasingmass ratio. Table B.2.
Example of an output table from the online data ( f pwl at b = ε φ ini σ . −
4) 0.00000 0.000 0.000 0.0002 . −
2) 0.00000 2.500(-3) 0.000 0.0005 . −
2) 8 . −
4) 0.000 0.000 0.0000.1287 9 . −
2) 3.042(-2) 8.750(-3) 2.917(-2)0.2288 0.3888 9.417(-2) 3.042(-2) 1.500(-2)0.9153 0.9575 0.6033 0.2438 0.11583.661 1.000 1.000 1.000 0.72718.238 1.000 1.000 1.000 1.00014.65 1.000 1.000 1.000 1.000However, not all numerical experiments could be fittedequally well. In fact, we had to compromise. It is likely that abetter fit involves more parameters, e.g., the elongation of theaggregates or their fractal dimension. Here, we have adopted ap-proximate fits that follow the qualitative picture in both the CCA( V = V ) and the PCA ( V ≪ V ) limit. Remark, finally, thatfor the molecular cloud environment the hit-and-stick regime isonly relevant in the initial stages of coagulation at densities of n ≥ cm − or grain sizes a < . µ m. B.3. The collisiontables
Given the level of complexity, it is not feasible to provide sim-ple analytical expressions for the collision outcome (in terms ofthe parameters listed in Table 3.3) as function of the collision pa-rameters ( E , φ σ , ˜ b , N / N ). Therefore, like in Paszun & Dominik(2009), the results are expressed in a tabular format. In total72 tables are provided. They describe six output quantities (seeTable 3.3) for six impact parameters b and for both the local andthe global recipes. Since listing all these tables here is impracti-cal, we will provide them in the digital form as online material accompanying this manuscript. We present two examples to il-lustrate the format.Each table lists one output quantity as function of the dimen-sionless energy parameter ε and the initial filling factor of ag-gregates φ σ . The only exception concerns the fraction of missedcollisions, f miss . This quantity provides a correction to the col-lision cross-section of particles, in our case calculated from theouter radius of an aggregate a out (cf. Appendix C). The fillingfactor φ σ is not an appropriate quantity to use here, because it isambiguous where it concerns the structure of particles. For ex-ample, low φ σ could mean either a very fractal structure (andcorrespondingly high number of missing collisions) or a porousbut homogeneous structure (and low number of missing colli-sions). Therefore, it is more appropriate to relate the probabilityof a collision miss to the radii with which the particle is charac-terized. Thus, f miss is provided as a function of the ratio of theouter radius over the projected surface equivalent radius, a out / a σ .Each table is preceded by a header that specifies: the corre-sponding recipe (keyword: global or local ), the correspondingimpact parameter b , and the quantity listed in the table (key-words are: fmiss , Nf , Sf , fpwl , q , Csig ). In the case of Table B.2the header is
Therefore, Table B.2 presents the fraction of mass in the power-law component, f pwl , for the global recipe and for head-on colli-sion. Table B.3.
Example of an output table from the online data ( f pwl at b = ε φ ini σ ε and the initial filling factor φ ini σ (or the ratio of the outer over the geometrical radii a out / a σ in thecase of f miss ), respectively. Here, ε denotes the collision energyscaled by a normalization constant that involves (i) the breaking or rolling energy and (ii) the reduced or total number of par-ticles, see Sect. 3.2 and Table 3.3. In the case of Table B.2 thescaling parameter is ε = E / E br N tot .Table B.3 is the second example. It is taken from the localrecipe and it presents the f pwl quantity for the head-on collision.The dimensionless energy parameter ε has fewer entries in thelocal recipe tables than in the global recipe. In Table B.3 the en-ergy is scaled by reduced number of monomers N µ (local recipescaling) and by the breaking energy E br (erosion scaling) as in-dicated in Table 3.3. The header in this case is Note that in the local recipe the filling factors are lower. In thiscase larger aggregates are used to model collisions at large massratio, N / N = . The fractal structure of these aggregatesresults in a lower filling factor. B.4. Limitationsof thecollisionrecipe
The main limitation of the collision recipe is that, due to com-putational constraints, the binary aggregate simulations can onlysimulate aggregates of a mass N . . For the recipe to be-come applicable for large aggregates scaling of the results of thecollision experiments is required. This is a critical point of therecipe for which suitable dimensionless quantities had to be de-termined. However, the extrapolation assumes that the collisionphysics that determines the outcomes of collisions at low- N alsoholds at large scales. This is a crucial assumption in which colli-sional outcomes like bouncing are a priori not possible becausethese do not take place at the low- N part of the simulations.Bouncing of aggregates is observed in laboratory exper-iments (Blum & M¨unch 1993; Blum 2006; Langkowski et al.2008; Weidling et al. 2009), whereas it does not occur in oursimulations. For silicates, bouncing occurs at sizes above ap-proximately 100 µ m (i.e., N > particles) and is not fullyunderstood from a microphysical perspective. In the case of ice-coated silicate grains, which provide stronger adhesion forces,our simulations show that growth proceeds to ∼ µ m sizes.In this case, therefore, bouncing might slow down the growthearlier than observed in our experiments, especially when theinternal structure has already re-adjusted to a compact state.However, it is presently unclear how these laboratory experi-ments apply to ice aggregates and hence whether and to whatextent the results would be a ff ected by bouncing. We recognizethat this may, potentially, present a limitation to the growth of aggregates in molecular clouds, but also emphasize it will not af-fect the main conclusions from this study as in only a few modelsaggregates grow to sizes ≫ µ m.Another assumption of the collision model is that the grainshave a spherical geometry. Again, computational constraints ruleout numerical modeling of randomly shaped particles. Whethererratically shaped grains would help or harm the sticking orbouncing is unclear. Because the strength of an aggregate is de-termined by the amount of contact area between two grains, thestrength of irregularly shaped monomers depends on the localradius of curvature. Therefore, highly irregular grains are heldby contacts of much smaller size, because they are connectedby surface asperities. On the other hand, irregular grains mayform more than one contact. However, the geometry of the grainsdoes not necessarily pose a bottleneck to the validity of the col-lision model. Instead, like the size distribution, the consequenceof irregularly shaped monomers is reflected in a di ff erent energyscaling. Appendix C: The Monte Carlo program
The advantage of a Monte Carlo (MC) approach to the calcula-tion of the collisional evolution is that collisions are modeledindividually and that they, therefore, bear a direct correspon-dence to the collision model. Furthermore, in a MC approachstructural parameters (like φ σ ) can be easily included and thecollisional outcome can be quantified in detail. From the twoparticle properties ( N and φ σ ) the collisional quantities are de-rived, e.g., the relative velocities ∆ v between the aggregates (seeAppendix A.1). Then, from ∆ v and the particle’s outer radii wecalculate the collision rates between all particles present in theMC simulation. After the MC model has selected the collisionpartners, the collision recipe is implemented. First, the particleproperties ( m , φ σ ) and the collision properties ( ∆ v ) are turnedinto a collision ‘grid point’ given by the dimensionless ε, φ σ and˜ b . The six collision quantities (Table 3.3) are then taken from theappropriate entries from the recipe tables. Finally, these quanti-ties specify the change to the initial particle properties ( m , φ σ )and also describe the properties of the collision fragments.By virtue of the scaling relations discussed in Sect. 3.2 theMC coagulation model is able to treat much larger aggregatesthan the binary collision experiments. A broad size distribu-tions, which may, e.g., result from injection of small particlesdue to fragmentation, can, however, become problematic for aMC approach, since the high dynamic range required consumescomputational resources. To overcome this problem we use thegrouping method outlined by Ormel & Spaans (2008). In thismethod the 1-1 correspondence between a simulation particleand a physical particle is dropped; instead, the simulated par-ticles are represented by groups of identical physical particles.The group’s mass is determined by the peak of the m f ( m ) massdistribution – denoted m p – and particles of smaller mass ‘travel’together in groups of total mass m p . Grouping entails that a largeparticle can collide with many small particles simultaneously –a necessary approximation of the collision process.Below, we present the way in which we have imple-mented the collision recipe with the grouping method ofOrmel & Spaans (2008). C.1. Collisionrates
The cycle starts with the calculation (or update) of the collisionrates between the groups of the simulation. The individual col-lision rate between two particles i and j is C i j = K i j / V (units: .W. Ormel et al.: Dust Coagulation and Fragmentation in Molecular Clouds I. 23 s − ), where V is the simulation volume and K the collision ker-nel. For grouped collisions C i j is larger because many particlesare involved in the collision. The collision kernel K is definedas K i j = σ C i j ∆ v i j with σ C i j = π ( a out , + a out , ) the collisionalcross section (uncorrected for missing collisions) and ∆ v i j theaverage root-mean-square relative velocity (See Appendix A.1).Thus, to calculate the collision rates we need the relative veloci-ties and the relation between the geometrical and the outer radius(Appendix B.2.2). C.2. Determinationofcollisionpartners
Random numbers determine which two groups collide and thenumber of particles that are involved from the i and j groups, η i and η j . Then, each i -particle collides with η j /η i j -particles. Thegrouping method implicitly assumes that collision rates do notchange significantly during the collision process. To enforce theplausibility of this assumption the grouping method limits thetotal mass of the j -particles colliding with the i -particle to be atmost 1% of the mass of an i -particle, i.e., η j m j /η i m i . f ε = − .Therefore, grouped collisions occur only in the local recipe. Forerosion or sticking this procedure works as intended. However,in collisions that result in breakage the grouping assumption ispotentially problematic, since the particle properties – and hencethe collision rates – then clearly change significantly over a sin-gle collision. Fortunately, in the local recipe breakage is rela-tively unimportant. Catastrophic disruptions (shattering) is prob-lematic for the same reasons, because when it occurs, there is no‘large’ aggregate left. However, for energetic reasons we expectthat shattering occurs mainly when two equal size particles areinvolved, for which the global recipe would apply (and no group-ing). In the following we continue with a collision of η t = η j /η i j -particles colliding with a single i -particle. C.3. Determiningthe collisionquantitiesingroupedcollisions
When the collision is in the ‘hit-and-stick’ regime the propertiesof the new particles are easily found by adding the masses ofthe j -particles to the i -particle and calculating their filling fac-tor using Eq. (B.2). We therefore concentrate here on the localor global recipe. The collision is then characterized by the threedimensionless parameters: normalized collision energy ε , fillingfactor φ σ and impact parameters ˜ b (Sect. 3.2). These three pa-rameters constitute an arbitrary point in the 3D ( ε, φ σ , b )-space,and will in general be confined by eight grid points ( k ) whichcorrespond to the parameters at which results from the collisionexperiments are available, see Fig. C.1. We next distribute the η t collisions over the grid point in which the weight of a gridpoint is inversely proportional to the ‘distance’ to ( ε, φ σ , b ) asexplained in Fig. C.1. Taking account of the collisions that re-sult in a miss, we have η t = η miss + X k = η k ; η miss = X k = η miss , k , (C.1)where η miss , k ≃ η t P k f miss , k denotes the number of collisions atthe grid point resulting in a miss. Here, P k denotes the weight ofthe grid point ( P k P k = f miss the fraction of missed collisionsat the grid point and the ≃ sign indicates this number is roundedto integer values. Similarly, the number of ‘hits’ at a grid point isgiven by η k ≃ η t P k (1 − f miss , k ). Not all of these grid points haveto be occupied (i.e., η k can be zero). In the special case withoutgrouping η t = ( ε, φ, b ) b φ ε P ( ε ) = ε − εε − ε P ( ε ) = ε − ε ε − ε ε ε ε − ε ε − ε Fig. C.1.
Illustration of the picking of the grid points. The colli-sion takes place at ( ε, φ σ , ˜ b ): a point that is generally surroundedby eight grid points (here corresponding to the nodes of the cube)at which results from the binary collision simulations are avail-able. Each node is then assigned a probability inversely propor-tional to the distance to the grid point. Thus, the probability thatthe energy parameter ε = ε is picked (corresponding to four ofthe eight grid nodes) is P = ( ε − ε ) / ( ε − ε ). The procedure isidentical for the other grid points.We continue with the general case of multiple occupied gridpoints. First, we consider the mass that is eroded, given by the f pwl , k quantities. The mass eroded at one grid point per collisionis given by M pwl , k = f pwl , k ( m i + m j ) (global recipe) or M pwl , k = f pwl , k m i m j / ( m i + m j ) (local recipe). Then, the total mass erodedby the group collision is M pwl = X k = M pwl , k η k . (C.2)If this quantity is more than m i , clearly there is no large fragmentcomponent. Otherwise, the mass of the large fragment compo-nent is M large = m i + ( η t − η miss ) m j − M pwl . Each M pwl , k quan-tity is distributed as a power-law with the exponent provided bythe slope q k of the grid point (see below). Concerning the large-fragment component, there is a probability that it will break,given by the N f , k and S f , k quantities. As argued before, breakagewithin the context of the grouping algorithm cannot be consis-tently modeled. Notwithstanding these concerns, we choose toimplement it in the grouping method. Because its probability issmall, we assume it happens at most only once during the groupcollision. The probability that it occurs is then P = − Y k = (1 − P , k ) η k , (C.3)where P , k is the probability that breakage occurs at a grid pointand follows from the S f and N f quantities. If breakage occurs, Recall that in grouped collisions ( η t >
1) this implies that the group-ing method is not fully accurate as the change in mass is of the order ofthe mass itself; but the procedure is always fine if η t = the masses M pwl are removed first and we divide the remainingmass M large in two.The last quantity to determine is the change in the fillingfactor of the large aggregate, denoted by the C φ symbol for onecollision. Like Eq. (C.3) we multiply the changes in C φ at theindividual grid nodes, φ σ, large = h φ σ i m Y k = C η k φ, k , (C.4)This completes the implementation of the collisional outcomewithin the framework of the grouping mechanism. That is, wehave the masses and the filling factor of the large fragment com-ponent ( M large , φ σ, large ), and have computed the distribution ofthe power-law component in terms of mass. Recall, finally, thatall these results are per i -particle, and that the multiplicity of theresults is η i . C.4. Pickingof the power-lawcomponentmasses
The final part of the MC cycle is to pick particles according tothe power-law distribution, under the constraints of a total mass η k M pwl , k and slope q k at each grid point k . The general formulafor picking the mass of the small fragments is m small = m + r m rem m ! + q − / (1 + q ) , (C.5)where m rem is the remaining mass of the distribution (it startsat m rem = η k M pwl , k and decreases every step by m small ) and r a random number between 0 and 1. From the definition of thepower-law component m small cannot be more than 25% of thetotal mass. In the MC program the number of distinct fragmentsthat can be produced is limited to a few per grid point. This is toprevent an influx of a very large number of species (non-identicalparticles; in this case, particles of di ff erent mass), which wouldlead to severe computational problems, filling-up the statevectorarray (see below). Therefore, if the same mass m small is pickedagain it is considered to be the same species, and the multiplicityof this species is increased by one. After we have obtained amaximum of η dis distinct species, we redistribute the mass M pow over the species. In this way the fragment distribution is onlysampled at a few discrete points. C.5. Merging/Duplication
The final part of the MC program consist of an inventory, andpossible adjustment, of the amount of groups and species ( N s )present in the program (Ormel & Spaans 2008). To combinea su ffi ciently high resolution with an e ffi cient computation interms of speed is one of the virtues of the grouping method.One key parameter, determining the resolution of the simula-tion, is the N ∗ s parameter (the target number of species in a sim-ulation). In order to obtain a su ffi cient resolution we require thata total mass of m p ( t ) N ∗ s is present in the simulation at all times,where m p ( t ) is the mass peak of the distribution, m p = h m i / h m i .Particles are duplicated to fulfill this criterion, adding mass tothe system. To prevent a pileup of species we followed the ‘equalmass method’ (Ormel & Spaans 2008). However, we found thatdue to the fragmentation many species were created at anyrate – too many, in fact ( N s > N ∗ s ) which would severely af-fect the e ffi ciency of the program. Therefore, when N s = N ∗ s was reached we used the ‘merging algorithm’ (Ormel & Spaans 2008) to combine neighboring species, a process that averagesover their (structural) parameters but conserves mass. This sig-nificantly improved the e ffi ciency (i.e., speed) of the simulation,although the many fragments created by the collisions (all con-tributing to a higher N s ) can be regarded as a redundancy, be-cause it requires a lot of subsequent regrouping. The alterna-tive would be to produce only 1 new species per collision event(see Zsom & Dullemond 2008). Here, we prefer to stick with amore detailed representation of each collision event by creating arange of particles, but we acknowledge that this amount of detailis to some extent lost by the subsequent merging. .W. Ormel et al.: Dust Coagulation and Fragmentation in Molecular Clouds I. 25 Appendix D: List of symbols
Symbol Description E ∗ Reduced modulus of elasticity (Eq. (9)) R gd Gas-to-dust ratio by mass ∆ v Relative velocity (Appendix A) γ Surface energy density (Eq. (9)) η Number of particles or groups (Appendix C) φ σ Geometrical filling factor (Eq. (7)) µ Molecular mass (Sect. 2) ν m , ν t Molecular, turbulent viscosity (Sect. 2.1 ξ crit Critical displacement for irreversible rolling (Eq. (9)) ρ s Material density, ρ s = .
65 g cm − (silicates) ρ g Gas density, ρ g = µ nm H σ Average projected surface area (Sect. 3) σ C12
Collisional cross section (Sect. 3) τ f Friction time (Eq. (A.1)) C φ Change in geometrical filling factor, C φ = φ σ /φ ini σ (Sect. 3.3) D f Fractal dimension E Collision energy, E = m µ ( ∆ v ) E roll Rolling energy (Eq. (9)) E br Breaking energy (Eq. (9)) N Number of grains in aggregate (dimensionless measureof mass) N µ Reduced number of monomers in collision N µ = N N / ( N + N ) N f Number of big fragments N tot Total number of monomers in collision, N tot = N + N Re Reynolds number (Eq. (5)) S f Spread in number of fragments of big component(Sect. 3.3)St Particle Stokes number (Appendix A) T Temperature (Sect. 2.1) a Monomer radius a out Aggregate outer radius (Fig. 1) a σ Aggregate geometrical radius (projected surface equiva-lent radius) (Fig. 1) a µ Reduced radius (Eq. (9)) b Impact parameter b max Sum of the outer radii, b max = a , out + a , out ˜ b Normalized impact parameter, ˜ b = b / b max c g Sound speed (gas) f miss Fraction of collision misses (Sect. 3.3) f pwl Fraction of mass in power-law component (Sect. 3.3) n Particle density (gas) m Particle mass m µ Reduced mass m H Hydrogen mass q Power-law exponent (size distribution) (Sect. 3.3) r Random deviate t ad Ambipolar di ff usion time (Eq. (2)) t coll , Initial collision time (Eq. (A.4)) t ff Free-fall time (Eq. (1)) t s Inner (Kolmogorov) eddy turn-over time (Eq. (6a)) v L Large eddy turn-over velocity (Sect. 2.1) v ss