Turning the knobs on dust evolution: Comparing codes, parameters and their effects on planet formation and disc observables

Protoplanetary discs contain a wide range of dust sizes that strongly influence their thermal structure and planet formation processes such as planetesimal formation and pebble accretion. Dust evolution models are therefore essential for both planet formation simulations and the interpretation of disc observations. Several open-source dust evolution codes are available, each adopting different methods and assumptions. We present a systematic comparison of 1D radial simulations using DustPy, TriPoD, and two-pop-py, and 2D radial-vertical simulations with TriPoD, mcdust, and cuDisc. The comparison includes dust size distributions, dust disc masses, planetary gap structures, millimetre fluxes and disc sizes from synthetic observations, planetesimal formation regions, and planetary growth via pebble accretion. We also perform a parameter study to assess how key dust-evolution parameters influence disc evolution, planet formation, and code agreement. In 1D, two-pop-py depletes dust masses faster and produces higher dust concentrations outside planetary gaps than DustPy or TriPoD. The latter two generally agree well, except when size distributions deviate strongly from a power law. While the calculated millimetre fluxes and disc radii typically agree well, planetesimal formation locations and pebble accretion rates vary significantly between codes. In 2D, we compare cuDisc, mcdust, and TriPoD in simulations of turbulence- and sedimentation-driven coagulation. The dust size distributions agree well, despite the completely different numerical approaches used to model dust coagulation. The largest differences arise in the upper atmosphere, where mcdust suffers from low mass resolution and TriPoD fails to reproduce the exact shape of size distributions that deviate from a power-law.

💡 Research Summary

This paper presents a systematic benchmark of six open‑source dust‑evolution codes that are widely used in the protoplanetary‑disc community. The authors compare three one‑dimensional (radial) codes—DustPy, TriPoD, and two‑pop‑py—and three two‑dimensional (radial–vertical) codes—cuDisc, mcdust, and TriPoD—using identical initial disc structures, gas temperature profiles, turbulence parameter (α), fragmentation velocity (v_frag), and grain internal density (ρ_int). In addition, a 0.1 M_Jup planet is inserted at 30 au to generate a gap and test how each code treats dust trapping at pressure bumps.

1D Results
DustPy and TriPoD solve the Smoluchowski coagulation equation on a discrete mass grid (≥7 bins per mass decade) and couple dust dynamics to gas via Stokes‑number‑based radial drift and α‑turbulent diffusion. Both codes include perfect sticking, fragmentation, and erosion, and they treat the Epstein and Stokes‑I drag regimes identically. Consequently, their predictions for dust surface density evolution, maximum grain size (a_max), millimetre (mm) fluxes, and synthetic disc radii agree within ~5 % for most of the simulation time. However, when the size distribution deviates strongly from a simple power law—particularly near a planetary pressure bump—TriPoD’s semi‑analytic reconstruction of the size distribution (based on a fixed power‑law exponent q) leads to differences of up to 30 % in local dust mass compared with DustPy, which directly integrates the full collision kernel.

Two‑pop‑py adopts a drastically simplified two‑population approach. It assumes fixed mass ratios between a small, well‑mixed grain population and a larger, drift‑dominated population, using “fudge factors” (f_f, f_d) to approximate the growth, drift, and fragmentation limits. This yields a dramatic speed‑up but at the cost of physical fidelity. The code depletes the total dust mass fastest of the three, and it predicts a pronounced dust pile‑up just outside the planetary gap because the small grains are assumed to be perfectly coupled to the gas while the large grains rapidly reach the drift‑fragmentation barrier. Consequently, the pebble (pebble‑accretion) flux in two‑pop‑py is the highest, leading to faster core growth in planet‑formation models. DustPy and TriPoD, by contrast, supply a more modest pebble flux, delaying core formation by 20–30 % under identical initial conditions.

