© 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 so-called “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 Poynting-Robertson 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 power-law 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 mass-dependence of the collisional cross-section if the model of collisional outcome is self-similar (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 gas-rich 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 power-law. 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 cm-sizes 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 steady-state emerges. In this paper, this situation is referred to as a fragmentation-coagulation 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 Poynting-Robertson 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 mm-size 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 steady-state 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 mm-sized 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 steady-state 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 steady-state 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 state-of-the-art 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. Power-law solutions for a dust coagulation-fragmentation equilibrium

In this section, we begin by summarizing some of the previous work on analytical self-similar 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 power-laws, we will solve for the size distribution in coagulation-fragmentation equilibrium. For simplicity, we consider a single monomer size only of mass m₀. 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 similar-sized 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 power-law 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 out-coming 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 mf since collision velocities for particles  > mf exceed the fragmentation threshold velocity.

2.1. The growth cascade

The fundamental quantity that governs the time-evolution 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₁,m₁ + dm₁]  and  [m₂,m₂ + dm₂] , 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₁ and m₂.

We use the same Ansatz as Tanaka et al. (1996), assuming that the collision kernel is given in the self-similar form (3)Here, h is any function which depends only on the masses through the ratio of m₂/m₁. By definition, the kernel C_m1,m2 has to be symmetric, therefore, Eq. (3) implies that h(m₁,m₂) 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 power-law, (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 non-local 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₁ < m with any other mass m₂ such that m₁ + m₂ > 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₁ to start from the mass of monomers m₀ (instead of 0) and the upper bound of the integral over m₂ 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₁ = m₁/m, x₂ = m₂/m, x₀ = m₀/m, and x_f = m_f/m one obtains (7)where K approaches a constant value in the limit of m ≫ m₀ 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 self-similar (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 power-law 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₀ 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₀ 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 steady-state 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₀.

• If ξ < 2, the m-dependence dominates the term in brackets in Eq. (13) and postulation of a steady-state, (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₂/m₁ in the h-function 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 similar-sized particles. This corresponds to case C in Fig. 1.

To determine the resulting power-law 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₂ is larger or smaller than m₁.

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 sweep-up 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 steady-state distribution by using large time steps. In this way, the time evolution is not resolved, but the steady-state 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 steady-state 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 steady-state 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 distributions2 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 steady-state size distribution between coagulation and fragmentation. The outcome of these simulations for a constant kernel (i.e., ν = 0) are power-law 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 “pile-up” is the following: grains typically grow mostly through collisions with similar-sized 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 power-law 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, dash-dotted 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 fragment-dominated 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 mid-plane. However, turbulent mixing counteracts this systematic motion. The vertical distribution of dust in a settling-mixing equilibrium can (close to the mid-plane) be estimated by a Gaussian distribution with a size-dependent 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³,

Fig. 4 Simulation result (solid line) and a “two power-law 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 mass-dependent dust scale height causes the number density distribution n(m) to depend on the vertical height, z. In disk-like 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 disk-like 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₁ and H₂ are the dust scale heights whereas h(m₂/m₁) 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 power-laws as can be seen in Fig. 4.

The fact that the distribution will locally follow a power-law 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 non-local 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. Non-constant kernels

We performed the same tests as in Sect. 3.1 also for non-constant 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 re-distribution 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 steady-state distributions are no longer power-law 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, dash-dotted 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 power-law.

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 coagulation-fragmentation equilibrium. We will therefore use a coagulation fragmentation code to find the steady-state solutions and to show how the steady-state 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 mid-plane temperature of the disk. Ormel & Cuzzi (2007) have derived closed form expressions for particle collision velocities induced by turbulence. They also provided easy-to-use 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₁ + a₂)², it is straight-forward 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.

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 break-up 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₂ 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 power-laws 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 mid-plane by (38)Here, σ_H2 is the cross section of molecular hydrogen (taken to be 2  ×  10^-15 cm²) 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 turn-over 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 turn-over 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 power-law index ξ, threshold fragmentation velocity uf, turbulence parameters αt, surface density Σg, particle internal density ρs, and mid-plane 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 steady-state 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 turn-over 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 steady-state). 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 steady-state criterion (“pile-up 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 two-fold: 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₁₂, 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 dust-to-gas ratio). It can be seen that not only the total mass is decreased due to the fixed dust-to-gas 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 steady-state would be reached faster and the new size distribution would be a scaled-up 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₁₂, 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 mid-plane temperature T (as in the case of the turbulence parameter) whereas the transition between the different turbulent regimes a₁₂ does not. Therefore, increasing the temperature decreases a_max in the same way as decreasing Σ_g does. However a₁₂ 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 steady-state 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³ (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 high-velocity impacts and a steady state will never be reached since particles can grow unhindered over the meter-size barrier.

There are also restrictions to a very small a_max: if a_max is close to or even smaller than a₁₂, 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₁₂, 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₀ are the Stokes numbers of a_i and monomers (a₀ = 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 aP derived from the fitting formula is about 25%.

  Table 3

  Power-law exponents of the distribution n(m)·m·a.

  Table 4

  Definition of the variables used in this paper.

  Fig. 8 Step-by-step construction of the fit distribution for the following parameters: Σg = 20 g cm-2, Σd = 0.2 g cm-2, αt = 1 × 10-4, uf = 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 power-law, the distribution is continuous apart from a jump at a12. 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: aBT (solid), asett (dash-dotted), a12 (dashed) and aP (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² and y_a = 1.6.

• 5.

  The power-law indices of the mass distribution δ_i for each interval between the regime boundaries (a_BT, a₁₂, 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 power-laws () 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₁₂. An example of this first version is shown in the top left panel of Fig. 8.

• 6.

  Mimic the cut-off 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 un-normalized) 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 power-law distribution n_f(m) ∝ m^−ξ, we derived analytical steady-state solutions for self-similar 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 power-law 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 ready-to-use 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.

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

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₀ 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₂ > m₁ or m₂ < m₁, (A.1)where I₁, I₂, and I₃ correspond (from left to right) to the three integrals defined in Eq. (22). We will now evaluate these integrals in the limits of m₀ ≪ m ≪ m_f, using the limiting behavior of h(m₂/m₁) as given in Eq. (21).

A.1. First integral

Carrying out the integration over m₂ in I₁ 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₀ ≪ m, we derive (A.5)

A.2. Second integral

We rewrite I₂ using the dimensionless variables x₁ = m₁/m and x₂ = m₂/m which yields (A.6)By integrating over x₂, we derive (A.7)The term in square brackets can be split into a sum of integrals, the first of which is straight-forward 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₂, are real. The numerical value of I₂/D for a range of values of a and b are shown in Fig. A.1.

Fig. A.1 Integral I2/D as function of a and b.

A.3. Third integral

Similarly, we can rewrite I₃ as (A.12)If Eq. (A.10) holds, then from m_f ≫ m follows (A.13)

A.4. Deriving the steady-state distribution

The first and the third integrand I₁ and I₃ show the same mass dependence and can be summed up to (A.14)while I₂ 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₁ + I₃ are much larger than I₂, 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 steady-state 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₁ + a₂)² and the Brownian motion relative velocities (cf. Eq. (32)): (B.1)where θ = m₂/m₁. Thus, the m₁ 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₁ − a₂| (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₁,a₂), 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₁ to a₂. 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₂ < a_sett) while t_IIb includes the fact that settling of particles to the mid-plane increases the dust density at the mid-plane (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₀ 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₁ > 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₁₂ ≤ 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 mf since collision velocities for particles  > mf exceed the fragmentation threshold velocity.

In the text

	Fig. 2 Grain size distributions2 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, dash-dotted line).
In the text

	Fig. 4 Simulation result (solid line) and a “two power-law 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, dash-dotted 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 power-law.

In the text

Fig. 6 Fiducial model and variations of the most important parameters: fragmentation power-law index ξ, threshold fragmentation velocity uf, turbulence parameters αt, surface density Σg, particle internal density ρs, and mid-plane 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 aP derived from the fitting formula is about 25%.
In the text

Fig. 8 Step-by-step construction of the fit distribution for the following parameters: Σg = 20 g cm-2, Σd = 0.2 g cm-2, αt = 1 × 10-4, uf = 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 power-law, the distribution is continuous apart from a jump at a12. 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: aBT (solid), asett (dash-dotted), a12 (dashed) and aP (dotted).

In the text

	Fig. A.1 Integral I2/D as function of a and b.
In the text