Issue 
A&A
Volume 525, January 2011



Article Number  A11  
Number of page(s)  16  
Section  Interstellar and circumstellar matter  
DOI  https://doi.org/10.1051/00046361/201015228  
Published online  26 November 2010 
Dust size distributions in coagulation/fragmentation equilibrium: numerical solutions and analytical fits
MaxPlanckInstitut für Astronomie, Königstuhl 17, 69117
Heidelberg,
Germany
email: birnstiel@mpia.de
Received:
17
June
2010
Accepted:
9
September
2010
Context. Grains in circumstellar disks are believed to grow by mutual collisions and subsequent sticking due to surface forces. Results of many fields of research involving circumstellar disks, such as radiative transfer calculations, disk chemistry, magnetohydrodynamic simulations largely depend on the unknown grain size distribution.
Aims. As detailed calculations of grain growth and fragmentation are both numerically challenging and computationally expensive, we aim to find simple recipes and analytical solutions for the grain size distribution in circumstellar disks for a scenario in which grain growth is limited by fragmentation and radial drift can be neglected.
Methods. We generalize previous analytical work on selfsimilar steadystate grain distributions. Numerical simulations are carried out to identify under which conditions the grain size distributions can be understood in terms of a combination of powerlaw distributions. A physically motivated fitting formula for grain size distributions is derived using our analytical predictions and numerical simulations.
Results. We find good agreement between analytical results and numerical solutions of the Smoluchowski equation for simple shapes of the kernel function. The results for more complicated and realistic cases can be fitted with a physically motivated “black box” recipe presented in this paper. Our results show that the shape of the dust distribution is mostly dominated by the gas surface density (not the dusttogas ratio), the turbulence strength and the temperature and does not obey an MRN type distribution.
Key words: accretion, accretion disks / stars: premainsequence / circumstellar matter / planets and satellites: formation / protoplanetary disks
© ESO, 2010
1. Introduction
Dust size distributions are a fundamental ingredient for many astrophysical models in the context of circumstellar disks and planet formation: whenever dust is present, it dominates the opacity of the disk, thereby influencing the temperature and consequently also the vertical structure of the disk (e.g., Dullemond & Dominik 2004). Small grains effectively sweep up electrons and therefore strongly affect the chemistry and the ionization fraction (also via grain surface reactions, see Vasyunin et al., in prep.) and thereby also the angular momentum transfer of the disk (e.g., Wardle & Ng 1999; Sano et al. 2000).
Today, it is well established that the dust distributions in asteroid belts and debris disks are governed by a socalled “collision cascade” (see Williams & Wetherill 1994): larger bodies in a gas free environment exhibit high velocity collisions ( ≳ km s^{1}), far beyond their critical fragmentation threshold, which lead to cratering or even complete shattering of these objects. The resulting fragments in turn suffer the same fate, thus producing ever smaller grains down to sizes below about a few micrometers where PoyntingRobertson drag removes the dust particles (e.g., Wyatt et al. 1999). The grain number density distribution in the case of such a fragmentation cascade has been derived by Dohnanyi (1969) and Williams & Wetherill (1994) and was found to follow a powerlaw number density distribution n(m) ∝ m^{−α} with index (which is equivalent to n(a) ∝ a^{3.5}), with very weak dependence on the mechanical parameters of the fragmentation process. Tanaka et al. (1996), Makino et al. (1998) and Kobayashi & Tanaka (2010) showed that this result is exactly independent of the adopted collision model and that the resulting slope α is only determined by the massdependence of the collisional crosssection if the model of collisional outcome is selfsimilar (in the context of fluid dynamics, the same result was independently obtained by Hunt 1982; and Pushkin & Aref 2002). The value of agrees well with the size distributions of asteroids (see Dohnanyi 1969) and of grains in the interstellar medium (MRN distribution, see Mathis et al. 1977; Pollack et al. 1985) and is thus widely applied, even at the gasrich stage of circumstellar disks.
However, in protoplanetary disks gas drag damps the motions of particles. Very small particles are tied to the gas and, as a result, have a relative velocity low enough to make sticking feasible. The size distribution therefore deviates from the MRN powerlaw. Theoretical models of grain growth indicate that particles can grow to sizes much larger than a few μm (see Nakagawa et al. 1981; Weidenschilling 1980, 1984, 1997; Dullemond & Dominik 2005; Tanaka et al. 2005; Brauer et al. 2008a; Birnstiel et al. 2009). Indeed, the observational evidence suggests that growth up to cmsizes is possible (e.g., Testi et al. 2003; Natta et al. 2004; Rodmann et al. 2006; Ricci et al. 2010).
However, as particles grow, they become more loosely coupled to the gas. This results in an increase in the relative velocity between the particles, a common feature of most sources of particles’ relative velocity (i.e., turbulence, radial drift, and settling motions). Therefore, we expect that the assumption of perfect sticking will break down at some point and other collisional outcomes (bouncing, erosion, catastrophic disruption) become possible (see Blum & Wurm 2008). It is expected, then, that at a certain point growth will cease for the largest particles in the distribution. Collisions involving these particles result in fragmentation, thus replenishing the small grains. On the other hand mutual collisions among small particles still result in coagulation. As a result, a steadystate emerges. In this paper, this situation is referred to as a fragmentationcoagulation equilibrium.
The situation in protoplanetary disks differs, therefore, from that in debris disks. In the latter only fragmentation operates. The mass distribution still proceeds towards a steady state but, ultimately, mass is removed from the system due to, radiation pressure or PoyntingRobertson drag.
In this paper, we consider the situation that the total mass budget in the system is conserved. For simplicity, we will ignore motions due to radial drift in this study. This mechanism effectively removes dust particles from the disk as well as providing particles with a large relative motion. However, the derived presence of mmsize dust particles in protoplanetary disks is somewhat at odds with the usual (laminar) prescriptions for the radial drift rate (Weidenschilling 1977). Recently, it was shown that mm observations of protoplanetary disks can be explained by steadystate size distributions if radial drift is inefficient (Birnstiel et al. 2010b). On the other hand, if radial drift would operate as effectively as the laminar theory predicts, then the observed populations of mmsized particles at large disk radii cannot be sustained (Brauer et al. 2007). Possible solutions to reduce the drift rate include bumps in the radial pressure profile (see Kretke & Lin 2007; Brauer et al. 2008b; Cossins et al. 2009) or zonal flows (see Johansen et al. 2009).
In this work we analytically derive steadystate distributions of grains in the presence of both coagulation and fragmentation. The analytical predictions are compared to numerical simulations and applied to grain size distributions in turbulent circumstellar disks. Both the theoretical and numerical results presented in this work are used to derive a fitting formula for steadystate grain size distributions in circumstellar disks.
The paper is outlined as follows: in Sect. 2, we briefly summarize and then generalize previous results by Tanaka et al. (1996) and Makino et al. (1998). In Sect. 3, we test our theoretical predictions and their limitations by a stateoftheart grain evolution code (see Birnstiel et al. 2010a). Grain size distributions in circumstellar disks are discussed in Sect. 4. In Sect. 5, we present a fitting recipe for these distributions that can easily be used in models where grain properties are important. Our findings are summarized in Sect. 6.
2. Powerlaw solutions for a dust coagulationfragmentation equilibrium
In this section, we begin by summarizing some of the previous work on analytical selfsimilar grain size distributions on which our subsequent analysis is based. We will then extend this to include both coagulation and fragmentation processes in a common framework. Under the assumption that the relevant quantities, i.e., the collisional probability between particles, the distribution of fragments, and the size distribution, behave like powerlaws, we will solve for the size distribution in coagulationfragmentation equilibrium. For simplicity, we consider a single monomer size only of mass m_{0}. This is therefore the smallest mass in the distribution.
Another key assumption of our analytical model is that we assume the existence of a sharp threshold mass m_{f}, above which collisions always result in fragmentation and below which collision always result in coagulation. As explained above, the physical motivation for this choice is that relative velocities increase with mass. This also means that in our theoretical model we neglect collision outcomes like bouncing or erosion (erosion is included in the simulations). Even though, small particles will in reality cause cratering/growth instead of complete fragmentation of the larger particle, we find that this assumption is often justified because fragmenting similarsized collisions prevent any growth beyond m_{f}. Thus, collisions with much smaller particles (m ≪ m_{f}) do not have an important influence on the maximum size, however, they can significantly change the amount of small particles in the stationary distribution (see Sect. 5). We further assume a constant porosity of the particles, which relates the mass and size according to (1)where ρ_{s} is the internal density of a dust aggregate.
In our analysis, we will identify three different regimes, which are symbolically illustrated in Fig. 1:

