DustPy: Protoplanetary Dust Evolution Code

• DustPy is a Python package for simulating dust evolution in protoplanetary disks, offering a comprehensive 1D framework that couples gas dynamics, turbulent diffusion, and collisional dust processes.
• It models key mechanisms such as viscous advection, radial drift, coagulation, and fragmentation via the Smoluchowski equation, ensuring realistic dust size distributions.
• Its modular design enables extensions for pressure traps, photoevaporation, and planetesimal formation, making it a benchmark for code-comparison studies in dust evolution.

Searching arXiv for recent and foundational DustPy papers to support the article. Searching arXiv for DustPy foundations, applications, and code-comparison papers. DustPy is a Python package for simulating dust evolution in protoplanetary disks. In its standard usage, it operates in a one-dimensional, radially resolved, vertically integrated framework and solves gas and dust transport together with collisional dust evolution, including viscous advection, diffusion, radial drift, coagulation, and fragmentation (Stammler et al., 2022). Across the recent literature, DustPy functions both as a general-purpose dust-evolution code and as a forward model that links microphysical assumptions about grain growth and transport to observables such as continuum radial profiles, spectral indices, disk sizes, and inferred pebble fluxes (Jiang et al., 2023). Its modular structure has made it a common numerical backbone for studies of turbulence constraints, pressure-trap leakage, compact disks, photoevaporative dispersal, planetesimal formation, and code-comparison benchmarks (Stammler et al., 2022).

1. Definition and scope

DustPy is designed to evolve the dust size distribution in protoplanetary disks while coupling that evolution to the gas disk and to aerodynamic transport processes (Stammler et al., 2022). In the usage documented across the cited studies, its core role is to compute the spatiotemporal evolution of $\Sigma_{\rm g}(r,t)$ and $\Sigma_{\rm d}(r,a,t)$ under viscous gas evolution, dust advection, turbulent diffusion, vertical settling, and the Smoluchowski coagulation–fragmentation equation (Jiang et al., 2023).

The code is routinely described as a “full physics” baseline for one-dimensional dust-evolution calculations because it evolves a polydisperse dust population directly rather than reducing the solids to a small number of representative sizes (Eriksson et al., 23 Mar 2026). This has made it a standard reference in comparisons with reduced models such as two-pop-py and TriPoD, especially in regimes where the size distribution departs strongly from a truncated power law or becomes bimodal (Eriksson et al., 23 Mar 2026, Pfeil et al., 2024).

A recurring theme in the literature is that DustPy is not merely a grain-growth calculator. It is used as a forward model that connects assumptions about fragmentation thresholds, turbulent transport, trapping, and disk thermodynamics to quantities inferred from ALMA, VLA, and synthetic radiative-transfer analyses, including maximum grain sizes, continuum radii, spectral indices, and radial brightness profiles (Jiang et al., 2023, Kurtovic et al., 12 Jun 2025).

2. Core equations and physical model

DustPy solves the coupled radial evolution of gas and dust in a viscously evolving disk while tracking the dust size distribution via the Smoluchowski coagulation equation (Jiang et al., 2023). In the studies summarized here, the gas is commonly evolved with an $\alpha$-viscosity prescription, $\nu=\alpha c_s H$, and the dust is transported with a size-dependent advection–diffusion equation that includes both drag-driven drift and turbulent diffusion (Jiang et al., 2023, Kurtovic et al., 12 Jun 2025).

The dust radial velocity is typically written in the Nakagawa–Takeuchi–Lin form,

$$v_{r,d} \;=\; \frac{v_{r,g}}{1+\mathrm{St}^2} \;-\; \frac{2\,\mathrm{St}}{1+\mathrm{St}^2}\,\eta\,v_K,$$

with $\eta$ set by the local pressure gradient and $\mathrm{St}$ the Stokes number (Jiang et al., 2023, Kurtovic et al., 12 Jun 2025). In the Epstein regime, many of the applications use

$$\mathrm{St} \;=\; \frac{\pi}{2}\,\frac{a\,\rho_\bullet}{\Sigma_g},$$

or an equivalent piecewise Epstein/Stokes expression when larger particles are present (Jiang et al., 2023, Kurtovic et al., 12 Jun 2025).

Turbulent diffusion is generally modeled as

$D_d \;=\; \frac{\delta\,c_s H}{1+\mathrm{St}^2},$

with several studies setting the radial and vertical dust diffusivities equal to the adopted turbulence parameter, while others explicitly decouple dust diffusivity from the gas accretion viscosity (Jiang et al., 2023, Houge et al., 13 Apr 2026). Vertical settling is usually represented through an equilibrium dust scale height such as

