0. Algorithm

TriPoDPy runs a parametric dust evolution model initially proposed in Pfeil et al (2024).

In the following, we will briefly introduce the algorithm, which represents the dust size distribution at each radius using merely three parameters. Namely, these are the dust surface densities of the small and large dust population \(\Sigma_0\) and \(\Sigma_1\) as well as the maximum particle size \(s_{max}\). For a detailed description of the algorithm please also have a look at Pfeil et al (2024).

Derivation

First, we will derive an expression for the dust surface density in a particle size interval. Using the particle size-dependency of the particle density \(n\left(a\right)\propto a^{\xi}\), we arrive at \(\sigma\left( a\right) \propto n\left(a\right) m_p \propto a^{\xi+3}\) for the surface density. Then, normalization to the total dust surface density yields:

\begin{align*} \Sigma_{tot} &= \int\limits_{s_{min}}^{s_{max}} \sigma\left(a\right) \mathrm{d}a = C \int\limits_{s_{min}}^{s_{max}} a^{\xi+3} \mathrm{d}a = \left\{ \begin{aligned} \frac{C}{\xi+4} \left(s_{max}^{\xi+4} - s_{min}^{\xi+4}\right) \textrm{for }\xi\neq -4\\ C \left[\log(s_{max}/s_{min})\right] \textrm{for }\xi = -4 \end{aligned} \right. \\ &\Rightarrow \sigma\left(a\right) = \left\{ \begin{aligned} \frac{\Sigma_{tot}(\xi+4)}{s_{max}^{\xi+4} - s_{min}^{\xi+4}} a^{\xi+3}\textrm{for }\xi\neq -4 \\ \frac{\Sigma_{tot}}{\log(s_{max}/s_{min})}a^{\xi+3} \textrm{for }\xi = -4 \end{aligned} \right. \end{align*}

\begin{align*} \Rightarrow \Sigma_{[a_0, a_1]} = \int\limits_{a_0}^{a_1} \sigma\left(a\right) \mathrm{d}a = \left\{ \begin{aligned} \Sigma_{tot} \frac{a_1^{\xi+4}-a_0^{\xi+4}}{s_{max}^{\xi+4}-s_{min}^{\xi+4}}\textrm{for }\xi\neq -4 \\ \Sigma_{tot} \frac{\log(a_1/a_0)}{\log(s_{max}/s_{min})} \textrm{for }\xi = -4 \end{aligned} \right. \end{align*}

Given the minimum and maximum particle size of the density distribution \(s_{min}\) and \(s_{max}\), we can calculate the intermediate particle size on a logarithmic scale (geometric mean) via \(s_{int}=\sqrt{s_{min}s_{max}}\). This motivates a split of the density distribution into two populations and their respective dust surface densities \(\Sigma_0\) (small population) and \(\Sigma_1\) (large population). These can be determined by:

\begin{align*} \Sigma_0 = \int\limits_{s_{min}}^{s_{int}} \sigma\left(a\right) \mathrm{d}a = \left\{ \begin{aligned} \Sigma_{tot} \frac{s_{int}^{\xi+4}-s_{min}^{\xi+4}}{s_{max}^{\xi+4}-s_{min}^{\xi+4}}\textrm{for }\xi\neq -4 \\ \Sigma_{tot} \frac{\log(s_{int}/s_{min})}{\log(s_{max}/s_{min})} \textrm{for }\xi = -4 \end{aligned} \right. \end{align*}

and

\begin{align*} \Sigma_1 = \int\limits_{s_{int}}^{s_{max}} \sigma\left(a\right) \mathrm{d}a = \left\{ \begin{aligned} \Sigma_{tot} \frac{s_{max}^{\xi+4}-s_{int}^{\xi+4}}{s_{max}^{\xi+4}-s_{min}^{\xi+4}}\textrm{for }\xi\neq -4 \\ \Sigma_{tot} \frac{\log(s_{max}/s_{int})}{\log(s_{max}/s_{min})} \textrm{for }\xi = -4 \end{aligned} \right. \end{align*}

We can now derive the relationship between the characteristic particle sizes of the density distribution, the surface densities of the two dust populations and the exponent of the density distribution

\begin{align*} \log\left(\frac{\Sigma_1}{\Sigma_0}\right) &= \log\left(\frac{s_{max}^{\xi+4}-s_{int}^{\xi+4}}{s_{int}^{\xi+4}-s_{min}^{\xi+4}}\right) = \left(\xi+4\right) \log\left(\frac{s_{max}}{s_{int}}\right) + \log \left[\frac{1-\left(\frac{s_{int}}{s_{max}}\right)^{\xi+4}}{1-\left(\frac{s_{min}}{s_{int}}\right)^{\xi+4}}\right] \\ &= \left(\xi+4\right) \log\left(\frac{s_{max}}{s_{int}}\right) + \log \left[\frac{1-\sqrt{\frac{s_{min}}{s_{max}}}^{\xi+4}}{1-\sqrt{\frac{s_{min}}{s_{max}}}^{\xi+4}}\right] = \left(\xi+4\right) \log\left(\frac{s_{max}}{s_{int}}\right) \\ &\Rightarrow \xi = \frac{\log\left(\Sigma_1/\Sigma_0\right)}{\log\left(s_{max}/s_{int}\right)} - 4 \end{align*}

Evolution

With the parametric model of how to describe the particle size distrbition we now have to describe how it evolves due to its interactions with the gas disk and due to mutual collisions. To model the growth of the maximal particle size and its stagnation at the fragmetation barrier \(s_{max}\) evolves according to the following equation: \begin{align*} \dot{s}_{max} = \frac{\Sigma_{1} \Delta v_{\max}}{\rho_{m} \sqrt{2 \pi} H_{1}}\left(\frac{1 - \left( \frac{v_{\mathrm{frag}}}{\Delta v_{\max}} \right)^{s}}{1 + \left( \frac{v_{\mathrm{frag}}}{\Delta v_{\max}}\right)^{s}}\right), \end{align*}

where \(\Delta v_{max} = \Delta v_{0.4 s_{max}, s_{max}}\), is the relative velocity between paricles of size \(s_{max}\) and pativles that are 0.4 times thier size. \(s>0\) is an number derived from calibration (\(s=3\)). The relative velocities between the particles take into acount contributions form the turbulence, settling, relative drift motion and brownian motion as decribed in Pfeil et al (2024).

To calculate the radial gas transport, we use the average particle sizes for the small and large dust population: \begin{align*} a_0 = \frac{\int\limits_{s_{min}}^{s_{int}}a \sigma\left(a\right) \mathrm{d}a}{\int\limits_{s_{min}}^{s_{int}} \sigma\left(a\right) \mathrm{d}a} = \frac{\int\limits_{s_{min}}^{s_{int}}a^{\xi+4} \mathrm{d}a}{\int\limits_{s_{min}}^{s_{int}} a^{\xi+3} \mathrm{d}a} = \left\{ \begin{aligned} \frac{s_{int}s_{min}}{s_{int}-s_{min}} \log\left(s_{int}/s_{min}\right) \textrm{ for }\xi = -5 \\ \frac{s_{int}-s_{min}}{\log\left(s_{int}/s_{min}\right)} \textrm{ for }\xi = -4 \\ \frac{\xi+4}{\xi+5} \frac{s_{int}^{\xi+5}-s_{min}^{\xi+5}}{s_{int}^{\xi+4}-s_{min}^{\xi+4}} \textrm{ for }\xi \neq -4, -5\end{aligned} \right. \end{align*}

\begin{align*} a_1 = \frac{\int\limits_{s_{int}}^{s_{max}}a \sigma\left(a\right) \mathrm{d}a}{\int\limits_{s_{int}}^{s_{max}} \sigma\left(a\right) \mathrm{d}a} = \frac{\int\limits_{s_{int}}^{s_{max}}a^{\xi+4} \mathrm{d}a}{\int\limits_{s_{int}}^{s_{max}} a^{\xi+3} \mathrm{d}a} = \left\{ \begin{aligned} \frac{s_{max}s_{int}}{s_{max}-s_{int}} \log\left(s_{max}/s_{int}\right) \textrm{ for }\xi = -5 \\ \frac{s_{max}-s_{int}}{\log\left(s_{max}/s_{int}\right)} \textrm{ for }\xi = -4 \\ \frac{\xi+4}{\xi+5} \frac{s_{max}^{\xi+5}-s_{int}^{\xi+5}}{s_{max}^{\xi+4}-s_{int}^{\xi+4}} \textrm{ for }\xi \neq -4, -5\end{aligned} \right. \end{align*}

The solution to the transport of our surface densitys is then given by the advection-diffusion equation using the velocities of \(a_0\) and \(a_1\) respectively: \begin{align*} \frac{\partial\Sigma_i }{\partial t} + \vec{\nabla} \left[ \Sigma_i \left(\vec{v}_g+\vec{v}_{dr,i}\right)\right] = 0 \end{align*}

Furthermore, we implement the drift limit by co-advection of the maximum particle size with the larger dust population described by the advection diffusion equation \begin{align*} \frac{\partial\left(\Sigma_1 s_{max}\right)}{\partial t} + \vec{\nabla} \left[ \Sigma_1 s_{max} \left(\vec{v}_g+\vec{v}_{dr,1}\right)\right] = \dot{s}_{max} \end{align*}

To model the collisions of particles we need to exchange mass between the small and large population due to sweep-up of small particles and fragmentation of large particles. The core idea is that the large grains shift mass to the small grain on their own collisional timescale due to fragmentations, adversly the large grains sweep up mass from the small gains on the small-large collison timescale. As the dust density is given by \(\rho_d(z)=\frac{\Sigma_d}{\sqrt{2\pi}H_d} \exp\left(-\frac{z^2}{2 H_d^2}\right)\), these two processes can be described via

\begin{align*} \dot{\Sigma}_{0 \rightarrow 1} &= \frac{\Sigma_0 \Sigma_1 \sigma_{01} \Delta v_{01}}{m_1 2 \pi H_0 H_1} \int\limits_{-\infty}^\infty \exp\left[-\frac{z^2}{2}\frac{H_0^2 + H_1^2}{H_0^2 H_1^2}\right] \mathrm{d}z = \frac{\Sigma_0 \Sigma_1 \sigma_{01} \Delta v_{01}}{m_1 \sqrt{2 \pi \left(H_0^2 + H_1^2\right)}} \\ \dot{\Sigma}_{1 \rightarrow 0} &= \frac{\Sigma_1^2 \sigma_{11} \Delta v_{11}}{m_1 2 \pi H_1^2} \mathcal{F} \int \limits_{-\infty}^\infty \exp \left[-\frac{z^2}{H_1^2}\right]\mathrm{d}z = \frac{\Sigma_1^2 \sigma_{11} \Delta v_{11}}{m_1 2 \sqrt{\pi} H_1} \mathcal{F} \end{align*},

where \(\sigma_{mn}\) and \(\Delta v_{mn}\) denote the collision cross-section and relative velocity between the populations \(m\) and \(n\). The factor \(\mathcal{F}\) modifies the relative effectiveness between sweep-up and fragmentation. Equating \(\dot{\Sigma}_{d, 0 \rightarrow 1}\) and \(\dot{\Sigma}_{d, 1 \rightarrow 0}\) in the steady state yields \begin{align*} \frac{\Sigma_1}{\Sigma_0} = \sqrt{\frac{2H_1^2}{H_0^2 + H_1^2}} \frac{\sigma_{01}\Delta v_{01}}{\sigma_{11}\Delta v_{11}} \mathcal{F}^{-1}. \end{align*}

Using \(\xi = \frac{\log\left(\Sigma_1/\Sigma_0\right)}{\log\left(s_{max}/s_{int}\right)} - 4\) from above, we arrive at \begin{align*} \mathcal{F} = \sqrt{\frac{2H_1^2}{H_0^2 + H_1^2}} \frac{\sigma_{01}\Delta v_{01}}{\sigma_{11}\Delta v_{11}}\left(\frac{s_{max}}{s_{min}}\right)^{-\left(\xi'+4\right)}, \end{align*} where we use \(\xi' = \Pi_{frag}\xi_{frag} + \Pi_{stick}\xi_{stick}\) to steer the model towards the desired density distribution exponents in the respective regimes. To this end, we require \(\xi_{frag} = -3.5\) and \(\xi_{stick} = -3.0\). In order to smoothly transition between the fragmentation and drift regime, we assign a transition function with \begin{align*} f\left(\Delta v_{11} / v_{frag}\right) = \left\{ \begin{aligned} \rightarrow 1 \textrm{ for } \Delta v_{11} \geq v_{frag} \\ \rightarrow 0 \textrm{ for } \Delta v_{11} \ll v_{frag} \end{aligned} \right. \end{align*}

to \(\Pi_{frag}\). Obviously, we also require \(\Pi_{stick} = 1-\Pi_{frag}\). Having formulated the expressions for population transfer via fragmentation and sweep-up, we can now write down the source terms for the dust population surface densities

\begin{equation} \dot{\Sigma}_0 = \dot{\Sigma}_{1 \rightarrow 0} - \dot{\Sigma}_{0 \rightarrow 1} = - \dot{\Sigma}_1, \end{equation}

which satisfy conservation of mass.

In the fragmentation regime the target power law exponent can change depending on the velocity regime:

\begin{equation} q_{\mathrm{frag}} = \begin{cases} q_{\mathrm{turb.1}} = -3.75, & \text{small particle turb. regime} \\ q_{\mathrm{turb.2}} = -3.5, & \text{fully intermediate turb. regime} \\ q_{\mathrm{drift\text{-}frag}} = -3.75, & \text{drift-dominated regime} \end{cases} \end{equation}

where the appropriate velocity thresholds and the trasition function can be found in Pfeil et al (2024).

Step-by-step Summary

1. Calculate density distribution exponent

\begin{align*} \xi = \frac{\log\left(\Sigma_1/\Sigma_0\right)}{\log\left(s_{max}/s_{int}\right)} - 4 \end{align*}

2. Determine characteristic particle sizes

\begin{align*} a_0 = \left\{ \begin{aligned} \frac{s_{int}s_{min}}{s_{int}-s_{min}} \log\left(s_{int}/s_{min}\right) \textrm{ for }\xi = -5 \\ \frac{s_{int}-s_{min}}{\log\left(s_{int}/s_{min}\right)} \textrm{ for }\xi = -4 \\ \frac{\xi+4}{\xi+5} \frac{s_{int}^{\xi+5}-s_{min}^{\xi+5}}{s_{int}^{\xi+4}-s_{min}^{\xi+4}} \textrm{ for }\xi \neq -4, -5\end{aligned} \right. \end{align*}

\begin{align*} a_1 = \left\{ \begin{aligned} \frac{s_{max}s_{int}}{s_{max}-s_{int}} \log\left(s_{max}/s_{int}\right) \textrm{ for }\xi = -5 \\ \frac{s_{max}-s_{int}}{\log\left(s_{max}/s_{int}\right)} \textrm{ for }\xi = -4 \\ \frac{\xi+4}{\xi+5} \frac{s_{max}^{\xi+5}-s_{int}^{\xi+5}}{s_{max}^{\xi+4}-s_{int}^{\xi+4}} \textrm{ for }\xi \neq -4, -5\end{aligned} \right. \end{align*}

\begin{align*} \langle a \rangle_m = \left\{ \begin{aligned} \frac{s_{max}s_{min}}{s_{max}-s_{min}} \log\left(s_{max}/s_{min}\right) \textrm{ for }\xi = -5 \\ \frac{s_{max}-s_{min}}{\log\left(s_{max}/s_{min}\right)} \textrm{ for }\xi = -4 \\ \frac{\xi+4}{\xi+5} \frac{s_{max}^{\xi+5}-s_{min}^{\xi+5}}{s_{max}^{\xi+4}-s_{min}^{\xi+4}} \textrm{ for }\xi \neq -4, -5\end{aligned} \right. \end{align*}

3. Derive source term for growth of maximum particle size

\begin{align*} \dot{s}_{max} = \frac{\Sigma_1}{\rho_s}\Delta v \left[\frac{1-\left(\frac{v_{frag}}{\Delta v}\right)^s}{1+\left(\frac{v_{frag}}{\Delta v}\right)^s}\right] \end{align*}

with \begin{align*} \Delta v = \Delta v_{0.4 s_{max}, s_{max}} \end{align*} and \(s=3\)

4. Update parameters

E.g. particle masses, Stokes numbers, scale heights, dust midplane densities, backreaction constants, relative velocities, diffusivities, etc.

5. Calculate source terms for dust population surface densities

\begin{align*} \dot{\Sigma}_{0 \rightarrow 1} = \frac{\Sigma_0 \Sigma_1 \sigma_{01} \Delta v_{01}}{m_1 2 \pi H_0 H_1} \int\limits_{-\infty}^\infty \exp\left[-\frac{z^2}{2}\frac{H_0^2 + H_1^2}{H_0^2 H_1^2}\right] \mathrm{d}z = \frac{\Sigma_0 \Sigma_1 \sigma_{01} \Delta v_{01}}{m_1 \sqrt{2 \pi \left(H_0^2 + H_1^2\right)}} \end{align*}

\begin{align*} \dot{\Sigma}_{1 \rightarrow 0} = \frac{\Sigma_1^2 \sigma_{11} \Delta v_{11}}{m_1 2 \pi H_1^2} \mathcal{F} \int \limits_{-\infty}^\infty \exp \left[-\frac{z^2}{H_1^2}\right]\mathrm{d}z = \frac{\Sigma_1^2 \sigma_{11} \Delta v_{11}}{m_1 2 \sqrt{\pi} H_1} \mathcal{F} \end{align*}

with \begin{align*} \mathcal{F} = \sqrt{\frac{2H_1^2}{H_0^2 + H_1^2}} \frac{\sigma_{01}\Delta v_{01}}{\sigma_{11}\Delta v_{11}}\left(\frac{s_{max}}{s_{min}}\right)^{-\left(\xi'+4\right)} \end{align*}

6. Integrate \(\Sigma_0, \Sigma_1\) and \(s_{max}\)

\begin{align*} \frac{\partial\Sigma_i }{\partial t} + \vec{\nabla} \left[ \Sigma_i \left(\vec{v}_g+\vec{v}_{dr,i}\right)\right] = \frac{\partial\Sigma_i }{\partial t} \vert_{coag} \end{align*}

with:

\begin{align*} \frac{\partial\Sigma_0 }{\partial t}\vert_{coag} = \dot{\Sigma}_{1 \rightarrow 0} - \dot{\Sigma}_{0 \rightarrow 1} = - \frac{\partial\Sigma_1 }{\partial t}\vert_{coag} \end{align*}

\begin{align*} \frac{\partial\left(\Sigma_1s_{max}\right)}{\partial t} + \vec{\nabla} \left[ \Sigma_1 s_{max} \left(\vec{v}_g+\vec{v}_{dr,1}\right)\right] = \dot{s}_{max} \end{align*}

with: \begin{align*} \dot{s}_{max} = \frac{\Sigma_1}{\rho_s}\Delta v \left[\frac{1-\left(\frac{v_{frag}}{\Delta v}\right)^s}{1+\left(\frac{v_{frag}}{\Delta v}\right)^s}\right] \end{align*}