Regime A represents a case where grains grow sequentially(i.e. hierarchically) by collisions with similarsized grains until they reach an upper size limit and fragment backto the smallest sizes. The emerging powerlaw slope of the sizedistribution depends only on the shape of the collisional kernel.

Regime B is similar to regime A, however in this case the fragmented mass is redistributed over a range of sizes and, thus, influencing the outcoming distribution.

In Regime C, the upper end of the distribution dominates grain growth at all sizes. Smaller particles are swept up by the upper end of the distribution and are replenished mostly by redistributed fragments of the largest particles. The resulting distribution function depends strongly on how the fragmented mass is distributed after a disruptive collision.
For each of these regimes, we will derive the parameter ranges for which they apply and the slopes of the resulting grain size distribution.
Fig. 1 Illustration of the three different regimes for which analytical solutions have been derived. Case A represents the growth cascade discussed in Sect. 2.1, case B the intermediate regime (see Sect. 2.2) and case C the fragmentation dominated regime which is discussed in Sect. 2.3. Particles are shattered once they reach the fragmentation barrier m_{f} since collision velocities for particles > m_{f} exceed the fragmentation threshold velocity. 
2.1. The growth cascade
The fundamental quantity that governs the timeevolution of the dust size distribution is the collision kernel C_{m1,m2}. It is defined such that (2)gives the number of collisions per unit time per unit volume between particles in mass interval [m_{1},m_{1} + dm_{1}] and [m_{2},m_{2} + dm_{2}] , where n(m) is the number density distribution. Once specified it determines the collisional evolution of the system. In the case that the number density n(m) does not depend on position, C_{m1,m2} is simply the product of the collision cross section and the relative velocity of two particles with the masses m_{1} and m_{2}.
We use the same Ansatz as Tanaka et al. (1996), assuming that the collision kernel is given in the selfsimilar form (3)Here, h is any function which depends only on the masses through the ratio of m_{2}/m_{1}. By definition, the kernel C_{m1,m2} has to be symmetric, therefore, Eq. (3) implies that h(m_{1},m_{2}) is not symmetric (see Eq. (21)). ν is called the index of the kernel or the degree of homogeneity. As we will see in the following, ν is one of the most important parameters determining the resulting size distribution. Different physical environments are represented by different values of ν. Examples include the constant kernel (i.e., ν = 0, mass independent), the geometrical kernel (i.e., ν = 2/3, velocity independent) or the linear kernel (i.e., ν = 1, as for grains suspended in turbulent gas).
It is further assumed that the number density distribution of dust particles follows a powerlaw, (4)The time evolution of the mass distribution obeys the equation of mass conservation, (5)where the flux F(m) does not represent a flux in a typical continuous way since coagulation is nonlocal in mass space (each mass can interact with each other mass) but rather an integration of all growth processes which produce a particle with mass greater than m out of a particle that was smaller than m (i.e. collisions of m_{1} < m with any other mass m_{2} such that m_{1} + m_{2} > m). The flux in the case of pure coagulation was derived by Tanaka et al. (1996)(6)where we changed the lower bound of the integration over m_{1} to start from the mass of monomers m_{0} (instead of 0) and the upper bound of the integral over m_{2} to a finite upper end m_{f} (instead of infinity), compared to the definition used by Tanaka et al. (1996).
Substituting the definitions of above and using the dimensionless variables x_{1} = m_{1}/m, x_{2} = m_{2}/m, x_{0} = m_{0}/m, and x_{f} = m_{f}/m one obtains (7)where K approaches a constant value in the limit of m ≫ m_{0} and m ≪ m_{f}.
Postulation of a steady state (i.e. setting the time derivative in Eq. (5) to zero), leads to the condition (8)from which it follows that the slope of the distribution is (9)This result was already derived for the case of fragmentation by Tanaka et al. (1996) and Dohnanyi (1969) and for the coagulation by Klett (1975) and Camacho (2001). The physical interpretation of this is a “reversed” fragmentation cascade: instead of a resupply of large particles which produces ever smaller bodies, this represents a constant resupply of monomers which produce ever larger grains (cf. case A in Fig. 1).
2.2. Coagulation fragmentation equilibrium
As Tanaka et al. (1996) and Makino et al. (1998) pointed out, the previous result is independent of the model of collisional outcomes as long as this model is selfsimilar (Eq. (3)). However this is no longer the case if we consider both coagulation and fragmentation processes happening at the same time.
We will now consider the case with a constant resupply of matter, due to the particles that fragment above m_{f}. We assume these fragments obey a powerlaw mass distribution and are produced at a rate (10)where ξ reflects the shape of the fragment distribution and N is a constant.
If there is a constant flux of particles F(m) as defined above, then the flux of fragmenting particles (i.e., the flux produced by particles that are growing over the fragmentation threshold) is given by (11)where K is the integral defined in Eq. (7).
The resulting (downward) flux of fragments F_{f}(m) can then be derived by inserting Eq. (10) into the equation of mass conservation (Eq. (5)), (12)Integration from monomer size m_{0} to m yields (13)The normalization factor N can be determined from the equilibrium condition that the net flux vanishes, (14)and was found to be (15)In Eq. (13), we can distinguish two cases:

If ξ > 2, the contribution of m_{0} dominates the term in brackets. This means that most of the fragment mass is redistributed to monomer sizes and the situation is the same as in the pure coagulation case (cf. case A in Fig. 1). The steadystate condition F(m) + F_{f}(m) = 0 (i.e., the net flux is zero) yields that (16)is constant, which leads to the same result as Eq. (9). Intuitively, this is clear since the majority of the redistributed mass ends up at m ~ m_{0}.