$$H_d \;\approx\; H_g\sqrt{\frac{\alpha}{\alpha+\mathrm{St}}},$$

or a closely related closure (Eriksson et al., 23 Mar 2026, Tong et al., 29 Sep 2025).

Collisional evolution is solved with the Smoluchowski equation on a logarithmic mass grid, with relative velocities assembled from Brownian motion, turbulence, radial drift, azimuthal drift, and settling (Jiang et al., 2023, Eriksson et al., 23 Mar 2026). Standard analytic barriers are repeatedly used to interpret DustPy outputs. For fragmentation-limited growth,

$\Sigma_{\rm d}(r,a,t)$0

and

$\Sigma_{\rm d}(r,a,t)$1

while drift-limited estimates are used to diagnose when radial drift rather than fragmentation truncates the size distribution (Jiang et al., 2023, Houge et al., 13 Apr 2026).

This framework makes DustPy especially suitable for problems in which observables depend sensitively on the full dust size spectrum rather than on a single representative grain size. A plausible implication is that its strongest advantage over reduced models emerges when fragmentation-fed small grains, re-coagulation, or multi-modal size distributions control transport or emission.

3. Numerical architecture and configurability

The published applications consistently use DustPy in a one-dimensional radial configuration with logarithmic radial grids and logarithmic mass or size binning (Jiang et al., 2023, Eriksson et al., 23 Mar 2026). However, the exact numerical setup varies substantially by problem. Examples include radial domains from 1 to $\Sigma_{\rm d}(r,a,t)$2 au or beyond, grid sizes from 100 to 500 radial cells, and mass grids spanning from sub-micron grains to decimeter- or meter-scale solids (Tong et al., 6 Feb 2025, Schöll et al., 29 Apr 2026, Ying et al., 26 Jan 2026).

The code’s modularity is a central theme in both the package description and the application literature (Stammler et al., 2022). Studies modify or extend DustPy in several ways:

• by imposing pressure traps through local viscosity modulations or prescribed gap profiles (Jiang et al., 2023, Kurtovic et al., 12 Jun 2025);
• by coupling it to radiative-transfer post-processing through RADMC-3D and, in some cases, DustPyLib (Tong et al., 6 Feb 2025, Tong et al., 29 Sep 2025);
• by adding external photoevaporation modules or stellar X-ray photoevaporation source terms (Delussu et al., 7 Feb 2026, Ying et al., 26 Jan 2026);
• by introducing porosity or fractal-dimension prescriptions as custom mappings from mass to size, bulk density, and cross section (Schöll et al., 29 Apr 2026);
• by implementing a bouncing barrier through custom collision-outcome probabilities (Dominik et al., 2023);
• by adding infall source terms, thermal solvers, and gravitoturbulent transport for embedded Class 0/I disks (Carrera et al., 17 Apr 2025).

The literature also shows that DustPy is frequently embedded in larger workflows rather than used in isolation. In observationally oriented studies it is combined with FRANK-derived radial profiles, opacity models such as DSHARP or Ricci et al. (2010), and synthetic imaging pipelines that produce Band 3, 6, or 7 continuum predictions (Luo et al., 2 Mar 2026, Kurtovic et al., 12 Jun 2025, Tong et al., 29 Sep 2025). In planet-formation studies it is further coupled to parameterized prescriptions for streaming-instability planetesimal formation or to pebble-accretion calculations that use the evolved size distribution as input (Eriksson et al., 23 Mar 2026, Ying et al., 26 Jan 2026).

4. Collision physics, growth barriers, and extensions beyond the default model

A large fraction of DustPy-based work concerns how collision microphysics modifies disk evolution. In many studies, collisions fragment when the relative speed exceeds a prescribed threshold velocity $\Sigma_{\rm d}(r,a,t)$3, with commonly explored values including $\Sigma_{\rm d}(r,a,t)$4, $\Sigma_{\rm d}(r,a,t)$5, $\Sigma_{\rm d}(r,a,t)$6, and $\Sigma_{\rm d}(r,a,t)$7 depending on the assumed material properties (Jiang et al., 2023, Tong et al., 6 Feb 2025, Ying et al., 26 Jan 2026). This parameter controls whether growth remains fragmentation-limited, whether pebbles survive in pressure maxima, and how much small dust is replenished by destructive collisions.