All three 1D codes produce mm fluxes and observable disc radii that differ by less than 10 % across the parameter space explored, suggesting that for low‑resolution observational diagnostics the choice of code may be less critical. However, the differences in pebble fluxes and dust‑mass depletion rates translate into substantial uncertainties for planet‑formation timelines.

2D Results
The two‑dimensional comparison focuses on the vertical structure of dust, turbulence‑driven coagulation, and sedimentation. cuDisc solves the Smoluchowski equation on a high‑resolution mass grid using a GPU‑accelerated finite‑volume Godunov scheme for advection–diffusion, while simultaneously evolving the gas surface density via a viscous diffusion equation and assuming vertical hydrostatic equilibrium. mcdust implements a Monte‑Carlo representative‑particle method with an adaptive grid to keep a roughly constant number of particles per cell; collisions are sampled stochastically, and turbulent kicks are added as random velocity perturbations. The 3D version of TriPoD (based on Athena++) uses a semi‑analytic prescription for a_max and a two‑fluid (small/large) representation of the dust, with vertical settling velocities derived from a prescribed Δv set.

Across the mid‑plane and moderate heights (z/H ≈ 0.5) the three codes produce remarkably consistent grain‑size distributions and growth rates, confirming that the underlying physics (Brownian motion, turbulent relative velocities, and gas drag) is robust against numerical implementation. In the upper atmosphere (z/H > 1) the differences become noticeable. mcdust’s limited number of representative particles leads to poor mass resolution, under‑representing the tail of the small‑grain distribution and reducing the atmospheric opacity by 5–10 %. TriPoD’s reliance on a power‑law reconstruction fails when the true distribution deviates from a single power law, causing an artificial suppression of the smallest fragments produced by fragmentation events. cuDisc, with its full mass‑grid solution, captures these fine details and therefore provides the most accurate atmospheric size distribution.

The resulting synthetic mm images differ by 5–12 % in total flux and by a comparable amount in inferred disc radius, indicating that even for high‑resolution ALMA observations, the choice of dust‑evolution code can introduce a non‑negligible systematic error.

Parameter Study
The authors systematically vary four key parameters: the turbulent α‑parameter (10⁻⁴–10⁻²), the fragmentation velocity v_frag (1–10 m s⁻¹), the grain internal density ρ_int (0.1–3 g cm⁻³), and the two‑pop‑py fudge factors (±20 %). The main trends are:

• Increasing α enhances vertical mixing, populating the disc atmosphere with larger grains and boosting the mm flux by ~15 %.
• Raising v_frag pushes the fragmentation barrier outward, allowing grains to grow larger, slowing dust mass loss, and increasing the pebble flux.
• Lower ρ_int reduces the Stokes number for a given size, accelerating radial drift and leading to faster dust depletion.
• Two‑pop‑py is highly sensitive to its fudge factors; a 20 % change can alter the total dust mass by up to 30 %, underscoring the semi‑analytic model’s dependence on calibration.

DustPy and TriPoD exhibit more linear and predictable responses to these parameter variations, with differences between the two codes becoming significant only when the parameters push the system into regimes where the size distribution strongly deviates from a simple power law.

Conclusions and Implications
The study demonstrates that the selection of a dust‑evolution code is not a neutral technical choice; it directly influences the inferred disc mass, the location and efficiency of pebble accretion, and the synthetic observables used to compare with ALMA data. Two‑pop‑py offers rapid parameter sweeps but can over‑estimate pebble fluxes and dust pile‑ups, potentially biasing planet‑formation timescales. DustPy and cuDisc, while computationally demanding, provide the most physically accurate treatment of the full size distribution and are preferable for studies that require precise modeling of dust opacity, vertical structure, or planetesimal formation thresholds.

The authors suggest future work should explore hybrid approaches—e.g., embedding DustPy‑style mass grids into a faster semi‑analytic framework—or Bayesian inference pipelines that marginalize over code‑specific systematic uncertainties. Such strategies would help reconcile differences and improve the reliability of both theoretical planet‑formation models and the interpretation of high‑resolution disc observations.