If ξ < 2, the mdependence dominates the term in brackets in Eq. (13) and postulation of a steadystate, (17)leads to an exponent (18)less than Eq. (9). In this case, the slope of the fragment distribution matters. This scenario is represented as case B in Fig. 1.
2.3. Fragment dominated regime
The result obtained in the previous section may seem to be quite general. However, it does not hold for low ξvalues as we will show in the following.
In our case, the integrals do not diverge due to the finite integration bounds. However Makino et al. (1998) used 0 and ∞ as lower and upper bounds for the integration and thus needed to investigate the convergence of the integral. They derived the following conditions for convergence: where γ gives the dependence of m_{2}/m_{1} in the hfunction of Eq. (21) (see Makino et al. 1998): (21)The first condition (Eq. (19)) considers the divergence towards the upper masses, whereas the second condition pertains the lower masses. We assume that Eq. (20) is satisfied and consider the case of decreasing ξ for α given by Eq. (18) which results in a steeper size distribution (where the mass is concentrated close to the upper end of the distribution). We see that for ξ < 1 + ν − 2γ, Eq. (19) is no longer fulfilled. The behavior of the flux integral changes qualitatively: growth is no longer hierarchical, but it becomes dominated by contributions by the upper end of the integration bounds.
Physically, this means that the total number of collisions of any grain is determined by the largest grains in the upper end of the distribution. Hence, smaller sized particles are predominantly refilled by fragmentation events of larger bodies instead of coagulation events from smaller bodies and they are predominantly removed by coagulation events with big particles (near the threshold m_{f}) instead of similarsized particles. This corresponds to case C in Fig. 1.
To determine the resulting powerlaw distribution, we again focus on Eq. (6). The double integral of the flux F(m) can now be split into three separate integrals, (22)according to whether m_{2} is larger or smaller than m_{1}.
It can be derived (see Appendix A) that if the condition above (Eq. (19)) is violated and if m_{f} ≫ m, then the first and the third integral in Eq. (22) dominate the flux due to the integration until m_{f}. In this cases, the flux F(m) is proportional to m^{2 + γ−α}.
A stationary state in the presence of fragmentation (Eq. (13)) is reached if the fluxes cancel out, which leads to the condition (23)This is the sweepup regime where small particles are cleaned out by big ones (cf. case C in Fig. 1).
2.4. Summary of the regimes
Summarizing these findings, we find that the resulting distribution is described by three scenarios (depicted in Fig. 1), depending on the slope of the fragment distribution: (24)
3. Simulation results interpreted
In this section, we will test the analytical predictions of the previous section by a coagulation/fragmentation code (Birnstiel et al. 2010a; see also Brauer et al. 2008a). The code solves for the time evolution of the grain size distribution using an implicit integration scheme. This enables us to find the steadystate distribution by using large time steps. In this way, the time evolution is not resolved, but the steadystate distribution is reliably and very quickly derived.
We start out with the simplest case of a constant kernel and then – step by step – approach a more realistic scenario (in the context of a protoplanetary disk). In Sect. 4, we will then consider a kernel taking into account relative velocities of Brownian motion and turbulent velocities and also a fragmentation probability which depends on the masses and the relative velocities of the colliding particles.
The following results are only steadystate solutions, whether or not this state is reached depends on several conditions. Firstly, particles need to fragment. If there is no upper boundary for growth, it will proceed unlimited and a steady state will never be reached. Secondly, radial motion needs to be slow enough to allow for a steady state. If quantities like the surface density or the temperature vary smoothly between neighbouring regions in the disk, the steady state solutions will also be similar and radial transport will happen on top of a steadystate grain distribution. However if radial drift is acting strongly on particles of the distribution, the steady state will not be reached. Thirdly, a distribution of initially subμm sized grains will need some time to get into an equilibrium state. This time scale can be as small as <1000 years inside of a few AU, while it can be of the order of a million years at 100 AU. We provide a rough estimate of this time scale in Appendix C.
3.1. Constant kernel
In the following section, we consider the case of a constant kernel and will include fragmentation above particle sizes of 1 mm because this represents an instructive test case.
Fig. 2 Grain size distributions^{2} for a constant kernel (i.e., ν = 0) and different distributions of fragments. The peak towards the upper end of the distribution is due to the fragmentation barrier (explanation in the text). The slope of the mass distribution corresponds to 6 − 3α. 
We iteratively solve for the steadystate size distribution between coagulation and fragmentation. The outcome of these simulations for a constant kernel (i.e., ν = 0) are powerlaw distributions where the slope of the distribution depends on the fragmentation law (the slope ξ, see Eq. (10)). Figure 2 shows the corresponding size distributions for some of the different fragment distributions: the steepest distribution corresponds to the case of ξ = 0.5. For larger ξvalues, the slope of the mass distribution flattens. In all cases, a bump develops towards the upper end of the distributions. The reason for this “pileup” is the following: grains typically grow mostly through collisions with similarsized or larger particles. Since the distribution is truncated at the upper end (defined as a_{max}, see also Eq. (43)), particles close to the upper end lack larger collision partners, the growth rate at these sizes would decrease if the distribution keeps its powerlaw nature. This, in turn means that the flux could not be constant below a_{max}. To keep a steady state, the number of particles at that point has to increase in order to replace the missing collision partners at larger sizes.
Fig. 3 Exponent of grain size distributions for a constant kernel (i.e. ν = 0) and different distributions of fragments (solid line) and the analytic solution for the growth cascade (case A, dotted line), the intermediate regime (case B, dashed line), and the fragment dominated regime (case C, dashdotted line). 
Figure 3 shows how the slope of the resulting distribution depends on the fragmentation slope ξ, where the three previously discussed regimes can be identified:

Case A, growth cascade: as predicted, this scenario holds forvalues of ξ ≳ 2 where most of the fragmenting mass is redistributed to fragments. This case corresponds to a “reversed” collisional cascade (just the direction in mass space is reversed as collisions mostly lead to growth instead of shattering).

Case B, the intermediate regime: when the fragment mass is more or less equally distributed over all sizes, mass gain by redistribution of fragments and the mass loss due to coagulation have to cancel each other.

