FROM DUST TO PLANETESIMALS: AN IMPROVED MODEL FOR COLLISIONAL GROWTH IN PROTOPLANETARY DISKS

Protoplanetary disks are universally thought to be the birthplace of planets. Furthermore, the ubiquity of terrestrial to gas-giant exoplanets discovered around other stars (Cassan et al. 2012; Howard et al. 2012; see also the compilation at http://www.circumstellardisks.org) strongly suggests that the planet formation process must be robust and efficient. This simple fact, however, has been until now rather difficult to reconcile with a number of observational properties of protoplanetary disks.

Consider for instance a scenario in which planets are created by the gradual coagulation of small dust grains into progressively larger particles, eventually leading to the formation of planetesimals and then planets. This idea turns out to be surprisingly hard to model within the constraints set by disk observations. Indeed, signatures of ongoing gas accretion onto the host star typically disappear after 10 Myr (Calvet et al. 2000). In the core-accretion scenario (Pollack et al. 1996), coagulation must be fast enough to drive growth from submicron-sized dust grains to Earth-sized rocky planets while gas is still present in sufficient quantities to form gas giants.

More crucially, particles as they grow must necessarily pass through a critical phase where their radial velocity (induced by gas drag) approaches a small but significant fraction of their orbital speed (Weidenschilling 1977). They are prone to drift into the central star within 10 to 100 orbits unless growth through that phase is rapid enough. In short, coagulation must be very efficient for planets to exist.

However, protoplanetary disks are usually detected via a characteristic near- and mid-infrared emission that is in excess of the expected pure photospheric flux of their host stars (Cohen & Kuhi 1979). Objects showing such a feature in their spectral energy distribution are known as low-mass T-Tauri or intermediate-mass Herbig AeBe stars (Gillett & Stein 1971; Strom et al. 1972). The main contribution to the excess emission comes from hot μm-sized dust grains located within a few AU of the star. The coagulation efficiency of these tiny dust grains, upon collision, is usually thought to be much higher than that of larger bodies (Dullemond & Dominik 2005). This implies that the growth process would quickly render them invisible at these wavelengths, shifting the excess emission to longer wavelengths, in contrast to what is observed (Haisch et al. 2001; Cieza et al. 2007; Su et al. 2009).

A very similar problem was raised with the observational discovery of "large" grains at large orbital radii away from the central star—more precisely, millimeter(mm)- to centimeter(cm)-sized grains, at radii at least as far out as a few tens of AU (Holland et al. 1998; Wilner et al. 2005; Rodmann et al. 2006; Lommen et al. 2009). Forming such grains at these orbital distances, where the dynamical timescale is more than a hundred times longer than around 1 AU, is just as difficult as forming meter-sized objects in the inner stellar system. Furthermore, mm- and cm-size grains at 30 AU are also subject to rapid inward drift. Again, coagulation has to be very efficient for growth to cm-size to happen despite these two problems, but cannot a priori be too efficient otherwise the outer (and inner) disk would disappear entirely on a much shorter timescale than the one observed (Garaud 2007; Hernández et al. 2008).

A possible way out of this impasse was proposed by Dullemond & Dominik (2005), who introduced the possibility of fragmentation in addition to coagulation. Fragmentation is a common outcome of high-velocity collisions (see review by Blum & Wurm 2008). As long as the fragmentation rate is high enough, and the size distribution of the fragments is sufficiently steep to favor the creation of small grains upon impact, then it is possible to maintain a healthy μm-size grain population at all times in the disk (Dullemond & Dominik 2005; Brauer et al. 2008). Unfortunately, while fragmentation solves one problem it introduces another—particle growth beyond a new "fragmentation barrier," corresponding to cm-size at a few AU and down to less than mm-size at 30 AU, seems to be impossible unless additional physics are taken into account (Brauer et al. 2008; Johansen et al. 2008).

Upon realization of the true complexity of this problem, much effort was dedicated to modeling the collisional dynamics of particles in protoplanetary disks in more detail. The ensemble approach traditionally used to study aerosol growth in planetary atmospheres (Kovetz & Olund 1969; Podolak & Podolak 1980) and grain growth in molecular clouds (Rossi et al. 1991; Ossenkopf 1993), in which one models the evolution of the particle size distribution function with the Smoluchowski coagulation/fragmentation equations (Smoluchowski 1916; Melzak 1957), has become one of the most commonly used techniques applied to the problem. We shall adopt it here too.

The Smoluchowski equations are integro-differential equations. Information on the collisional dynamics of particles must be provided in the form of "collisional kernels," which summarize all aspects of the collisions, including the relative velocities of the colliding particles, their cross sections, their sticking and/or fragmentation probabilities, the fragment distribution, etc. Formally speaking, the coagulation kernel of two particles "i" and "j" (the latter being two indices that reference the particle masses m_i and m_j for instance) should be written as

$$\begin{eqnarray} K_{ij} & =& \int _{{\cal S}_i} d{\bf S}_i \int _{{\cal S}_j} d {\bf S}_j \int _{\cal I} d{\bf I} \int _0^\infty d \Delta _{ij} \int _0^1 d \epsilon _{ij}^s \nonumber \\ && \times \big\lbrace \big[ \Delta _{ij} \epsilon ^s_{ij} \big] p_\epsilon \big(\epsilon ^s_{ij}| {\bf S}_i,{\bf S}_j,\Delta _{ij}, \delta _{ij};{\bf D} \big) p_\Delta (\Delta _{ij}|{\bf S}_i,{\bf S}_j;{\bf D}) \nonumber\\ && \times p_I({\bf I}|{\bf S}_i,{\bf S}_j;{\bf D}) p_{S_i}({\bf S}_i;{\bf D}) p_{S_j}({\bf S}_j;{\bf D}) \big\rbrace \end{eqnarray} \tag{ 1 }$$

or, in other words, a high-dimension multivariate integral over the following.

• 1.  
  All possible variations in each particle's internal properties that result in the same particle mass, loosely summarized here as the vector S, which may include shape, composition, porosity, or other, and is drawn from a multivariate probability distribution function (pdf hereafter) p_S. Note that if the masses of the two particles are different, then the pdfs for each of the particles may also be different. Since the grain composition could depend on the local temperature, the pdfs should also depend on the local disk properties, loosely described here as the vector D. Finally, p_S may also be a function of time, if the repeated collisions lead to internal changes in the particles—this has not been explicitly written here for simplicity.
• 2.  
  All possible impact configurations, loosely summarized here as the vector I, which includes impact parameter and solid angle, and drawn from the pdf p_I. If collisional outcomes are independent of I, then the integral separates and yields the mean area cross section of the collision $\bar{a}_{ij}$.
• 3.  
  All possible relative velocities Δ_ij (drawn from the pdf $p_\Delta$) for the collisions of the two particles. Note that since particle motion is induced by collisional/frictional interaction with the gas, this pdf depends on the respective structure of both particles, and on the disk model considered.
• 4.  
  All possible sticking efficiencies ^s_ij for the collisions of the two particles, drawn from the pdf p. This pdf depends on everything else, including the local disk properties.

The integrand in the square brackets is the product of the collisional relative velocity Δ_ij, and the sticking efficiency ^s_ij, discussed above.

Clearly, Equation (1) is far too complicated to be of practical use. Tanaka et al. (2005) and Dullemond & Dominik (2005) dramatically simplified the kernel expression to a product of three means, which may themselves only depend on mean quantities (and on the disk model):

$$\begin{equation} K_{ij} = \bar{a}_{ij}(\bar{\bf S}_i,\bar{\bf S}_j;{\bf D}) \bar{\Delta }_{ij}(\bar{\bf S}_i,\bar{\bf S}_j;{\bf D}) \bar{\epsilon }^s_{ij} (\bar{\bf S}_i,\bar{\bf S}_j,\bar{\Delta }_{ij};{\bf D}) \mbox{, } \end{equation} \tag{ 2 }$$

where $\bar{a}_{ij}$ is the mean collisional cross section of a pair of particles of mass m_i and m_j, assuming mean internal properties $\bar{\bf S}_i$ and $\bar{\bf S}_j$; $\bar{\Delta }_{ij}$ is the mean relative velocity and $\bar{\epsilon }^s_{ij}$ is the mean sticking efficiency of the collisions. A similar equation is used to approximate the fragmentation kernel.

Having generally adopted this much simpler formulation, studies of the evolution of the particle size distribution function then focused on refining parametric prescriptions for the mean quantities involved. Progress was made on all three fronts—in characterizing the collisional velocities, the particle structure (which determines the collisional cross section) and collisional outcomes.

Thanks to the vast advances in supercomputing, it is now possible to study the dynamics of a vast number of inertial particles interacting with the disk turbulence. Using this information, one can then infer statistical properties of individual particle velocities (Nelson & Gressel 2010; Carballido et al. 2011), and of pairwise relative velocities (Carballido et al. 2008, 2010). The latter have helped begin to validate a theoretical prescription for the rms relative particle velocity of colliding particles proposed by Ormel & Cuzzi (2007), which is commonly used in dust coagulation models to construct the quantity $\bar{\Delta }_{ij}$ used in Equation (2). However, the study of the pdfs of collisional velocities is still very much in its infancy (see Carballido et al. 2010 for a first attempt at the problem using direct simulations of turbulence, and Hubbard 2012 using shell models for turbulence).

In parallel, much work has also been done to improve our understanding of the outcome of particle collisions, and typical particle structures (the latter often depending on the collisional history) using both experimental and numerical methods. It is generally accepted that collisions involving only μm-size particles have very high sticking efficiency and that successive collisions result in the formation of fractal aggregates (see Ossenkopf 1993; also see the review by Blum & Wurm 2008). The collisions of these aggregates with one another (or with monomers) then lead to further sticking (with very high probability), but also to compaction (Dominik & Tielens 1997; Blum & Wurm 2000). The resulting objects are roughly mm- or sub-mm-sized grains that have varying degrees of porosity (Blum et al. 2006). In short, for any collision involving only sub-mm particles, one can usually adopt a mean sticking efficiency $\bar{\epsilon }_{ij}^s \simeq 1$; the calculation of the area cross section as a function of the particle mass, by contrast, is more complicated if one wishes to take into account all aforementioned effects (Ossenkopf 1993; Ormel et al. 2007; Suyama et al. 2008; Okuzumi et al. 2009; Zsom et al. 2010).

Collisions involving larger particles have much more varied outcomes, depending on their internal properties, relative velocities, and mass ratio. Given the sheer dimensionality of parameter space, and the 35 orders of magnitude mass range between mm-sized grains and Earth-sized objects, a comprehensive study of the problem is impossible. It is generally agreed that sticking continues to occur for low-velocity (alternatively low-energy) collisions, that intermediate-velocity (energy) collisions may lead to bouncing with compaction (Blum & Münch 1993; Weidling et al. 2009), and that the outcome of high-velocity collisions very much depends on the mass ratio of the two particles—with equal-mass collisions being much more likely to lead to complete fragmentation than unequal-mass collisions (see the review by Blum & Wurm 2008). In fact, high-velocity collisions in high-mass-ratio events can also lead to partial or complete sticking (Teiser & Wurm 2009; Kothe et al. 2010). If the fraction of the projectile mass that sticks to the target is larger than the total amount of mass ejected away on impact—a phenomenon called "mass transfer" hereafter—then high-mass-ratio collisions can lead to the growth of the target.

However, the details of the transitions between each of these regimes are highly dependent on the particle masses and their internal structure, on the impact parameter, and on the collisional energy. Güttler et al. (2010) recently proposed a complex model for sticking, bouncing, and fragmentation outcomes that summarizes a suite of laboratory experiments. Numerical simulations of larger, cm-sized, dust aggregates show that the fragmentation velocities may be larger than the previously assumed 100 cm s⁻¹ (Geretshauser et al. 2011a) and that the outcome may be dependent on the inhomogeneity of the aggregate (Geretshauser et al. 2011b), the target and projectile sizes, as well as the porosity of the aggregates (Meru et al. 2013).

A number of papers have been published since 2005 which gradually take into account more and more of the aforementioned effects, as well as additional disk physics (Birnstiel et al. 2010), while keeping the kernel formulation given in Equation (2).

Without giving an exhaustive list, we note in particular the work of Brauer et al. (2008) and Windmark et al. (2012a). Brauer et al. (2008) incorporate the collisional dust dynamics with a whole-disk model to study simultaneously the evolution of the particle size distribution function through coagulation and fragmentation as well as mean motions. They model the mean sticking and fragmentation probabilities using piecewise linear or piecewise parabolic functions, and consider sticking and fragmentation, but no bouncing. They also include effects such as "cratering" (or equivalently, erosion⁸), in which high-mass-ratio collisions only lead to the partial rather than the complete fragmentation of the larger body. Later on, Windmark et al. (2012a) discussed even more complex models, which include sticking, fragmentation, cratering, and bouncing, taking into account the mass dependence of the regime boundaries (Weidling et al. 2012). Furthermore, they also considered the possibility of "mass transfer," where high-mass-ratio collisions lead to the net growth of the target. Windmark et al. (2012a) showed that the latter could lead to the growth of particles well beyond the fragmentation barrier, but only if they are introduced—somewhat artificially—beyond it to start with. Nevertheless, despite all these added physics, most models fail to robustly produce planetesimals at 1 AU and cm-sized particles far out in the disk, for reasonable disk assumptions, while keeping a μm-sized dust grain population consistent with observations.

By contrast, the manner in which the coagulation and fragmentation kernels are constructed has not really been revisited since 2005. However, the statement that they can be approximated using products of three separately taken averages is only technically correct if the quantities considered are mutually independent—which is certainly not the case. We demonstrate in this paper that the failure of previous models to reconcile disk observations with exoplanetary observations simply stems from the use of Equation (2), which vastly oversimplifies the problem and results in misleading conclusions concerning particle growth in protoplanetary disks.

Our goal here is not to propose a complete new model for the construction of collisional kernels—certainly, none as complicated as the one suggested in Equation (1). Instead, we aim to demonstrate that adding even just one element of stochasticity to the model, taking into account the full pdf of relative velocities instead of their mean only, resolves a number of the aforementioned problems, at least qualitatively speaking.

To our knowledge, Okuzumi et al. (2011) were the first to take into account velocity pdfs explicitly in collisional dust growth models with application to protoplanetary disks. They focused on studying very small particles only, which interact via Brownian motion and mean drift. Here, we also consider larger particles, and include the effect of turbulence in driving stochastic motions. We loosely follow here the same general methodology first introduced by Okuzumi et al. (2011) and later in the ISIMA (International Summer Institute for Modeling in Astrophysics) report by Galvagni et al. (2011) and Windmark et al. (2012b). We construct a plausible pdf for the relative velocities of two particles interacting with the gaseous disk. We then use the latter to construct new collisional kernels, and study the resulting collisional particle growth and fragmentation. Our work nevertheless differs notably from previously published ideas, in several ways. First, we explicitly take into account the effects of turbulence in constructing the velocity pdf, in contrast with Okuzumi et al. (2011). Second, while Galvagni et al. (2011) and Windmark et al. (2012b) did take turbulence into account, they modeled the relative velocities pdf as a Maxwellian, which neglects the difference between chaotically driven particle motion (e.g., Brownian motion and turbulence) and regular particle motion (radial drift, azimuthal drift, and vertical settling induced by gas drag). Moreover, the model proposed here does take drift into account more carefully, and is much more generally applicable. Finally, we point out and correct a mathematical error in the work of Windmark et al. (2012b), which affects their coagulation and fragmentation kernels.

In what follows, we present in Section 2 the general framework associated with the use and solution of the Smoluchowski equations. The construction of the kernels is discussed in Section 3, and contrasts the traditional approach with our new proposed one. Section 4 presents a sample disk model, and Section 5 shows how the evolution of the particle size distribution function in this disk dramatically changes between the two approaches. The principal finding is that particle growth beyond the traditional bouncing and fragmentation barriers is possible, and is a robust feature of all models that include relative velocity pdfs in the kernel construction and consider the possibility of mass transfer in high-mass-ratio collisions. The particle size distribution function rapidly evolves into two populations, one composed of small dust grains and one composed of much larger particles. The quantitative details of these two populations, however, are strongly model dependent. This is discussed in Section 6, together with recommendations for future work directions.

Assuming that the masses of the dust particles can take any values among a discrete "mass-mesh" {m₁, m₂, ..., m_I}, we define N_i as the number density of particles of mass m_i. The coupled nonlinear evolution equations for the functions N_i(t) are the discrete form of the Smoluchowski coagulation–fragmentation equations (Smoluchowski 1916; Melzak 1957), in which, for i = 1...I,

$$\begin{eqnarray} \frac{d N_i}{dt} &=& \frac{1}{2} \sum _{j,k=1}^I C_{jki} K_{jk} N_j N_k - \sum _{j=1}^I K_{ij} N_i N_j \nonumber \\ && + \frac{1}{2}\sum _{j,k=1}^I F_{jk} N_j N_k N_{ijk}^f - \sum _{j=1}^I F_{ij} N_i N_j \mbox{. } \end{eqnarray} \tag{ 3 }$$

In each of these equations, the first two terms on the right-hand-side model coagulation, while the last two terms model fragmentation. The quantities K_ij and F_ij are the coagulation and fragmentation kernels, respectively, which contain crucial information on the relative velocities of particles of mass m_i and m_j, their collisional cross section, and their probability of sticking and fragmenting as discussed in Section 1. The construction of these kernels is discussed in Sections 2.2 and 3.

The first term on the right-hand side of Equation (3) models the increase in the number density of particles of mass m_i resulting from the coagulation of particles of mass m_j and m_k such that m_j + m_k ≃ m_i. If the mass-mesh is linearly spaced, then a mass-point m_i exists which is exactly equal to m_j + m_k, namely, the point for which i = j + k. In that case, this coagulation term reduces to the more commonly used $\frac{1}{2} \sum _{j=1}^{i} K_{j,i-j} N_j N_{i-j}$.

In order to follow the growth of particles across many orders of magnitude in mass, however, one cannot realistically use a linearly spaced mass-mesh—it is much more common to use a logarithmically spaced one instead. With a non-uniform mesh, however, m_j + m_k usually falls in between two existing consecutive mass-points m_{i −} and m_{i +}, with m_{i −} < m_j + m_k < m_{i +}. This leads to the introduction of the third-rank tensor C_jki, which assigns different probabilities for the coagulation of m_j and m_k into either m_{i −} or m_{i +}. The details of how to construct C_jki in such a way as to conserve mass and coagulation rates were first discussed by Kovetz & Olund (1969), and are very nicely summarized in Brauer et al. (2008). Following their work, we take

$$\begin{eqnarray*} && C_{jki} = p_{jk} \mbox{ if } i = i- \nonumber \\ && C_{jki} = (1-p_{jk}) \mbox{ if } i = i+ \nonumber \\ && C_{jki} = 0 \mbox{ otherwise,} \nonumber \end{eqnarray*}$$

where

$$\begin{equation} p_{jk} = \frac{m_{i+} - (m_j + m_k) }{m_{i+}- m_{i-}} \mbox{. } \end{equation} \tag{ 4 }$$

The second and fourth terms on the right-hand side of Equation (3) describe how particles of mass m_i disappear when they collide with particles of mass m_j, and either coagulate (second term) or fragment (fourth term).

Finally, the third term on the right-hand side of Equation (3) describes how N^f_ijk particles of mass m_i are created when particles m_j and m_k (with m_j + m_k > m_i) collide and fragment. In what follows, we assume that in any fragmentation event the number density of fragments generated is a power law with fixed index −ξ, where ξ = 11/6 = 1.83 as in Brauer et al. (2008) for instance. Since the mass bins are logarithmically distributed, we have

$$\begin{equation} N^f_{ijk} = A_{jk} m_i^{1-\xi } H(L-i) \mbox{, } \end{equation} \tag{ 5 }$$

where H(L − i) is a discrete Heaviside function (defined so that H(L − i) = 0 if i > L), L is the index of the mass-point m_L, which is the mass of the largest fragment generated, and the normalization A_jk is uniquely defined by conservation of mass during the collision:

$$\begin{equation} m_j + m_k = A_{jk} \sum _{i=1}^{L} m_i^{2-\xi } \mbox{. } \end{equation} \tag{ 6 }$$

We assume here for simplicity that m_L is either equal to the largest mass-point still satisfying m_L < m_j + m_k, or, in the occasional very high mass collisions events with m_j + m_k ⩾ m_I (where m_I is the largest mass-point in our mesh), m_L is chosen to be equal to m_I.

It can be shown that this set of equations conserves the total mass of particles in the system exactly, with

$$\begin{equation} \frac{d}{dt} \sum _{i=1}^I m_i N_i =0 \mbox{. } \end{equation} \tag{ 7 }$$

Numerical truncation errors in the evaluation of the sums, however, can cause the total mass to drift (see below).

The coagulation–fragmentation equations are intrinsically nonlinear, and—aside from very specific kernels and initial conditions—must be evolved forward in time numerically. This system of equations is also inherently very stiff, so it is crucial to use an implicit time-stepping method. Here, we have chosen to use a simple Euler first-order time-centered scheme, and calculate the Jacobian of the right-hand side explicitly. We use an automatic adaptive time-stepping scheme to ensure the required level of precision. The algorithm was successfully tested against well-known analytical solutions of these equations, reviewed by Aldous (1999).

We continuously check for conservation of the total mass. As in Brauer et al. (2008), when the ratio of the smallest to the largest particle masses considered falls below the numerical precision of the system (which is the case in some of the simulations presented below), mass conservation becomes difficult to enforce. When the total mass drifts by more than an acceptable margin (typically, one part in a million) away from the initial mass, the difference is added back to the mass distribution function in the form of monomers at the smallest size considered. Comparison of the results using this simple fix, with those obtained using a quad-precision real arithmetic (which yields reliable results up to much larger particle sizes) shows that the difference is minimal. Since quad-precision real arithmetic requires about 10 times longer integration times, we use the fix in all simulations shown below.

Finally, let us first discuss a subtle but rather important point concerning the discrete version of the Smoluchowski equations. In the derivations presented so far, we have defined N_i to be the number density of particles of mass m_i. This definition is useful for numerical purposes, but can be troublesome when interpreting the results. Indeed, consider for instance a uniformly distributed mass-mesh with a total of 100 mass-points, and a number density of N₀ particles per cm³. If the particle mass distribution function is itself uniform, then N_i = N₀/100 for all i. By contrast, had we selected 200 mass-points instead, i.e., a finer mass-mesh, then N_i = N₀/200 for all i. This simple example illustrates that N_i is a resolution-dependent quantity.

A much better approach is to consider the continuous particle mass distribution function dn/dm, defined so that, for any two masses m₁ and m₂, $\int _{m_1}^{m_2} (dn/dm) dm$ is the number density of particles of mass between m₁ and m₂. This definition is now entirely independent of the mass discretization used. In order to relate dn/dm to the numerically determined N_is, simply note that, as long as there are many mass-points between m₁ and m₂, $\int _{m_1}^{m_2} (dn/dm) dm \simeq \sum _{i = i_1}^{i_2} N_i$, where i₁ is the index of the closest mass-point to m₁, and similarly for i₂. Looking at the integral as a Riemann sum, one can then straightforwardly identify

$$\begin{equation} N_i = \left. \frac{dn}{dm}\right|_{m=m_i} dm_i = \left. \frac{dn}{dm}\right|_{m=m_i} (m_i - m_{i-1}) \mbox{. } \end{equation} \tag{ 8 }$$

In Section 5, we present all our results in terms of dn/dm; the latter are calculated from the discretized N_i using Equation (8).

While the solution of the Smoluchowski coagulation–fragmentation equation poses interesting mathematical and computational problems, most of the physical complexity of the system studied manifests itself in the construction of the coagulation and fragmentation kernels, K_ij and F_ij. The latter must be constructed for each (i, j) pair (i.e., for all collisions involving one particle of mass m_i and one particle of mass m_j). As discussed in Section 1, they have traditionally been simplified into

$$\begin{equation} K_{ij} = \bar{a}_{ij} \bar{\Delta }_{ij} \bar{\epsilon }^s_{ij} \mbox{ and } F_{ij} = \bar{a}_{ij} \bar{\Delta }_{ij} \bar{\epsilon }^f_{ij} \mbox{, } \end{equation} \tag{ 9 }$$

where $\bar{a}_{ij}$ is the mean collisional cross section of (i, j) particle pairs (averaged over all possible individual internal structures resulting in the same mass), $\bar{\Delta }_{ij}$ is an estimate of the mean relative velocity of all (i, j) pairs (averaged over all possible realizations), and $\bar{\epsilon }^s_{ij}$ is the mean sticking probability of a collisional event, given this mean structure and at that mean relative velocity (and similarly for $\bar{\epsilon }^f_{ij}$). This formulation is equivalent to assuming that all (i, j) pairs have the same collisional cross section (neglecting the possibility of variations in density through compaction, and variations in shape), collide with the same relative velocity $\bar{\Delta }_{ij}$ (ignoring the chaotic nature of particle motions), and that all these collisions have the same sticking and fragmentation probabilities $\bar{\epsilon }_{ij}^s$ and $\bar{\epsilon }_{ij}^f$ (neglecting variations with impact parameter and inherent variability in collisional outcomes).

Often used with specific oversimplified prescriptions for the fragmentation and sticking probabilities $\bar{\epsilon }^s_{ij}(\bar{\Delta }_{ij})$ and $\bar{\epsilon }^f_{ij}(\bar{\Delta }_{ij})$, this formulation has, over the years, led to the introduction and discussion of arguably artificial "bouncing" and "fragmentation" barriers. As we shall demonstrate, however, taking into account even a single stochastic component—namely, the chaotic nature of the particle motion—in the description of the collisions overcomes these barriers, and can explain the growth of large particles while retaining a population of small particles quite naturally. In what follows, we therefore still adopt a simple collisional cross-section model and a simple model for the fragmentation and sticking probabilities in individual collisions, but look at the effects of modeling the stochastic nature of the particle's relative velocities. We now describe these ingredients in more detail.

Collisional cross section. The collisional cross section of two particles i and j is typically taken to be their area cross section, although electrostatic or gravitational effects could be important for tiny dust particles and planetesimals, respectively. The calculation of an area cross section for two given particles requires knowledge of their shape and internal density, which both depend on their collisional history. Taking these effects into account properly can only be done via Monte Carlo simulations, as in Zsom & Dullemond (2008) for instance. In ensemble evolution models (such as the one presented here) one is forced instead to make simplifying assumptions. Following most previous work on the subject (Tanaka et al. 2005; Dullemond & Dominik 2005) we will consider only spherical particles of solid density ρ_s, and use

$$\begin{equation} a_{ij} = \bar{a}_{ij} = \pi (s_i+ s_j)^2 \mbox{, } \end{equation} \tag{ 10 }$$

and s_i is the assumed radius of a particle of mass m_i. The two are related via $m = \frac{4}{3} \pi \rho _s s^3$, where ρ_s = 1 g cm⁻³ as in Windmark et al. (2012b).

Sticking, bouncing, and fragmenting. The complex dependence of possible collisional outcomes on the collisional velocity and particle properties was discussed in Section 1, and is still very much the subject of ongoing investigations. However, since our purpose is not to give quantitatively accurate predictions for the evolution of the particle size distribution function, but rather to study qualitatively the impact of the introduction of our new model, we use the simplest possible prescription for the fragmentation and sticking probabilities. As in Windmark et al. (2012b) we assume that, in any specific collisional event, particles stick for Δ_ij < v_b and fragment if Δ_ij > v_f. Between v_b and v_f lies the "bouncing" region. Hence we select, for individual collisions,

$$\begin{eqnarray} &&\epsilon _{ij}^f = H(\Delta _{ij} - v_f) \mbox{, } \nonumber \\ &&\epsilon _{ij}^s = H(v_b - \Delta _{ij}) \mbox{, } \end{eqnarray} \tag{ 11 }$$

where H is a Heaviside function. In all that follows, we use for the sake of comparison with the work of Windmark et al. (2012b), the values v_b = 5 cm s⁻¹ (unless otherwise specified) and v_f = 100 cm s⁻¹.

Relative velocity. The various possible sources of relative motion for solid particles in protoplanetary disks, as reviewed by Weidenschilling & Cuzzi (1993) for instance, can be divided into two categories: regular motion caused by frictional interaction with the mean component of the gas velocity (usually divided into radial drift, azimuthal drift, and vertical settling), and chaotic motion caused by collisions with gas molecules (Brownian motion) and interaction with the turbulent component of the gas velocity.

Let us begin by considering the motion of particles in a gas at rest (i.e., supposing that the macroscopic gas velocity is zero in the frame of reference considered). Each particle undergoes Brownian motion via collisions with the gas molecules, and thereby acquires a random velocity, which is isotropic and has an amplitude drawn from a Maxwellian pdf. Since the velocities of two particles both undergoing Brownian motion are uncorrelated, the pdf of their relative velocity Δ_ij is also a Maxwellian with

$$\begin{equation} p_{\rm B}(\Delta _{ij}) = \sqrt{\frac{2}{\pi }} \frac{\Delta _{ij}^2}{\big(\sigma ^{\rm B}_{ij} \big)^3} \exp \left(- \frac{\Delta _{ij}^2}{2 \big(\sigma ^{\rm B}_{ij} \big)^2 }\right) \mbox{, } \end{equation} \tag{ 12 }$$

where (σ^B_ij)² = kT(m_i + m_j)/m_im_j, k is the Boltzmann constant, and T is the local gas temperature. The mean value of the relative collisional velocities is

$$\begin{equation} \bar{\Delta }_{ij}^B = \sqrt{ \frac{8 kT (m_i + m_j)}{\pi m_i m_j }} = \sqrt{\frac{8}{\pi }} \sigma ^{\rm B}_{ij}\mbox{. } \end{equation} \tag{ 13 }$$

Note that, as defined, σ^B_ij is exactly $1/\sqrt{3}$ times the rms collisional velocity and $\bar{\Delta }_{ij}^B$ is $\sqrt{8/3\pi }$ times the rms collisional velocity. This implies that for a Maxwellian, to a good approximation, the mean and rms velocities are equal to one another—a well-known result that we shall use later.

In protoplanetary disks, however, particles can have significant mean velocities with respect to the gas. This is particularly true for the larger particles, which orbit around the star at near-Keplerian speeds, while the azimuthal gas velocity is notably sub-Keplerian (since the gas is pressure-supported while the particles are not). As a result of this differential motion, collisions with gas molecules transfer net momentum to the particles. This effect is often modeled as an added drag term in the particle's equation of motion, and leads to radial and azimuthal drift within the disk (Whipple 1972; Weidenschilling 1977; Nakagawa et al. 1986), as well as net vertical settling toward the mid-plane over time (Goldreich & Ward 1973; Dubrulle et al. 1995; Garaud et al. 2004).

Let us consider a frame that is rotating with the local Keplerian angular velocity Ω_K. In that frame, assuming that the surface density of solids is much smaller than the surface density of the gas, particles of mass m_i have mean radial and azimuthal velocities $\bar{u}_i$ and $\bar{v}_i$, respectively, with

$$\begin{eqnarray} && \bar{u}_i = \frac{1}{S^2_i + 1} (u_g - 2 S_i\eta v_{\rm K} ) \mbox{, } \nonumber \\ && \bar{v}_i = \frac{1}{2S_i }(\bar{u}_i - u_g) \mbox{, } \end{eqnarray} \tag{ 14 }$$

where u_g is the radial gas velocity, and η is related to the deviation of the gas azimuthal velocity from the local azimuthal Keplerian velocity v_K = rΩ_K. All these quantities depend on the selected disk model, and are presented in more detail in Section 4. Finally, S_i is the Stokes number of particles of mass m_i, defined as

$$\begin{equation} S_i = \tau _i \Omega _{\rm K} = \frac{s_i \rho _s}{\rho _m c} \Omega _{\rm K} = \frac{ \sqrt{2\pi } s_i \rho _s }{\Sigma } \mbox{, } \end{equation} \tag{ 15 }$$

where τ_i is the particle stopping time, ρ_m is the mid-plane disk density, c is the local sound speed, and Σ is the local surface density of the gas (see Section 4). Note that we have assumed here that particles remain in the Epstein gas drag regime, which is true for small particles but may not be accurate above a certain size. This assumption should be dropped should one require a quantitatively more accurate prediction for the evolution of the particle size distribution function, but is satisfactory within our stated qualitative goals (see Section 1).

The mean settling velocity of a particle of mass m_i is calculated as in Birnstiel et al. (2010):

$$\begin{equation} \bar{w}_i = - h_i \Omega _{\rm K} \min \left(S_i, \frac{1}{2} \right) \mbox{, } \end{equation} \tag{ 16 }$$

where h_i is the vertical scale height of particles of mass m_i, and where the minimum function guarantees that the vertical velocity is at most equal to that of the epicyclic motion. The particle disk scale height h_i can be obtained by solving an advection–diffusion equation (Dubrulle et al. 1995; Garaud 2007), which yields

$$\begin{equation} \frac{1}{h_i^2} = \frac{ S_i \Omega _K }{D_i} + \frac{1}{h^2}\mbox{, } \end{equation} \tag{ 17 }$$

where h is the local gas disk scale height, $D_i = \frac{\nu }{1+S_i}$ is the reduced turbulent diffusivity for the particles of size s_i, and ν is the local turbulent viscosity of the gas (see Section 4).

As should be obvious from these derivations, particles of different sizes have different velocities relative to that of the gas. As a result, they also acquire relative velocities with respect to each other. This time, however, no random effects are involved so the mean value of the relative velocities of particles i and j, both induced by gas drag, is

$$\begin{equation} \bar{\Delta }^D_{ij} = \sqrt{ (\bar{u}_i - \bar{u}_j)^2 + (\bar{v}_i - \bar{v}_j)^2 + (\bar{w}_i - \bar{w}_j)^2}\mbox{, } \end{equation} \tag{ 18 }$$

and the pdf is well approximated by a δ-function centered around the mean.

However, the effect of gas drag described above is only strictly valid in a "non-turbulent accretion disk"—somewhat of an oxymoron. As first described by Voelk et al. (1980), particles interacting with turbulent eddies have inherently stochastic motions. Their velocity pdf is non-isotropic (when the turbulence itself is anisotropic), and depends sensitively on the particle stopping time compared with the eddy turnover time at the energy dissipation and injection scales (Voelk et al. 1980; Ormel & Cuzzi 2007). To add to the complexity of the problem, particle trajectories are no longer necessarily independent of one another—two small particles trapped in the same eddy have strongly correlated motion. Recent progress has been made to characterize the problem using full numerical simulations (Carballido et al. 2010) and simplified vortex gas models (Rast & Pinton 2011), although the results still have limited applicability in both cases. A theoretical model of the rms relative velocities in turbulence, which we assume is very similar to their mean relative velocity $\bar{\Delta }_{ij}^T$ by analogy with the case of Brownian motion, was proposed by Ormel & Cuzzi (2007, see their Equations (27)–(29)):

$$\begin{eqnarray} \bar{\Delta }_{ij}^{T} &=& \sqrt{\frac{S_i - S_j}{S_i + S_j}}\left[ \frac{S_i^2}{S_i + {\rm Re}^{-1/2}} - \frac{S_j^2}{S_j + {\rm Re}^{-1/2}} \right]^{1/2} \nonumber\\ && \times v_e \mbox{ if } S_{i,j} \le {\rm Re}^{-1/2} \mbox{, } \nonumber \\ \bar{\Delta }_{ij}^{T}&=& \left[ 2.2 S_i - S_j + \frac{2 S^2_i}{S_i + S_j} \left(\frac{1}{2.6} + \frac{S_j^3}{1.6 S_i^3 + S_i^2 S_j} \right) \right] ^{1/2} \nonumber\\ && \times v_e \mbox{ if } {\rm Re}^{-1/2} \le S_i \le 1 \mbox{ for all } j<i \mbox{, } \nonumber \\ \bar{\Delta }_{ij}^T &=& \left[ \frac{1}{1 + S_i} + \frac{1}{1 + S_j} \right]^{1/2} v_e \mbox{, } \mbox{ if } 1 < S_i \mbox{ for all } j<i, \end{eqnarray} \tag{ 19 }$$

where Re is the local Reynolds number, which we will take to be 10⁸ (see Section 4), $v_e = \sqrt{\alpha } c$ is the typical eddy velocity, and α is the turbulent intensity parameter (see Section 4). Note that in these equations we have assumed that S_i ⩾ S_j. If the opposite is true, then the indices should be switched.

As discussed in the previous section, the collisional and fragmentation kernels are usually given by Equation (9). The problem then reduces to calculating a single value for the mean relative velocities of the two particles $\bar{\Delta }_{ij}$, combining information from all the possible sources of motion described in Section 2.2. Traditionally, $\bar{\Delta }_{ij}$ is constructed as in Tanaka et al. (2005), with

$$\begin{equation} \bar{\Delta }_{ij} = \sqrt{ \big(\bar{\Delta }^B_{ij} \big)^2 + \big(\bar{\Delta }^D_{ij} \big)^2 + \big(\bar{\Delta }^T_{ij} \big)^2 } \mbox{, } \end{equation} \tag{ 20 }$$

and the fragmentation and coagulation probabilities $\bar{\epsilon }^s_{ij}$ and $\bar{\epsilon }^f_{ij}$ are simply taken to be

$$\begin{equation} \bar{\epsilon }_{ij}^s = \epsilon _{ij}^s(\bar{\Delta }_{ij}) \mbox{ and } \bar{\epsilon }_{ij}^f = \epsilon _{ij}^f(\bar{\Delta }_{ij}) \mbox{, } \end{equation} \tag{ 21 }$$

where the functions ^s_ij and ^f_ij are constructed (for instance⁹) as in Equation (11).

In order to take into account the stochastic nature of the particle velocities, collisional and fragmentation kernels need to be rewritten¹⁰ as

$$\begin{eqnarray} && K_{ij} = \bar{a}_{ij} \int _0^\infty \Delta _{ij} p(\Delta _{ij}) \epsilon ^s_{ij}(\Delta _{ij}) d \Delta _{ij} \mbox{, } \nonumber \\ && F_{ij} = \bar{a}_{ij} \int _0^\infty \Delta _{ij} p(\Delta _{ij}) \epsilon ^f_{ij}(\Delta _{ij}) d\Delta _{ij} \mbox{, } \end{eqnarray} \tag{ 22 }$$

which one may recognize as a much-simplified version of Equation (1), where p(Δ_ij) is a single pdf for the relative velocities of the two particles (dropping the suffix Δ on the p for clarity of notation), which ought to combine information from all the possible sources of motion described in Section 2.2. The main difficulty here comes from the fact that, while the pdfs of particles undergoing Brownian motion or mean drift are known, we still do not have any knowledge of the shape of the pdf of relative velocities for particles undergoing collisions via turbulent motions. One is left to choose the latter somewhat arbitrarily.

Galvagni et al. (2011) and Windmark et al. (2012b) both proposed that Δ_ij should be distributed as a Maxwellian. Galvagni et al. (2011) considered all possible regimes simultaneously, and constructed the Maxwellian so that its mean is given by $\bar{\Delta }_{ij}$ as defined in Equation (20). Windmark et al. (2012b) only considered the turbulent regime, and required that the mean of the Maxwellian be $\bar{\Delta }^T_{ij}$ instead. Both models naturally yield the same pdf in a regime that is dominated by turbulent motions (i.e., when $\bar{\Delta }^T_{ij} \simeq \bar{\Delta }_{ij} \gg \bar{\Delta }^D_{ij}$) but both models fail in the limit where the system is dominated by regular motion. The proposed model by Galvagni et al. (2011) overestimates the particle dispersion in that limit, while Windmark et al. (2012b) do not address the question.

In what follows, we build on the work of Okuzumi et al. (2011) and propose an alternative, more rigorous approach for the construction of p(Δ_ij), which takes into account the fact that individual particles have deterministic velocities (induced by radial drift, azimuthal drift, and vertical settling), in addition to stochastic motions induced by Brownian motion and turbulence. For simplicity, in all that follows we assume that the stochastic motions are isotropic.¹¹

Let us focus first on a given direction (taking here, for the sake of illustration, the radial direction). We model the one-dimensional velocity pdf of particle i as a Gaussian with mean $\bar{u}_i$ (the mean radial drift velocity discussed earlier) and standard deviation σ_i (specified later). Hence,

$$\begin{equation} p_r(u_i) = \frac{1}{\sqrt{2\pi } \sigma _i } \exp \left(-\frac{(u_i-\bar{u}_i)^2}{2\sigma _i^2}\right) \mbox{, } \end{equation} \tag{ 23 }$$

where the index r denotes the radial direction. Note that, in this context, u_i can be positive or negative. Given a second particle indexed by j, with a similarly defined radial velocity pdf, a well-known (although non-trivial) result from probability theory is that the pdf of their relative radial velocities is

$$\begin{eqnarray} p_r (u_i - u_j) &=& p_r(\Delta _{r,ij})\qquad\nonumber\\ & =& \frac{1}{\sqrt{2 \pi } \sigma _{ij}} \exp \left[ - \frac{(\Delta _{r,ij} -(\bar{u}_{i} - \bar{u}_{j}))^2}{2 \sigma _{ij}^2} \right]\mbox{, }\qquad \end{eqnarray} \tag{ 24 }$$

where σ²_ij = σ_i² + σ²_j. Note that Δ_{r, ij} = u_i − u_j can be positive or negative. Furthermore, it changes sign upon permutation of i and j.

Similar expressions can be obtained for the pdfs of relative velocities in the azimuthal (φ) and vertical (z) directions. Assuming independence of the distributions in each respective direction, the pdf of the three-dimensional relative velocity is simply the product of the three derived pdfs:

$$\begin{equation} p_{{\rm 3D}}({\mathbf {\Delta }}_{ij}) d {\mathbf {\Delta }}_{ij} = p_r(\Delta _{r,ij}) p_\varphi (\Delta _{\varphi ,ij}) p_z(\Delta _{z,ij}) d \Delta _{r,ij} d \Delta _{\varphi ,ij} d \Delta _{z,ij} \mbox{. } \end{equation} \tag{ 25 }$$

From this expression, with a little bit of algebra (see Appendix A) and assuming isotropy in the particle stochastic motion (i.e., assuming that σ_i and σ_j are independent of the direction selected), one can show that the pdf of the amplitude of the relative velocity Δ_ij = |Δ_ij| is

$$\begin{eqnarray} p(\Delta _{ij}) &=& \frac{1}{\sqrt{2\pi } \sigma _{ij}} \frac{\Delta _{ij}}{\bar{\Delta }^D_{ij}} \Bigg[ \exp \left(-\frac{ \big(\Delta _{ij} - \bar{\Delta }^D_{ij} \big)^2 }{ 2\sigma _{ij}^2} \right) \nonumber\\ && - \exp \left(-\frac{ \big(\Delta _{ij} + \bar{\Delta }^D_{ij} \big)^2 }{ 2\sigma _{ij}^2} \right) \Bigg] \mbox{, } \end{eqnarray} \tag{ 26 }$$

where $\bar{\Delta }^D_{ij}$ is given by Equation (18). This time, both Δ_ij and $\bar{\Delta }^D_{ij}$ are always positive by construction, and invariant under permutation of i and j. It can be shown that this pdf is always positive, and that the integral over all possible values of Δ_ij is indeed 1, as required.

In the limit where $\bar{\Delta }^D_{ij} \ll \sigma _{ij}$, which corresponds to cases where the deterministic relative velocities of the particles are much smaller than their rms velocities (i.e., in the case where relative motion is dominated by Brownian motion or turbulent motion), p(Δ_ij) reduces to a Maxwellian distribution (to lowest order in a Taylor expansion in $\bar{\Delta }^D_{ij}$), with

$$\begin{equation} p(\Delta _{ij}) \simeq \sqrt{\frac{2}{\pi }} \frac{\Delta _{ij}^2}{\sigma _{ij}^3} \exp \left(- \frac{\Delta _{ij}^2}{2\sigma _{ij}^2} \right) \mbox{, } \end{equation} \tag{ 27 }$$

which has a mean

$$\begin{equation} \bar{\Delta }_{ij} = \sqrt{\frac{8}{\pi }} \sigma _{ij}\mbox{. } \end{equation} \tag{ 28 }$$

In order to ensure that the mean velocity of that distribution is the same as that obtained from Brownian motion or turbulent motion when these dominate the dynamics of the system, we take

$$\begin{equation} \sigma ^2_{ij} = \frac{\pi }{8} \big[ \big(\bar{\Delta }^B_{ij} \big)^2 + \big(\bar{\Delta }^T_{ij} \big)^2 \big] \mbox{, } \end{equation} \tag{ 29 }$$

where $\bar{\Delta }^B_{ij}$ and $\bar{\Delta }^T_{ij}$ are defined in Equations (13) and (19), respectively. Our formalism recovers the idea proposed by Galvagni et al. (2011) and by Windmark et al. (2012b) in the limit where turbulence is the dominant factor. Note that, as defined, σ_ij is not the rms velocity of the distribution defined in Equation (26), although it is related to its variance. In the limit where the mean drift velocity is close to zero, σ_ij is, as discussed earlier, $1/\sqrt{3}$ times the rms velocity. Our derivation is so far very similar to that of Okuzumi et al. (2011), and recovers their relative velocity pdf in the limit where the particle dispersion is dominated by Brownian motion.

When $\bar{\Delta }^D_{ij} \gg \sigma _{ij}$, which correspond to cases where the deterministic relative velocities of the particles are much larger than their rms velocities, the second Gaussian term in Equation (26) is negligible, and the relative velocity pdf reduces to

$$\begin{equation} p(\Delta _{ij}) \simeq \frac{1}{\sqrt{2\pi } \sigma _{ij}} \frac{\Delta _{ij}}{\bar{\Delta }^D_{ij}} \exp \left(-\frac{\big(\Delta _{ij} - \bar{\Delta }^D_{ij} \big)^2 }{ 2\sigma _{ij}^2} \right) \mbox{. } \end{equation} \tag{ 30 }$$

It is interesting to note that this is not a Gaussian, but has a somewhat more significant tail for larger values of Δ_ij.

It is worth noting that the mean relative drift velocity, $\bar{\Delta }^D_{ij}$, is not equal to the actual mean relative velocity $\bar{\Delta }_{ij}$; the latter is given by

$$\begin{eqnarray} \bar{\Delta }_{ij} &=& \int _0^\infty \Delta _{ij} p(\Delta _{ij}) d \Delta _{ij} = \frac{ \sigma _{ij}^2 + \big(\bar{\Delta }^D_{ij} \big)^2}{\bar{\Delta }^D_{ij}} {\rm erf}\left(\frac{\bar{\Delta }^D_{ij}}{\sqrt{2} \sigma _{ij}}\right) \nonumber\\ && +\sqrt{\frac{2}{\pi } } \sigma _{ij} \exp \left(- \frac{\big(\bar{\Delta }^D_{ij} \big)^2}{2 \sigma _{ij}^2} \right) \mbox{. } \end{eqnarray} \tag{ 31 }$$

The relationship between $\bar{\Delta }^D_{ij}$ and $\bar{\Delta }_{ij}$ is shown in Figure 1. We see that $\bar{\Delta }_{ij}$ tends to $\bar{\Delta }^D_{ij}$ when $\bar{\Delta }^D_{ij} \gg \sigma _{ij}$, i.e., when the mean particle velocity is large compared with the rms particle velocity, and to $\sqrt{8/\pi } \sigma _{ij}$ when $\bar{\Delta }^D_{ij} \ll \sigma _{ij}$, i.e., when the mean particle velocity is small compared with the rms particle velocity. This behavior is consistent with expectations.

Zoom In Zoom Out Reset image size

Finally, we can now calculate the new coagulation and fragmentation kernels, following Equation (22), by convolving the velocity pdf constructed in Equation (26) with the individual collision sticking and fragmentation probabilities. With the simple expressions for ^s_ij and ^f_ij given by Equation (11), it is possible to calculate these kernels analytically¹²

$$\begin{eqnarray} F_{ij} &=& \bar{a}_{ij} \int _0^\infty \Delta _{ij} p(\Delta _{ij}) \epsilon _{ij}^f(\Delta _{ij}) d\Delta _{ij} = \bar{a}_{ij} \int _{v_f}^\infty \Delta _{ij} p(\Delta _{ij}) d\Delta _{ij} \nonumber \\ & =& \bar{a}_{ij} \left\lbrace \frac{\big(\bar{\Delta }^D_{ij} \big)^2 + \sigma _{ij}^2}{2 \bar{\Delta }^D_{ij}} A(v_f) + \frac{\sigma _{ij}}{\sqrt{2\pi } \bar{\Delta }^D_{ij}} B(v_f) \right\rbrace \mbox{, } \end{eqnarray} \tag{ 32 }$$

and

$$\begin{eqnarray} K_{ij} &=& \bar{a}_{ij} \int _0^{\infty } \Delta _{ij} p(\Delta _{ij}) \epsilon _{ij}^s(\Delta _{ij}) d \Delta _{ij} = \bar{a}_{ij} \int _0^{v_b} \Delta _{ij} p(\Delta _{ij}) d \Delta _{ij} \nonumber\\ & =& \bar{a}_{ij} \Bigg\lbrace \frac{\big(\bar{\Delta }^D_{ij} \big)^2 + \sigma _{ij}^2}{2 \bar{\Delta }^D_{ij}} \left[ 2 {\rm erf} \left(\frac{ \bar{\Delta }^D_{ij} }{\sqrt{2} \sigma _{ij}} \right) - A(v_b) \right] \nonumber\\ && + \frac{\sigma _{ij}}{\sqrt{2\pi } \bar{\Delta }^D_{ij}} \left[ 2 \bar{\Delta }^D_{ij} \exp \left(- \frac{ \big(\bar{\Delta }^D_{ij} \big)^2}{2\sigma _{ij}^2} \right) - B(v_b) \right] \Bigg\rbrace \mbox{, } \end{eqnarray} \tag{ 33 }$$

where for simplicity of notation we have defined

$$\begin{eqnarray} A(v) &=& {\rm erf} \left(\frac{v + \bar{\Delta }^D_{ij}}{\sqrt{2} \sigma _{ij}} \right) - {\rm erf} \left(\frac{v - \bar{\Delta }^D_{ij}}{\sqrt{2} \sigma _{ij}} \right) \mbox{, } \nonumber \\ B(v) &=& \big(v + \bar{\Delta }^D_{ij} \big) \exp \left(- \frac{ \big(v- \bar{\Delta }^D_{ij} \big)^2}{2\sigma _{ij}^2}\right) \nonumber\\ && - \big(v- \bar{\Delta }^D_{ij} \big) \exp \left(- \frac{ \big(v + \bar{\Delta }^D_{ij} \big)^2}{2\sigma _{ij}^2}\right) \mbox{. } \end{eqnarray} \tag{ 34 }$$

Note that if $\bar{\Delta }^D_{ij} \ll \sigma _{ij}$ then (to lowest order in a Taylor expansion in $\bar{\Delta }^D_{ij}$)

$$\begin{eqnarray} F_{ij} &=& \bar{a}_{ij} \bar{\Delta }_{ij} \frac{2\sigma ^2_{ij} + v_f^2}{2\sigma ^2_{ij}} \exp \left(- \frac{ v_f^2}{2\sigma _{ij}^2}\right) \mbox{, } \nonumber \\ K_{ij} &=& \bar{a}_{ij} \bar{\Delta }_{ij} \left[ 1 -\frac{2\sigma ^2_{ij} + v_b^2}{2\sigma ^2_{ij}} \exp \left(- \frac{ v_b^2}{2\sigma _{ij}^2}\right) \right] \mbox{. } \end{eqnarray} \tag{ 35 }$$

If desired, one can thereby construct "mean" sticking and fragmentation probabilities

$$\begin{eqnarray} && \bar{\epsilon }_{ij}^s = \frac{K_{ij}}{\bar{a}_{ij} \bar{\Delta }_{ij}} = \frac{ \int _0^\infty \Delta _{ij} p(\Delta _{ij}) \epsilon _{ij}^s(\Delta _{ij}) d \Delta _{ij} }{ \int _0^\infty \Delta _{ij} p(\Delta _{ij}) d \Delta _{ij} } \mbox{, } \nonumber \\ && \bar{\epsilon }_{ij}^f = \frac{F_{ij}}{\bar{a}_{ij} \bar{\Delta }_{ij}} = \frac{ \int _0^\infty \Delta _{ij} p(\Delta _{ij}) \epsilon _{ij}^f(\Delta _{ij}) d \Delta _{ij} }{ \int _0^\infty \Delta _{ij} p(\Delta _{ij}) d \Delta _{ij} }\mbox{, } \end{eqnarray} \tag{ 36 }$$

so that, using these expressions, we have an exact analog for Equation (9). Note that the second term in each of these expressions is a rather general formula for the mean sticking and fragmentation probabilities, which is valid for any assumed velocity pdf, and any assumed individual collision and fragmentation probabilities, but does assume that the collisional outcomes are independent of the particle density and shape.

It can easily be verified that when ^s_ij(Δ_ij) + ^f_ij(Δ_ij) = 1 for individual particles, then we also have $\bar{\epsilon }_{ij}^s + \bar{\epsilon }_{ij}^f = 1$, or equivalently,

$$\begin{equation} K_{ij} + F_{ij} = \bar{a}_{ij} \bar{\Delta }_{ij}\mbox{. } \end{equation} \tag{ 37 }$$

This happens for instance when v_b = v_f (i.e., in the absence of bouncing) in the simple Heaviside model given by Equation (11), and can be verified directly by adding Equations (32) and (33).

Note that the coagulation kernel derived in Equation (33) looks like the one presented by Okuzumi et al. (2011) in their Equation (12), but is slightly different because the sticking probability of individual collisions in their paper is taken to be a piecewise linear function of the collisional energy rather than a Heaviside function of the collisional velocity. Our results would of course recover one another had we chosen the same sticking probability function, and in the limit where the particle dispersion is not affected by turbulence (i.e., in the limit where $\bar{\Delta }^T_{ij} = 0$). By contrast, Windmark et al. (2012b) would not recover the same formula in any limit.

Indeed, note that Windmark et al. (2012b) discuss the role of the total sticking and fragmentation probabilities (see their Equations (4) and (5)), in our notation written as

$$\begin{equation} \big[\epsilon _{ij}^s \big] = \int _0^{v_b} p_M(\Delta _{ij}) d\Delta _{ij} \quad \mbox{and} \quad \big[\epsilon _{ij}^f \big] = \int _{v_f}^\infty p_M(\Delta _{ij}) d\Delta _{ij} \end{equation} \tag{ 38 }$$

(where p_M denotes a Maxwellian distribution) rather than the mean sticking and fragmentation velocities $\bar{\epsilon }_{ij}^s$ and $\bar{\epsilon }_{ij}^f$. However, as should be clear from the derivation presented above, these "total" probabilities play no role in the calculation of the kernels—the means must be used instead. Their results should therefore be considered incorrect even though they look qualitatively similar to ours (see Section 5).

Figure 2 shows $\bar{\epsilon }_{ij}^s$ and $\bar{\epsilon }_{ij}^f$, as calculated in Equation (36), as a function of σ_ij and $\bar{\Delta }_{ij}^D$, for v_b = 5 cm s⁻¹ and v_f = 100 cm s⁻¹ (see Section 2.2). We see that when both σ_ij and $\bar{\Delta }_{ij}^D$ are small compared with the bouncing threshold, sticking is very probable, while when either σ_ij or $\bar{\Delta }_{ij}^D$ are large compared with v_f, fragmentation is very probable. The transition between the two regimes, however, is much smoother when σ_ij is large than when it is small; this behavior is, again, consistent with expectations.

Zoom In Zoom Out Reset image size

It is important to note that the effective sticking and fragmentation probabilities $\bar{\epsilon }_{ij}^s$ and $\bar{\epsilon }_{ij}^f$ are never zero, by contrast with the individual collision functions ^s_ij and ^f_ij. This result is very general, and stems from the fact that $\bar{\epsilon }_{ij}^s$ and $\bar{\epsilon }_{ij}^f$ are convolution integrals of ^s_ij and ^f_ij with the assumed particle velocity pdfs times the particle relative velocity. As a result, fragmentation is always possible (albeit unlikely) even when the mean relative velocity of two particles is very small, and sticking is always possible (albeit unlikely) even when the mean relative velocity of two particles is very large. This has a number of important consequences, as shown in Section 5.

We assume that the surface density of the gas follows a truncated power law:

$$\begin{equation} \Sigma (r) = \frac{M}{2\pi R r} e^{-r/R}\mbox{, } \end{equation} \tag{ 39 }$$

where M is the disk mass, r is the orbital radius, and R is the truncation radius. As discussed by Lynden-Bell & Pringle (1974; see also Hartmann et al. 1998; Garaud 2007), this model has the desirable property of being an attracting self-similar solution of the equations describing the "viscous spreading" of the disk, as long as the turbulent viscosity ν varies linearly with radius:

$$\begin{equation} \nu (r) = \nu _{\rm AU} \frac{r}{1\,{\rm AU}} \mbox{, } \end{equation} \tag{ 40 }$$

where the subscript AU refers to a quantity measured at 1 AU. Conservation of angular momentum implies that the radial velocity of the gas is

$$\begin{equation} u_g(r) = -\frac{3}{r^{1/2} \Sigma (r)} \frac{\partial }{\partial r} (r^{1/2} \nu (r) \Sigma (r) ) = - 3 \frac{\nu _{\rm AU} }{1\,{\rm AU}} \left(\frac{1}{2} - \frac{r}{R} \right) \mbox{, } \end{equation} \tag{ 41 }$$

while the radial force balance near the mid-plane implies that the azimuthal velocity of the gas is

$$\begin{equation} v_g(r) = v_{\rm K}(r) + \frac{1}{2\rho _m(r) \Omega _K(r) } \frac{\partial p_m}{\partial r} = (1-\eta (r)) v_{\rm K}(r) \mbox{, } \end{equation} \tag{ 42 }$$

where p_m is the mid-plane gas pressure, and where

$$\begin{equation} \eta (r) = - \frac{1}{2r \rho _m(r) \Omega ^2_{\rm K}(r)} \frac{\partial p_m}{\partial r} \mbox{. } \end{equation} \tag{ 43 }$$

As illustrated by Equations (41) and (42), the local gas velocity depends on the local thermodynamical structure of the disk. We now turn to describing the latter in more detail.

Assuming that the disk is thin, vertically isothermal and in hydrostatic equilibrium implies that the density profile is

$$\begin{equation} \rho (r,z) = \rho _m(r) \exp \left(- \frac{z^2}{2 h^2(r)} \right) \mbox{, } \end{equation} \tag{ 44 }$$

where the disk scale height h(r) is related to the local sound speed c(r) via

$$\begin{equation} h(r) = \frac{c(r) }{\Omega _K(r) }\mbox{. } \end{equation} \tag{ 45 }$$

Integrating the density profile across the disk yields

$$\begin{equation} \rho _m(r) = \frac{\Sigma (r)}{\sqrt{2\pi } h(r) } \mbox{. } \end{equation} \tag{ 46 }$$

Meanwhile, the mid-plane gas pressure is given by the equation of state (which we assume here to be a perfect gas), so

$$\begin{equation} p_m(r) = \frac{{\cal R} \rho _m(r) T(r) }{\mu } = c^2(r) \rho _m(r) \mbox{, } \end{equation} \tag{ 47 }$$

where μ is the mean molecular weight of the gas and ${\cal R}$ is the gas constant.

In both cases, we need to know the local sound speed to evaluate ρ_m and p_m. Following the standard α model for turbulent viscosity, we use

$$\begin{equation} \nu (r)= \alpha c(r) h(r) \mbox{, } \end{equation} \tag{ 48 }$$

where α is usually assumed to be constant in the entire disk. Combining Equations (40) and (48) we have

$$\begin{eqnarray} c^2(r) &=& c^2_{\rm AU} \left(\frac{r}{1\,{\rm AU}} \right)^{-1/2} \mbox{ where } c^2_{\rm AU} = \frac{\nu _{\rm AU}}{\alpha } \Omega _{\rm AU} \mbox{, } \nonumber \\ h(r) &=& h_{\rm AU} \left(\frac{r}{1\,{\rm AU}} \right)^{5/4} \mbox{ where } h_{\rm AU} = \frac{c_{\rm AU}}{\Omega _{\rm AU}} \mbox{. } \end{eqnarray} \tag{ 49 }$$

Having laid out the governing equations for the disk, we need to choose its actual parameters. For the sake of comparison with the work of Windmark et al. (2012b), we select the parameters for our disk to have the same properties at 1 AU (see their Table 1). This immediately yields α = 10⁻⁴, c_AU = 10⁵ cm s⁻¹, and the Reynolds number Re = 10⁸. Using a mean molecular weight of μ = 2.3, the temperature at 1 AU is 276 K, close to the value reported by Windmark et al. (2012b). They choose a gas surface density of Σ = 1700 g cm⁻², which we can recover to a good approximation by selecting (for instance¹³), a central star of mass M_⋆ = 0.75 M_☉, surrounded by a disk of mass M_d = 0.05 M_⋆, with a R = 30 AU truncation radius. Using these parameters we find that the viscous timescale at 1 AU is about 400,000 years, so the gas disk would not evolve significantly during the simulation.

At 30 AU, our disk model has a local gas surface density of 21 g cm⁻² and a temperature of 50 K. At this radius, the viscous evolution timescale is larger than 10 Myr, so again, the gas disk would not evolve significantly during the simulation. Finally, in all that follows we select a dust-to-gas mass ratio of Z = 0.01. This ensures that the dust dynamics cannot influence the gas dynamics, as assumed throughout this work.

In what follows, we shall be interested in two representative disk regions, located close to and far from the central star, respectively. The close-in region is selected to be at 1 AU, for ease of comparison with the work of Windmark et al. (2012b), and the far-out region is selected to be at 30 AU, to address the problem of particle growth beyond cm-size discussed in Section 1.

As shown in Figure 3, relative particle motions are dominated by turbulence for sizes up to about 10 cm at 1 AU. Collisions induced by differential radial drift dominate for larger particles. Hence, we expect our model to become more relevant than the one proposed by Galvagni et al. (2011) and Windmark et al. (2012b) if growth occurs beyond 10 cm at 1 AU.¹⁴ At 30 AU, this transition happens for mm- to cm-sized particles.

Zoom In Zoom Out Reset image size

Figure 3 also shows the logarithm of the mean relative velocity $\bar{\Delta }_{ij}$ of all (i, j) particle pairs, and contours of the bouncing and fragmentation thresholds v_b = 5 cm s⁻¹ and v_f = 100 cm s⁻¹, respectively. In what follows, we shall often refer to these contours as the "original" bouncing and fragmentation barriers, as they often play that role in simulations which use the traditional method for kernel construction (see Section 3.1). In any case, we find that bouncing is expected to affect particle growth beyond mm-size at 1 AU, and 10 μm size at 30 AU, for the selected value of v_b. Fragmentation is expected to be important for particles around cm-size at 1 AU, and mm-size at 30 AU.

Let us first look at particle growth at 1 AU, in the case where the possibility of bouncing is ignored (taking v_b = v_f = 100 cm s⁻¹ in all simulations). We compare the three possible evolution scenarios laid out above. We find that in all cases the particle mass distribution function rapidly evolves to a steady state (within less than 10,000 years) shown in Figure 4. In each case, the steady-state distribution dn/dm is well approximated by a truncated power law, with the same slope as that of the fragment mass distribution (see Section 2.1). The truncation size, however, depends on the model considered.

Zoom In Zoom Out Reset image size

In Model O, particles grow up to a few centimeters in size. Predicting the maximum mass/size achievable in this model is in fact very easy—one can simply read it at the intersection of the horizontal axis and the $\bar{\Delta }_{ij} = v_f$ contour¹⁵ in Figure 3. Indeed, since smaller particles are much more numerous than larger ones, high-mass-ratio collisions are much more common than equal-mass collisions. The growth of a large particle "i" is thus stalled when collisions with much smaller ones begin to lead to fragmentation instead of sticking. In Model O, where this transition occurs through a Heaviside function, the maximum particle size s_i achievable is then simply given by $\bar{\Delta }_{ij} = v_f$, taking j = 1 to identify collisions with the smallest particles in the simulation.

As expected from the discussion in Section 3.2, since particles are not seen to grow up to a size where radial drift becomes important, the "Maxwellian" (M) model and the "New" (N) model predict very similar outcomes. In both cases, the maximum size achieved is significantly smaller than in Model O. This is quite different from the results of Windmark et al. (2012b), and could be due to their possible use of "total" probabilities instead of the mean fragmentation and sticking probabilities (see Section 3). Furthermore, we also find that the maximum particle size achievable is sensitive to the minimum particle size considered in the simulation. This is illustrated in Figure 4, which shows, in addition to the three models described above, the outcome of a different run using Model N, with exactly the same parameters and same resolution, but simply decreasing the minimum size considered by two orders of magnitude. The maximum particle size achieved in this case is three times smaller than before.

After further investigation, we found that this effect is in fact generic to any model for which the fragmentation probability does not drop exactly to zero below a certain threshold velocity (which is the case in Models M and N, but not in Model O)—and is in fact very easy to understand from a physical point of view. Indeed, as discussed above, collisions are much more frequent the smaller the projectile. In the very crude collisional model considered in this section, it is possible (albeit unlikely) for a μm-sized particle to collide with a much larger one and cause its fragmentation. Even though the fragmentation probability is very low, the sheer number of collisions that take place ensure that fragmentation does happen on a regular basis. Since the number density of particles is a steep function of their mass (or size), the new maximum particle size achievable thus sensitively depends on the smallest particle considered in the numerical integration, as well as the rate at which the fragmentation probability drops to zero when $\bar{\Delta }_{ij}$ drops below v_f. This can be shown analytically, and will be the subject of a separate publication.

We note, however, that this trend disappears when considering models in which high-mass-ratio collisions do not lead to complete fragmentation, see Section 5.3 for details. In this sense, our comments on the effect of the minimum particle size are somewhat academic.

We now consider the effect of bouncing, by setting v_b = 5 cm s⁻¹. Again, we compare the three scenarios discussed above. We set the minimum particle size in the simulation to be s_min = 0.01 cm, for ease of comparison with Windmark et al. (2012b) and with the results of the previous section, bearing in mind the potential effect of that choice on the outcome. In all cases, the simulation is integrated for 30,000 years, and the resulting particle size/mass distribution functions are shown in Figure 5.

Zoom In Zoom Out Reset image size

At that point in time, Model O is still evolving. By construction, it will continue to do so until all particles have grown to a size above which any collisions they may undergo only result in bouncing. Little by little, all small particles disappear from the system. In many ways, this is the worst case scenario in terms of comparison with observations: small particles are entirely depleted, but growth to sizes larger than the bouncing barrier is also impossible. Note that the outcome in this case is strongly dependent on the initial conditions chosen. Suppose for instance that all particles are initialized within the bouncing region—then, the particle size distribution function is in fact in a steady state. Here, we have selected to initialize all particles within the sticking region, to emphasize the effect of the bouncing barrier.

In contrast with Model O, Models M and N have already reached a steady state by 30,000 years. As in the previous section, the two are very similar to one another, as expected from the fact that particles do not appear to grow beyond a size for which radial drift becomes important. As discussed by Windmark et al. (2012b), the maximum particle size achieved is, this time, larger than in Model O since occasional low-velocity collisions are still possible even when the mean relative velocity of the two particles is already beyond the bouncing barrier. More importantly, although fragmentation is rare (since the mean fragmentation probability, at the sizes achieved, is very small), it is nevertheless sufficient to replenish the small-particle population. In other words, models M and N provide more growth than Model O, and maintain a significant population of small grains, but the largest particle size achieved at steady state is still much smaller than what is needed to form planetesimals.

As discussed in Section 1, high-mass-ratio collisions are unlikely to result in the complete destruction of the larger body, as assumed in Sections 5.1 and 5.2. Instead, the smaller particle (dubbed the "projectile" hereafter) is more likely to excavate a small crater from the larger one (dubbed "the target" hereafter) and/or could itself partially or completely stick to the target (Langkowski et al. 2008; Teiser & Wurm 2009; Kothe et al. 2010), thus resulting in mass transfer.

Unfortunately, our understanding of the conditions under which mass transfer occurs, and in particular its dependence on the particle masses, porosity, and velocities, is still in its infancy. For this reason, we use here a very simplified model of the process, in which we assume that the projectile completely sticks to the target if the mass ratio of two particles is larger than a certain critical value ϕ (which we take to be 50, by analogy with Windmark et al. (2012b)), and if the collision would have otherwise led to fragmentation. To implement this numerically, one simply needs to use the following algorithm. (1) Calculate the kernels K_ij and F_ij in the absence of "mass transfer." (2) If the mass ratio of the two particles is greater than ϕ, first redefine K_ij as K_ij := K_ij + F_ij, and then set F_ij := 0.

This prescription differs from that of Windmark et al. (2012b), who assume that only 10% of the mass of the projectile sticks to the target, and the other 90% fragments. We recognize that our own choice is likely to overestimate the actual amount of mass transfer occurring during a collision, but since the actual value of ϕ is very poorly constrained anyway, the error made is within the same order of the approximation. We also show in Section 5.5 that the results presented below still hold qualitatively with lower mass-transfer rates. The advantage of this approach is that the mathematical structure of the effective fragmentation/sticking probabilities is now very easy to visualize. For low mass ratio, they are as depicted in Figure 2. For high mass ratio, $\bar{\epsilon }_{ij}^f = 0$ while $\bar{\epsilon }_{ij}^s$ now looks like the sum of the previously calculated sticking and fragmentation probability functions (i.e., the sum of the two surfaces shown in Figure 2). In other words, $\bar{\epsilon }_{ij}^s$ is close to unity for low collisional velocities and for high collisional velocities, with a substantial "gap" (or perhaps a "moat" may be more appropriate) in between that corresponds to the bouncing region. The latter is delimited by the original bouncing and fragmentation barriers (see Section 4.2).

Mass transfer has fundamental implications for particle growth: as long as $\bar{\epsilon }_{ij}^s$ never strictly drops to zero anywhere, which is the case in Models M and N as discussed in Section 3.2, there is always a possibility for growth across the "moat" in high-mass-ratio collisions. Indeed, while bouncing in this regime may be the most likely outcome of a collision, there is always a possibility of very low and very high velocity events that result in coagulation and growth. Hence, once particles are large enough, they may always continue to grow by sweeping up smaller ones. As long as fragmentation via low-mass-ratio collisions is rare enough, the net effect is growth beyond the original barriers.

Figure 6 presents the results of simulations that include the simple mass-transfer model discussed above. We compare simulations using Models O, M, and N, this time after 10,000 years only. The effect we are interested in is already very clear, and only becomes more pronounced beyond that time. One immediately notices that mass transfer has no effect on Model O. This was already pointed out by Windmark et al. (2012b), and is an obvious result in the light of the discussion above. Indeed, in Model O, both $\bar{\epsilon }_{ij}^s$ and $\bar{\epsilon }_{ij}^f$ strictly drop to zero when $\bar{\Delta }_{ij} \in [v_b,v_f]$. This means that, even for high-mass-ratio collisions, $\bar{\epsilon }_{ij}^s$ strictly drops to zero for $\bar{\Delta }_{ij} > v_b$, preventing any possibility for growth beyond the original bouncing barrier.

Zoom In Zoom Out Reset image size

For Models M and N, however, we see a dramatic increase in particle growth. Furthermore, since the largest particles are now in a regime that is dominated by radial drift rather than turbulence, the predicted evolution of the size distribution functions for models M and N differ significantly. After only 10,000 years, the maximum particle size has increased by more than 4 orders of magnitude for Model M (corresponding to more than 12 orders of magnitude increase in mass, as found by Windmark et al. 2012b). The effect is even more pronounced for Model N, where the maximum particle size is yet another order of magnitude larger. Furthermore, the surface density of larger particles is three to four orders of magnitude larger in model N than in Model M. In short, correctly modeling the difference between regular and stochastic processes when constructing the particle relative velocity distribution function turns out to be highly beneficial to growth. Finally, note that even though fragmentation no longer completely suppresses the growth of large particles, it is still amply sufficient to replenish the small grain population.

Using Model N, we can now study other regions of the disk. Here, we consider the disk at 30 AU, and include both bouncing and mass transfer with the same parameters as in the previous sections: v_b = 5 cm s⁻¹, v_f = 100 cm s⁻¹, and ϕ = 50. As seen in Figure 3, the bouncing threshold now occurs for much smaller particles than at 1 AU. For this reason, we shift the mass-mesh so that s ∈ [10⁻⁴, 10⁵] cm instead, and keep the same resolution.

In order to present complementary yet comparable information to that of the previous sections, we show the actual temporal evolution of the particle distribution function. It is clear from Figure 7 that the system now contains two interacting but distinct particle populations. The small-particle population is distributed as a truncated power law, with a power index dictated by the assumed fragmentation law of the model (here, with dn/dm ∼ m^−1.83, see Section 2.1), and a maximum size around 0.1 mm. At early times, the bouncing barrier induces a pileup just above that size, giving the appearance of a fairly mono-disperse population. However, the pileup gradually disappears later on, as the fragmentation of an increasing number of larger particles replenishes the small-particle population. Meanwhile, the larger particles steadily sweep up the smaller ones, resulting in the gradual emergence of a large-particle population.

Zoom In Zoom Out Reset image size

After about 500,000 years (in this particular simulation), the large-particle population has become very significant indeed. Its distribution, especially at later times, can also be viewed as a truncated power law. The power index is shallower than that of the small particles, with dn/dm ∼ m^−1.2. The maximum particle size within the large-particle population initially grows quickly with time, increasing by three orders of magnitude in the first 100,000 yr. Later on, the growth of the maximum particle size slows down in favor of a steady increase in the number density of the large particles while keeping the same overall shape of the distribution.

It is clear from this simulation that the effects described here have the potential of solving simultaneously the two important puzzles raised by observations of protoplanetary disks—the persistence of a small-particle population for millions of years, and the inferred presence of rather large particles (>cm size) far out in the disk. Furthermore, the emergence of two particle populations with different mass distribution functions, rather than a single continuous power law, is consistent with the observations of Wilner et al. (2005). We will return to these points in Section 6.2. At this point, however, it is time to look into the dependence of the model results on the input parameters.

In all previous simulations, we used the same values for the fragmentation threshold v_f, the bouncing threshold v_b, and the mass ratio above which mass transfer happens. We now vary the latter two to see their effects on evolution of the particle distribution function. The results are presented in Figure 8.

Zoom In Zoom Out Reset image size

We consider here the same disk at 30 AU (very similar results apply at 1 AU). In all simulations, we fix v_f = 100 cm s⁻¹, but vary v_b between 5 cm s⁻¹ and 20 cm s⁻¹, and ϕ between 50 and 200. We evolve the system of Equations (3) for 500,000 years. In all cases we observe the same qualitative evolution of the particle size distribution function into the two particle populations described in the previous section (small particles and large particles). However, we also see that the transfer of mass from the small- to the large-particle population, which in turn controls the large-particle population growth rate, depends quite sensitively on ϕ, and on the width of the bouncing region (i.e., on v_f − v_b).

First, we find that the larger ϕ is, the slower the mass flux from small to large particles. This effect is intuitive: a given target mass m_i will effectively sweep all particles of mass m_i/ϕ or smaller. The total amount of mass available for growth can be evaluated from ∫^m_i/ϕ₀m(dn/dm)dm, and is naturally smaller when ϕ increases (the details of exactly how much smaller depends on the specific shape of the mass distribution function for small particles).

This mass flux is also dependent on the width of the bouncing gap. To understand this, let us consider for instance the tiniest particles in the system (taking i = 1), and look at their sticking probability with other particles around. In a mass-transfer scenario, collisions with anything of mass m_j > ϕm₁ do not lead to fragmentation. For ϕ = 50 and s₁ = 10⁻⁴ cm, this corresponds to s_j > 4 × 10⁻⁴ cm—in other words, very few of the collisions involving s₁ lead to fragmentation. However, the mass-transfer model considered here only leads to sticking for collisions that would otherwise lead to fragmentation. Bouncing remains a barrier to growth, and thus the larger the bouncing region, the stronger this barrier is. Figure 9 illustrates this clearly, by showing the mean sticking probability $\bar{\epsilon }_{1j}^s$ of a particle of index i = 1 colliding with a particle of size s_j. The wider the gap, the deeper the minimum in the function $\bar{\epsilon }_{1j}^s$. This minimum is the growth bottleneck of the system, and thus controls the flux of mass from the small-particle population to the large-particle population.

Zoom In Zoom Out Reset image size

The results of Section 5 suggest that growth beyond the traditional "bouncing" and "fragmentation" barriers, and the maintenance of a small-particle population, merely requires two ingredients: (1) a threshold above which high-mass-ratio collisions lead to sticking rather than to fragmentation and (2) mean sticking and fragmentation probabilities that never strictly drop to zero.

In fact, these ingredients are always naturally expected to arise in any "physically motivated" model for the evolution of the particle size distribution function. Their absence from previous work can only be viewed as unfortunate and misleading oversimplifications of the problem.

Let us begin by considering the case of high-mass-ratio collisions. In the original work of Dullemond & Dominik (2005), for instance, the fragmentation and sticking probabilities were constructed by considering the total kinetic energy of the collision in comparison with the binding energy of the particle. In these kinds of models, high-mass-ratio collisions are much less likely to lead to the fragmentation of the larger body, and sticking becomes a more plausible outcome. In some sense, such prescriptions already incorporate the idea of mass transfer, without any need to define it artificially. Unfortunately, much of the work published between 2005 and 2011 dropped this energy-dependent view of the fragmentation/sticking thresholds in favor of a velocity-dependent one (such as the one used, for the sake of illustration, in this paper). We advocate that future work should return to using the energy-based approach in which the notion of mass transfer naturally emerges, in conjunction with the use of velocity pdfs for the construction of the kernels.

Mass transfer on its own, however, is not sufficient: the second condition is equally important. This was illustrated in Section 5.3, where we showed that even with mass transfer Model O does not lead to growth beyond the fragmentation barrier. In other words, in any model where a given particle size s_i exists for which the mean sticking efficiency $\bar{\epsilon }_{ij}^s =0$ for all possible particle pairs (i, j), growth beyond s_i is naturally impossible. This was a rather common outcome of models that used the "traditional approach" in conjunction with piecewise defined functions $\bar{\epsilon }_{ij}^s$ and $\bar{\epsilon }_{ij}^f$ that drop strictly to zero across a particular velocity threshold, as in the work of Brauer et al. (2008) and the piecewise linear model used by Birnstiel et al. (2010).

In reality, however, the mean sticking probability for a given particle size is never expected to drop to zero entirely. Even when the mean collisional velocity/energy of a particle pair (i, j) is high, low-velocity/energy collisions are always possible. This is explicitly taken into account when using pdfs of relative velocities in the calculation of the mean sticking and fragmentation probabilities (see Section 3.2). Indeed, since the latter can be written as convolution products of the relative velocity pdf and the sticking and fragmentation probabilities of individual collisions, they are always strictly positive. Again, we advocate that such a model should always be used in the future.

As long as it is using the two ingredients listed above, we expect that any local model will yield answers that are similar to the ones shown in Section 5 for Model N, and will reveal the emergence of two populations of particles: a small-particle population constantly replenished by the fragmentation of larger bodies (with a size distribution function controlled by a collisional fragmentation cascade) and a large-particle population that grows by sweeping the smaller particles (with a size distribution controlled by a coagulation/fragmentation balance). In this sense, anyone's preferred prescription for the relative velocity pdfs, for the mass-transfer scenario, for the effect of porosity on the collisions, should yield qualitatively similar results—the latter are certainly not model-dependent.

However, we also showed that the details of the mass flux between the two populations are, unfortunately, very much dependent on the model considered, and within a given choice of model, on the selected parameters. We illustrated this in Section 5.5 by considering a range of parameters within our own toy model, and showed that the predicted surface density of large particles in the system can vary by 5–10 orders of magnitude simply by changing the parameters by a few. This somewhat unsatisfactory result can be understood, as shown earlier, by the fact that the mass flux between the small- and large-particle populations is ultimately controlled by the sticking efficiency bottleneck, as well as the available surface density of small particles that a larger one can "sweep," which are themselves strongly dependent on the parameters.

Furthermore, we showed in Appendix B that the results are, unfortunately, also sensitively dependent on the numerical resolution used to discretize the mass-mesh, if the latter is too small. In particular, a low resolution appears to shorten the growth timescale artificially, by causing an artificial diffusion in mass-space of the particle size distribution function. Interestingly, however, we find that higher resolution runs eventually do reach very similar states as the lower resolution ones, but take longer to do so. In this sense, the emergence of two particle populations is robust. We also find that the typical particle sizes achievable in both populations is very similar in high- and low-resolution runs, demonstrating that this is again a generic property of the system.

As discussed in Section 1, the model presented in this work was, by and large, selected for simplicity rather than physical realism. Comparing results directly with observations will require much more sophistication in the physics included. The particle structures (porosities, composition) should be taken into account. The sticking and fragmentation probabilities for a given collision are more likely to depend, in a complex manner, on the collisional energy and the material strength of the two particles. This should also be taken into account instead of using Equation (11). In addition, the relative velocity pdfs of two particles undergoing turbulent motion is not likely to be a Maxwellian, as we had to assume here. Indeed, the velocity pdf of a single particle interacting with turbulent eddies is known to be better represented by Levy-type distributions than by Gaussian/Maxwellian distributions. The collisional relative velocities pdfs could therefore differ significantly from the one given in Equation (26). Levy-type distributions are defined by much more slowly decaying tails; this is likely to affect the particle size distribution evolution significantly (although quantitatively rather than qualitatively).

In addition, we must remember that the models studied in Section 5 are local models, in the sense that they neglect the radial flux of particles in and out of the domain considered. This approximation is adequate as long as the mass flux in and out of a mass bin caused by radial drift is small compared with that caused by coagulation and/or fragmentation. As we saw, including mass transfer and velocity pdfs is always sufficient, in a local model, to trigger the growth of a large-particle population. In a global disk model, however, this effect will have to be sufficiently strong to overcome the effect of radial drift on the large particle growth. Interestingly, however, we now understand what controls the growth rate beyond the fragmentation barrier, so the observations of the ubiquity of exoplanets, as well as of the presence of a cm-sized particle population at large radii, are still reconcilable with our type of model, and will help place constraints on assumed sticking and fragmentation parameterization and on the particle velocity pdf selected.

Finally, note that taking into account the effect of radial drift in our velocity pdfs enhanced particle growth significantly compared with the similarly local Maxwellian-only type of model proposed by Windmark et al. (2012b), in particular in the cases with mass transfer. This is because, for high-enough mass ratio between the two particles, high-velocity collisions (which are more frequent when radial drift is taken into account) result in sticking. As such, it is much more likely to yield growth to large sizes despite drift than theirs, and should therefore be used preferentially.

Various observational studies of protoplanetary disks have inferred the presence of a wide distribution of particle sizes from submicron-sized dust grains to cm-sized pebbles (Weintraub et al. 1989; Dutrey et al. 1996; Holland et al. 1998; Beckwith 1999; Lada et al. 2000; Wilner et al. 2005; Sicilia-Aguilar et al. 2006; Andrews & Williams 2008; Lommen et al. 2009). Intriguingly, no conclusive evidence has been found for a correlation between the age of primordial disks and the properties of their population of small dust grains (as revealed by silicate features; Kessler-Silacci et al. 2006; Furlan et al. 2009; Oliveira et al. 2010; Ricci et al. 2010). Williams & Cieza (2011) attribute this to a balance between grain growth and destruction on the one hand, and between crystallization and amorphization on the other, concluding that this balance seems to persist throughout the duration of the primordial disk stage (at least in the surface layers of the disk). The consequence for the scenario in which planet formation proceeds via dust growth is the coexistence of particles growing to planetary sizes with a population of small grains that are resupplied by continuous fragmentation. Future observations of planets embedded in a young protoplanetary disk or the confirmation of candidate planets in transition disks, such as T Cha (Huélamo et al. 2011) or LkCa15 (Kraus & Ireland 2012), would support this picture.

In the meantime, other proxies have been sought for evidence of ongoing planetary formation. Unfortunately, the failure of theoretical collisional models to produce results in which large and small particles coexist has, until now, hindered these efforts. The new model presented in this paper paves the way for resolving this problem. Our physically motivated coagulation and fragmentation kernels are able to simultaneously establish a sustained growth of particles and to preserve a μm-sized dust population in a protoplanetary disk. The particle size distribution that naturally emerges is one in which the small- and large-particle populations have notably different power laws: one that is dominated by fragmentation, and one that is dominated by coagulation via sweeping. Since the upper size cutoffs and total mass in each population are strongly dependent on the model considered, observations may help to rule out certain aspects of the parameter space that affect the theoretical results, enabling the model to be constrained.

From a qualitative point of view, the existence of two particle populations far from the central star agrees well with detailed modeling and cm observations of the TW Hydrae disk by Calvet et al. (2002) and Wilner et al. (2005). Furthermore, we expect large and small particles to coexist at different radial positions in the disk, albeit with a strong variation in the maximum size achievable by the larger particles—in other words, one may anticipate strong radial gradients in inferred particle growth beyond mm-size. This result is interesting in the light of the work of Guilloteau et al. (2011), who reported observational evidence for a radial dependence of the grain size in protoplanetary disks. We expect that more detailed future observations probing different parts of a disk will be able to determine whether or not two particle populations do indeed coexist at multiple radii in a disk, as suggested by our model. In any case, in light of our results, observers should always consider the possibility that two populations of particles exist when fitting spectral energy distributions of disks, rather than a single continuous population from small to large sizes.