The literature shows that DustPy has been used extensively to explore fragile-aggregate scenarios. In ringed disks, runs with $\Sigma_{\rm d}(r,a,t)$8 and low turbulence naturally produce radially smooth maximum grain-size profiles and favor fragmentation-limited growth even inside pressure maxima (Jiang et al., 2023). In compact-disk models, moderate turbulence beyond a dead zone plus fragile dust prevents mm grains from surviving in the active outer disk, confining mm emission to radii set by the dead-zone edge (Tong et al., 6 Feb 2025).

Bouncing is not part of the simplest sticking–fragmentation interpretation used in many papers, but it has been implemented explicitly in DustPy. The bouncing extension adds a third collision outcome between sticking and fragmentation, using a Maxwell–Boltzmann distribution of relative velocities and a rolling-force-based threshold (Dominik et al., 2023). In those models, bouncing drives the dust population toward a narrow, almost mono-disperse distribution, removes most micrometer-sized grains, modifies settling, and changes the conditions for both streaming instability and pebble accretion (Dominik et al., 2023). A common misconception is that adding bouncing merely lowers the maximum grain size. The DustPy results indicate a stronger effect: the entire shape of the size distribution changes, with consequences for midplane concentration and observables (Dominik et al., 2023).

Porosity and fractal growth have also been incorporated as custom DustPy extensions. One study prescribes the fractal dimension $\Sigma_{\rm d}(r,a,t)$9 as a free parameter and modifies the mapping between particle mass, radius, filling factor, and aerodynamic area while keeping the collision thresholds independent of porosity (Schöll et al., 29 Apr 2026). Under those assumptions, lower $\alpha$0 yields larger masses but does not raise the maximum Stokes number in fragmentation-limited growth; in bouncing-limited growth it lowers the maximum $\alpha$1, making streaming instability less favorable (Schöll et al., 29 Apr 2026). This suggests that DustPy can be used to isolate dynamical consequences of porosity even when the collision kernel itself remains simplified.

5. Pressure traps, leakage, and substructured disks

Pressure maxima are among the most common DustPy use cases. They are imposed in multiple ways: via localized depressions or enhancements in $\alpha$2, via analytic gap profiles based on Kanagawa et al. fits, or through torque-injection schemes designed to enforce target gap structures (Jiang et al., 2023, Stammler et al., 2023, Delussu et al., 7 Feb 2026). In all cases, the purpose is similar: modify the local pressure gradient so that large drifting pebbles are trapped while smaller, more gas-coupled grains may continue to cross the gap.

DustPy has been used to show that such traps are often leaky. In a solar-nebula context, simulations of an early Jupiter-induced gap found that particles trapped at the pressure bump fragment into smaller grains that diffuse through the gap, contaminating the inner disk on timescales inconsistent with the meteoritic CC/NC dichotomy (Stammler et al., 2023). Even a $\alpha$3 gap leaked a large fraction of the outer-disk dust mass in those runs, and low turbulent diffusivity reduced but did not eliminate the problem in time-dependent proto-Jupiter growth scenarios (Stammler et al., 2023).

A later large parameter study generalized this result and argued that most outer traps are leakier than previously thought when the full size distribution is evolved with DustPy (Houge et al., 13 Apr 2026). Leakage is quantified by comparing the cumulative mass crossing the gap in simulations with and without a planet-induced trap. Across much of the explored parameter space, especially for $\alpha$4 or $\alpha$5, the blocking efficiency remains low, because trapped pebbles continually fragment into small grains that diffuse and advect across the gap (Houge et al., 13 Apr 2026). Highly blocking traps occur only under low viscosity and weak turbulence, and in that regime the strongest suppression of the inner pebble flux comes primarily from planetesimal formation in the trap rather than from the trap alone (Houge et al., 13 Apr 2026).

DustPy has also been used to investigate more global consequences of substructures. In the AGE-PRO modeling, weak or strong traps are favored over smooth disks for reproducing the observed 1.3 mm fluxes, effective radii, and spectral-index distribution, while smooth disks drift-deplete too rapidly (Kurtovic et al., 12 Jun 2025). In a separate photoevaporation study, primordial substructures retained enough dust to keep dispersing disks millimeter-bright, whereas initially smooth disks evolved toward the millimeter-faint transition-disk population (Gárate et al., 2023).

6. Applications to turbulence, photoevaporation, and planetesimal formation

DustPy has become a major tool for inferring turbulence and testing planetesimal-formation pathways. One direct inversion method uses observed maximum grain sizes in rings, together with the assumption that fragmentation rather than drift limits growth inside pressure maxima, to estimate the midplane turbulence coefficient $\alpha$6 (Jiang et al., 2023). Applied to seven disks, this framework found low turbulence coefficients, typically $\alpha$7 in five systems when $\alpha$8, with IM Lup standing out at $\alpha$9 and HL Tau showing an increasing $\nu=\alpha c_s H$0 toward larger radii (Jiang et al., 2023).