Case C: most of the fragmented mass remains at large sizes. Therefore the mass distribution is dominated by the largest particles. In other words, growth is not hierarchical anymore. In this test case γ = 0, and from Eq. (19) it follows that the transition between the intermediate and fragmentdominated regime lies at ξ = 1, which can be seen in Fig. 3.
The measured slopes of the size distribution for m ≪ m_{f} are in excellent agreement with the model outlined in Sect. 2.
3.2. Including settling effects
The simplest addition to this model is settling: as grains become larger, they start to settle towards the midplane. However, turbulent mixing counteracts this systematic motion. The vertical distribution of dust in a settlingmixing equilibrium can (close to the midplane) be estimated by a Gaussian distribution with a sizedependent dust scale height H_{d}. Smaller particles are well enough coupled to the gas to have the same scale height as the gas, H_{g}, while larger particles decouple and their scale height is decreasing with grain size, (25)(see e.g., Dubrulle et al. 1995; Schräpler & Henning 2004; Youdin & Lithwick 2007) where the Stokes number St is a dimensionless quantity which describes the dynamic properties of a suspended particle. Very small particles have a small Stokes number and are therefore well coupled to the gas. Particles which have different properties (e.g., size or porosity) but the same St behave aerodynamically the same. In our prescription of turbulence, St is defined as the product of the particles stopping time τ_{st} and the orbital frequency Ω_{k}. We focus on the Epstein regime where the Stokes number can be approximated by (26)with Σ_{g} being the gas surface density and ρ_{s} being the internal density of the particles which relates mass and size via m = 4π/3ρ_{s} a^{3},
Fig. 4 Simulation result (solid line) and a “two powerlaw fit” to the vertically integrated dust distribution for a constant kernel with settling included. The dotted, vertical line denotes the grain size above which grains are affected by settling. 
Settling starts to play a role as soon as the Stokes number becomes larger than the turbulence parameter α_{t} which can be related to a certain size (27)Equation (27) only holds within about one gas pressure scale height because the dust vertical structure higher up in the disk deviates from the a Gaussian profile (see also Dubrulle et al. 1995; Schräpler & Henning 2004; Dullemond & Dominik 2004).
The massdependent dust scale height causes the number density distribution n(m) to depend on the vertical height, z. In disklike configurations, it is therefore customary to consider the column density, (28)Similar to Eq. (2), we can write the vertically integrated number of collisions as (29)which gives the total number of collisions that take place over the entire column of the disks.
The dependence of the collisional probability on scale height, is now reflected in the modified kernel (see Birnstiel et al. 2010a, Appendix A for derivation): (30)The point to realize here is that, due to the symmetry between Eqs. (2) and (29), the analysis in Sect. 2 holds also for disklike configuration, if the kernel is now replaced by . The resulting exponent α then concerns the column density dependence (N(m) ∝ m^{ − α}).
If we consider the case that St > α_{t} and substitute Eqs. (25) and (26) into (30), we find that (31)has an index ν = 1/6 for grain sizes larger than a_{sett} and ν = 0 otherwise (it should be noted that H_{1} and H_{2} are the dust scale heights whereas h(m_{2}/m_{1}) represents the function defined in Eq. (3)).
The theory described in Sect. 2 is strictly speaking only valid for a constant νindex, but if this index is constant over a significant range of masses, then the local slope of the distribution will still adapt to this index. In the case of settling, we can therefore describe the distribution with two powerlaws as can be seen in Fig. 4.
The fact that the distribution will locally follow a powerlaw is an important requirement for being able to construct fitting formulas which reproduce the simulated grain size distributions. It allows us in some cases to explain the simulation outcomes with the local kernel index (although coagulation and fragmentation are nonlocal processes in mass space, since each mass may interact with each other mass). A physically motivated recipe to fit the numerically derived distribution functions for the special case of ξ = 11/6 is presented in Sect. 5.
3.3. Nonconstant kernels
We performed the same tests as in Sect. 3.1 also for nonconstant kernels, i.e. kernels with ν > 0. The measured slopes for the cases of (ν = 5/6,γ = 0) (corresponding to the second turbulent regime, see Sect. 4.1), (ν = 1/3,γ = 1/3), and (ν = 1/6,γ = −1/2) (i.e., Brownian motion, see Sect. 4.1) are shown in Fig. 5. Similar to Fig. 3, the distribution always follows the minimum index of the three different regimes (cases A − C). It should be noted that for all indices ξ of the fragmentation law, all processes – coagulation, fragmentation, and redistribution of fragments – take place, but the relative importance of them is what determines the resulting slope.
In Fig. 5, it can be seen that case B (defined in Sect. 2.2) vanishes for the kernel in the upper panel, while it is present for a large range of ξ for the kernel in the middle panel. This can be explained by the definitions of the three regimes which were summarized in Eq. (24): with (ν = 5/6,γ = 0) (cf. upper panel in Fig. 5), case B is confined between ξ = 11/6 and 2. The grain size distribution, therefore, switches almost immediately from being growth dominated (case A) to fragmentation dominated (case C). In the case of a kernel with ν = 1/3 and γ = 1/3, this range spans from 2/3 to 2, as can be seen in the central panel of Fig. 5.
The lower panel shows the distribution for a Brownian motion kernel (i.e., ν = 1/6 and γ = −1/2). The grey shaded area highlights the range of ξ values where our predictions do not apply, that is, Eq. (20) is no longer fulfilled and thus, the resulting steadystate distributions are no longer powerlaw distributions.
Fig. 5 Exponents of the grain size distributions for three different kernels as function of the fragment distribution index. The plot shows the simulation results (solid lines), the analytic solution for the growth cascade (case A, dotted lines), the intermediate regime (case B, dashed lines), and the fragment dominated regime (case C, dashdotted lines). The upper panel was calculated with a (ν = 5/6,γ = 0)kernel (i.e. turbulent relative velocities, see Sect. 4.1), the middle panel with a (ν = 1/3,γ = 1/3)kernel and the lower panel with a (ν = 1/6,γ = − 1/2)kernel (i.e. Brownian motion relative velocities). In the grey shaded area, Eq. (20) is not fulfilled, and the size distribution is not a powerlaw. 
4. Grain size distributions in circumstellar disks
In this section, we will leave the previous “clean” kernels and focus on the grain size distribution in circumstellar disks including relative velocities due to Brownian and turbulent motion of the particles, settling effects, and a fragmentation probability as function of particle mass and impact velocity.
Combining all these effects makes it impossible to find simple analytical solutions in the case of a coagulationfragmentation equilibrium. We will therefore use a coagulation fragmentation code to find the steadystate solutions and to show how the steadystate distributions in circumstellar disks depend on our input parameters.
We will first discuss the model “ingredients”, i.e. the relative velocities and the prescription for fragmentation/sticking. Afterwards, we will define characteristic sizes at which the shape of the size distribution changes due to the underlying physics. In the last subsection, we will then show the simulation results and discuss the influence of different parameters.
4.1. Relative velocities
We will now include the effects of relative velocities due to Brownian motion, and due to turbulent mixing. The Brownian motion relative velocities are given by (32)where k_{b} is the Boltzmann constant and T the midplane temperature of the disk. Ormel & Cuzzi (2007) have derived closed form expressions for particle collision velocities induced by turbulence. They also provided easytouse approximations for the different particle size regimes which we will use in the following. Small particles (i.e. stopping time of particle ≪ eddy crossing time) belong to the first regime of Ormel & Cuzzi (2007) with velocities proportional to (33)The relative velocities of particles in the second turbulent regime of Ormel & Cuzzi (2007) are given by (34)where St_{max} is the larger of the particles Stokes numbers. Velocities in this regime show also a weak dependence on the ratio of the Stokes numbers which we will neglect in the following discussion.
Together with the geometrical cross section σ_{geo} = π(a_{1} + a_{2})^{2}, it is straightforward to estimate the indices of the kernel, ν and γ, as defined in Eqs. (3) and (21) for all these regimes, without settling (assuming that only Brownian motion or turbulent motion dominates the relative velocities). The indices for these three sources of relative particle motion are summarized in Table 1 (for a derivation, see Appendix B). If settling is to be included, then the ν index for particle sizes above a_{sett} has to be increased by (see Sect. 3.2) while γ remains the same. The ν and γ indices for all three regimes can be found in Table 1.
Kernel indices for the different regimes without settling.
4.2. Fragmentation and cratering
We introduce fragmentation and cratering according to the recipe of Birnstiel et al. (2010a): whether a collision leads to sticking or to fragmentation/cratering is determined by the fragmentation probability (35)which means that all impacts with velocities above the critical breakup threshold u_{f} lead to fragmentation or cratering while impacts with velocities below u_{f} − δu lead to sticking. The width of the linear transition region δu is chosen to be 0.2 u_{f}, since laboratory experiments suggest that there is no sharp fragmentation velocity (Blum & Muench 1993; Güttler 2010, priv. comm.). Our simulation results do not show a strong dependence on this value, which was tested by varying δu by a factor of 2.
The mass ratio of the two particles determines whether the impact completely fragments the larger body (masses within one order of magnitude) or if the smaller particle excavates mass from the larger body (masses differ by more than one order of magnitude). This distinction between an erosion and a shattering regime follows the numerical studies of Paszun & Dominik (2009) and Ormel et al. (2009) and experimental studies of Güttler et al. (2010). These works do not precisely constrain the mass ratio which distinguishes between both regimes, however our simulation results do only weakly depend on it.
In the case of fragmentation, the whole mass of both collision partners is redistributed to all masses smaller than the larger body according to Eq. (10).
In the case of cratering, it is assumed that the smaller particle (with mass m_{imp}) excavates its own mass from the larger body. The mass of the impactor as well as the excavated mass is then redistributed to masses smaller than the impactor mass according to Eq. (10). Thus, the total redistributed mass equals (36)and the mass of the larger body is reduced by 2 m_{imp}.
Most parameters such as the fragmentation velocity or the amount of excavated material during cratering are not yet well enough constrained (for the most recent experimental work, see Blum & Wurm 2008; Güttler et al. 2010, and references therein). Experiments suggest fragmentation velocities of a few m s^{1} and fragment distributions with ξ between 1 and 2 (Güttler et al. 2010 find ξ values between 1.07 and 1.37 for SiO_{2} grains). Simulations of silicate grain growth around 1 AU show that also bouncing (i.e. collisions without sticking or fragmentation) can play an important role (see Weidling et al. 2009; Zsom et al. 2010). However, changes in material composition such as organic or ice mantles or the monomer size are expected to change this picture. As there is still a large parameter space to be explored, we continue with the rather simple recipe of sticking, fragmentation and cratering outlined above.
4.3. Regime boundaries
From Eq. (24), we can calculate the slope of the distribution in the different regimes if we assume that the slope of the distribution at a given grain size always follows from the kernel index ν. To construct a whole distribution consisting of several powerlaws for each regime, we need to know where each of the different relative velocity regime applies.
It is important to note that relative velocities due to Brownian motion decrease with particle size whereas the relative velocities induced by turbulent motion increase with particle size (up to St = 1). Therefore, Brownian motion dominates the relative velocities for small particles, while for larger particles, turbulence dominates. From numerical simulations, we found that at those sizes where the highest turbulent relative velocities (i.e. collisions with the smallest grains) start to exceed the smallest Brownian motion relative velocities (i.e. collisions with similar sized grains), the slope of the distribution starts to be determined by the turbulent kernel slope. By equating the approximate relative velocity of Ormel & Cuzzi (2007) and Eq. (32), the according grain size can be estimated to be (37)where we approximate the Reynolds number (i.e., the ratio of turbulent viscosity ν_{t} = αc_{s} H_{p} over molecular viscosity) near the disk midplane by (38)Here, σ_{H2} is the cross section of molecular hydrogen (taken to be 2 × 10^{15} cm^{2}) and μ = 2.3 is the mean molecular weight in proton masses m_{p}.
Turbulent relative velocities strongly increase for grains with a stopping time that is larger or about the turnover time of the smallest eddies. More specifically, the Stokes number of the particles at this change in the relative velocity is (39)where t_{L} = 1/Ω_{k} and are the turnover times of the largest and the smallest eddies, respectively, and Ω_{k} is the Kepler frequency. Ormel & Cuzzi (2007) approximated the factor y_{a} to be about 1.6. The corresponding grain size in the Epstein regime is therefore given by (40)As mentioned above, the Brownian motion relative velocities of small grains decrease with their size. For larger sizes, the relative velocities due to turbulent motion are gaining importance, which are increasing with size until a Stokes number of unity. For typical values of the sound speed (41)and the turbulence parameter α_{t}, the largest turbulent relative velocity exceed the critical collision velocity of the grains (which is of the order of a few m s^{1}) and therefore leads to fragmentation of the dust particles. In the case of very quiescent environments and/or larger critical collision velocities, particles do not experience this fragmentation barrier and can continue to grow. Hence, a steady state is never reached. The work presented here focuses on the former case where (42)and grain growth is always limited by fragmentation.
Fig. 6 Fiducial model and variations of the most important parameters: fragmentation powerlaw index ξ, threshold fragmentation velocity u_{f}, turbulence parameters α_{t}, surface density Σ_{g}, particle internal density ρ_{s}, and midplane temperature T. The shape of the vertically integrated size distributions does not depend on the stellar mass or the distance to the star (only via the radial dependence of the parameters above). 
As relative turbulent velocities are (in our case) increasing with grain size, we can relate the maximum turbulent relative velocity and the critical collision velocity to derive the approximate maximum grain size which particles can reach (see Birnstiel et al. 2009) (43)
4.4. Resulting steadystate distributions
The parameter space is too large to even nearly discuss all possible outcomes of steady state grain size distributions. We will therefore focus on a few examples and rather explain the basic features and the most general results only. For this purpose, we will adopt a fiducial model and consider the influences of several parameters on the resulting grain size distribution: ξ, u_{f}, α_{t}, Σ_{g}, ρ_{s}, and T (see Fig. 6).
The fiducial model (see the solid black line in Fig. 6) shows the following features: steep increase from the smaller sizes until a few tenth of a micrometer. This relates to the regime dominated by Brownian motion relative velocities. The upper end of this regime can be approximated by Eq. (37). The flatter part of the distribution is caused by a different kernel index ν in the parts of the distribution which are dominated by turbulent relative velocities. The dip at about 60 μm (cf. Eq. (40)) is due to the jump in relative velocities as the stopping time of particles above this size exceeds the shortest eddy turnover time (see Ormel & Cuzzi 2007).
The upper end of the distribution is approximately at a_{max}. The increased slope of the distribution and the bump close to the upper end are caused by two processes. Firstly, a boundary effect: grains mostly grow by collisions with similar or larger sized particles. Grains near the upper end of the distribution lack larger sized collision partners and therefore the number density needs to increase in order to keep the flux constant with mass (i.e. in order to keep a steadystate). Secondly, the bump is caused by cratering: impacts of small grains onto the largest grains do not cause growth or complete destruction of the larger bodies, instead they erode them. Growth of these larger bodies is therefore slowed down and, similar to the former case, the mass distribution needs to increase in order to fulfill the steadystate criterion (“pileup effect”).
The upper left panel in Fig. 6 shows the influence of the distribution of fragments after a collision event: larger values of ξ mean that more of the fragmented mass is redistributed to smaller sizes. Consequently, the mass distribution at smaller sizes increases relative to the values of smaller ξ values.
The strong influence of the fragmentation threshold velocity u_{f} can be seen in the upper right panel in Fig. 6: according to Eq. (43), an order of magnitude higher u_{f} leads to a 100 times larger maximum grain size.
The grain size distributions for different levels of turbulence are shown in the middle left panel of Fig. 6. The effects are twofold: firstly, an increased α_{t} leads to increased turbulent relative velocities, thus, moving the fragmentation barrier a_{max} to smaller sizes (cf. Eq. (43)). Secondly, a larger α_{t} shifts the transition within the turbulent regime, a_{12}, to smaller sizes. Consequently, the second turbulent regime gains importance as α_{t} is increased since its upper end lower boundary extend ever further.
The middle right panel in Fig. 6 displays the influence of an decreased gas surface density Σ_{g} (assuming a fixed dusttogas ratio). It can be seen that not only the total mass is decreased due to the fixed dusttogas ratio but also the upper size of the distribution a_{max} decreases. This is due to the coupling of the dust to the gas: with larger gas surface density, the dust is better coupled to the gas. This is described by a decreased Stokes number (see Eq. (26)) which, in turn, leads to smaller relative velocities and hence a larger a_{max}. Interestingly, the shape of the grain size distribution does not depend on the total dust mass, but on the total gas mass, as long as gas is dynamically dominating (i.e., Σ_{g} ≫ Σ_{d}). If more dust were to be present, grains would collide more often, thus, a steadystate would be reached faster and the new size distribution would be a scaledup version of the former one. However the velocities at which grains collide are determined by the properties of the underlying gas disk. In this way, the dust grain size distribution is not only a measure of the dust properties, but also a measure of the gas disk physics like the gas density and the amount of turbulence.
The shape of the size distribution for different grain volume densities ρ_{s} does not change significantly. However most regime boundaries (a_{sett}, a_{12}, and a_{max}) are inversely proportional to ρ_{s} (because of the coupling to the gas, described by the Stokes number). A decrease (increase) in ρ_{s} therefore shifts the whole distribution to larger (smaller) sizes as can be seen in the lower left panel in Fig. 6.
The upper end of the distribution, a_{max}, is inversely proportional to the midplane temperature T (as in the case of the turbulence parameter) whereas the transition between the different turbulent regimes a_{12} does not. Therefore, increasing the temperature decreases a_{max} in the same way as decreasing Σ_{g} does. However a_{12} is not influenced by temperature changes, therefore the shape of the size distribution changes in a different way than in the case of changing Σ_{g} as can be seen by comparing the middle right and the lower right panels of Fig. 6.
5. Fitting formula for steadystate distributions
In this section, we will describe a simple recipe which allows us to construct vertically integrated grain size distributions which fit reasonably well to the simulation results presented in the previous section.
The recipe does not directly depend on the radial distance to the star or on the stellar mass. A radial dependence only enters via radial changes of the input parameters listed in Table 4. This recipe has been tested for a large grid of parameter values^{3} (shown in Table 2), however there are some restrictions.
5.1. Limitations
These fits strictly apply only for the case of ξ = 11/6. In this case, the slopes of the distribution agree well with the predictions of the intermediate regime (case B, defined in Sect. 2.2). For smaller values of ξ, the slopes do not strictly follow the analytical predictions. This is due to the fact that we include cratering, which is not covered by our theory. Erosion is therefore an important mode of fragmentation: it dominates over complete disruption through the high number of small particles (see also Kobayashi & Tanaka 2010) and it is able to redistribute significant amounts of mass to the smallest particle sizes.
One important restriction for this recipe is the upper size of the particles a_{max}: it needs to obey the condition of Eq. (42), since otherwise, particles will not experience the fragmenting highvelocity impacts and a steady state will never be reached since particles can grow unhindered over the metersize barrier.
There are also restrictions to a very small a_{max}: if a_{max} is close to or even smaller than a_{12}, then the fit will not represent the true simulation outcome very well. In this case, the upper end of the distribution, and in more extreme cases the whole distribution will look much different. Thus, the sizes should obey the condition (44)where each inequality should be within a factor of a few.
5.2. Recipe
The following recipe calculates the vertically integrated mass distribution of dust grains in a turbulent circumstellar disk within the the above mentioned limitations. The recipe should be applied on a logarithmic grain size grid a_{i} with a lower size limit of 0.025 μm and a fine enough size resolution (a_{i + 1}/a_{i} ≲ 1.12). For convenience, all variables are summarized in Table 4. The steps to be performed are as follows:

1.
Calculate the grain sizes which represent the regime boundaries a_{BT}, a_{12}, a_{sett} which are given by Eqs. (37), (40), and (27).

2.
Calculate the turbulent relative velocities for each grain size. For this, we approximate the equations which are given by Ormel & Cuzzi (2007). Collision velocities with monomers are given by (45) where Re is the Reynolds number (see Eq. (38)), and St_{i} and St_{0} are the Stokes numbers of a_{i} and monomers (a_{0} = 0.025 μm), respectively (cf. Eq. (26)). u_{gas} is given by (46)and the interpolation parameter ϵ is defined as (47)and collisions with similar sized bodies are approximated as (48)A comparison between these approximations and the formulas of Ormel & Cuzzi (2007) is shown in Fig. 7.
Fig. 7 Comparison of the turbulent relative velocities of Ormel & Cuzzi (2007) to the fitting formula used in the recipe (see Eqs. (48) and (45)) for grain size distributions. The largest error in the resulting upper grain size a_{P} derived from the fitting formula is about 25%.
Table 3Powerlaw exponents of the distribution n(m)·m·a.
Table 4Definition of the variables used in this paper.
Fig. 8 Stepbystep construction of the fit distribution for the following parameters: Σ_{g} = 20 g cm^{2}, Σ_{d} = 0.2 g cm^{2}, α_{t} = 1 × 10^{4}, u_{f} = 1 m s^{1}, ξ = 1.833, T = 50 K, ρ_{s} = 1.6 g cm^{3}. The upper left panel shows the distribution after step 5: each interval obeys a different powerlaw, the distribution is continuous apart from a jump at a_{12}. The upper right panel displays the fit including an increase at the upper end, according to step 6. The bump caused by cratering (cf. Eq. (53)) is shown in the bottom left panel while the bottom right panel compares the final fit distribution (solid curve) to the simulation result (dashed curve). The vertical lines correspond to the regime boundaries: a_{BT} (solid), a_{sett} (dashdotted), a_{12} (dashed) and a_{P} (dotted).

