Context. Dust coagulation and fragmentation impact the structure and evolution of protoplanetary disks and set the initial conditions for planet formation. Dust grains dominate the opacities they determine the cooling times of the gas via thermal accommodation in collisions they influence the ionization state of the gas and the available grain surface area is an important parameter for the chemistry in protoplanetary disks. Therefore dust evolution is an effect that should not be ignored in numerical studies of protoplanetary disks. Available dust coagulation models are however too computationally expensive to be implemented in largescale hydrodynamic simulations. This limits detailed numerical studies of protoplanetary disks including these effects mostly to onedimensional models. Aims. We aim to develop a simple  yet accurate  dust coagulation model that can be easily implemented in hydrodynamic simulations of protoplanetary disks. Our model shall not significantly increase the computational cost of simulations and provide information about the local grain size distribution. Methods. The local dust size distributions are assumed to be truncated power laws. Such distributions can be fully characterized by only two dust fluids large and small grains and a maximum particle size truncating the power law. We compare our model to state oftheart dust coagulation simulations and calibrate it to achieve a good fit with these sophisticated numerical methods. Results. Running various parameter studies we achieved a good fit between our simplified threeparameter model and DustPy a stateoftheart dust coagulation software. Conclusions. We present TriPoD a subgrid dust coagulation model for the PLUTO code. With TriPoD we can perform twodimensional vertically integrated dust coagulation simulations on top of a hydrodynamic simulation. Studying the dust distributions in twodimensional vortices and planetdisk systems is thus made possible.
['methods: numerical', 'planets and satellites: formation', 'protoplanetary disks', 'stars: protostars', 'dust', 'extinction', 'Astrophysics - Earth and Planetary Astrophysics']
TriPoD TriPopulation size distributions for Dust evolution Coagulation in vertically integrated hydrodynamic simulations of protoplanetary disks
{'Coagulation in vertically integrated hydrodynamic simulations of protoplanetary disks.': '<!-- image -->  \nThomas Pfeil  \n123  \n, Til Birnstiel  \n14  \n, and Hubert Klahr  \n2  \n- 1 University Observatory, Faculty of Physics, Ludwig-Maximilians-Universität München, Scheinerstr. 1, D-81679 Munich, Germany\n- 3 Center for Computational Astrophysics, Flatiron Institute, 162 Fifth Avenue, New York, NY 10010, USA e-mail: tpfeil@flatironinstitute.org\n- 2 Max-Planck-Institut für Astronomie, Königstuhl 17, D-69117 Heidelberg, Germany\n- 4 Exzellenzcluster ORIGINS, Boltzmannstr. 2, D-85748 Garching, Germany  \nNovember 1, 2024', 'ABSTRACT': 'Context. Dust coagulation and fragmentation impact the structure and evolution of protoplanetary disks and set the initial conditions for planet formation. Dust grains dominate the opacities, they determine the cooling times of the gas via thermal accommodation in collisions, they influence the ionization state of the gas, and the available grain surface area is an important parameter for the chemistry in protoplanetary disks. Therefore, dust evolution is an e ff ect that should not be ignored in numerical studies of protoplanetary disks. Available dust coagulation models are, however, too computationally expensive to be implemented in large-scale hydrodynamic simulations. This limits detailed numerical studies of protoplanetary disks, including these e ff ects, mostly to one-dimensional models. Aims. We aim to develop a simple-yet accurate-dust coagulation model that can be easily implemented in hydrodynamic simulations of protoplanetary disks. Our model shall not significantly increase the computational cost of simulations and provide information about the local grain size distribution.  \nMethods. The local dust size distributions are assumed to be truncated power laws. Such distributions can be fully characterized by only two dust fluids (large and small grains) and a maximum particle size, truncating the power law. We compare our model to state-of-the-art dust coagulation simulations and calibrate it to achieve a good fit with these sophisticated numerical methods.  \nResults. Running various parameter studies, we achieved a good fit between our simplified three-parameter model and DustPy , a state-of-the-art dust coagulation software.  \nConclusions. We present TriPoD , a sub-grid dust coagulation model for the PLUTO code. With TriPoD , we can perform twodimensional, vertically integrated dust coagulation simulations on top of a hydrodynamic simulation. Studying the dust distributions in two-dimensional vortices and planet-disk systems is thus made possible.  \nKey words. protoplanetary disks - dust evolution - hydrodynamics - methods: numerical', '1. Introduction': "Models of dust coagulation in protoplanetary disks are required to understand the formation of cm-sized pebbles (Brauer et al. 2008; Birnstiel et al. 2009) and km-sized planetesimals (Wetherill & Stewart 1989; Schlichting & Sari 2011; Kobayashi et al. 2016; Lau et al. 2022; Dr˛a˙zkowska et al. 2023); they are indispensable for the interpretation of observational data (Birnstiel et al. 2018; Dullemond et al. 2018) and necessary to simulate the assembly of whole planetary systems (Lichtenberg et al. 2021; Emsenhuber et al. 2021). The size of dust grains also determines their aerodynamic properties and thus sets the timescales at which grains drift towards the central star or collect in local pressure maxima (Whipple 1972; Weidenschilling 1977). Furthermore, dust is the dominating source of opacity in circumstellar disks, which means the size distribution of the dust grains has a strong influence on the disks' thermal structure (Muley et al. 2023) and hydrodynamics (Lesur et al. 2023), as well as on the interpretation of observations (Birnstiel et al. 2018; Leiendecker et al. 2022; Bergez-Casalou et al. 2022; Antonellini et al. 2023). In addition, the presence of small grains sets limits to the  \ndisks' ionization and is thus also important for studies of magnetohydrodynamic mechanisms like the MRI (Balbus & Hawley 1991) and magnetized disk winds (Blandford & Payne 1982) with non-ideal MHD e ff ects (Guillet et al. 2020; Pascucci et al. 2023; Tsukamoto & Okuzumi 2022). In 2018, observation with the Atacama Large Millimeter / submillimeter Array (ALMA) revealed that a broad variety of substructures exist in the spatial distribution of dust in protoplanetary disks (Andrews et al. 2018). Numerous gaps, spirals, and vortices have since been observed in the dust continuum emission (e.g., Pérez et al. 2018; Baruteau et al. 2019; Tsukagoshi et al. 2022) and also in molecular line observations (see Öberg et al. 2021). The existence of these structures raises questions regarding their origins and how they impact the formation and composition of planetesimals and planets within the disks-which makes hydrodynamic simulations of protoplanetary disks, including dust coagulation models, necessary (Birnstiel et al. 2018; Dr˛a˙zkowska et al. 2019).  \nModeling the evolution of the solid content of protoplanetary disks has proven to be an expensive task in computational astrophysics. Many dust evolution models (e.g. Nakagawa et al. 1981; Weidenschilling 1980; Brauer et al. 2008; Birnstiel et al. 2010;  \nFig. 1: Comparison between the full coagulation model DustPy (upper row, 171 dust fluids) and our new three-parameter power-law prescription TriPoD , which we implemented in the PLUTO code (two dust fluids). We show two snapshots of the one-dimensional disk models in the first two panels in each row. The third panel shows the local dust size distribution of the respective model at 10 au for both snapshots. The light grid in the background represents size distribution power laws with n ( a ) ∝ a -3 . 5 and n ( a ) ∝ a -3 . 0 .  \n<!-- image -->  \nmax  \n(  \nΣ  \nd  \n(  \na  \n))  \nCharnoz & Taillifet 2012; Dr˛a˙zkowska et al. 2019), numerically solve a discretized form of the Smoluchowski coagulation equation (Smoluchowski 1916). The coagulation equation's integrodi ff erential nature makes solving it, however, numerically costly. This procedure utilizes a grid of grain sizes, meaning that dozens or hundreds of dust fluids have to be modeled in a single simulation, each representing grains of a di ff erent size that interact with every other grain size via collisions. Studies of dust coagulation are often carried out in one-dimensional disk models in either the radial direction (e.g., Lenz et al. 2020; Dr˛a˙zkowska et al. 2021; Pinilla et al. 2021; Gárate et al. 2021; Burn et al. 2022), the vertical direction (e.g, Zsom et al. 2011; Krijt & Ciesla 2016), or limited two-dimensional studies (Dr˛a˙zkowska et al. 2019).  \nTherefore, e ff orts are pursued to solve the coagulation equations more e ffi ciently, such as the use of new numerical methods (e.g., Estrada & Cuzzi 2008; Lombart & Laibe 2021; Lombart et al. 2022, 2024; Marchand et al. 2021, 2022, 2023), or by applying far-reaching simplifications to the physics of dust coagulation (Birnstiel et al. 2012) that make it possible to implement dust coagulation as a subgrid model in hydrodynamic simulations (Tamfal et al. 2018; Vorobyov et al. 2018). Machinelearning-aided techniques also promise a fast, yet simplified approach to model dust coagulation on top of a hydrodynamic simulation (Pfeil et al. 2022).  \nHere, we present a semi-analytic model of dust coagulation, which is based on a two-population approach, originally developed by (Birnstiel et al. 2012) and employed in various forms by others (Tamfal et al. 2018; Vorobyov et al. 2019; Vorobyov & Elbakyan 2019; Vorobyov et al. 2020). This model, however, has the critical disadvantage of not evolving the full dust size distribution but only the maximum particle size. With our new model, we are now able to conduct two-dimensional, vertically integrated hydrodynamics simulations of protoplanetary disks with an evolving dust size distribution at low computational cost. Figure 1 shows an example comparison between a full coagulation simulation, conducted with DustPy (Stammler & Birnstiel 2022)-a full-fledged dust coagulation softwareand a one-dimensional hydrodynamic simulation with the PLUTO  \ncode equipped with our new method. While one-dimensional hydrodynamics simulations do not allow for a direct comparison of performance due to the di ff erent methods to handle the transport ( DustPy utilizes an implicit integration scheme and does not solve the equations of hydrodynamics, but an advectiondi ff usion equation for a Keplerian disk), they allow for detailed tests of the accuracy of our new model.  \nThis article is structured as follows: In Section 2, we briefly review the physics of dust dynamics and coagulation. We also give a short description of two-pop-py (Birnstiel et al. 2012), the progenitor of our new model. Section 3 introduces our new three-parameter dust coagulation model TriPoD and how it is integrated in the PLUTO hydrodynamics code (Mignone et al. 2007). As TriPoD is a highly simplified prescription of dust coagulation, we have to calibrate it to achieve a good fit with full-coagulation models. The respective calibration runs are presented in Section 4. In Section 5, we present test simulations that demonstrate the accuracy of TriPoD in comparison with DustPy simulations. We also give an example of a two-dimensional simulation of a planet-disk system in which we compare the outcome of our new model to the old two-pop-py model. We discuss the limitations of our approach in Section 6 and summarize in Section 7.", '2.1. Dust-gas relative motion': "Dust particles moving in the gaseous protoplanetary disk can experience drag forces due to di ff erences between the equilibrium velocities of dust and gas particles. The gas in a PPD is radially stratified, which means the disk has a radial pressure gradient. Consequently, gas moves on a slightly sub-Keplerian orbit, where hydrostatic equilibrium is given by  \nΩ 2 = 1 R ρ g ∂ P ∂ R + Ω 2 K , (1)  \nwhere P and ρ g are the gas pressure and density respectively, R is the cylindrical, stellocentric radius, and Ω 2 K = GM ∗ / R 3 is the Kep-  \nn angular frequency (see Appendix A for a list of symbols). Conversely, radial pressure forces are negligible for the solid particles, which means their equilibrium orbits would be Keplerian (ignoring the gas drag and additional e ff ects like radiation pressure, etc.). Gas, at velocity u g and dust particles, at velocity u d are aerodynamically coupled via a friction force density  \nf fric = ρ d u d -u g t fric , (2)  \nwhere the strength of the coupling can be characterized by the stopping timescale t fric. In the Epstein regime (Epstein 1924), the friction time can be written  \nt fric = r π 8 ρ m a ρ g c s , (3)  \nwhere ρ m is the particles' material density, a denotes the particle radius, and c s is the sound speed. A useful dimensionless measure of the strength of the coupling between gas and dust particles is the Stokes number St B t fric Ω K. For St ≪ 1, gas and dust particles are well-coupled and the particles quickly adjust to changes in the gas velocity. For St ≫ 1 however, particles are decoupled from the gas and the friction force is no longer significant for the trajectories of the dust grains.  \nThe stopping time is the timescale on which the particles and the gas approach a steady state. The respective terminal velocities of the grains in force equilibrium were derived by Nakagawa et al. (1986). Ignoring additional velocity components of the gas, for example, due to viscous evolution, the respective relative radial velocity between the grains and the gas then follows as  \nv d-g , R = St(1 + ε ) St 2 + (1 + ε ) 2 1 Ω K ρ g ∂ P ∂ R (4) ≈ St St 2 + 1 1 Ω K ρ g ∂ P ∂ R for ε ≪ 1 , (5)  \nwhere ε is the dust-to-gas density ratio. Thus particles drift towards pressure maxima and reach their maximum terminal velocity at a Stokes number of one.", '2.2. Dust-dust relative motion': "Relative velocities between the gas and dust depend on the aerodynamic properties of the dust. Di ff erently sized grains therefore experience relative velocities due to gas drag. For small grains, Brownian motion is of importance. Additionally, turbulence causes random variations in the gas velocities that act on the dust particles according to their aerodynamic coupling to differently sized eddies in the gas.  \nBrownian motion Random molecular motion of particles leads to relative velocities that depend on the respective particles' masses (Brauer et al. 2008)  \n∆ v 01 Brown = r 8 k B T ( m 0 + m 1) π m 0 m 1 , (6)  \nwhere m 0 and m 1 denote the particles' masses, T is the gas temperature, and k B is the Boltzmann constant. This e ff ect is only relevant for the smallest particles on micrometer scales.  \nRelative drift velocities For two particles with Stokes numbers St0 and St1, the relative drift velocities in the case of low dustto-gas ratio, are given by Eq. 5  \n∆ v 01 drift = GLYPH<12> GLYPH<12> GLYPH<12> GLYPH<12> GLYPH<12> GLYPH<12>       St0 St 2 0 + 1 -St1 St 2 1 + 1       1 Ω K ρ g ∂ P ∂ R GLYPH<12> GLYPH<12> GLYPH<12> GLYPH<12> GLYPH<12> GLYPH<12> . (7)  \nRelative settling velocities Dubrulle et al. (1995) derived the vertical dust distribution in a disk with turbulent di ff usion as  \nρ d( z ) = ρ d , mid exp -z 2 2 H d ! , (8)  \nwhere H d is the scale height of the dust, given by  \nH d = H q 1 + St δ , (9)  \n(see also Fromang & Nelson 2009; Binkert 2023). Here, H refers to the gas scale height, and δ denotes the turbulent di ff usivity parameter, which for now is assumed to be equal to the turbulent gas viscosity parameter α . Thus, particles of di ff erent sizes are, on average, also found at di ff erent heights, and thus have di ff erent terminal velocities. The average relative velocities between the two particle populations can then be approximated as  \n∆ v 01 set = | H d1St1 -H d0St0 | Ω K . (10)  \nRelative velocities due to turbulence Ormel & Cuzzi (2007) derived closed-form expressions for the relative particle velocities in di ff erent turbulence regimes, which depend on the Stokes numbers, the turbulent gas velocity, and the local Reynolds number  \nRe = ν turb ν mol ≈ α c s H c s λ mfp = α H ρ g σ H2 µ m p , (11)  \nwhere λ mfp is the mean free path of the gas molecules, µ m p is the mean molecular mass of the gas, and σ H2 is the collisional cross-section of two gas molecules (here H2). The respective derivations assume a Kolmogorov turbulent energy cascade (Kolmogorov 1941). We use an implementation of these velocities identical to the one utilized in the full dust coagulation code DustPy .", '2.3. Dust coagulation': 'Dust particles in protoplanetary disks undergo collision since they experience di ff erential velocities due to their interaction with the gas. Surface forces can lead to sticking in such collisions and thus facilitate the growth of dust particles. If collision velocities are too high, they can lead to fragmentation. The Smoluchowski equation (Smoluchowski 1916) describes the evolution of continuous mass (or size) distributions of grains n ( m ), as a consequence of these processes. In this work, however, we are not numerically solving the coagulation equation as in full-fledged coagulation models like DustPy (Stammler & Birnstiel 2022). Instead, we use the results gained with such elaborate numerical methods to construct a simplified, semi-analytic prescription for dust coagulation. For this, it is instructive to have a look at some main results obtained with full coagulation models and simple analytic derivations. One of these simplifying assumptions is a monodisperse size distribution. Kornet et al.  \n(2001) have shown that for such a case, the particle growth rate can be written as  \n˙ a = ρ d ρ m ∆ v, (12)  \nwhere ∆ v denotes the relative velocity between the grains. If the relative velocities are assumed to be caused by gas turbulence in the fully intermediate regime, one finds that the growth of the particles is occurring on a timescale  \nt grow = 1 ε Ω K . (13)  \nHowever, dust grains can only grow in size as long as their relative velocities due to di ff erent aerodynamic properties are not above a critical fragmentation velocity v frag, for which collisions would lead to the destruction of the particles. Another possible outcome of grain collisions is bouncing, as shown by laboratory experiments (Güttler et al. 2010) and studied in numerical models (Zsom et al. 2010; Dominik & Dullemond 2024). It is possible to derive analytic estimates for the maximum reachable particle size, given a certain fragmentation velocity. Birnstiel et al. (2012) derived the respective maximum particle size in the turbulent fragmentation limit as  \na turb-frag = r 8 π ρ g 3 ρ m v 2 frag α c s Ω K . (14)  \nFurthermore, dust grains undergo radial drift, which results in relative velocities between grains of di ff erent sizes (see Eq. 5). The resulting collisions can also lead to a drift-fragmentation limit, which is given by  \na drift-frag = r 8 π ρ g ρ m v frag c s P 1 - N GLYPH<12> GLYPH<12> GLYPH<12> GLYPH<12> GLYPH<12> ∂ P ∂ r GLYPH<12> GLYPH<12> GLYPH<12> GLYPH<12> GLYPH<12> -1 , (15)  \nwhere the constant N is approximately 0.5 (Birnstiel et al. 2012). Dust particles can thus not reach sizes larger than  \na frag = min( a drift-frag , a turb-frag) . (16)  \nIn these cases, the typically reached size distributions are approximately power laws with characteristic exponents. Birnstiel et al. (2011) used coagulation simulations and analytical calculations to find these fragmentation-limited size distributions. They derived analytic expressions for the resulting power-law exponent of a mass distribution n ( m ) ∝ m -q m in three di ff erent regimes, translating to the size distribution n ( a ) ∝ a q = a -3 q m + 2 . If coagulation and fragmentation happen simultaneously, the power-law exponent can be written as  \nq = 2 -3 2 ( ν + ξ + 1) , (17)  \nwhere ν is a kernel index that depends on how the relative velocities change with particle size. The parameter ξ determines the typical distribution of fragments in a destructive collision. It is usually set to the canonical value of 11 / 6. Birnstiel et al. (2011) determined that in the first regime of turbulence, as derived by Ormel & Cuzzi (2007), where the relative velocities scale linearly with particle size, ν = 1 and thus q = -3 . 75. The same is true for relative drift velocities, which also linearly depend on the di ff erence in Stokes numbers between the colliding particles (as long as St ≪ 1). This means that typical size distributions in coagulation-fragmentation equilibrium, with collisions driven by  \nFig. 2: Typical size distributions in regions of a protoplanetary disk model that are either turbulence-dominated (dark red lines) or drift-dominated (blue lines). Both regimes result in distinct power-law exponents of the distribution that we overplot as dashed lines. Assuming a power law as the overall shape of the distribution makes it possible to approximately describe it with only three parameters: The cuto ff particle size a max and the two densities Σ 0 (contained in the size interval [ a min , √ a max a min) and Σ 1 (contained in the size interval [ √ a max a min , a max]). This approximation is the basis of our new dust coagulation model.  \n<!-- image -->  \ndi ff erential drift or the first regime of turbulence, can be approximated as power laws n ( a ) ∝ a -3 . 75 , as illustrated in Fig. 2. This case is thus relevant whenever the particle fragmentation velocity is low and only small particles exist, or if the drift velocities dominate over the turbulent velocities, as in the outer regions of protoplanetary disks.  \nFinally, we have the stages of dust growth in which the fragmentation barrier is not yet reached and the largest particles are undergoing a sweep-up growth. In this case, typical size distributions are steeper and we assume q = -3 . 0 in this case (Simon et al. 2022; Birnstiel 2023).  \nThe other prominent case is relevant whenever turbulence is causing relative velocities between particles in the so-called fully intermediate regime. Then, ν = 5 / 6, and one finds that q = -3 . 5, which is equal to the Mathis, Rumpl, & Nordsieck (1977) (MRN) distribution. This particular case is relevant in the inner parts of protoplanetary disks, where the particles grow to the largest sizes and where radial drift is less relevant.  \nRadial drift itself also sets a limit to the maximum particle size, which is reached when the radial drift time scale becomes equal to the local growth time scale  \na drift ˙ a ! = R | ˙ R | a drift = r 8 π ρ d ρ m v K c s γ , (18)  \nwhere v K is the Keplerian velocity and γ = GLYPH<12> GLYPH<12> GLYPH<12> GLYPH<12> d log P d log R GLYPH<12> GLYPH<12> GLYPH<12> GLYPH<12> denotes the radial double-logarithmic pressure gradient. The drift limit is relevant in the outer regions of protoplanetary disks, where drift  \ncan become rapid for large particles. In the drift limit, grains do not undergo fragmentation. Therefore, typical size distributions are sweep-up dominated and accordingly steep. We assume q = -3 . 0, which is typically a good fit to full coagulation simulations.', '2.4. The two-pop-py model by (Birnstiel et al. 2012)': 'With two-pop-py , Birnstiel et al. (2012) introduced a strongly simplified and very fast method to model the e ff ects of dust coagulation in protoplanetary disks. In this method, dust is realized as a single fluid that drifts relative to the gas. The flux calculation, however, considers two dust species. The small population represents the monomers with fixed size a mon and is assumed to move along with the gas. Larger, grown grains make up the population of size a gr( t ), which is evolving in time, following the monodisperse dust growth rate, which is limited by the abovediscussed growth barriers  \na gr( t ) = min " min( a drift , a turb-frag , a drift-frag) , a mon exp t -t 0 t grow !# , (19)  \nto simulate the initial phase of coagulation and the growth limits. Both species are associated with a drift velocity according to their Stokes numbers Stmon / gr = a mon / gr ρ m π / 2 Σ g, as  \nv mon / gr = v gas 1 + St 2 mon / gr + Stmon / gr 1 + St 2 mon / gr 1 ρ g Ω K ∂ P ∂ r , (20)  \nwhere the first term takes into account gas velocities arising from viscous disk evolution and the second term is the radial drift velocity derived by Nakagawa et al. (1986) in the limit of small dust-to-gas ratios. The total dust flux velocity follows as a mass average of both populations, via  \n¯ v = (1 -f m) v mon + f m v gr . (21)  \nThe ratio f m is dependent on the limiting factor of grain growth: fragmentation-limited, drift-fragmentation-limited, drift-limited, or neither when the dust is still in the growth phase  \nf m = ( 0 . 97 if a drift < a frag , a drift-frag 0 . 75 otherwise. (22)  \nGiven the resulting velocity, a flux can be calculated that is used to evolve the dust surface density in time. This method is fast and can be easily implemented in a hydrodynamics code. It has however some serious drawbacks:  \n- -Dust growth is always assumed to be limited by an equilibrium of either fragmentation and coagulation or drift and coagulation. This is mostly true if no substructure is present. In some situations like planetary gaps, however, fluxes and grain sizes are no longer determined by coagulation in the gap but by the supply of small grains that di ff use into the gap and the e ffi cient removal of larger grains. In such cases, grain sizes can be underestimated by two-pop-py , which would assume the drift limit.\n- -Although the maximum grain size is known, no information on the actual size distribution is provided beyond the knowledge of whether the distribution is drift-limited or fragmentation-limited.\n- -The di ff erences between the size distribution in the fragmentation limit and the drift-fragmentation limit are not taken into account. The mass fraction f m only considers whether the maximum particle size is drift-limited.\n- -As growth is always assumed to be driven by collisions in the fully intermediate regime of turbulence, the actual growth timescale can be underestimated (Powell et al. 2019).\n- -The fragmentation limit is reached instantaneously and not gradually. This is especially problematic if the timescale for dust advection is short.\n- -The model is only calibrated to reproduce the dust size evolution in protoplanetary disks around 1 M ⊙ stars.\n- -The model overestimates the concentration of dust in pressure bumps due to the treatment of dust being transported as one fluid. Although the flux is calculated in a mass-averaged way that has been determined experimentally, intermediately sized grains are neglected, which would lead to a wider dust distribution in pressure bumps.  \nWith our new model, TriPoD , we aim to mitigate these problems.', '3. The TriPoD model': 'Our new three-parameter dust evolution model TriPoD makes use of the knowledge gained from full-fledged coagulation models that are discussed in the previous section. In particular, the power-law prescriptions of the dust size distribution are the basis on which we build our method. TriPoD describes the powerlaw size distribution with only three parameters: the dust column densities of a small population Σ 0 and the large population Σ 1, as well as a maximum particle size a max at which the distribution is truncated (see Appendix A for a list of symbols used in the context of TriPoD ). The populations are separated at the size a int = √ a min a max, which then defines the power-law exponent of a full size distribution via  \nq = log( Σ 1 / Σ 0) log( a max / a int) -4 . (23)  \nNote, that this formula is in general independent of whether we describe a size distribution in the sense of column densities or volume densities. Although we develop the TriPoD method in this paper for use in vertically integrated size distributions, we could also do so for volume densities. In the following, we will thus oftentimes use the dust-to-gas ratio ε , which can be interpreted as either the vertically integrated, or the local version. Assuming the dust size distribution to follow a truncated power law n ( a ) ∝ a q , we can write the dust-to-gas ratio size distribution as ε ( a ) ∝ a q m ∝ a q a 3 . Normalizing this to the total dust-to-gas ratio ε tot, we get  \nε ( a ) =                  ε tot( q + 4) a q + 4 max -a q + 4 min a q + 3 for q , -4 ε tot log( a max) -log( a min) 1 a for q = -4 , (24)  \nwhere a min is the minimum particle size and a max is the maximum particle size of the truncated distribution. From this, the dust-to-gas ratio within a given size interval [ a I , a II] with a min ≤  \na I < a II ≤ a max, can be calculated as  \nε a II a I B Z a II a I ε ( a ) d a =                    ε tot a q + 4 II -a q + 4 I a q + 4 max -a q + 4 min for q , -4 ε tot log( a II) -log( a I) log( a max) -log( a min) for q = -4 . (25)  \nSimilarly, the mass-averaged particle size of this distribution is defined as  \n⟨ a ⟩ a II a I B R a II a I ε ( a ) a d a R a II a I ε ( a ) d a =                                  a II a I a II -a I log a II a I ! for q = -5 q + 4 q + 5 a q + 5 II -a q + 5 I a q + 4 II -a q + 4 I for q , -5 , -4 a II -a I log( a II) -log( a I) for q = -4 . (26)  \nIn the TriPoD model, we define two particle populations that together contain the entire dust density of the distribution  \nε 0 = Z a int a min ε ( a ) d a ε 1 = Z a max a int ε ( a ) d a . (27)  \nThe mass-averaged particle sizes of both populations are then given by a 0 B ⟨ a ⟩ a int a min and a 1 B ⟨ a ⟩ a max a int (see Eq. 26). It can be shown that the two populations ε 0 and ε 1 exactly represent the power-law distribution if a int is defined as the geometric mean of the maximum and minimum size a int = √ a max a min. The powerlaw exponent q is then given by Eq. 23. Knowing only ε 0 and ε 1 and the maximum size a max thus allows us to reconstruct the entire size distribution.  \nOn a size grid with N cells ai , spanning from a min to aN ≥ a max, we can write the mass in a single size bin, and likewise the entire size distribution as  \nε i =                    ε tot a q + 4 i + 1 / 2 -a q + 4 i -1 / 2 a q + 4 max -a q + 4 min θ ( a ) for q , -4 ε tot log( ai + 1 / 2) -log( ai -1 / 2) log( a max) -log( a min) θ ( a ) for q = -4 , (28)  \nwhere ai -1 / 2 and ai + 1 / 2 denote the cell interfaces on the size grid, ε tot = P N i = 0 ε i is the total dust-to-gas density ratio and θ ( a ) = Θ ( ai -a min) Θ ( ai -a max) represents two Heaviside step functions that cut o ff the distribution at the minimum and maximum particle sizes. This means we can directly compare our model results with full dust coagulation models like DustPy , which evolve a large grid of sizes instead of just two fluids in our case. We illustrate this in Fig. 2, where we overplot the detailed size distributions, obtained with a full coagulation model, with their respective three-parameter size distribution representation. The respective dust column densities Σ 0 and Σ 1 of the two populations are shown as the horizontal bars spanning the respective size ranges.  \nIn the following, we describe how we evolve the threeparameter size distribution (i.e., ε 0 , ε 1 and a max) in time, using a semi-analytic description of dust coagulation. The prescriptions given in this paper represent the first iteration of our new TriPoD model that is derived and calibrated for vertically integrated disk models. Therefore, all calculations include gas and dust column densities instead of volume densities.', '3.1. Particle growth (evolution of a max )': 'We model particle growth within the monodisperse approximation (Eq. 12). The growth limits are realized by comparing a given fragmentation velocity v frag with the velocities between large grains ∆ v max, given by turbulence, di ff erential settling and drift, and Brownian motion. For this, we modify our growth rate by a sigmoid-like function, which leads to growth for ∆ v max < v frag, and decay for ∆ v max > v frag, resulting in  \n˙ a max = Σ 1 ∆ v max ρ m √ 2 π H 1          GLYPH<16> v frag ∆ v max GLYPH<17> s -1 1 + GLYPH<16> v frag ∆ v max GLYPH<17> s          , (29)  \nwith s being a parameter controlling the steepness of the transition from growth to fragmentation, H 1 being the scale height of large dust grains. The relative grain velocity ∆ v max is defined between grains of size a max and f ∆ v a max, where f ∆ v ∈ (0 , 1) is a model parameter. Determining f ∆ v and s is the main task during the calibration of our model with respect to the full coagulation code DustPy (see Section 4.1 and the red numbers in Table 3).', '3.2. Fragmentation and sweep-up (evolution of Σ 0 and Σ 1 )': "In our model, two e ff ects account for the evolution of the threeparameter size distribution and the interaction between the populations; fragmentation transfers mass from the large population to the small population, while collisions between larger and smaller particles lead to sweep-up, and thus mass transfer from the small bin to the large bin. Erosion of large particles due to collisions with small grains is not accounted for in our model. Wedescribe the sweep-up process via the collision rates between large and small particles  \n˙ ρ d , 0 → 1 = ρ 2 g ε 0 ε 1 m 0 m 1 σ 01 ∆ v 01 m 0 , (30)  \nwhere m 0 and m 1 are the representative particle masses, ∆ v 01 is the representative relative velocity between the large and small particles, and σ 01 is the representative collision cross section. These quantities are derived from the respective population's mass-averaged particle sizes (see Eq. 26).  \nFragmentation predominantly occurs in collisions between two large grains. The corresponding transfer rate is thus determined by the collision rates between large grains  \n˙ ρ d , 1 → 0 = ρ 2 g ε 2 1 m 2 1 σ 11 ∆ v 11 m 1 F . (31)  \nHere F represents a function that regulates the relative e ff ectiveness of sweep-up and fragmentation. It is a function of the grain size, the desired power-law exponent of the size distribution, and the relative velocities ∆ v 11, defined for collisions between particles of size a 1 and f ∆ v a 1. As we cannot model the microphysics of collisions between grains in our simplified framework for dust evolution, we base the functional form of F on the  \nwell-understood results of full coagulation models, which treat the evolution of the size distribution as the result of collisions between grains of all present sizes.  \nFor this first version of our three-parameter model, we are only considering vertically integrated disk models, which means our mass transfer rates are given by  \n˙ Σ 0 → 1 = Σ 0 Σ 1 σ 01 ∆ v 01 m 1 q 2 π ( H 2 0 + H 2 1 ) (32)  \n˙ Σ 1 → 0 = Σ 2 1 σ 11 ∆ v 11 m 1 q 4 π H 2 1 ˜ F , (33)  \n(see Appendix B). In order to determine the vertically integrated version of F , named ˜ F , we consider the steady state between fragmentation and sweep-up. In this equilibrium, a steady size distribution would be reached. Given ˙ Σ d , 0 → 1 = ˙ Σ d , 1 → 0, we arrive at  \n˜ F = s 2 H 2 1 H 2 0 + H 2 1 σ 01 σ 11 ∆ v 01 ∆ v 11 a max a int ! -( q + 4) , (34)  \nwhere q is the desired power-law exponent of the distribution which will be reached on the dominating collisional timescale.  \nMass redistribution due to sweep-up and fragmentation is realized by defining the source terms of both populations in a totalmass-conserving manner as  \n˙ Σ 0 = ˙ Σ 1 → 0 -˙ Σ 0 → 1 ˙ Σ 1 = -˙ Σ 0 . (35)  \nDetermining the size distribution power-law exponent q  \nWe have summed up the typical particle size distribution exponents in Section 2.3, which were determined by Birnstiel et al. (2011) for distributions in coagulation-fragmentation equilibrium. These are given by  \nq frag =            q turb . 1 = -3 . 75 small particle turb. regime q turb . 2 = -3 . 5 fully intermediate turb. regime q drift-frag = -3 . 75 drift-dominated regime  \nq sweep = -3 . 0 not in equilibrium (see Birnstiel 2023)  \nOur task is now to find a way to smoothly switch between these regimes in our three-parameter coagulation model depending on which regime is prevailing under the given conditions.  \nFirstly, we can determine whether to apply the equilibrium size distributions q frag or whether the dust has not yet reached coagulation-fragmentation equilibrium and collisions lead predominantly to sticking which results in q sweep. For this, we define a transition function  \np frag B ( → 1 for ∆ v max > v frag → 0 for ∆ v max < v frag p sweep = 1 -p frag ⇒ q = q frag · p frag + q sweep · p sweep . (36)  \nWe now have to determine the equilibrium size distribution exponent p frag. For this, we can again define transition functions. We determine whether the small particle turbulence regime  \n(turb.1) or the fully intermediate regime (turb.2) dominates the relative turbulent velocities  \np turb.1 B ( → 1 for ∆ v turb.1 > ∆ v turb.2 → 0 for ∆ v turb.1 < ∆ v turb.2 p turb.2 = 1 -p turb.1 ⇒ q turb-frag = q turb.1 · p turb.1 + q turb.2 · p turb.2 . (37)  \nLastly, we must determine whether we are in the turbulencedominated regime or in the drift-dominated regime. Similar to before, we define  \np drift B ( → 1 for ∆ v drift > ∆ v turb → 0 for ∆ v drift < ∆ v turb p turb = 1 -p drift ⇒ q frag = q drift-frag · p drift + q turb-frag · p turb , (38)  \nwhere q turb-frag comes from Eq. 37. With this, we have everything we need to calculate q from Eq. 36 and we can determine the mass exchange rates from Eq. 35. The exact form of the transition functions is not of great importance as long as the transition is su ffi ciently fast but still smooth enough to not cause issues during numerical integration. The choices that worked best in our subsequent tests are listed in Appendix C.", '3.3. Passive dust fluids (in the PLUTO code)': "Weuse the PLUTO 1 code to solve the equations of hydrodynamics in our calibration and test simulations with TriPoD . The Euler equations, solved by PLUTO , read  \n∂ρ g ∂ t + ∇ · ( ρ g u ) = 0 (39)  \n∂ρ g u ∂ t + ∇ · ( ρ g u ⊗ u ) = -∇ P -ρ g ∇ Φ+ ∇ · T , (40)  \nwhere u is the gas velocity vector, Φ is the gravitational potential, and T is the viscous stress tensor. The ideal equation of state is used as a closure relation  \nP = k B T µ m p ρ g . (41)  \nThe PLUTO code allows for the treatment of passive tracer fluids, which are simply advected with the gas following  \n∂ ( ρ g ε ) ∂ t + ∇ · ( ε ρ g u g) = 0 . (42)  \nIn this work, we consider vertically integrated protoplanetary disks and the advected quantities in our TriPoD model are thus the local dust-to-gas ratios of our two dust populations ε 0 = Σ 0 / Σ g and ε 1 = Σ 1 / Σ g. The maximum particle size is defined as a tracer of the large dust population, meaning our third tracer fluid is a max ε 1.  \nThe respective tracer fluxes are modified to simulate a dust fluid that is aerodynamically coupled to the gas and thus undergoes radial and azimuthal drift (in the terminal velocity approximation). To achieve this within PLUTO 's tracer prescription, we add a flux component corresponding to the relative velocity between dust and gas  \nu Σ 0 / 1 = St Σ 0 / 1 St 2 Σ 0 / 1 + 1 1 Ω K ρ g ∇ P , (43)  \nwhere St Σ 0 / 1 is the mass-averaged Stokes number of the respective population. The drift velocities are limited to a fraction of the soundspeed. The third tracer ( a max ε 1) is given the same drift velocity as the large dust population ( ε 1). The tracer fluxes are calculated with the upstream dust density and density-weighted maximum particle size based on the drift velocities at the respective cell interface as  \nF drift Σ 0 , i + 1 / 2 =Σ 0 , i max(0 , v Σ 0 , i + 1 / 2) + Σ 0 , i + 1 min( v Σ 0 , i + 1 / 2 , 0) (44)  \nF drift Σ 1 , i + 1 / 2 =Σ 1 , i max(0 , v Σ 1 , i + 1 / 2) + Σ 1 , i + 1 min( v Σ 1 , i + 1 / 2 , 0) (45)  \nF drift a max Σ 1 , i + 1 / 2 = a max , i Σ 1 , i max(0 , v Σ 1 , i + 1 / 2) + a max , i + 1 Σ 1 , i + 1 min( v Σ 1 , i + 1 / 2 , 0) . (46)  \nDust di ff usion is implemented in a flux-limited manner, where the transport velocity is limited to the turbulent gas velocity v max = √ δ c s 1 + St 2 . The di ff usion fluxes are given by  \nF di ff Σ 0 , i + 1 / 2 = -λ lim δ c s H 1 + St 2 0 Σ g , i + 1 / 2 ∇ ε 0 (47)  \nF di ff Σ 1 , i + 1 / 2 = -λ lim δ c s H 1 + St 2 1 Σ g , i + 1 / 2 ∇ ε 1 (48)  \nF di ff a max Σ 1 , i + 1 / 2 = -λ lim δ c s H 1 + St 2 1 Σ g , i + 1 / 2 ∇ ( a max ε 1) , (49)  \nwhere λ lim is the respective flux limiter function (see Appendix D). All quantities are interpolated to the cell interface. The drift and di ff usion fluxes are added to the advective flux of the dust fluids stemming from the gas motion. We have used this approach in Pfeil et al. (2023), where we also presented some simple test cases of the method. Note that this approach makes use of the terminal velocity approximation and is thus strictly speaking only valid for St ≪ 1. The TriPoD model itself, being a local model, is not bound to this form of dust advection scheme and could be implemented in any multi-fluid-capable hydrodynamics code.", '3.4. Complete right-hand side of the conservation equations (in the PLUTO code)': "All modifications made for our three-parameter dust evolution model can be applied within the framework of PLUTO and without changing the underlying reconstruct-solve-average scheme of the code.  \nSource terms are added to the right-hand side of the conservative hydrodynamics equations describing the evolution of the three-parameter dust size distribution in the framework of PLUTO . For each dimension and evolving variable, we have the right-hand side of the conservation equation  \nR i = -∆ t ∆ V i GLYPH<2> ( A F ) i + 1 / 2 -( A F ) i -1 / 2 GLYPH<3> + ∆ t S i , (50)  \nwhere ∆ t is the timestep, ∆ V , is the cell volume, A the respective cell surface, F is the flux through the interfaces, determined from advection with the gas plus relative terminal velocity and di ff usion fluxes, and S i is the source term, given by fragmentation, sweep-up, and growth.  \nFor the two dust fluids, the fluxes are calculated via Eq. 44 to Eq. 49, and the source term of the respective dust fluid is calculated from Eq. 35  \nF Σ 0 / 1 , i ± 1 / 2 = F di ff Σ 0 / 1 , i ± 1 / 2 + F drift Σ 0 / 1 , i ± 1 / 2 + F adv Σ 0 / 1 , i ± 1 / 2 (51) S Σ 0 / 1 , i = ˙ Σ 0 / 1 , i , (52)  \nArticle number, page 8 of 22  \nwhere F adv Σ 0 / 1 , i ± 1 / 2 is the flux component that is calculated by PLUTO for the passive advection of the tracers. The maximum particle size is advected together with a large population as a tracer. The fluxes and growth rate are  \nFa max Σ 1 , i ± 1 / 2 = F di ff a max Σ 1 , i ± 1 / 2 + F drift a max Σ 1 , i ± 1 / 2 + F adv a max Σ 1 , i ± 1 / 2 (53)  \nS a max Σ 1 , i = a max , i ˙ Σ 1 , i + Σ 1 , i ˙ a max , i . (54)  \nHere, F adv a max Σ 1 , i ± 1 / 2 is the flux component that is calculated by PLUTO for the passive advection of the tracer and ˙ a max is the particle growth rate, given by Eq. 29. The evolution equations are evolved in time with PLUTO 's standard third-order Runge-Kutta routine. The timestep is limited by the CFL condition, which we don't have to adjust since the dust drift velocities are limited to a fraction of the speed of sound and the typical growth timescale is much longer than an orbital timescale.", '4. Calibration': "Our model has several free parameters, which have to be calibrated in comparison to full dust coagulation simulations. For the growth rate of the dust (Eq. 29), the main parameters are s (determining the steepness of the transition from growth to fragmentation), and the parameter f ∆ v , which determines the relative velocities ∆ v 11 and ∆ v max through the size ratio between the grains.  \nWe ran a series of one-dimensional simulations for calibration of the model against the full-coagulation code DustPy . We set up a 1D disk model similar to the standard DustPy model. We excluded viscous evolution since we are only interested in comparing the dust evolution in both models for the time being. The radial disk structure follows  \nΣ g( R ) = (2 + β Σ ) M disk 2 π R 2 c R R c ! β Σ exp        -R R c ! 2 + β Σ        B Σ g , 0 R R 0 ! β Σ exp        -R R c ! 2 + β Σ        , (55)  \nwhere we employed the code units Σ g , 0 = 733 . 28 g cm -2 (derived from a DustPy setup with M disk = 0 . 05 M ⊙ ), β Σ = -0 . 85, β T = -0 . 5, and R 0 = 1 au. In PLUTO , the radial temperature structure is expressed in terms of the speed of sound, which is given by  \nc 2 s = T 0 k B µ m p = 0 . 05 L ∗ 4 π R 2 σ SB ! 1 4 k B µ m p B v 2 0 H 0 R 0 ! 2 R R 0 ! β T , (56)  \nwith v 0 = √ GM ∗ / R 0 . In our simulations, this was parameterized by the disk aspect ratio H 0 / R 0 at reference radius R 0.  \nThe simulations were initialized with a total dust-to-gas ratio of 1%, and an initial maximum particle size of 1 µ m. Particles larger than the initial drift limit were excluded from the initial dust profile to avoid an inwards drifting wave of large dust at the beginning of the simulation. The initial dust size distribution in both DustPy and TriPoD follows a power law with q = -3 . 5 . The stellar and disk structure parameters are shown in Table 1 and the parameters for the di ff erent calibration runs are shown in Table 3.  \nFig. 3: Parameter study for the factor f ∆ v , which is the most important parameter for the growth rate calibration of our model. For this, we calculated the deviation between the mass-averaged particle size's growth rate in DustPy and TriPoD . In the top row, we show models for di ff erent global values of f ∆ v . It can be seen that f ∆ v = 0 . 4 seems to fit the dust growth rate best.  \n<!-- image -->  \nFig. 4: Parameter study for the factor s , which determines the steepness of the transition from sweep-up growth to coagulationfragmentation equilibrium in Eq. 29. Heat maps show the deviation of the growth rate in TriPoD from the growth rate in DustPy . The lower row shows the average over the size dimension. The influence of a variation in s can be seen around the fragmentation limit, where a larger s corresponds to a steeper transition (fast growth close to the fragmentation limit), while a small s corresponds to a slower transition where growth rates decrease earlier and more slowly. The case of s = 3 has the overall best fit with the DustPy model.  \n<!-- image -->  \nTable 1: Stellar parameters and disk parameters for the simulations presented in this work.  \nNotes. Calibration runs are performed for the one-solar-mass case.  \nTable 2: Dust properties for the calibration runs (Section 4.1, Section 4.2) and the test simulations with di ff erent stellar masses Section 5.1.", '4.1. Calibrating the growth rate': 'Determining the appropriate relative grain velocities to reproduce the full coagulation model with our simplified prescription is the topic of this section.  \nThe particle growth rate in our model can be adjusted by varying the relative velocity ∆ v max in Eq. 29. This is controlled  \nTable 3: Parameters of the one-dimensional simulations performed to calibrate the three-parameter model to the DustPy simulations.  \nNotes. Each simulation in one block is compared to the same DustPy simulation with identical α and δ . The bold red numbers are our final calibrated values used in the following test simulations.  \nover the parameter f ∆ v ∈ (0 , 1), which sets the size ratio between the colliding particles a max and f ∆ v a max, ( ∆ v max) as well as a 1 and f ∆ v a 1, ( ∆ v 11). We have thus determined which particle size ratio best reproduces the relative velocity that determines the overall growth rate. Furthermore, the parameter s in Eq. 29 can be used to adjust the growth rate around the transition from growth to fragmentation, where a small s corresponds to a reduced growth rate close to the fragmentation limit, and a large s corresponds to a steep transition from particle growth to fragmentation-coagulation equilibrium.  \nTo characterize the e ff ects in comparison with the full coagulation model DustPy , we measured the rate of change of the mass-averaged particle size in TriPoD and DustPy and determined the respective deviations for the given parameters. We used the mass-averaged size and not the maximum particle size because it is also a measure of the shape of the size distribution itself and not just the upper cuto ff . To calibrate the growth rate we ran DustPy and TriPoD models without transport (all fluxes are zero). We were, however, still considering the relative drift and sedimentation velocities in the calculation of ∆ v 01 , ∆ v 11 and ∆ v max. We set up the disk models for a solar-mass pre-mainsequence star with a 0 . 05 M ⊙ gas disk and a dust-to-gas ratio of 1 % (see Table 1 and Table 2 for details).  \nFirst, we ran a parameter study for the factor f ∆ v . This parameter can be physically interpreted as the particle size ratio in the  \ncollisions that dominate the growth process. Since we know that growth is dominated by particle collisions where the size difference between projectile and target is not too large, we chose values between 0.2 and 0.6. As can be seen in Fig. 3, a factor of f ∆ v = 0 . 4 leads to the overall best agreement between the particle growth rate in TriPoD and DustPy .  \nThe second parameter setting the growth rate is s , which determines the steepness of the transition from growth to fragmentation in Eq. 29. We know that this transition is usually not leading to an instantaneous cuto ff of the growth rate. In reality, the particles have collision velocities following a MaxwellBoltzmann distribution (Stammler & Birnstiel 2022). Very large values of s would thus lead to a time evolution that does not agree well with DustPy . We ran models with values from 2 to 6, which we compared in Fig. 4. As can be seen, a value of s = 3 leads to the best fit between the growth rates close to the fragmentation limit.', '4.2. Calibrating the dust transport': "Dust transport is realized by modifying the dust tracer flux (see Section 3.3). For this, the mass-averaged particle size of each dust population is calculated to determine the upstream flux through the respective cell interface. This means that the steepness of the size distribution, which determines the amount of mass in each bin, also partially determines the radial mass flux. As discussed in the introduction, the power law's shape in fragmentation-coagulation equilibrium is well known. In the drift limit, however, it can not be easily deduced from analytical arguments. Numerical dust coagulation simulations with DustPy show that the power law is generally more top-heavy in the drift limit / sweep-up phase and probably best approximated by a value between -2.5 and -3. We have therefore conducted TriPoD simulations covering this parameter range to determine what value best reproduces the late-stage mass evolution of the disk. The mass evolution is generally well reproduced for q sweep = -3 . 0, which is why this value is used for all following simulations.  \nIn comparison with a more realistic size distribution, the power-law prescription lacks the gradual decrease in mass, shortly before the maximum particle size (visible in the right panels of Fig. 1). This means that, even if our power-law size distribution reproduces the real size distribution very well, the mass-averaged sizes of both models will not be identical. For this reason, our dust fluxes will also be slightly di ff erent from the total flux in DustPy . This is generally not a significant e ff ect. We nonetheless accounted for it by using a calibration factor f drift that is multiplied to the calculated mass-averaged particle sizes before the calculation of the dust flux, similar to the fudge factor in the two-pop-py model (Birnstiel et al. 2012). For this, we ran models for di ff erent values of this calibration factor, including dust transport. We neglected any gas transport in these simulations and focused on the dust fluxes, where we conducted one set of simulations with a dust di ff usivity of δ = 10 -3 and one with δ = 0 (no di ff usion). Except for the inclusion of dust transport, the simulations are identical to the setups used in the growth rate calibration. We applied the best-fitting parameters from the growth rate calibration. We use 150 log-spaced radial grid cells to resolve our simulation domain, which was defined between 2 au and 250 au. This is the same resolution as used in the DustPy runs. For comparison, we present a time series of the dust column densities in these simulations runs for f drift = 0 . 6, 0 . 7, 0 . 8, 0 . 9 and 1 . 0 in the upper row of Fig. 5 (Fig. E.1 depicts the case without di ff usion). The lower two panels depict the dust  \n2  \nFig. 5: Comparison between DustPy and our model in a setup with dust di ff usion ( δ = 10 -3 ) and with di ff erent drift calibration factors f drift. Solid lines show the results of our TriPoD calibration runs and dashed lines show the respective DustPy simulation, to which we calibrate our model. The upper row shows a timeseries of the dust column density evolution in three snapshots. In the lower row, we show the mass evolution and the errors with respect to the full coagulation model DustPy . For a factor of f drift = 0 . 8, the mass evolution of the full coagulation model is well reproduced by our three-parameter model.  \n<!-- image -->  \nmass evolution (lower left) and the respective deviations from the full coagulation model DustPy (lower right). As can be seen, the choice of f drift determines the mass flux throughout the disk's dust evolution. Smaller values of f drift correspond to slower dust velocities and lower fluxes. All values result in deviations of ≲ 40 % from the full coagulation model at all times. The value of f drift = 0 . 8, however, shows the smallest mass error. For this value, the disk's dust mass only deviates by 10 % from DustPy during most of the mass evolution. The absolute error only increases to ≲ 30 % after 90 % of the mass has already drifted out of the simulation domain. The overall trend in the mass evolution with f drift seems to be independent of the di ff usion parameter, as can be seen by comparison of Fig. 5 and Fig. E.1. We therefore chose a value of f drift = 0 . 8 for all following simulations.", '4.3. Treatment of planetary gaps': "Since the particle size in TriPoD can only change due to coagulation / fragmentation, we have to introduce an additional source term that takes into account the particle size reduction due to the depletion of large particles via transport. The classic example of such a case is a gap carved by a planet. In this case, dust diffusion into the gap from the outer disk, removal of large grains within the gap, and coagulation determine what particle sizes are present. Due to the pressure maximum at the gap's outer edge, larger grains can no longer drift into the co-orbital region of a planet. However, smaller particles are still di ff used through the pressure maximum by turbulence and thus slightly increase the solid content of the mostly depleted gap. As the largest particles are removed from the gap, we reduce the maximum particle size on the respective dust depletion timescale in TriPoD . For this, we set a lower limit for the fraction of large particles at which  \nwe begin the size reduction process. We define a hypothetical source term for ε 1, which would set a lower limit to ε 1 in terms of a critical fraction of the total dust density d ε 1 d t = f crit( ε 0 + ε 1) -ε 1 ∆ t where ∆ t is the current simulation timestep. Since our goal is to reduce the maximum particle size on the depletion timescale, we set τ depletion = ε 1 / d ε 1 d t and define the respective size reduction rate  \nd a max d t = a max τ depletion 1 -a max a lim ! , (57)  \nwhere a lim = 1 µ m is a minimum size that we define to limit the shrinking to a reasonable value. To retain a meaningful powerlaw exponent for the size distribution during the shrinking process, we set the mass change rate of the large population to  \nd ε 1 d t = ∂ε 1 ∂ a max d a max d t . (58)  \nHere, ∂ε 1 /∂ a max follows analytically from Eq. 25. Equation 57 and Eq. 58 are added as source terms to the respective right-hand side of the conservation equations if the condition for size reduction is fulfilled ( ε 1 < f crit ε tot). In that way we achieve the following:  \n- -Whenever transport is largely dominating over coagulation, thus removing the large particles, we reduce the maximum particle size on the depletion timescale of the large grain population.\n- -We thus set a lower limit for the power-law exponent and therefore retain a physically meaningful size distribution, even in planetary gaps.  \nDetermining the critical value of ε 1 at which we begin the size reduction process is an experimental task. The size distribution  \nFig. 6: Test for the size reduction in planetary gaps. The dashed lines correspond to the maximum grain size of the respective model (see color bars for the model parameters), while the solid lines show the dust column density for each model. Thick black lines correspond to the respective DustPy simulation that we want to reproduce. The particle size is reduced on the dust depletion timescale when the large dust population makes up less than the critical fraction f crit of the total dust density. We run these simulations to determine which value leads to the best agreement with the DustPy simulation.  \n<!-- image -->  \nwithin a planetary gap is typically not top-heavy due to the e ffi -cient removal of the largest grains (e.g., Dr˛a˙zkowska et al. 2019). The large mass bin should thus contain ≲ 50 % of the total dust density in the gap.  \nWe set up a one-dimensional simulation with a gap corresponding to a Jupiter-mass planet, following the description of Du ff ell (2020) in TriPoD and DustPy . The gap is not evolving in time, but pre-defined in the disk's gas structure and the initial dust density structure. We ran simulations for several values of the critical mass fraction f crit to determine which values give us the best fit between TriPoD and DustPy in a planetary gap. For the limiting particle size, we set either 1 µ m or 10 µ m. The resulting dust column density profiles are shown in Fig. 6. As can be seen, we reached a reasonably good fit for a value of f crit = 0 . 425 and a limiting particle size of a lim = 1 µ m. We have therefore chosen these values for all following simulations.", '5. Test Simulations': 'Here, we present a series of test cases of our TriPoD model for comparison with DustPy . Our fiducial model is depicted in Fig. 1. We set up a smooth protoplanetary disk around a solarmass star for this. The disk structure is identical to the simulations in the calibration runs. The first two panels in each row show the snapshots of the dust size distributions in DustPy (top row) and TriPoD after 20 500 yr and 3 × 10 6 yr of evolution. For this run, gas evolution has been turned o ff to get a more direct comparison between the dust evolution and transport models. The shape of the size distribution is well reproduced in the drift-limited and fragmentation-limited cases.', '5.1. Different stellar masses': 'We tested the accuracy of our model in simulations of systems with di ff erent stellar masses and disk masses. For this, we took the stellar properties of stars with M ∗ = 0 . 1, 0 . 3, 0 . 5, 0 . 7 and 0 . 9 M ⊙ from the pre-main sequence evolution tracks of Bara ff e et al. (2015) at ∼ 10 5 yr of evolution (their first recorded snapshot). We kept the overall disk structure constant for all setups with the values from Table 1 and set the disk mass to 0 . 05 M ∗ . We ran the simulations for all stellar masses once with dust diffusion and a di ff usion parameter of δ = 10 -3 and once without dust di ff usion. Gas evolution was turned o ff in these simulations, as we were only interested in the accuracy of our dust evolution model. We present three of the simulations with di ff usion in Fig. 7, where the first two panels in each row show snapshots of the simulations. The respective DustPy model is shown on top and the TriPoD results below. We found that the overall size distribution evolution is well reproduced for most of the simulation domain and runtime. Our model is also able to capture the shape of the size distributions, as can be seen in the right panels. In the early stages of dust evolution (left panel of each row in Fig. 7), the dust size distribution is in coagulation-fragmentation equilibrium, which results in a typical power-law shape of the distribution with exponent p = -3 . 5. In the case of the 0 . 1 M ⊙ star system, the dust is still in its initial growth phase at the time of the earlier snapshot. Therefore, the distribution has not yet fully reached the equilibrium state and is still more top-heavywhich is typical for the growth phase in which smaller particles are swept up by bigger particles. In this phase, the distribution is not so well reproduced by a power law. Our three-parameter model nonetheless captures the steeper slope well around the size distribution peak, which contains most of the mass.  \nWe present a more detailed look at the column density and mass evolution in Fig. 8 (simulations without di ff usion can be found in Fig. E.2). The relative mass error with respect to the full coagulation model is always ≲ 20 % until ∼ 90 % of the mass has been accreted. The small errors at the beginning of the simulations are seen to add up in the late stages of disk evolution, where the absolute relative mass error in simulations is ≲ 35 % at the end of the simulation when more than 99 % of the mass has drifted out of the simulation domain after 3 × 10 6 yr.  \nFor stellar masses above 0 . 3 M ⊙ , the errors are even smaller. Here, we determine a deviation of less than 10 % until ∼ 90 % of the dust mass has drifted out of the simulation domain. Afterward, in the final stages of dust mass evolution, absolute errors increase again and are generally < 30 %. Only for the system with the least massive star, we measured a maximum deviation of ∼ 35 % at the very end of the dust mass evolution. One reason for this could be the smaller disk size, which means the dust growth front reaches the outer edge of the disk earlier. This area of strong radial dust-to-gas ratio gradients seems to be the origin of most of the deviations, as can be seen in the upper right panel of Fig. 8. Overall, the mass evolution in disks around stars of various masses is very well reproduced.', '5.2. Varying fragmentation velocity': 'To assess the accuracy of our model for di ff erent fragmentation velocities we chose a simple setup with a temperature-dependent  \nFig. 7: Comparison of the dust size distributions in three simulations with di ff erent stellar masses and di ff erent disk masses. The first row in each group depicts a DustPy simulation (full treatment of coagulation). The lower rows show the respective TriPoD simulations with the PLUTO code. Panels on the right depict the local dust size distributions at 10 au at both snapshots. The gray lines in the right panel follow power laws with q = -3 . 5 and q = -3 . 0.  \n<!-- image -->  \nFig. 8: Comparison between DustPy and our model in setups with di ff erent stellar masses. Solid lines show the results of our TriPoD simulations and dashed lines show the respective DustPy simulations. The upper row shows a time series of the dust column density evolution in three snapshots. In the lower row, we show the mass evolution and the errors with respect to the full coagulation model DustPy .  \n<!-- image -->  \nfragmentation velocity. Following 2  \nv frag =                  100 cm s -1 for T > 150 K (ices evaporated) 1000 cm s -1 for 150 K > T > 80 K (H2O ice) 700 cm s -1 for 80 K > T > 44 K (NH3 ice) 100 cm s -1 for 44 K > T (CO2 ice) , (59)  \nwe assumed particles covered in water ice to be the stickiest, followed by NH3 ice and bare silicates or CO2 ice. We used this highly simplified model only to illustrate that TriPoD accurately reproduces the DustPy results for setups with di ff erent fragmentation velocities. We display the resulting dust size distributions for two snapshots of the simulations in Fig. 9. In the disk regions within ∼ 2 au, temperatures exceed 150 K and the fragmentation velocity is accordingly low since all ices are assumed to be evaporated. Thus particles only grow up to millimeter sizes. In the regions beyond the water snow line, we find the two regions with larger fragmentation velocities; first the region with water ice, then, outside of ∼ 8 au the region with ammonia ice which we assigned a slightly lower fragmentation velocity. The outer regions of the disk beyond the CO2 snow line have again been given the lowest fragmentation velocity and have thus also lower maximum particle sizes. As can be seen in Fig. 9, the overall evolution of the disk and the local dust size distributions in the distinct regions of our DustPy and TriPoD models show excellent agreements in dust column densities and maximum particle sizes.', '5.3. Planetary gap': "As a next test case, we conducted a simulation with a planetinduced gap. For the gap in our one-dimensional DustPy and TriPoD simulations we used the same gap profile as in Section 4.3. Instead of letting the disk viscously evolve, we again turned o ff gas evolution for these simulations and imposed the gap profile immediately on the disk's column density structure. Wecompare the resulting size distributions in Fig. 10 with a similar DustPy model in the top row. As can be seen, dust particles collect in the pressure maximum outside of the planetary gap. The gap itself becomes depleted of dust, as the particles drift out and cannot be replenished due to the e ff ect of the pressure bump. As a result, the maximum particle size in the full coagulation model is reduced within the gap. As large grains are removed, the maximum particle size is also reduced in TriPoD , and thus a reasonable size distribution is retained. In that way, a good fit with the particle size in the gap of the full coagulation model is achieved. The size distribution in the full coagulation model is no longer top-heavy inside the gap but declines toward smaller sizes. The same is achieved in TriPoD with our size reduction method, as seen in the right panels of Fig. 10. Particles collecting in the pressure bump outside of the gap are fragmentation-limited and can thus grow to larger sizes than in the outer regions, which are drift-limited. Most mass from the outer disk collects around the pressure maximum and does not reach the inner disk.  \nThe size distribution in the regions inside of the gap becomes drift-limited and dust densities decrease significantly due to the decreased inflow from the outer regions. Only di ff usion of small particles through the gap maintains some particle flow in the inner disk. For comparison, we also show the dust column densities and maximum particle sizes calculated with the old two-pop-py method in Fig. 10. As can be seen, the new  \nFig. 9: Comparison between a DustPy simulation and a PLUTO simulation with our TriPoD dust coagulation model of a protoplanetary disk with radially varying fragmentation velocity (see Eq. 59). The left and center columns depict the global dust distribution at two di ff erent snapshots of the simulations, with the DustPy simulation in the top row and the TriPoD simulation at the bottom. The right plot shows the local dust size distributions in the second snapshot at the locations of the vertical lines of the same color in the central panel. The gray lines in the right panel follow power laws with q = -3 . 5 and q = -3 . 0.  \n<!-- image -->  \nmethod agrees much better with DustPy . In the two-pop-py simulation, we find a much narrower peak in the dust densities around the pressure maximum due to the lack of a small particle population. Grain sizes are also underestimated in the gap, as two-pop-py sets the maximum particle size to the drift limit. This assumes an equilibrium between transport and coagulation which is not given in the gap, where the present particles are di ff using in from the outer edge and large particles are quickly removed by drift.", '5.4. Two-dimensional planet-disk simulation': "As a last test case, we ran a two-dimensional simulation of a protoplanetary disk with a Jupiter-mass planet. We set up a simulation in polar coordinates. The domain is spanning 4-34 au radially and full 2 π azimuthally at a resolution of 1024 cells in the radial and azimuthal direction. The planet is represented by an additional gravitational potential following  \nΦ p = -GM p p d 2 + r 2 sm , (60)  \nwhere M p is the planets mass, d is the distance to the planet, and r sm = 0 . 7 H is the gravitational smoothing length. The full gravitational potential in our simulation domain is then given as  \nΦ tot = Φ ∗ + Φ p , (61)  \nwhere Φ ∗ is the gravitational potential of a Solar-mass star. The disk was set up as a radial power law in column density and temperature with an isothermal equation of state. Details on the disk structure and simulation setup can be found in Table 4. We employed a viscosity ν = α c s H , with α = 10 -3 . The hydrodynamic equations were solved with the HLL Riemann solver, using the third-order accurate Runge-Kutta scheme for the time integration and piece-wise-polynomial spatial reconstruction scheme to the fifth-order.  \nFor the gas velocities, we used a zero-gradient boundary condition at the inner boundary, where we kept the gas density and  \nthe azimuthal velocity fixed to the initial values in the ghost cells. Similar boundary conditions were applied at the outer domain edge, with the di ff erence that we allowed for outgoing velocities but not for inflow. The dust densities at the outer boundary were also kept at the initial value for 1000 orbital time scales after which we began to decrease them exponentially on a timescale of 1000 orbits in order to simulate a reduced dust inflow.  \nWecompare the results to a simulation with the two-pop-py model, which only has one dust fluid. The comparison is presented in Fig. 11, where the top row shows the TriPoD simulation, and the middle row shows two-pop-py after 2000 planetary orbits of evolution. In the lower panels, we present the respective quantities' azimuthal averages, for which we have masked out the region around the planet that has been marked with a white circle.  \nThe di ff erent grain sizes and the redistribution of mass because of fragmentation and coagulation between the two populations lead to di ff erent dust density distributions in the simulation with TriPoD compared to the simulation with two-pop-py . The grains in the two-pop-py simulation collect in two narrowly peaked overdensities at the pressure maximum outside of the planetary gap and a weak second pressure perturbation. The gap is more depleted than in the simulation with the TriPoD model. Furthermore, dust densities in two-pop-py are strongly enhanced in the weak second pressure perturbation that has almost no visible e ff ect in the simulation with TriPoD . The reason for this is the absence of a smaller, separately advected grain population in two-pop-py for which the e ff ect of trapping would be much weaker. This is accounted for in the TriPoD simulation, where the smaller dust population broadens the dust peaks significantly due to the weaker trapping of small grains. This also has the e ff ect that more dust can di ff use through the gap, which means that the densities inside and in the gap itself are slightly higher in TriPoD than in two-pop-py .  \nInside the gap, TriPoD limits the size distribution exponent to a minimum value of q ≈ 4 due to the applied size reduction rates which means the distribution is no longer fragmentation-  \nFig. 10: Comparison between a DustPy simulation and a PLUTO simulation with our TriPoD dust coagulation model in a protoplanetary disk with a pre-defined gap. We also show the particle sizes and dust column density profiles calculated with the old two-pop-py model as dotted lines. The plot on the right-hand side depicts the local dust size distributions at 5 au (within the gap) and at 8 . 5 au (at the pressure maximum), for two di ff erent snapshots. The gray lines in the right panel follow power laws with q = -3 . 5 and q = -3 . 0.  \n<!-- image -->  \nmited but dominated by transport e ff ects. The particles remaining in the gap are accordingly small and the size distributions are no longer top-heavy, as seen in the work by Dr˛a˙zkowska et al. (2019). The two-pop-py model instead assumes the drift limit in the gap and inside of it. The size distribution, assumed by two-pop-py is thus vastly di ff erent from the almost flat or declining distribution seen in TriPoD and in full coagulation models. Similar to the results of Dr˛a˙zkowska et al. (2019), we find that the two-pop-py method still approximates the dust column densities inside of the gap fairly well. The Dust grain sizes are, however, underestimated inside and around the gap and the density profile at the pressure trap is not well reproduced by two-pop-py .  \nOur results are thus qualitatively similar to the conclusions by Dr˛a˙zkowska et al. (2019). Coagulation makes it possible for small grains to pass through the gap via di ff usion and coagulate again inside of the planet's orbit. The dust accumulates in the pressure bump, but the peak in the dust density is broader than in a simulation without coagulation, due to the presence of small grains, which are constantly produced in our simulation as a result of fragmentation.", '6. Discussion': 'In its current form in the PLUTO code, our model operates in the terminal velocity approximation. The dust fluids are implemented as passive tracers with modified fluxes. Thus, feedback from the dust on the gas is not included in our two-dimensional test simulation. However, as shown by Dr˛a˙zkowska et al. (2019), the dust feedback in their two-dimensional planet-disk simula-  \nTable 4: Initial conditions for our two-dimensional PLUTO simulation with TriPoD .Notes. Identical to the initial conditions in Dr˛a˙zkowska et al. (2019).  \non only had a minor e ff ect on the simulation outcome. E ff ects like the streaming instability (Youdin & Goodman 2005), cannot be simulated with our current version of TriPoD due to this limitation of our very simple dust fluid implementation. Unfortunately, PLUTO does currently not have a multi-fluid feature. A  \n]  \n2  \n-  \ncm  \ng  \n[  \nj  \n〉  \ng  \nΣ  \n〈  \nFig. 11: Comparison between two-dimensional PLUTO simulations with the TriPoD coagulation model and with the two-pop-py model. Both simulations have the exact same gas density structure and are only distinguished by the applied dust coagulation model. The lower panels show the azimuthally averaged quantities. For the azimuthal averages, we have masked the region around the planet that has been marked with the white circle. The red line in the left panel shows the radial pressure gradient parameter  \n<!-- image -->  \nη = 1 2 GLYPH<16> H R GLYPH<17> 2 ∂ log P ∂ log R . It can be seen that in our simulation a weak second pressure maximum appears at ∼ 20 au, which leads to strong trapping in the two-pop-py model (Birnstiel et al. 2012). This e ff ect is much less pronounced if a polydisperse size distribution is considered as done in our TriPoD simulation, where the small dust population is less e ffi ciently trapped.  \nfuture iteration of the model could be combined with a code that supports dust fluids that are treated more self-consistently, like FARGO3D (Benítez-Llambay & Masset 2016; Benítez-Llambay et al. 2019) or Athena++ (Stone et al. 2020; Huang & Bai 2022). The terminal velocity approach also means that e ff ects occurring on timescales shorter than the friction timescale of the particles cannot be modeled in the current form of TriPoD . Possible cases could be the e ff ects of shock waves, or spiral density waves.  \nAnother shortcoming arises from the way we calculate the dust fluxes. As the maximum particle size changes throughout a simulation domain, neighboring cells usually do not have the  \nsame particle sizes. Therefore, the bin interfaces separating Σ 0 and Σ 1 are also di ff erent between neighboring cells. Our advection scheme does not take this into account in its current form, where transport is only occurring between the same size bin and not across bins. Due to the smooth variation of a max throughout simulations of protoplanetary disks, neighboring cells will typically still have very similar maximum particle sizes. The effect of this inaccuracy in the advection scheme is therefore likely small. However, one could construct extreme cases in which the error would be large, for example, if two neighboring cells had  \nvastly di ff erent maximum particle sizes. We will have to address this issue in a future version of the model.  \nDetails of the grain size distribution that are modeled in fullfledged coagulation models like DustPy cannot be reproduced in our simplified prescription. To still achieve a good fit with such models, our model has several calibration factors that have been adjusted to reach a good fit with DustPy . However, this treatment of dust coagulation cannot reproduce the finer details of the coagulation process.  \nOur model omits the e ff ects of bouncing (see Dominik & Dullemond 2024). This e ff ect can lead to a steep, top-heavy size distribution and smaller maximum particle sizes. Bouncing could be implemented in a future version of TriPoD .  \nTriPoD , as well as DustPy are designed to simulate the dust evolution in vertically integrated models of protoplanetary disks, assuming settling-mixing equilibrium at all times. E ff ects like sedimentation-driven coagulation (Zsom et al. 2011; Krijt & Ciesla 2016) can thus not be modeled by DustPy or TriPoD in its current form. This is however not a fundamental limitation of TriPoD , which could be adapted to work in three-dimensional setups in a future version.', '7. Summary and conclusion': "We present TriPoD , an accurate and computationally inexpensive sub-grid model for dust coagulation in vertically integrated hydrodynamic simulations of protoplanetary disks. TriPoD only requires two dust fluids and a tracer fluid representing the maximum particle size to model the evolution of a polydisperse dust size distribution. This makes it possible to run inexpensive simulations of planet-disk systems, protoplanetary disks with vortices, etc., and deduce the particle properties and dust densities.  \nWe have implemented our model in the PLUTO code. The workflow of the model during one timestep of a hydrodynamic simulation is as follows:  \n- 1. PLUTO calculates the gas fluxes across the cell interfaces with one of the standard Riemann solvers and reconstruction methods that are provided with the code.\n- 2. The dust fluxes and the particle size flux are calculated (Eq. 44 to Eq. 49). Flux velocities are calculated from interpolated interface values. The drift flux is then taken upstream depending on the limited dust drift velocity. For dust di ff usion, we calculate the limited di ff usion flux in a manner identical to the DustPy model. The di ff usion flux is then added to the total dust flux. (modifications in the PLUTO code made in file adv\\_flux.c )\n- 3. The relative particle velocity ∆ v max is calculated for the particle sizes a max and f ∆ v a max. It is used to determine the particle growth rate (Eq. 29).\n- 4. The source terms for the large and the small dust fluids, are calculated. These fragmentation and sweep-up rates (Eqs. 32 and 33) are then added to the source terms. If the large dust population is depleted to less than f crit ε tot, we calculate the respective dust depletion time and add the size reduction terms to the source terms in a way that conserves the current power law (Eq. 57). Mass is shifted accordingly from the small to the large population (Eq. 58). The source terms of the two dust populations are calculated according to (Eq. 35) and added to the equations' right-hand side (modification in PLUTO made in file rhs\\_source.c ).\n- 5. The array of conserved quantities is advanced by one timestep by PLUTO using a standard time integration scheme provided with the code.  \nAlthough TriPoD is a highly simplified model, it can predict the dust mass evolution accurately for millions of years of evolution and simulate the e ff ects of disk sub-structures on the dust size distribution, as demonstrated in our tests and calibration runs. Applications of this first version of TriPoD could include:  \n- -More accurate studies of chemical networks in protoplanetary disks, which are highly dependent on the available grain surface area.\n- -Parameter studies of planet-disk systems with dust coagulation, that were so far infeasible due to the high computational cost of full coagulation models.\n- -Better radiative transfer post-processing of simulations, given the knowledge of the grain size distribution.\n- -Simulations with self-consistently calculated thermal relaxation times.  \nIn the future, we plan to extend this first version of the model to three-dimensional simulations of protoplanetary disks.", 'Acknowledgments': "T.P., H.K., and T.B. acknowledge the support of the German Science Foundation (DFG) priority program SPP 1992 'Exploring the Diversity of Extrasolar Planets' under grant Nos. BI 1816 / 72 and KL 1469 / 16-1 / 2. T.B. acknowledges funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program under grant agreement No 714769 and funding by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under grants 361140270, 325594231, and Germany's Excellence Strategy EXC-2094 - 390783311. Computations were conducted on the computational facilities of the University of Munich (LMU). We thank the anonymous referee for their constructive comments, which helped us to improve the quality of this article.", 'References': 'Andrews, S. M., Huang, J., Pérez, L. M., et al. 2018, ApJ, 869, L41  \nAntonellini, S., Kamp, I., & Waters, L. B. 2023, A&A, 672, A92  \nBalbus, S. A. & Hawley, J. F. 1991, ApJ, 376, 214  \n- Bara ff e, I., Homeier, D., Allard, F., & Chabrier, G. 2015, A&A, 577  \nBaruteau, C., Barraza, M., Pérez, S., et al. 2019, MNRAS, 486, 304  \n- Benítez-Llambay, P., Krapp, L., & Pessah, M. E. 2019, Astrophys. J. Suppl. Ser., 241, 25\n- Benítez-Llambay, P. & Masset, F. S. 2016, Astrophys. J. Suppl. Ser., 223, 11 Bergez-Casalou, C., Bitsch, B., Kurtovic, N. T., & Pinilla, P. 2022, A&A, 659, A6\n- Binkert, F. 2023, MNRAS, 000, 1\n- Birnstiel, T. 2023, Annu. Rev. Astron. Astrophys. [ arXiv:2312.13287 ]\n- Birnstiel, T., Dullemond, C. P., & Brauer, F. 2009, A&A, 503, 5\n- Birnstiel, T., Dullemond, C. P., & Brauer, F. 2010, A&A, 513, 79  \nBirnstiel, T., Dullemond, C. P., Zhu, Z., et al. 2018, ApJ, 869, L45  \nBirnstiel, T., Klahr, H., & Ercolano, B. 2012, A&A, 539, 148  \nBirnstiel, T., Ormel, C. W., & Dullemond, C. P. 2011, A&A, 525, 11  \nBlandford, R. D. & Payne, D. G. 1982, MNRAS, 199, 883  \n- Brauer, F., Dullemond, P., & Henning, T. 2008, A&A, 480, 859  \nBurn, R., Emsenhuber, A., Weder, J., et al. 2022, A&A, 666  \nCharnoz, S. & Taillifet, E. 2012, ApJ, 753, 119  \n- Dominik, C. & Dullemond, C. P. 2024, A&A, 682, A144\n- Dr˛a˙zkowska, J., Bitsch, B., Lambrechts, M., et al. 2023, in Protostars Planets VII, ed. S.-i. Inutsuka, Y. Aikawa, T. Muto, K. Tomida, & M. Tamura, Vol. 534 (San Francisco: ASPCS), 717-759\n- Dr˛a˙zkowska, J., Li, S., Birnstiel, T., Stammler, S. M., & Li, H. 2019, ApJ, 885, 91\n- Dr˛a˙zkowska, J., Stammler, S. M., & Birnstiel, T. 2021, A&A, 647, 15  \nDubrulle, B., Morfill, G., & Sterzik, M. 1995, Icarus, 114, 237  \nDu ff ell, P. C. 2020, ApJ, 889, 16  \n- Dullemond, C. P., Birnstiel, T., Huang, J., et al. 2018, ApJ, 869, L46\n- Emsenhuber, A., Mordasini, C., Burn, R., et al. 2021, A&A, 656, 69  \nEpstein, P. S. 1924, Phys. Rev., 23, 710  \nEstrada, P. R. & Cuzzi, J. N. 2008, ApJ, 682, 515 Fromang, S. & Nelson, R. P. 2009, A&A, 496, 597 Gárate, M., Delage, T. N., Stadler, J., et al. 2021, A&A, 655, 18 Guillet, V., Hennebelle, P., Pineau Des Forêts, G., et al. 2020, A&A, 643, 17 Güttler, C., Blum, J., Zsom, A., Ormel, C. W., & Dullemond, C. P. 2010, A&A, 513, A56  \nHuang, P. & Bai, X.-N. 2022, Astrophys. J. Suppl. Ser., 262, 11  \nKobayashi, H., Tanaka, H., & Okuzumi, S. 2016, ApJ, 817, 105  \nKolmogorov, A. N. 1941, Dokl. AN SSSR, 30, 299  \nKornet, K., Stepinski, T. F., & Rózyczka, M. 2001, A&A, 378, 180  \nKrijt, S. & Ciesla, F. J. 2016, ApJ, 822, 111  \n- Lau, T. C. H., Dr˛a˙zkowska, J., Stammler, S. M., Birnstiel, T., & Dullemond, C. P. 2022, A&A, 668, 170  \nLeiendecker, H., Jang-Condell, H., Turner, N. J., & Myers, A. D. 2022, ApJ, 941, 172  \nLenz, C. T., Klahr, H., Birnstiel, T., Kretke, K., & Stammler, S. 2020, A&A, 640 Lesur, G., Ercolano, B., Flock, M., et al. 2023, in Protostars Planets VII, ed.  \nS.-i. Inutsuka, Y. Aikawa, T. Muto, K. Tomida, & M. Tamura, Vol. 534 (San  \nFrancisco: ASPSC), 465-500  \nLevermore, G. D. & Pomraning, G. C. 1981, ApJ, 248, 321  \nLichtenberg, T., Dr˛a˙zkowska, J., Schönbächler, M., Golabek, G. J., & Hands,  \nT. O. 2021, Science (80-. )., 371, 365  \nLombart, M., Bréhier, C.-E., Hutchison, M., & Lee, Y.-N. 2024, MN-  \narXiv:2404.11851  \nRAS[  \n]  \nLombart, M., Hutchison, M., & Lee, Y. N. 2022, MNRAS, 517, 2012  \nLombart, M. & Laibe, G. 2021, MNRAS, 501, 4298  \n- Marchand, P., Guillet, V., Lebreuilly, U., & Mac Low, M. M. 2021, A&A, 649, 50\n- Marchand, P., Guillet, V., Lebreuilly, U., & Mac Low, M. M. 2022, A&A, 666, 27\n- Marchand, P., Lebreuilly, U., Mac Low, M. M., & Guillet, V. 2023, A&A, 670, 61\n- Mathis, J. S., Rumpl, W., & Nordsieck, K. H. 1977, ApJ, 217, 425\n- Mignone, A., Bodo, G., Massaglia, S., et al. 2007, Astrophys. J. Suppl. Ser., 170, 228  \nMuley, D., Fuksman, J. D. M., & Klahr, H. 2023, A&A, 678, A162  \nNakagawa, Y., Nakazawa, K., & Hayashi, C. 1981, Icarus, 45, 517  \nNakagawa, Y., Sekiya, M., & Hayashi, C. 1986, Icarus, 67, 375  \n- Öberg, K. I., Guzmán, V. V., Walsh, C., et al. 2021, Astrophys. J. Suppl. Ser., 257, 1  \nOrmel, C. W. & Cuzzi, J. N. 2007, A&A, 466, 413  \nOrmel, C. W., Cuzzi, J. N., & Tielens, A. G. G. M. 2008, ApJ, 679, 1588  \nPascucci, I., Cabrit, S., Edwards, S., et al. 2023, in Protostars Planets VII, ed.  \n- S.-i. Inutsuka, Y. Aikawa, T. Muto, K. Tomida, & M. Tamura (San Francisco: ASPCS), 567-604  \nPérez, L. M., Benisty, M., Andrews, S. M., et al. 2018, ApJ, 869, L50  \nPfeil, T., Birnstiel, T., & Klahr, H. 2023, ApJ, 959, 121  \n- Pfeil, T., Cranmer, M., Ho, S., et al. 2022, in Mach. Learn. Phys. Sci. Work. NeurIPS  \nPinilla, P., Lenz, C. T., & Stammler, S. M. 2021, A&A, 645, 70  \nPowell, D., Murray-Clay, R., Pérez, L. M., Schlichting, H. E., & Rosenthal, M. 2019, ApJ, 878, 116  \nSchlichting, H. E. & Sari, R. 2011, ApJ, 728, 68  \nSimon, J. B., Blum, J., Birnstiel, T., & Nesvorný, D. 2022, Accept. Publ. Comets III [ arXiv:2212.04509 ]  \nSmoluchowski, M. 1916, Phys. Z., 17, 557  \nStammler, S. M. & Birnstiel, T. 2022, ApJ, 935, 35  \n- Stone, J. M., Tomida, K., White, C. J., & Felker, K. G. 2020, Astrophys. J. Suppl. Ser., 249, 4  \nTamfal, T., Dr˛a˙zkowska, J., Mayer, L., & Surville, C. 2018, ApJ, 863, 97  \nTsukagoshi, T., Nomura, H., Muto, T., et al. 2022, ApJ, 928, 49  \nTsukamoto, Y. & Okuzumi, S. 2022, ApJ, 934, 88  \nVorobyov, E. I., Akimkin, V., Stoyanovskaya, O., Pavlyuchenkov, Y., & Liu, H. B. 2018, A&A, 614, 98  \nVorobyov, E. I. & Elbakyan, V. G. 2019, A&A, 631, A1  \nVorobyov, E. I., Matsukoba, R., Omukai, K., & Guedel, M. 2020, A&A, 638, A102  \nVorobyov, E. I., Skliarevskii, A. M., Elbakyan, V. G., et al. 2019, A&A, 627, A154  \nWeidenschilling, S. J. 1977, MNRAS, 180, 57  \nWeidenschilling, S. J. 1980, Icarus, 44, 172  \nWetherill, G. W. & Stewart, G. R. 1989, Icarus, 77, 330  \nWhipple, F. 1972, in From Plasma to Planet, ed. A. Evlius (New York: Wiley Interscience Division), 211  \nYoudin, A. N. & Goodman, J. 2005, ApJ, 620, 459  \nZsom, A., Ormel, C. W., Dullemond, C. P., & Henning, T. 2011, A&A, 534, A73 Zsom, A., Ormel, C. W., Güttler, C., Blum, J., & Dullemond, C. P. 2010, A&A, 513, A57', 'General symbols': 'Table A.1: Lists of symbols', 'Appendix B: Column density formulation': 'Assuming σ 01 , 11, ∆ v 01 , 11, m 0 , 1, and F to be vertically constant we deduce the vertically integrated mass exchange rates (Eq. 32 and Eq. 33)  \n˙ Σ d , 0 → 1 = Σ 0 Σ 1 σ 01 ∆ v 01 m 12 π H 0 H 1 Z ∞ -∞ exp       -z 2 2       H 2 0 + H 2 1 H 2 0 H 2 1             d z (B.1)  \n= Σ 0 Σ 1 σ 01 ∆ v 01 m 1 q 2 π ( H 2 0 + H 2 1 ) (B.2)  \n˙ Σ d , 1 → 0 = Σ 2 1 σ 11 ∆ v 11 m 12 π H 2 ˜ F Z ∞ -∞ exp       -z 2 H 2       d z  \n= Σ 1 σ 11 ∆ v 11 m 1 q 4 π H 2 1 ˜ F  \n1 1 (B.3) 2 (B.4)  \nThus, also F has a modified form in the vertically integrated setup  \n˜ F = s 2 H 2 1 H 2 0 + H 2 1 σ 01 σ 11 ∆ v 01 ∆ v 11 a max a int ! -( q + 4) . (B.5)', 'Appendix C: Transition functions': 'In order to determine the transition functions, we have taken into account the results of local dust coagulation simulations and theoretical models. The transitions should not be too rapid to not cause issues during the numerical integration.  \nFor the transition between the small regime of turbulence and the fully intermediate regime of turbulence we take the approximate transition criterion from Ormel et al. (2008)  \nf turb . 2 turb . 1 = 5 t s τ max ( > 1 small particle regime (turb.1) < 1 fully inter. regime (turb.2), (C.1)  \nwhere τ max is the friction time of the largest particles, and t s = Re -1 / 2 Ω -1 K is the small eddy turnover time. Using the same functional form as for the growth rate Eq. 29, but now from 0 to 1 instead of -1 to 1, we define  \np turb . 1 B 1 2           1 -GLYPH<16> f turb . 2 turb . 1 GLYPH<17> 4 -1 GLYPH<16> f turb . 2 turb . 1 GLYPH<17> 4 + 1           . (C.2)  \nThis expression is approaching 0 if τ max < 5 t s and 1 if τ max > 5 t s. We apply a very similar function to the transition from the drift-dominated regime to the turbulence-dominated regime. We define  \nf turb drift = ∆ v turb ∆ v drift ( > 1 turbulence dominated < 1 drift dominated , (C.3)  \nwhere ∆ v turb is the turbulent collision velocity. We define  \np drift B 1 2           1 -GLYPH<16> f turb drift GLYPH<17> 6 -1 GLYPH<16> f turb drift GLYPH<17> 6 + 1           . (C.4)  \nwhere we found that a slightly steeper transition fits better with the results of local coagulation simulations.  \nFinally, we have to introduce the transition from growth to fragmentation. For this, we chose a steeper transition, following  \np frag = exp        -5 min ∆ v tot v frag , 1 . 0 ! -1 . 0 !! 2        , (C.5)  \nwhich ensures that the equilibrium size distribution is reached on a fast enough timescale.', 'Appendix D: Flux-limited dust diffusion': 'In order to avoid unrealistically large di ff usion fluxes in the presence of strong gradients in the dust-to-gas ratio, we introduce a flux limiter that is conceptually identical to the one by Levermore & Pomraning (1981), which was used by Stammler & Birnstiel (2022). Without the flux limiter, the di ff usion flux is given by  \nF di ff = -D Σ g ∇ ε . (D.1)  \nThe transport velocity shall however not be larger than the turbulent velocities that drive the di ff usion, meaning the maximum allowed flux is  \nF di ff , max B v turb ε Σ g = √ δ c s 1 + St 2 ε Σ g . (D.2)  \nThe flux limiter is now defined via  \nλ lim = 1 + χ 1 + χ + χ 2 , (D.3)  \nwhere  \nχ = F di ff F di ff , max , (D.4)  \nis the ratio of the di ff usion flux and the maximum flux. The limited flux is then given by  \nF di ff , lim = λ lim F di ff . (D.5)  \nNote that although the flux is a vector (with components for each direction), all operations here are component-wise. In PLUTO , we define all quantities in the above equations at the cell interfaces.', 'Appendix E: Calibrations and test simulations without diffusion': 'Here, we present the same calibration and test simulations shown in the main part of this article, but now without dust di ff usion ( δ = 0).  \n2  \nFig. E.1: Comparison between DustPy and our model in a setup without dust di ff usion and with di ff erent drift calibration factors f drift. Solid lines show the results of our TriPoD calibration runs and dashed lines show the respective DustPy simulation, to which we calibrate our model. The upper row shows a timeseries of the dust column density evolution in three snapshots. In the lower row, we show the mass evolution and the errors with respect to the full coagulation model DustPy . For a factor of f drift = 0 . 8, the mass evolution of the full coagulation model is well reproduced by our three-parameter model.  \n<!-- image -->  \n2  \nFig. E.2: Comparison between DustPy and our model in setups with di ff erent stellar masses without dust di ff usion. Solid lines show the results of our TriPoD simulations and dashed lines show the respective DustPy simulations. The upper row shows a time series of the dust column density evolution in three snapshots. In the lower row, we show the mass evolution and the errors with respect to the full coagulation model DustPy .  \n<!-- image -->'}