In broader model grids, DustPy also links turbulence and dust fragility to multi-wavelength ALMA observables. One study found that only two families of parameter combinations reproduce the targeted ring morphologies: fragile dust with $\nu=\alpha c_s H$1–$\nu=\alpha c_s H$2 in disks with $\nu=\alpha c_s H$3, or more resilient dust with $\nu=\alpha c_s H$4–$\nu=\alpha c_s H$5 in disks with $\nu=\alpha c_s H$6 (Tong et al., 29 Sep 2025). In those models, successful cases required the observed rings to be optically thick at both 1.3 and 3 mm, and the inferred maximum grain sizes in the emitting layer were smaller than the true midplane maximum sizes because only small grains were lifted to the emitting surface (Tong et al., 29 Sep 2025).

Photoevaporation is another major application area. DustPy has been coupled to XEUV or X-ray photoevaporation to study how clearing fronts, cavity edges, and associated pressure maxima affect the solids (Gárate et al., 2023, Ying et al., 26 Jan 2026). In the transition-disk context, the simulations indicate that dust trapping determines the millimeter flux while photoevaporation controls cavity opening and expansion, effectively decoupling brightness from cavity size (Gárate et al., 2023). In a late-stage planetesimal-formation context, X-ray photoevaporation creates an expanding pressure maximum at the cavity edge that can trigger streaming-instability-driven planetesimal formation (Ying et al., 26 Jan 2026). The fiducial model in that study formed $\nu=\alpha c_s H$7 of planetesimals with a dust-to-planetesimal conversion efficiency of $\nu=\alpha c_s H$8, and larger disks, higher metallicities, lower viscosities, higher fragmentation thresholds, and stronger X-ray luminosities all increased the final planetesimal mass (Ying et al., 26 Jan 2026).

DustPy has also been extended to embedded Class 0/I disks with infall, heating and cooling, and snowline-dependent fragmentation thresholds (Carrera et al., 17 Apr 2025). Those calculations recover outward advection of grains and a water-snowline “advection-condensation-drift” loop, but under $\nu=\alpha c_s H$9 the midplane dust-to-gas ratio remains at least an order of magnitude below the streaming-instability threshold, even after including recent external-turbulence criteria (Carrera et al., 17 Apr 2025). This suggests that early planetesimal formation in such disks may require low-turbulence angular-momentum transport, such as winds, or a mechanism other than standard streaming instability.

7. Position within the dust-evolution ecosystem

DustPy occupies an important reference position among open-source dust-evolution codes. In recent code-comparison work, it serves as the benchmark “full coagulation–fragmentation” model against which reduced 1D and 2D approaches are calibrated (Eriksson et al., 23 Mar 2026, Pfeil et al., 2024). Compared with TriPoD, DustPy generally agrees well on dust masses and observables except when the true size distribution becomes strongly non-power-law (Eriksson et al., 23 Mar 2026). Compared with two-pop-py, DustPy usually predicts less extreme gap-edge trapping and slower dust depletion, because it does not collapse the solids into two representative populations with analytic barrier prescriptions (Eriksson et al., 23 Mar 2026).

The comparison literature also clarifies where DustPy is strongest. It is most reliable as a reference solver for one-dimensional, radially integrated dust evolution with a realistic polydisperse distribution, especially when the target quantity depends on fragmentation-fed small grains, leakage across traps, or the distinction between peak, average, and full-distribution Stokes numbers (Eriksson et al., 23 Mar 2026). At the same time, the same studies emphasize standard limitations: one-dimensional geometry, vertically integrated gas dynamics, no fully self-consistent gas–dust backreaction in many applications, and frequent reliance on prescribed pressure bumps rather than self-generated substructures (Eriksson et al., 23 Mar 2026, Kurtovic et al., 12 Jun 2025).

The broader DustPy literature therefore presents a coherent picture. DustPy is not a universal disk simulator, but a modular and comparatively high-fidelity framework for one-dimensional dust evolution that is especially valuable when the research question turns on the full evolving size distribution. Its extensive application to turbulence constraints, leakage problems, compact disks, photoevaporative clearing, and planetesimal formation shows that it has become a central numerical instrument in connecting dust microphysics to disk observables and early planet formation (Stammler et al., 2022, Jiang et al., 2023).