3.
Using Eqs. (45) and (48) for the relative velocities, and the transition width δu = 0.2u_{f}, find the grain sizes which correspond to the following conditions:

a_{L}: particles above this size experience impacts with monomers with velocities of (i.e., cratering starts to become important).

a_{P}: particles above this size experience impacts with equal sized grains with velocities of (i.e., fragmentation becomes important).

a_{R}: particles above this size experience impacts with similar sized grains with velocities of (i.e. every impact causes fragmentation/cratering).


4.
Calculate the factor J according to the recipe where μ = 2.3, σ_{H2} = 2 × 10^{15} cm^{2} and y_{a} = 1.6.

5.
The powerlaw indices of the mass distribution δ_{i} for each interval between the regime boundaries (a_{BT}, a_{12}, a_{sett}, a_{P}) are calculated according to the intermediate regime (cf. Eq. (18)). The slopes have to be chosen according to the kernel regime (Brownian motion, turbulent regime 1 or 2) and according to whether the regime is influenced by settling or not (i.e., if a is larger or smaller than a_{sett}). The resulting slopes of the mass distribution are given in Table 3. The first version of the fit f(a_{i}) (where a_{i} denotes the numerical grid point of the particles size array) is now constructed by using powerlaws () in between each of the regimes, up to a_{P}. The fit should be continuous at all regime boundaries except a drop of 1/J at the transition at a_{12}. An example of this first version is shown in the top left panel of Fig. 8.

6.
Mimic the cutoff effects which cause an increase in the distribution function close to the upper end: linearly increase the distribution function for all sizes a_{inc} < a < a_{P}: (51)with (52)The resulting fit after this step is shown in the upper right panel of Fig. 8.

7.
The bump due to cratering is mimicked by a Gaussian, (53)where σ is defined as (54)but should be limited to be at least (55)

8.
The fit ℱ(a_{i}) is now constructed by using the maximum of f(a_{i}) and b(a_{i}) in the following way: (56)

