DustPy: A Python Package for Dust Evolution in Protoplanetary Disks

The case of pure coagulation, i.e., perfect sticking of two particles forming a new larger body, has a number of computational challenges. First, the mass grid of DustPy is logarithmically spaced to cover a large dynamic range from submicrometer-sized particles to large boulders. This has the disadvantage that the resulting mass of two colliding particles m^coll = m_i + m_j will in general not fall exactly onto the mass grid itself, but in between two mass bins. We follow the approach of Brauer et al. (2008) to linearly distribute the newly formed particle between the two adjacent mass bins. Assuming the mass of the particle resulting from a sticking collision of particle m_i and m_j falls in between the two mass bins m_m ≤ m^coll < m_n , we split the coagulation rate linearly between both mass bins:

$$\begin{eqnarray}&&{K}_{{ijk}}=\left\{\begin{array}{l}\epsilon \,\mathrm{if}\,{{\rm{m}}}_{k}={{\rm{m}}}_{m}\\ 1-\epsilon \,\mathrm{if}\,{{\rm{m}}}_{k}={{\rm{m}}}_{n}\\ 0\,\mathrm{else},\end{array}\right.,\end{eqnarray} \tag{ 15 }$$

where is given by

$$\begin{eqnarray}&&\epsilon =\displaystyle \frac{{m}_{n}-{m}^{\mathrm{coll}}}{{m}_{n}-{m}_{m}}.\end{eqnarray} \tag{ 16 }$$

If m^coll = m_m ⇒ = 1 and the entire particle will be distributed into mass bin m_m . Since by definition m_n > m^coll, mass will always be distributed into a mass bin that is larger than the combined mass of the colliding particles, leading to artificial growth. The mass grid, therefore, needs to be fine enough to limit this numerical inaccuracy. This limitation is discussed in Section 4.1 in more detail.

Another challenge is purely computational and is caused by the limited precision of computers. For double-precision numbers, for example, we face the problem that m_i + m_j = m_i , if m_i is more than 15 orders of magnitudes more massive than m_j . This would prevent large particles from growing by sweeping up many small particles. This problem can be solved by rearranging the sums in the discrete Smoluchowski equation. We, again, follow the approach of Brauer et al. (2008). Figure 2 shows a sketch of three types of particle collision that are dealt with separately in this section.

Starting from Equation (14) we can separate the diagonals from the first term on the right-hand side, i.e., the equal-mass particle collisions, and combine it with the Kronecker δ of the negative term:

$$\begin{eqnarray}&&\begin{array}{l}{\dot{n}}_{k}^{{\rm{s}}}=\displaystyle \sum _{i=1}^{k}{R}_{{ii}}^{{\rm{s}}}{n}_{i}^{2}\left({K}_{{iik}}-{\delta }_{{ik}}\right)\\ \,+\,\displaystyle \sum _{i=1}^{k}\displaystyle \sum _{j=1}^{i-1}{K}_{{ijk}}{R}_{{ij}}^{{\rm{s}}}{n}_{i}{n}_{j}-{n}_{k}\displaystyle \sum _{j=1}^{{N}_{m}}{n}_{j}{R}_{{jk}}^{{\rm{s}}}.\end{array}\end{eqnarray} \tag{ 17 }$$

The superscript "s" denotes that these are the source terms and collisions rates for purely sticking collisions. The sum over i does not need to go all the way to N_m but only up to k, because sticking collisions involving particles larger than m_k can never positively contribute to n_k . We now look more closely at the second and third terms on the right-hand side and separate the case k = i from the second term, i.e., those collisions that can be affected by machine-precision errors, when a large particle with mass m_k is sweeping up a small particle m_j such that the resulting mass is in between m_k and m_k+1:

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \sum _{i=1}^{k}\displaystyle \sum _{j=1}^{i-1}{K}_{{ijk}}{R}_{{ij}}^{{\rm{s}}}{n}_{i}{n}_{j}-{n}_{k}\displaystyle \sum _{j=1}^{{N}_{m}}{n}_{j}{R}_{{jk}}^{{\rm{s}}}\\ =\,\displaystyle \sum _{i=1}^{k-1}\displaystyle \sum _{j=1}^{i-1}{K}_{{ijk}}{R}_{{ij}}^{{\rm{s}}}{n}_{i}{n}_{j}\\ \,+\,\displaystyle \sum _{j=1}^{k-1}{K}_{{kjk}}{R}_{{kj}}^{{\rm{s}}}{n}_{k}{n}_{j}-{n}_{k}\displaystyle \sum _{j=1}^{{N}_{m}}{n}_{j}{R}_{{jk}}^{{\rm{s}}}.\end{array}\end{eqnarray} \tag{ 18 }$$

The second term on the right-hand side with K_kjk represents these special collisions for which particles with masses m_k and m_j collide, but still have a positive contribution to n_k (type 1 in Figure 2). The term only describes particle collisions for which m_k + m_j ≤ m_k+1, otherwise the resulting mass of the collision would be too large to positively contribute to n_k . We therefore introduce a number, c, which is defined as the smallest integer for which the condition m_k + m_k+1−c ≤ m_k+1 is fulfilled. In general, c would depend on k. But since the mass grid of DustPy is regular logarithmic, c will be a constant as long as the mass grid does not change. We can now replace the upper boundary of the sum in the second term with k + 1 − c and combine it with the respective negative part of the sum in the third term:

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \sum _{j=1}^{k-1}{K}_{{kjk}}{R}_{{kj}}^{{\rm{s}}}{n}_{k}{n}_{j}-{n}_{k}\displaystyle \sum _{j=1}^{{N}_{m}}{n}_{j}{R}_{{jk}}^{{\rm{s}}}\\ =\,\displaystyle \sum _{j=1}^{k+1-c}{R}_{{jk}}^{{\rm{s}}}{n}_{j}{n}_{k}\left({K}_{{kjk}}-1\right)\\ \,-\,{n}_{k}\displaystyle \sum _{j=k+2-c}^{{N}_{m}}{n}_{j}{R}_{{jk}}^{{\rm{s}}}.\end{array}\end{eqnarray} \tag{ 19 }$$

Please note that c ≥ 2, since the equal-size collisions $\left(j=k\right)$ are already included in the first term in Equation (17). In any case, mass grids with K_kkk ≠ 0 that lead to m_k ≤ m_k + m_k < m_k+1 would have fewer than ${\mathrm{log}}_{2}10\approx 3.3$ mass bins per decade. Simulations should have at least seven mass bins per decade for simple collision models (Ohtsuki et al. 1990), and even more for complex collision models (Drażkowska et al. 2014). Further, note that the collision rates are symmetric, i.e., ${R}_{{jk}}^{{\rm{s}}}={R}_{{kj}}^{{\rm{s}}}$. Collisions of particle m_j with m_k occur at the same rate as collisions of particle m_k with m_j . Since K_kjk in the first term means that m_k ≤ m_k + m_j < m_k+1 we can set m_m = m_k and use Equation (15) to get

$$\begin{eqnarray}&&\begin{array}{l}{K}_{{kjk}}-1=\epsilon -1=\displaystyle \frac{{m}_{k+1}-\left({m}_{k}+{m}_{j}\right)}{{m}_{k+1}-{m}_{k}}-1\\ =\,-\displaystyle \frac{{m}_{j}}{{m}_{k+1}-{m}_{k}}.\end{array}\end{eqnarray} \tag{ 20 }$$

In collisions prone to machine-precision errors, a computer would falsely calculate m_k + m_j = m_k . Already manipulating K_kjk − 1 in advance in these collisions eliminates these errors. Using this for the first term and combining it with the second term, we get

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \sum _{j=1}^{k+1-c}\left({K}_{{kjk}}-1\right){R}_{{jk}}^{{\rm{s}}}{n}_{j}{n}_{k}-{n}_{k}\displaystyle \sum _{j=k+2-c}^{{N}_{m}}{n}_{j}{R}_{{jk}}^{{\rm{s}}}\\ =\,\displaystyle \sum _{j=1}^{{N}_{m}}{D}_{{jk}}{R}_{{jk}}^{{\rm{s}}}{n}_{j}{n}_{k},\end{array}\end{eqnarray} \tag{ 21 }$$

with

$$\begin{eqnarray}&&{D}_{{jk}}=\left\{\begin{array}{l}-\displaystyle \frac{{m}_{j}}{{m}_{k+1}-{m}_{k}}\,\mathrm{if}\,j\leqslant k+1-c\\ -1\,\mathrm{if}\,j\gt k+1-c\end{array}\right..\end{eqnarray} \tag{ 22 }$$

Up to this point the coagulation equation reads

$$\begin{eqnarray}&&\begin{array}{l}{\dot{n}}_{k}^{{\rm{s}}}=\displaystyle \sum _{i=1}^{k}{R}_{{ii}}^{{\rm{s}}}{n}_{i}^{2}\left({K}_{{iik}}-{\delta }_{{ik}}\right)+\displaystyle \sum _{i=1}^{k-1}\displaystyle \sum _{j=1}^{i-1}{K}_{{ijk}}{R}_{{ij}}^{{\rm{s}}}{n}_{i}{n}_{j}\\ \,+\,\displaystyle \sum _{j=1}^{{N}_{m}}{D}_{{jk}}{R}_{{jk}}^{{\rm{s}}}{n}_{j}{n}_{k}.\end{array}\end{eqnarray} \tag{ 23 }$$

Now we look at the second term on the right-hand side for the case i = k − 1:

$$\begin{eqnarray}&&\sum _{j=1}^{k-2}{K}_{k-1,{jk}}\ {R}_{k-1,j}^{{\rm{s}}}\ {n}_{k-1}\ {n}_{j}.\end{eqnarray} \tag{ 24 }$$

These are the other types of collision that can be affected by machine-precision errors. In this case particles with masses m_k−1 and m_j collide and have a positive contribution to n_k . We can distinguish two cases here. In the first case the resulting mass of the colliding particles falls in between m_k−1 ≤ m_k−1 + m_j < m_k (type 2a in Figure 2). This means m_k = m_n in Equation (15) and therefore

$$\begin{eqnarray}&&\begin{array}{l}{K}_{k-1,{jk}}=1-\epsilon =1-\displaystyle \frac{{m}_{k}-\left({m}_{k-1}+{m}_{j}\right)}{{m}_{k}-{m}_{k-1}}\\ =\,\displaystyle \frac{{m}_{j}}{{m}_{k}-{m}_{k-1}}.\end{array}\end{eqnarray} \tag{ 25 }$$

These collisions are identical to type 1 but look at the mass that is distributed into the larger mass bin. In the second case the resulting mass falls in between m_k ≤ m_k−1 + m_j < m_k+1 (type 2b in Figure 2). Here we have m_k = m_m in Equation (15) and therefore

$$\begin{eqnarray}&&\begin{array}{l}{K}_{k-1,{jk}}=\epsilon \\ =\,\displaystyle \frac{{m}_{k+1}-\left({m}_{k-1}+{m}_{j}\right)}{{m}_{k+1}-{m}_{k}}{\rm{\Theta }}\left({m}_{k+1}-{m}_{k-1}-{m}_{j}\right)\\ =\,\left[1-\displaystyle \frac{{m}_{j}+{m}_{k-1}-{m}_{k}}{{m}_{k+1}-{m}_{k}}\right]{\rm{\Theta }}\left({m}_{k+1}-{m}_{k-1}-{m}_{j}\right).\end{array}\end{eqnarray} \tag{ 26 }$$

Cases with m_k−1 + m_j > ≥ m_k+1 do not contribute positively toward n_k . The Heaviside step function takes care of these cases. In both cases either 1 − or itself can be affected by machine-precision errors. It is therefore advisable to manipulate directly these cases in advance, as shown above. We can now split the sum into both cases using the constant c that has been introduced earlier:

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \sum _{j=1}^{k-2}{K}_{k-1,{jk}}\ {R}_{k-1,j}^{{\rm{s}}}\ {n}_{k-1}\ {n}_{j}\\ =\,\displaystyle \sum _{j=1}^{k-2}{E}_{{jk}}\ {R}_{k-1,j}^{{\rm{s}}}\ {n}_{k-1}\ {n}_{j}.\end{array}\end{eqnarray} \tag{ 27 }$$

The matrix E is given by

$$\begin{eqnarray}\begin{array}{rcl}{E}_{{jk}} & = & \left\{\begin{array}{l}\displaystyle \frac{{m}_{j}}{{m}_{k}-{m}_{k-1}}\,{\rm{if}}\,{j}\leqslant {k}-{c}\\ \left[1-\displaystyle \frac{{m}_{j}-{m}_{k-1}-{m}_{k}}{{m}_{k+1}-{m}_{k}}\right]{\rm{\Theta }}\left({m}_{k+1}-{m}_{j}-{m}_{k-1}\right)\,{\rm{if}}\,{j}\gt {k}-{c}.\end{array}\right.\end{array}\end{eqnarray} \tag{ 28 }$$

The full coagulation equation now reads

$$\begin{eqnarray}&&\begin{array}{l}{\dot{n}}_{k}^{{\rm{s}}}=\displaystyle \sum _{i=1}^{k}\left({K}_{{iik}}-{\delta }_{{ik}}\right){R}_{{ii}}^{{\rm{s}}}{n}_{i}^{2}+\displaystyle \sum _{j=1}^{{N}_{m}}{D}_{{jk}}{R}_{{jk}}^{{\rm{s}}}{n}_{j}{n}_{k}\\ \,+\,\displaystyle \sum _{i=1}^{k-2}\displaystyle \sum _{j=1}^{i-1}{K}_{{ijk}}{R}_{{ij}}^{{\rm{s}}}{n}_{i}{n}_{j}+\displaystyle \sum _{j=1}^{k-2}{E}_{{jk}}\ {R}_{k-1,j}^{{\rm{s}}}\ {n}_{k-1}\ {n}_{j}.\end{array}\end{eqnarray} \tag{ 29 }$$

It is useful to bring the equation into a double sum form:

$$\begin{eqnarray}&&\begin{array}{l}{\dot{n}}_{k}^{{\rm{s}}}=\displaystyle \sum _{i=1}^{{N}_{m}}\displaystyle \sum _{j=1}^{{N}_{m}}\left({K}_{{ijk}}-{\delta }_{{ik}}\right){R}_{{ij}}^{{\rm{s}}}{n}_{i}{n}_{j}{\delta }_{{ij}}+\displaystyle \sum _{i=1}^{{N}_{m}}\displaystyle \sum _{j=1}^{{N}_{m}}{D}_{{ij}}{R}_{{ij}}^{{\rm{s}}}{n}_{i}{n}_{j}{\delta }_{{jk}}\\ \,+\,\displaystyle \sum _{i=1}^{{N}_{m}}\displaystyle \sum _{j=1}^{{N}_{m}}{K}_{{ijk}}{R}_{{ij}}^{{\rm{s}}}{n}_{i}{n}_{j}{\rm{\Theta }}\left(k-i-\displaystyle \frac{3}{2}\right){\rm{\Theta }}\left(i-j-\displaystyle \frac{1}{2}\right)\\ \,+\,\displaystyle \sum _{i=1}^{{N}_{m}}\displaystyle \sum _{j=1}^{{N}_{m}}{E}_{j,i+1}{R}_{{ij}}^{{\rm{s}}}{n}_{i}{n}_{j}{\delta }_{i,k-1}{\rm{\Theta }}\left(k-j-\displaystyle \frac{3}{2}\right),\end{array}\end{eqnarray} \tag{ 30 }$$

where the Kronecker δ and the Heaviside step function Θ are used to pick the correct values and ranges for i and j. In that way the coagulation equation can be written with one single double sum. Note that the inner sum has to go up until N_m , because the matrices D_ji and E_jk are not symmetric. Since i and j are integer numbers, the terms of $\tfrac{1}{2}$ and $\tfrac{3}{2}$ in the Heaviside step functions are used to avoid a potentially undefined behavior for ${\rm{\Theta }}\left(0\right)$.

Zoom In Zoom Out Reset image size

From a computational perspective it is beneficial to bring the equation into a symmetrical form to save half of the iterations. DustPy, therefore, solves the following equation:

$$\begin{eqnarray}&&{\dot{n}}_{k}^{{\rm{s}}}=\sum _{i=1}^{{N}_{m}}\sum _{j=1}^{i}{C}_{{ijk}}{R}_{{ij}}^{{\rm{s}}}{n}_{i}{n}_{j},\end{eqnarray} \tag{ 31 }$$

with

$$\begin{eqnarray}&&{C}_{{ijk}}=\left\{\begin{array}{l}{\tilde{C}}_{{ijk}}+{\tilde{C}}_{{jik}}\,\mathrm{if}\,i\geqslant j\\ 0\qquad \qquad \quad \mathrm{else}\end{array}\right.\end{eqnarray} \tag{ 32 }$$

and

$$\begin{eqnarray}&&\begin{array}{l}{\tilde{C}}_{{ijk}}=\displaystyle \frac{1}{2}{K}_{{ijk}}{\delta }_{{ij}}+{D}_{{ij}}{\delta }_{{jk}}\\ \,+\,{K}_{{ijk}}{\rm{\Theta }}\left(k-i-\displaystyle \frac{3}{2}\right){\rm{\Theta }}\left(i-j-\displaystyle \frac{1}{2}\right)\\ \,+\,{E}_{j,i+1}{\delta }_{i,k-1}{\rm{\Theta }}\left(k-j-\displaystyle \frac{3}{2}\right).\end{array}\end{eqnarray} \tag{ 33 }$$

Please note that diagonal entries $\left(i=j\right)$ must not be counted twice. All diagonals need therefore a factor of $\tfrac{1}{2}$. Only the first two terms of Equation (30) contain diagonals. The first term contains only diagonal entries, while the second term contains both diagonal and off-diagonal entries. However, since c ≥ 2, this means for the diagonals of the second term D_kk = −1. This is one particle that gets removed from the distribution in collisions with two equal-size particles. The other particle is removed by δ_ik in the first term of Equation (30). We can therefore omit the δ_ik in the first term and the factor of $\tfrac{1}{2}$ for the second term in Equation (33).

For any given collision of particles with masses m_i and m_j , C_ijk will either have zero, three, or four nonzero elements. C_ijk will always have two positive entries for k = m and k = n (the two mass bins in between which the resulting collisional mass falls) and negative entries for k = i and k = j (the colliding particles), leading to four nonzero entries. In special cases, when i = j, k = i, or k = j, this number can be reduced to three. If m_i + m_j is larger than the largest mass of the mass grid C_ijk will be set to zero to prevent mass loss through the upper boundary of the mass grid.

A peculiar property of C_ijk for any particle collision is

$$\begin{eqnarray}&&\sum _{k=1}^{{N}_{m}}{C}_{{ijk}}=\left\{\begin{array}{l}-1\\ 0\,\mathrm{if}\,{m}_{i}+{m}_{j}\gt {m}_{{N}_{m}}.\end{array}\right.\end{eqnarray} \tag{ 34 }$$

Since the coagulation equation described above works on number densities, the sum of C_ijk over k for any combination of i and j has to be −1 as long as the mass of the colliding particles is within the mass grid. Two particles collide, stick, and form a single larger particle. Therefore, for every sticking collision the total number of particles is reduced by one.

C_ijk only needs to be calculated once in the beginning of the simulation as long as the mass grid does not change. Because there are a maximum of four nonzero elements for any combination of i and j, the coagulation problem is of the order ${ \mathcal O }\left({N}_{m}^{2}\right)$.

Equation (31) can be written in matrix form:

$$\begin{eqnarray}&&\displaystyle \frac{\partial }{\partial t}{\boldsymbol{n}}={{\mathbb{J}}}^{s}\cdot {\boldsymbol{n}},\end{eqnarray} \tag{ 35 }$$

with the sticking Jacobian J^s being defined as

$$\begin{eqnarray}&&{J}_{{ki}}^{{\rm{s}}}=\sum _{j=1}^{i}{C}_{{ijk}}{R}_{{ij}}{n}_{j}.\end{eqnarray} \tag{ 36 }$$

A sketch of the structure of the Jacobian with the contributions of the four terms in the definition of $\tilde{C}$ in Equation (33) is shown in Figure 3 for a mass resolution of seven mass bins per decade. The last column is always empty, since collisions with particles of mass ${m}_{{N}_{m}}$ will always result in a particle exceeding the mass grid. In this setup the element J₂₁ is empty, because it represents collisions involving at least one particle of mass m₁ that have a positive contribution to n₂. However, the mass grid is fine enough, such that m₁ + m₁ > m₃, which means that n₂ cannot be filled from these types of collisions. The first term represents equal particle collisions, the second term contains the ${\mathbb{D}}$ matrix with the negative contributions to the distribution, the fourth term contains the contribution of the ${\mathbb{E}}$ matrix, and the third term contains the remaining collisions.

Zoom In Zoom Out Reset image size

If the relative velocity of the colliding particles exceeds the fragmentation velocity, particles fragment rather than stick and grow. DustPy distinguishes by default two types of fragmentation events: full fragmentation and erosion.

Full fragmentation means that both colliding particles fully fragment, leaving behind a fragment distribution that follows a power law:

$$\begin{eqnarray}&&n\left(m\right){\rm{d}}m\propto {m}^{\gamma }{\rm{d}}m.\end{eqnarray} \tag{ 37 }$$

The exponent γ has to be determined experimentally. DustPy uses by default $\gamma =-\tfrac{11}{6}$, taken from Dohnanyi (1969). Erosion, on the other hand, happens when both colliding particles differ significantly in mass. The smaller projectile particle then fully fragments while chipping off some mass from the larger target particle. The outcome of a erosive collision is a fragment distribution and a slightly less massive remnant target particle. In DustPy the transition between full fragmentation and erosion is by default at a particle mass ratio of 10.

To calculate the contribution of fragmenting collisions, DustPy uses the ${ \mathcal O }\left({N}_{m}^{2}\right)$ algorithm developed by Rafikov et al. (2020). We slightly modified the algorithm to make it strictly mass conserving and to account for the DustPy code units, where the dust quantities are integrated over the mass bin. For this purpose we define a normalized fragment distribution:

$$\begin{eqnarray}&&{\varphi }_{{ij}}=\displaystyle \frac{{m}_{j}^{1+\gamma }}{{\sum }_{j=1}^{i}{m}_{j}^{1+\gamma }}{\rm{\Theta }}\left(i-j+\displaystyle \frac{1}{2}\right)\times \displaystyle \frac{1}{{\rm{g}}\,{\mathrm{cm}}^{3}}.\end{eqnarray} \tag{ 38 }$$

φ_ij is the amount that gets added to n_j from a fragment distribution with a total mass 1 and a largest fragment mass of m_i . The Heaviside step function sets φ_ij to zero if the index j is greater than the largest mass bin of the fragment distribution. The exponent is 1 + γ instead of γ, because the quantity is integrated over the mass bin. As another quantity we define the total mass of fragments that is created in a single particle collision event:

$$\begin{eqnarray}&&{A}_{{ij}}=\left\{\begin{array}{l}{m}_{i}+{m}_{j}\,\mathrm{for}\,\mathrm{full}\,\mathrm{fragmentation}\\ \left(1+\chi \right){m}_{j}\,\mathrm{for}\,\mathrm{erosion}.\end{array}\right.\end{eqnarray} \tag{ 39 }$$

For fully fragmenting collisions the fragment mass is the total mass of the colliding particles. For erosive collisions the projectile particle chips off a fraction of χ of its own mass from the target particle. In DustPy χ = 1 by default. Without loss of generality, m_j is always the mass of the smaller projectile particle. Another quantity that is needed is the largest mass, m_k , of the fragment distribution. The index of the largest fragment is given by

$$\begin{eqnarray}&&{k}_{{ij}}^{\mathrm{lf}}=\left\{\begin{array}{l}i\,\mathrm{for}\,\mathrm{full}\,\mathrm{fragmentation}\\ j\,\mathrm{for}\,\mathrm{erosion}\end{array}\right..\end{eqnarray} \tag{ 40 }$$

In fully fragmenting collisions the fragment distribution goes all the way up to the largest particle. In erosive collisions, the largest particle of the fragment distribution has the mass of the projectile particle. With these quantities we can now sum up the contribution of all collisions to the total fragment distribution multiplied with their individual fragment mass and weighted by their collision rates and store them in a vector, ${A}_{k}^{* }$, at the position of the largest fragment:

$$\begin{eqnarray}&&{A}_{k}^{* }=\sum _{i=1}^{{N}_{m}}\sum _{j=1}^{i}{A}_{{ij}}{R}_{{ij}}^{{\rm{f}}}{n}_{i}{n}_{j}{\delta }_{k,{k}_{{ij}}^{\mathrm{lf}}},\end{eqnarray} \tag{ 41 }$$

where the superscript "f" denotes the collision rates for fragmenting collisions. The contribution from fragments of all collisions into n_k is then given by

$$\begin{eqnarray}&&{\dot{n}}_{k}^{\mathrm{fragments}}=\sum _{i=k}^{{N}_{m}}{A}_{i}^{* }{\varphi }_{{ik}}.\end{eqnarray} \tag{ 42 }$$

In both cases, full fragmentation and erosion, the smaller projectile particle will fully fragment and has to be removed from the particle distribution:

$$\begin{eqnarray}&&{\dot{n}}_{k}^{\mathrm{projectile}}=-\sum _{i=1}^{{N}_{m}}\sum _{j=1}^{i}{R}_{{ij}}^{{\rm{f}}}{n}_{i}{n}_{j}{\delta }_{{jk}}.\end{eqnarray} \tag{ 43 }$$

Similarly, the larger target particle has to be removed in fully fragmenting collisions. In erosive collisions, however, the target particle has to be removed from the distribution and then added, as a remnant particle, at another place in the distribution:

$$\begin{eqnarray}&&\begin{array}{l}{\dot{n}}_{k}^{\mathrm{target}}=-\displaystyle \sum _{i=1}^{{N}_{m}}\displaystyle \sum _{j=1}^{i}{R}_{{ij}}^{{\rm{f}}}{n}_{i}{n}_{j}{\delta }_{{ik}}\\ \,+\,\displaystyle \sum _{i=1}^{{N}_{m}}\displaystyle \sum _{j=1}^{i}{\tilde{H}}_{{ijk}}{R}_{{ij}}^{{\rm{f}}}{n}_{i}{n}_{j}.\end{array}\end{eqnarray} \tag{ 44 }$$

The matrix ${\tilde{H}}_{{ijk}}$ is similar to K_ijk in Equation (15) from the previous section and decides between which two mass bins the remnant particle has to be distributed, but with the mass of the remnant particle instead of the total mass of both collision partners. If the remnant particle has a mass m_i−1 ≤ m_i − χ m_j < m_i , mass would be removed and then added into the target particle's mass bin n_i . If that is the case, a similar manipulation as for the coagulation can be performed to avoid machine-precision errors. Since the mass grid of DustPy is logarithmically spaced and the transition between full fragmentation and erosion is defined by the mass ratio of the colliding particles, we can define a constant p, such that full fragmentation happens if j ≥ i − p:

$$\begin{eqnarray}&&\begin{array}{l}{\dot{n}}_{k}^{\mathrm{target}}=\displaystyle \sum _{i=1}^{{N}_{m}}\displaystyle \sum _{j=1}^{i-p-1}{H}_{{ijk}}{R}_{{ij}}^{{\rm{f}}}{n}_{i}{n}_{j}\\ \,-\,\displaystyle \sum _{i=1}^{{N}_{m}}\displaystyle \sum _{j=i-p}^{i}{R}_{{ij}}^{{\rm{f}}}{n}_{i}{n}_{j}{\delta }_{{ik}}.\end{array}\end{eqnarray} \tag{ 45 }$$

The first term now holds both the positive and the negative contribution for the target particle from erosive collisions. We can distinguish two cases. In the first case the positive and negative contributions can be combined, since they both affect the same mass bin:

1. m_i−1 ≤ m_i − χ m_j < m_i :

$$\begin{eqnarray}&&\begin{array}{l}{H}_{{ij},i-1}=\epsilon =\displaystyle \frac{{m}_{i}-\left({m}_{i}-\chi {m}_{j}\right)}{{m}_{i}-{m}_{i-1}}=\displaystyle \frac{\chi {m}_{j}}{{m}_{i}-{m}_{i-1}}\\ {H}_{{iji}}=1-\epsilon -1=-\epsilon =-\displaystyle \frac{\chi {m}_{j}}{{m}_{i}-{m}_{i-1}}.\end{array}\end{eqnarray} \tag{ 46 }$$

2. m_m ≤ m_i − χ m_j < m_n < m_i :

$$\begin{eqnarray}&&\begin{array}{l}{H}_{{ij},n-1}=\epsilon =\displaystyle \frac{{m}_{n}-\left({m}_{i}-\chi {m}_{j}\right)}{{m}_{n}-{m}_{m}}\\ {H}_{{ij},n}=1-\epsilon =1-\displaystyle \frac{{m}_{n}-\left({m}_{i}-\chi {m}_{j}\right)}{{m}_{n}-{m}_{m}}\\ {H}_{{iji}}=-1.\end{array}\end{eqnarray} \tag{ 47 }$$

The full equation for fragmentation and erosion is therefore the sum of all three contributions:

$$\begin{eqnarray}&&{\dot{n}}_{k}^{{\rm{f}}}={\dot{n}}_{k}^{\mathrm{fragments}}+{\dot{n}}_{k}^{\mathrm{projectile}}+{\dot{n}}_{k}^{\mathrm{target}}.\end{eqnarray} \tag{ 48 }$$

The source terms of fragmentation can also be written in matrix form. A sketch of the fragmentation Jacobian is shown in Figure 4 for a model where every collision leads to a fragmentation event. The fragmentation Jacobian is a simple upper triangular matrix.

Zoom In Zoom Out Reset image size

This fragmentation and erosion prescription is of the order ${ \mathcal O }\left({N}_{m}^{2}\right)$. The full dust growth equation can then be simply combined to

$$\begin{eqnarray}&&\displaystyle \frac{\partial }{\partial t}{n}_{k}={\dot{n}}_{k}^{{\rm{s}}}+{\dot{n}}_{k}^{{\rm{f}}},\end{eqnarray} \tag{ 49 }$$

with the dust Jacobian being the sum of the sticking and fragmentation Jacobians:

$$\begin{eqnarray}&&{{\mathbb{J}}}^{\mathrm{coag}}={{\mathbb{J}}}^{{\rm{s}}}+{{\mathbb{J}}}^{{\rm{f}}}.\end{eqnarray} \tag{ 50 }$$

The collision rates for sticking/fragmenting collisions are the product of the geometrical cross section, the relative velocities of the particles, and the sticking/fragmentation probabilities, respectively:

$$\begin{eqnarray}&&{R}^{{\rm{s}}/{\rm{f}}}=\frac{1}{1+{{\delta }}_{ij}}{{\sigma }}_{{\rm{g}}{\rm{e}}{\rm{o}}}{v}_{ij}^{{\rm{r}}{\rm{e}}{\rm{l}}}{p}^{{\rm{s}}/{\rm{f}}},\end{eqnarray} \tag{ 51 }$$

with ${\sigma }_{\mathrm{geo}}=\pi {\left({a}_{i}+{a}_{j}\right)}^{2}$. However, please note that the DustPy quantities are vertically integrated, which has been ignored so far. Vertical integration of the Smoluchowski equation introduces a correction factor, which will be incorporated into the collision rates. This is discussed in Section 3.2.6 in more detail. Further note, that a population of N equal-sized particles has $\frac{1}{2}N(N-1)$ possible collisions amongst each other, which is $\approx \frac{{N}^{2}}{2}$ in the limit of large N. Equal particle collisions, therefore, need a factor of $\frac{1}{2}$, which is incorporated via the ${\delta }_{ij}$ in the collision rates in Equation (51).

DustPy considers by default five different sources of relative velocities between dust particles: Brownian motion, radial and azimuthal drift, vertical settling, and turbulence.

Figure 5 shows all five contributions to the relative velocities in an example simulation of the default DustPy model at a distance 1 au from the star. Brownian motion is especially important as a driver of initial dust growth, when the particles are rather small.

Zoom In Zoom Out Reset image size

Brownian motion. The relative velocities due to Brownian motion are given by

$$\begin{eqnarray}&&{v}_{{ij}}^{\mathrm{rel},\mathrm{brown}}=\min \left(\sqrt{\displaystyle \frac{8{k}_{{\rm{B}}}T\left({m}_{i}+{m}_{j}\right)}{\pi {m}_{i}{m}_{j}}},{c}_{{\rm{s}}}\right).\end{eqnarray} \tag{ 52 }$$

Since this formula is diverging for very small particle masses, the relative velocities are limited to the sound speed c_s. However, one should note that for very small particles and high temperatures, the relative velocities due to Brownian motion can easily exceed typical values for the fragmentation velocity. In the simple collision model that is used by default in DustPy, there is no distinction on particle size when deciding between sticking and fragmentation. Even though these small particles would in reality still stick (or bounce) at these velocities (Chokshi et al. 1993; Blum & Wurm 2008), DustPy would treat those collisions as fragmentation events.

Azimuthal drift. Since dust particles of different sizes have different degrees of sub-Keplerian motion, this leads to a relative velocity in azimuthal direction, which is given by

$$\begin{eqnarray}&&{v}_{{ij}}^{\mathrm{rel},\mathrm{azi}}=\left|{v}_{\mathrm{drift}}^{\max }\left(\displaystyle \frac{1}{1+{\mathrm{St}}_{i}^{2}}-\displaystyle \frac{1}{1+{\mathrm{St}}_{j}^{2}}\right)\right|.\end{eqnarray} \tag{ 53 }$$

Dust particles of the same Stokes number do not experience any relative velocity due to azimuthal drift, because they drift at the same speed.

Radial drift. Dust particles of different sizes have different radial drift speeds. This induces relative velocities between dust particles. They are given by

$$\begin{eqnarray}&&{v}_{{ij}}^{\mathrm{rel},\mathrm{rad}}=\left|{v}_{{\rm{d}},i}-{v}_{{\rm{d}},j}\right|,\end{eqnarray} \tag{ 54 }$$

with the radial dust velocities from Equation (8).

Vertical settling. Dust particles of different sizes settle with different velocities toward the midplane. DustPy uses the descriptions of Dullemond & Dominik (2004) and Birnstiel et al. (2010) to account for this effect:

$$\begin{eqnarray}&&{v}_{{ij}}^{\mathrm{rel},\mathrm{sett}}=\left|{h}_{i}\min \left({\mathrm{St}}_{i},\displaystyle \frac{1}{2}\right)-{h}_{j}\min \left({\mathrm{St}}_{j},\displaystyle \frac{1}{2}\right)\right|{{\rm{\Omega }}}_{{\rm{K}}},\end{eqnarray} \tag{ 55 }$$

where h_i is the dust scale height, given by Dubrulle et al. (1995) as

$$\begin{eqnarray}&&{h}_{i}={H}_{{\rm{P}}}\displaystyle \frac{{\delta }_{z}}{{\delta }_{z}+{\mathrm{St}}_{i}}.\end{eqnarray} \tag{ 56 }$$

δ_z is the vertical settling parameter similar to the turbulent α parameter.

Turbulent motion. To calculate the relative velocities due to turbulent motion we follow the prescription of Ormel & Cuzzi (2007). Instead of the turbulent α parameter we use the δ_t parameter, which works in an identical way but allows us to disentangle both effects.

The total relative velocity is then the quadratic sum of all contributions:

$$\begin{eqnarray}&&{v}_{{ij}}^{\mathrm{rel}}=\sqrt{\sum _{k}{\left({v}_{{ij}}^{\mathrm{rel},k}\right)}^{2}}.\end{eqnarray} \tag{ 57 }$$

If the relative collision velocity exceeds the fragmentation velocity, particles start to fragment instead of growing to larger bodies. In the default DustPy model the fragmentation velocity is set to v_frag = 1 ms⁻¹.

Different particles, however, do not collide with a single relative velocity as described in Section 3.2.4. Instead, the relative velocities follow the Maxwell–Boltzmann distribution:

$$\begin{eqnarray}&&{ \mathcal P }\left({\rm{\Delta }}v;{v}_{\mathrm{rms}}\right)=\sqrt{\displaystyle \frac{54}{\pi }}\displaystyle \frac{{\rm{\Delta }}{v}^{2}}{{v}_{\mathrm{rms}}^{3}}\exp \left[-\displaystyle \frac{3}{2}{\left(\displaystyle \frac{{\rm{\Delta }}v}{{v}_{\mathrm{rms}}}\right)}^{2}\right],\end{eqnarray} \tag{ 58 }$$

with the rms velocity v_rms, which DustPy assumes to be the single velocity derived in Section 3.2.4.

The collision rate for fragmenting collisions described in Section 3.2.3 would therefore be an integral over all possible relative velocities in the Maxwell–Boltzmann distribution, which are above the fragmentation velocity:

$$\begin{eqnarray}&&\begin{array}{l}{R}_{{ij}}^{{\rm{f}}}={\displaystyle \int }_{{v}_{\mathrm{frag}}}^{\infty }{\sigma }_{\mathrm{geo}}{\rm{\Delta }}v{ \mathcal P }\left({\rm{\Delta }}v,{v}_{\mathrm{rms}}\right){\rm{d}}{\rm{\Delta }}v.\end{array}\end{eqnarray} \tag{ 59 }$$

This integral has an analytical solution, and therefore the fragmentation probability in Equation (51) can be written as

$$\begin{eqnarray}&&\begin{array}{l}{p}_{{ij}}^{{\rm{f}}}=\displaystyle \frac{{R}_{{ij}}^{{\rm{f}}}}{{\sigma }_{\mathrm{geo}}{\rm{\Delta }}\bar{v}}\\ =\,\left(\displaystyle \frac{3}{2}{\left(\displaystyle \frac{{v}_{\mathrm{frag}}}{{v}_{{ij}}^{\mathrm{rel}}}\right)}^{2}+1\right)\exp \left[-\displaystyle \frac{3}{2}{\left(\displaystyle \frac{{v}_{\mathrm{frag}}}{{v}_{{ij}}^{\mathrm{rel}}}\right)}^{2}\right],\end{array}\end{eqnarray} \tag{ 60 }$$

with the mean velocity of the Maxwell–Boltzmann distribution ${\rm{\Delta }}\bar{v}=\sqrt{\tfrac{8\pi }{3}}{v}_{{ij}}^{\mathrm{rel}}$. In that way, it is sufficient to only calculate one relative velocity per particle collision while accounting for the velocity distribution in the fragmentation probability.

The sticking probability is then given by

$$\begin{eqnarray}&&{p}_{{ij}}^{{\rm{s}}}=1-{p}_{{ij}}^{{\rm{f}}}.\end{eqnarray} \tag{ 61 }$$

Bouncing, i.e., neither sticking nor fragmentation, is not included in the default model of DustPy, but can be easily implemented if ${p}_{{ij}}^{{\rm{s}}}+{p}_{{ij}}^{{\rm{f}}}\lt 1$. Figure 6 shows the Maxwell–Boltzmann distribution in the case of an rms velocity equal to the fragmentation velocity and the sticking/fragmentation probabilities for different relative velocities.

Zoom In Zoom Out Reset image size

The same approach of a velocity distribution was used by Windmark et al. (2012), while Birnstiel et al. (2010) originally used a simple formula for the transition between sticking and fragmentation at the fragmentation velocity.

So far we have ignored the vertical dimension of the disk. The Smoluchowski equation discussed in previous chapters works on volume densities. DustPy, on the other hand, only has one spatial dimension, the distance from the star r. Here we describe the method of Birnstiel et al. (2010) in vertically integrating the Smoluchowski equation. We assume that the vertical dust distribution can be described with a Gaussian:

$$\begin{eqnarray}&&{n}_{k}\left(r,z\right)={n}_{k}\left(r,0\right)\cdot \exp \left[-\displaystyle \frac{1}{2}{\left(\displaystyle \frac{z}{{h}_{k}\left(r\right)}\right)}^{2}\right],\end{eqnarray} \tag{ 62 }$$

with the dust scale height h_k given by Equation (56). Integration over z leads to

$$\begin{eqnarray}&&\begin{array}{l}{N}_{k}\left(r\right)={\displaystyle \int }_{-\infty }^{\infty }{n}_{k}\left(r,z\right){\rm{d}}z\\ =\,\sqrt{2\pi }{h}_{k}\left(r\right){n}_{k}\left(r,0\right).\end{array}\end{eqnarray} \tag{ 63 }$$

Every single term in the collisional dust source terms introduced in Sections 3.2.1 and 3.2.2 contain the product of two densities with the collision rates R_ij n_i n_j . Integrating this term over z and assuming to zeroth order that the collision rates do not depend on z leads to

$$\begin{eqnarray}&&\begin{array}{l}{R}_{{ij}}{\displaystyle \int }_{-\infty }^{\infty }{n}_{i}\left(r,z\right){n}_{j}\left(r,z\right){\rm{d}}z\\ =\,{R}_{{ij}}{n}_{i}{n}_{j}{\displaystyle \int }_{-\infty }^{\infty }\exp \left[-\displaystyle \frac{1}{2}{z}^{2}\left(\displaystyle \frac{1}{{h}_{i}^{2}}+\displaystyle \frac{1}{{h}_{j}^{2}}\right)\right]{\rm{d}}z\\ =\,\displaystyle \frac{{R}_{{ij}}}{\sqrt{2\pi \left({h}_{i}^{2}+{h}_{j}^{2}\right)}}{N}_{i}{N}_{j}\equiv {\tilde{R}}_{{ij}}{N}_{i}{N}_{j}.\end{array}\end{eqnarray} \tag{ 64 }$$

Vertically integrating the Smoluchowski equation can therefore be achieved by simply replacing the midplane volume densities n_i with the vertically integrated number densities N_i and by multiplying the collision rates with a correction factor. The quantity ${\tilde{R}}_{{ij}}$ is stored as kernel in DustPy.

DustPy does not store the vertically integrated number densities but uses dust surface densities, which are simply given by

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{{\rm{d}},i}={m}_{i}{N}_{i}.\end{eqnarray} \tag{ 65 }$$

Note that by using surface densities instead of number densities, this introduces further mass factors in quantities like C_ijk , ${A}_{i}^{* }$, or H_ijk introduced above. This is straightforward to do, but increases the complexity of the equations presented here. We therefore refer the interested reader to the software for details on the implementation.

Further note, that only the dust densities have been vertically integrated. Other quantities, like the relative velocities needed for ${\tilde{R}}_{{ij}}$, are still calculated in the midplane. This is a valid simplification, since most of the mass will settle anyway rather quickly toward the midplane.

The solution of Equation (71) with the constant kernel is given by

$$\begin{eqnarray}&&\begin{array}{l}n\left(m,t\right)=\displaystyle \frac{{N}_{0}}{{m}_{0}}{\left(\displaystyle \frac{2}{\alpha {N}_{0}t}\right)}^{2}\\ \,\times \,\exp \left[\displaystyle \frac{2}{\alpha {N}_{0}t}\left(1-\displaystyle \frac{m}{{m}_{0}}\right)\right],\end{array}\end{eqnarray} \tag{ 72 }$$

where m₀ is the smallest possible mass and N₀ the initial total number density of particles:

$$\begin{eqnarray}&&{N}_{0}={\int }_{0}^{\infty }n\left(m,0\right){\rm{d}}m.\end{eqnarray} \tag{ 73 }$$

Note that the DustPy code units are the number densities integrated over the mass bin. We therefore initialize the simulation by setting the first mass bin to N₀ and all other mass bins to zero.

The result is shown in Figure 8 (left panel) for a simulation with the default mass resolution of DustPy with dotted lines and for a simulation with a four times higher mass resolution with solid lines for α = 1. The analytical solutions given by Equation (72) are plotted with black solid lines. There is a slight deviation visible at the upper mass tail of the distribution in the default resolution run. An explanation for this is given at the end of the next section.

Zoom In Zoom Out Reset image size

The solution of the linear kernel is given by

$$\begin{eqnarray}&&\begin{array}{l}n\left(m,t\right)=\displaystyle \frac{{N}_{0}}{2\sqrt{\pi }{m}_{0}^{2}}\displaystyle \frac{g}{{\left(1-g\right)}^{0.75}}\\ \,\times \,\exp \left[-\displaystyle \frac{m}{{m}_{0}}{\left(1-\sqrt{1-g}\right)}^{2}\right],\end{array}\end{eqnarray} \tag{ 74 }$$

with N₀ being the initial total number density as for the constant kernel and

$$\begin{eqnarray}&&g=\exp \left[-\alpha {N}_{0}{m}_{0}t\right].\end{eqnarray} \tag{ 75 }$$

Initially, only the first mass bin was filled with N₀ while all other mass bins were set to zero.

The result is shown in Figure 8 (right panel) for a simulation with the default mass resolution of DustPy with dotted lines and for a simulation with a four times higher mass resolution with solid lines for α = 1. The analytical solutions given by Equation (74) are plotted with black solid lines. As for the constant kernel there is a deviation at the upper mass end of the distribution that is worse the lower the mass resolution is.

The reason for this is the algorithm described in Section 3.2.1. Since the mass grid is logarithmically spaced, the combined mass of both colliding particles in a sticking collision will not directly fall onto the mass grid itself, but has to be distributed between the adjacent mass bins as described in Equation (15). This leads to artificial growth, because material will be added to a bin that is more massive than the combined mass of the collision partners. This causes the simulations to be more massive at the higher mass end and—due to mass conservation—less massive at the lower mass end compared to the analytic solutions. The situation is worse for the linear kernel, because the kernel is proportional to the colliding mass itself. An overestimation of mass will overestimate the kernel itself.

But since the computational time is highly sensitive to the number of mass bins, one has to find a compromise between accuracy and execution time. For the simple collision model in the default DustPy simulation, the default mass resolution should be sufficient, since growth will be eventually halted by the fragmentation or by the drift barrier. Only the growth timescale might be slightly underestimated. For more complex collision models, including, for example, mass transfer, a higher mass resolution might be crucial. For more details on this we refer to Drażkowska et al. (2014), which performed mass resolution tests for more complex collision models. In any case, we advise to always run selected simulations with a higher mass resolution to verify that the default mass resolution was sufficient.

Figure 9 shows the relative errors in mass for the two benchmark models of the constant and linear kernels with the default and the high-resolution runs. In all cases the errors are very close to machine-precision levels for double-precision floating point numbers. The coagulation algorithm of DustPy is therefore mass conserving.

Zoom In Zoom Out Reset image size