9.
Finally, the fit needs to be normalized to the dust surface density at the given radius. ℱ is a (yet unnormalized) vertically integrated mass distribution (shown in the bottom right panel of Fig. 8). To translate this mass distribution to a vertically integrated number density distribution N(a), we need to normalize it as (57)
6. Conclusions
In this work, we generalize the analytical findings of previous works to the case of grain size distributions in a coagulation/fragmentation equilibrium. Under the assumption that all grains above a certain size, a_{max}, fragment into a powerlaw distribution n_{f}(m) ∝ m^{−ξ}, we derived analytical steadystate solutions for selfsimilar kernels and determined three different cases (see Sect. 2.4). Cratering is not covered by our theory. However our simulations that include cratering agree with the theoretical predictions for a fragmentation law with ξ = 11/6.
Results show that dust size distributions in circumstellar disks do not necessarily follow the often adopted MRN powerlaw distribution of n(a) ∝ a^{3.5} (see Mathis et al. 1977; Dohnanyi 1969; Tanaka et al. 1996; Makino et al. 1998; Garaud 2007) when both coagulation and fragmentation events operate. We performed detailed simulations of grain growth and fragmentation to test the analytical predictions and found very good agreement between the theory and the simulation results.
We applied the theory to the gaseous environments of circumstellar disks. Unlike the models of Garaud (2007), the upper end of the size distribution is typically not limited by the growth time scale but by fragmentation because relative velocities increase with grain size and reach values large enough to fragment grains. The shape of the dust distribution is determined by the gaseous environment (e.g., gas surface density, level of turbulence, temperature and others) since the gas is dynamically dominant as long as the gas surface density significantly exceeds the dust surface density. The total dust mass merely provides the normalization of the distribution and the time scale in which an equilibrium is reached. The results presented in this work show that the physics of growth and fragmentation directly link the upper and the lower end of the dust distribution in circumstellar disks.
A readytouse recipe for deriving vertically integrated dust size distributions in circumstellar disks for a fixed value of ξ = 11/6 is presented in Sect. 5. Although the collision kernel in circumstellar disks is complicated, we found good agreement with our fitting recipe for a fragment distribution with ξ = 11/6. The recipe can readily be used for further modeling such as disk chemistry or radiative transfer calculations.
It should be noted that in this paper we will plot the distributions typically in terms of n(a)·m·a which is proportional to the distribution of mass. The advantage of plotting it this way instead of plotting n(a) is the following: when n(m) follows a powerlaw m^{−α}, then the grain size distribution n(a) describes a powerlaw with exponent 2−3α and the mass distribution exponent is 6−3α. For typical values of α, this is less steep and differences between a predicted and a real distribution are more prominent.
See also http://www.mpia.de/distributionfits
Acknowledgments
We like to thank Jürgen Blum and Carsten Güttler for useful discussions and Antonella Natta, Luca Ricci, Francesco Trotta, Taku Takeuchi, Gijs Mulders, Sean Andrews and the anonymous referee for helpful comments.
References
 Birnstiel, T., Dullemond, C. P., & Brauer, F. 2009, A&A, 503, L5 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Birnstiel, T., Dullemond, C. P., & Brauer, F. 2010a, A&A, 513, A79 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Birnstiel, T., Ricci, L., Trotta, F., et al. 2010b, A&A, 516, L14 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Blum, J., & Muench, M. 1993, Icarus, 106, 151 [NASA ADS] [CrossRef] [Google Scholar]
 Blum, J., & Wurm, G. 2008, ARA&A, 46, 21 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Brauer, F., Dullemond, C. P., Johansen, A., et al. 2007, A&A, 469, 1169 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Brauer, F., Dullemond, C. P., & Henning, T. 2008a, A&A, 480, 859 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Brauer, F., Henning, T., & Dullemond, C. P. 2008b, A&A, 487, L1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Camacho, J. 2001, Phys. Rev. E, 63, 46112 [NASA ADS] [CrossRef] [Google Scholar]
 Cossins, P., Lodato, G., & Clarke, C. J. 2009, MNRAS, 393, 1157 [Google Scholar]
 Dohnanyi, J. S. 1969, J. Geophys. Res., 74, 2531 [NASA ADS] [CrossRef] [Google Scholar]
 Dubrulle, B., Morfill, G. E., & Sterzik, M. F. 1995, Icarus, 114, 237 [NASA ADS] [CrossRef] [Google Scholar]
 Dullemond, C. P., & Dominik, C. 2004, A&A, 421, 1075 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Dullemond, C. P., & Dominik, C. 2005, A&A, 434, 971 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Garaud, P. 2007, ApJ, 671, 2091 [NASA ADS] [CrossRef] [Google Scholar]
 Güttler, C., Blum, J., Zsom, A., Ormel, C. W., & Dullemond, C. P. 2010, A&A, 513, A56 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hunt, J. R. 1982, J. Fluid Mech., 122, 169 [NASA ADS] [CrossRef] [Google Scholar]
 Johansen, A., Youdin, A., & Klahr, H. 2009, ApJ, 697, 1269 [NASA ADS] [CrossRef] [Google Scholar]
 Klett, J. D. 1975, J. Atmos. Sci., 32, 380 [NASA ADS] [CrossRef] [Google Scholar]
 Kobayashi, H., & Tanaka, H. 2010, Icarus, 206, 735 [NASA ADS] [CrossRef] [Google Scholar]
 Kretke, K. A., & Lin, D. N. C. 2007, ApJ, 664, L55 [NASA ADS] [CrossRef] [Google Scholar]
 Makino, J., Fukushige, T., Funato, Y., & Kokubo, E. 1998, New Astron., 3, 411 [NASA ADS] [CrossRef] [Google Scholar]
 Mathis, J. S., Rumpl, W., & Nordsieck, K. H. 1977, ApJ, 217, 425 [NASA ADS] [CrossRef] [Google Scholar]
 Nakagawa, Y., Nakazawa, K., & Hayashi, C. 1981, Icarus, 45, 517 [NASA ADS] [CrossRef] [Google Scholar]
 Natta, A., Testi, L., Neri, R., Shepherd, D. S., & Wilner, D. J. 2004, A&A, 416, 179 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ormel, C. W., & Cuzzi, J. N. 2007, A&A, 466, 413 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ormel, C. W., Paszun, D., Dominik, C., & Tielens, A. G. G. M. 2009, A&A, 502, 845 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Paszun, D., & Dominik, C. 2009, A&A, 507, 1023 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Pollack, J. B., McKay, C. P., & Christofferson, B. M. 1985, Icarus, 64, 471 [NASA ADS] [CrossRef] [Google Scholar]
 Pushkin, D., & Aref, H. 2002, Phys. Fluids, 14, 694 [NASA ADS] [CrossRef] [Google Scholar]
 Ricci, L., Testi, L., Natta, A., et al. 2010, A&A, 512, A15 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Rodmann, J., Henning, T., Chandler, C. J., Mundy, L. G., & Wilner, D. J. 2006, A&A, 446, 211 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Sano, T., Miyama, S. M., Umebayashi, T., & Nakano, T. 2000, ApJ, 543, 486 [NASA ADS] [CrossRef] [Google Scholar]
 Schräpler, R., & Henning, T. 2004, ApJ, 614, 960 [NASA ADS] [CrossRef] [Google Scholar]
 Tanaka, H., Inaba, S., & Nakazawa, K. 1996, Icarus, 123, 450 [Google Scholar]
 Tanaka, H., Himeno, Y., & Ida, S. 2005, ApJ, 625, 414 [NASA ADS] [CrossRef] [Google Scholar]
 Testi, L., Natta, A., Shepherd, D. S., & Wilner, D. J. 2003, A&A, 403, 323 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Wardle, M., & Ng, C. 1999, MNRAS, 303, 239 [NASA ADS] [CrossRef] [Google Scholar]
 Weidenschilling, S. J. 1977, MNRAS, 180, 57 [NASA ADS] [CrossRef] [Google Scholar]
 Weidenschilling, S. J. 1980, Icarus, 44, 172 [NASA ADS] [CrossRef] [Google Scholar]
 Weidenschilling, S. J. 1984, Icarus, 60, 553 [NASA ADS] [CrossRef] [Google Scholar]
 Weidenschilling, S. J. 1997, Icarus, 127, 290 [NASA ADS] [CrossRef] [Google Scholar]
 Weidling, R., Güttler, C., Blum, J., & Brauer, F. 2009, ApJ, 696, 2036 [NASA ADS] [CrossRef] [Google Scholar]
 Williams, D. R., & Wetherill, G. W. 1994, Icarus, 107, 117 [NASA ADS] [CrossRef] [Google Scholar]
 Wyatt, M. C., Dermott, S. F., Telesco, C. M., et al. 1999, ApJ, 527, 918 [NASA ADS] [CrossRef] [Google Scholar]
 Youdin, A. N., & Lithwick, Y. 2007, Icarus, 192, 588 [NASA ADS] [CrossRef] [Google Scholar]
 Zsom, A., Ormel, C. W., Güttler, C., Blum, J., & Dullemond, C. P. 2010, A&A, 513, A57 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Appendix A: Derivation of the fragment dominated size distribution
In this section, we will derive the slope of the size distribution which is dominated by the largest particles. Since in our scenario, the integrals are confined between the monomer mass m_{0} and the largest particles at the fragmentation barrier m_{f}, the integrals do not diverge as in the scenario of Tanaka et al. (1996) and Makino et al. (1998), who consider integration bounds of zero and infinity. However, if the first condition of Makino et al. (1998) is not fulfilled, then the mass flux (cf. Eq. (6)) is dominated by the upper bound of the integral. This is the case which we will consider in the following.
As noted in Sect. 2.3, the flux integral (Eq. (6)) can be split into three separate integrals, in order to distinguish cases of m_{2} > m_{1} or m_{2} < m_{1}, (A.1)where I_{1}, I_{2}, and I_{3} correspond (from left to right) to the three integrals defined in Eq. (22). We will now evaluate these integrals in the limits of m_{0} ≪ m ≪ m_{f}, using the limiting behavior of h(m_{2}/m_{1}) as given in Eq. (21).
A.1. First integral
Carrying out the integration over m_{2} in I_{1} leads to (A.2)from which we can derive Eq. (19), the first convergence criterion of Makino et al. (1998).
We consider the case where this condition does not hold, i.e., (A.3)Then, the m_{f} term in brackets dominates over the other term. Thus, the term in brackets is constant and carrying out the integration yields (A.4)Now, if the second condition, Eq. (20) holds, using m_{0} ≪ m, we derive (A.5)
A.2. Second integral
We rewrite I_{2} using the dimensionless variables x_{1} = m_{1}/m and x_{2} = m_{2}/m which yields (A.6)By integrating over x_{2}, we derive (A.7)The term in square brackets can be split into a sum of integrals, the first of which is straightforward to evaluate as (A.8)while the second can be identified as a sum of a Beta function B(a,b) and an incomplete Beta function , (A.9)The conditions from above, assure that the Beta functions, and thus I_{2}, are real. The numerical value of I_{2}/D for a range of values of a and b are shown in Fig. A.1.
Fig. A.1 Integral I_{2}/D as function of a and b. 
A.3. Third integral
Similarly, we can rewrite I_{3} as (A.12)If Eq. (A.10) holds, then from m_{f} ≫ m follows (A.13)
A.4. Deriving the steadystate distribution
The first and the third integrand I_{1} and I_{3} show the same mass dependence and can be summed up to (A.14)while I_{2} is proportional to (A.15)and therefore (A.16)Since the constants of proportionality are factors of order unity and m < m_{f}, the integrals I_{1} + I_{3} are much larger than I_{2}, therefore, the flux F(m) is proportional to m^{γ − α + 2}.
In the case of a steady state, this flux and the downward flux of fragments (which is proportional to m^{2 − ξ}) have to cancel each other. Therefore, the exponents of the mass dependence need to cancel out, i.e. (A.17)which is the slope of the steadystate condition in the fragmentation dominated regime (case C in Fig. 1).
Appendix B: Derivation of the degree of homogeneity
In the following, we will derive the degree of homogeneity (cf. Eq. (3)) as shown in Table 1. For simplicity, we will drop all constant factors and consider only the proportionalities. The Brownian motion kernel is a product of the geometrical cross section σ_{geo} = π (a_{1} + a_{2})^{2} and the Brownian motion relative velocities (cf. Eq. (32)): (B.1)where θ = m_{2}/m_{1}. Thus, the m_{1} dependence of the kernel gives . For θ ≪ 1, it follows that and by comparison to Eq. (21), we can derive .
The relative velocities in the first regime of turbulent relative velocities is proportional to a_{1} − a_{2} (cf. Eq. (33)), thus, the kernel can be written as (B.2)In this case, the function h(θ) is different: h(θ) ∝ 1 for θ ≪ 1, and we derive ν = 1 and γ = 0.
In the second turbulent regime, the relative velocities are proportional to the square root of the Stokes number of the larger particle (see Eq. (34)). By using a limit representation of max(a_{1},a_{2}), we can write (B.3)Thus, for the limit of small θ, we derive and γ = 0. If settling is included (see Eq. (31)), then ν increases by an additional factor of 1/6.
Appendix C: Estimating the equilibration time scale
As shown in the appendix of Birnstiel et al. (2010a), monodisperse growth (i.e. assuming all particles have the same size) provides a good estimate of growth time scales. In this approximation, the growth rate is given by (C.1)where τ_{coll} is the collision time scale. Integration of Eq. (C.1) yields the time a particle needs to grow from size a_{1} to a_{2}. For Brownian motion and the first two turbulent velocity regimes (see Eqs. (32) and (34)), we can derive the growth times Here the last two growth times belong to two distinct cases: t_{IIa} is the time in the second turbulent regime (see Eq. (34)) without settling effects (i.e., a_{2} < a_{sett}) while t_{IIb} includes the fact that settling of particles to the midplane increases the dust density at the midplane (denoted by ρ_{d}) and thus accelerates growth.
Typically, coagulation starts by Brownian motion growth from subμm sized particles to sizes where turbulent velocities become important (i.e., sizes larger than a_{BT}, see Eq. (37)). Therefore, we can estimate this time by (C.5)where a_{0} is the monomer size.
We will neglect the time needed by particles to grow until the second turbulent regime because it is typically much shorter than the other time scales involved. If particles in the second turbulent regime are already influenced by settling (a_{1} > a_{sett}), then growth proceeds according to Eq. (C.4) and the time is given by (C.6)If the size range is entirely below a_{sett}, then the timescale is given by (C.7)For the case that a_{12} ≤ a_{sett} ≤ a_{max}, we need to add both contributions (C.8)By comparison to our time evolving simulations of particle growth and fragmentation, we found that (C.9)estimates within a factor of a few the time at which the distribution reaches a steady state.
All Tables
All Figures
Fig. 1 Illustration of the three different regimes for which analytical solutions have been derived. Case A represents the growth cascade discussed in Sect. 2.1, case B the intermediate regime (see Sect. 2.2) and case C the fragmentation dominated regime which is discussed in Sect. 2.3. Particles are shattered once they reach the fragmentation barrier m_{f} since collision velocities for particles > m_{f} exceed the fragmentation threshold velocity. 

In the text 
Fig. 2 Grain size distributions^{2} for a constant kernel (i.e., ν = 0) and different distributions of fragments. The peak towards the upper end of the distribution is due to the fragmentation barrier (explanation in the text). The slope of the mass distribution corresponds to 6 − 3α. 

In the text 
Fig. 3 Exponent of grain size distributions for a constant kernel (i.e. ν = 0) and different distributions of fragments (solid line) and the analytic solution for the growth cascade (case A, dotted line), the intermediate regime (case B, dashed line), and the fragment dominated regime (case C, dashdotted line). 

In the text 
Fig. 4 Simulation result (solid line) and a “two powerlaw fit” to the vertically integrated dust distribution for a constant kernel with settling included. The dotted, vertical line denotes the grain size above which grains are affected by settling. 

In the text 
Fig. 5 Exponents of the grain size distributions for three different kernels as function of the fragment distribution index. The plot shows the simulation results (solid lines), the analytic solution for the growth cascade (case A, dotted lines), the intermediate regime (case B, dashed lines), and the fragment dominated regime (case C, dashdotted lines). The upper panel was calculated with a (ν = 5/6,γ = 0)kernel (i.e. turbulent relative velocities, see Sect. 4.1), the middle panel with a (ν = 1/3,γ = 1/3)kernel and the lower panel with a (ν = 1/6,γ = − 1/2)kernel (i.e. Brownian motion relative velocities). In the grey shaded area, Eq. (20) is not fulfilled, and the size distribution is not a powerlaw. 

In the text 
Fig. 6 Fiducial model and variations of the most important parameters: fragmentation powerlaw index ξ, threshold fragmentation velocity u_{f}, turbulence parameters α_{t}, surface density Σ_{g}, particle internal density ρ_{s}, and midplane temperature T. The shape of the vertically integrated size distributions does not depend on the stellar mass or the distance to the star (only via the radial dependence of the parameters above). 

In the text 
Fig. 7 Comparison of the turbulent relative velocities of Ormel & Cuzzi (2007) to the fitting formula used in the recipe (see Eqs. (48) and (45)) for grain size distributions. The largest error in the resulting upper grain size a_{P} derived from the fitting formula is about 25%. 

In the text 
Fig. 8 Stepbystep construction of the fit distribution for the following parameters: Σ_{g} = 20 g cm^{2}, Σ_{d} = 0.2 g cm^{2}, α_{t} = 1 × 10^{4}, u_{f} = 1 m s^{1}, ξ = 1.833, T = 50 K, ρ_{s} = 1.6 g cm^{3}. The upper left panel shows the distribution after step 5: each interval obeys a different powerlaw, the distribution is continuous apart from a jump at a_{12}. The upper right panel displays the fit including an increase at the upper end, according to step 6. The bump caused by cratering (cf. Eq. (53)) is shown in the bottom left panel while the bottom right panel compares the final fit distribution (solid curve) to the simulation result (dashed curve). The vertical lines correspond to the regime boundaries: a_{BT} (solid), a_{sett} (dashdotted), a_{12} (dashed) and a_{P} (dotted). 

In the text 
Fig. A.1 Integral I_{2}/D as function of a and b. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.