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

Abstract: 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.

• The paper compares six dust-evolution approaches and finds DustPy and TriPoD generally agree, while two-pop-py depletes dust too quickly and produces unrealistic gap-edge concentrations.
• The paper shows integrated millimetre fluxes and dust masses typically agree within about 50%, but similar fluxes can conceal major structural errors in synthetic disc images.
• The paper finds planet-formation outcomes are far more sensitive to code choice, representative Stokes number, and dust-size treatment, with pebble-accretion masses differing by factors of several.

Overview and motivation

Dust evolution models underpin both planet formation simulations and the interpretation of protoplanetary disc observations, yet the community relies on a heterogeneous set of open-source codes built on very different numerical approaches. Eriksson et al. (2603.22550) present a systematic comparison of three 1D radial dust evolution codes — DustPy (a full Smoluchowski solver), TriPoD (a semi-analytic two-fluid model), and two-pop-py (a fudge-factor prescription) — together with a smaller comparison of three 2D radial–vertical codes: cuDisc (GPU-based Smoluchowski solver), mcdust (Monte Carlo representative-particle method), and TriPoD in an Athena++-based implementation. Beyond code agreement, the study performs a one-parameter-at-a-time parameter study spanning stellar mass ($0.05$–$1.3\,M_\odot$), accretion rate, cut-off radius, dust-to-gas ratio $Z = 10^{-4}$–$10^{-1}$, turbulence $\alpha_{\rm t} = 10^{-5}$–$10^{-2}$, fragmentation velocity $v_{\rm frag} = 0.1$–$25\,\mathrm{m\,s^{-1}}$, internal grain density, initial grain size, mass-bin resolution, and treatment of initially drifting particles. The diagnostics compared include size distributions, dust disc masses, planetary gap structures, millimetre fluxes and disc radii from synthetic observations, streaming-instability (SI) planetesimal formation regions, and pebble-accretion growth tracks.

The central finding is a sharp dichotomy: bulk disc properties such as dust masses and millimetre fluxes are relatively robust to the choice of code, whereas planet formation outcomes — where and when planetesimals form, and how fast planets grow by pebble accretion — depend strongly on it. This has direct consequences for population synthesis studies that treat dust evolution as interchangeable between prescriptions.

Methodological setup

All 1D simulations use a viscous $\alpha$-disc with $\alpha = 5\times10^{-3}$, a passively irradiated temperature profile, a grid extending to 1000 au, and a nominal 3 Myr evolution with no disc dispersal mechanism. The nominal case uses $1.3\,M_\odot$0, $1.3\,M_\odot$1, $1.3\,M_\odot$2, and $1.3\,M_\odot$3. For TriPoD, the authors use TriPoDPy for radial runs and a modified Athena++ version for 2D, including a recalibration of the shrinkage parameter from $1.3\,M_\odot$4 to $1.3\,M_\odot$5 and a change in how the maximum particle size diffuses; this recalibration was tuned against cuDisc, which should be kept in mind when interpreting the good 2D agreement. Planet-induced gaps are imposed via Duffell gap profiles applied inversely to the viscosity rather than through explicit planetary potentials.

For planetesimal formation, the paper applies three SI criteria: YG05 ($1.3\,M_\odot$6), LI24 (forced turbulence), and LI25 (turbulence-free clumping threshold). A key methodological caveat is stated plainly: these criteria are derived for monodisperse particles at a fixed global pressure gradient of 0.05, while the simulated discs have evolving, spatially varying pressure gradients — particularly near planetary gaps. The choice of representative Stokes number (density-weighted average versus distribution peak) is shown to matter as much as the choice of criterion.

1D results: sizes, masses, and gaps

In unperturbed discs, DustPy and TriPoD agree well on total dust mass, while two-pop-py systematically depletes dust too rapidly. TriPoD deviates from DustPy precisely when the size distribution departs from a power law — at low gas densities or high fragmentation velocities, where drift-limited growth produces bimodal distributions. Parameter trends follow analytic expectations: particle sizes scale as $1.3\,M_\odot$7, higher turbulence slows depletion (smaller grains drift more slowly), and low-$1.3\,M_\odot$8 discs lose only about 60% of their dust over 3 Myr because growth never reaches either barrier.

The most consequential result concerns planetary gaps. two-pop-py produces nonphysically strong dust concentrations at gap edges, far exceeding DustPy and TriPoD, consistent with independent work by Houge et al. (in prep.) showing suppressed cross-gap transport in two-pop-py for $1.3\,M_\odot$9. The authors conclude explicitly that two-pop-py is unsuitable for discs with planetary gaps except where little pile-up is expected. Notably, the presence of a gap-opening planet does not change total dust mass evolution in the nominal setup, implying efficient cross-gap filtering for those parameters.

Observational diagnostics

Millimetre fluxes (at 1 mm, using DSHARP opacities) broadly track dust masses, and all three codes agree within roughly 50% across most of the parameter space, with deviations concentrated at low gas densities. However, RADMC-3D intensity maps of the $Z = 10^{-4}$0 case reveal that the similar total fluxes mask radically different morphologies: two-pop-py's over-bright ring at the pressure maximum compensates for a fainter outer disc. Integrated fluxes alone can therefore hide large structural errors — a caution for studies fitting only integrated ALMA photometry. Disc radii ($Z = 10^{-4}$1) agree well for TriPoD but show large time-variable errors for two-pop-py in some cases, traced to its size-independent growth timescale assumption in the outer disc.

Planetesimal formation via the streaming instability

Where SI criteria are satisfied varies significantly between codes, especially when criteria are marginally met. The YG05 criterion ($Z = 10^{-4}$2) is essentially never reached in smooth discs, and LI24 similarly requires strong trapping; LI25 is met more readily. Using the density-weighted average Stokes number requires significant trapping or supersolar metallicity regardless of other parameters, whereas the peak Stokes number yields substantially larger formation regions — meaning the representative-$Z = 10^{-4}$3 choice dominates over physical parameter variations. Among parameter trends, larger cut-off radii widen formation regions (outer-disc replenishment sustains inner-disc dust-to-gas ratios), higher turbulence suppresses formation, and increasing $Z = 10^{-4}$4 slightly reduces it. All three codes predict planetesimal formation at gap edges provided turbulence stays below $Z = 10^{-4}$5, though the authors note that accounting for the pressure-gradient dependence of the clumping threshold would make gap-edge formation easier than reported here.

Pebble accretion

Pebble accretion proves the most sensitive diagnostic. With polydisperse accretion, TriPoD final planet masses differ from DustPy by at most a factor of ~2 — often larger than the discrepancy in total dust mass itself, demonstrating that small differences in dust evolution amplify into larger differences in growth outcomes. Monodisperse assumptions bracket the polydisperse result in opposite directions: density-weighted-average $Z = 10^{-4}$6 yields lower-mass planets, peak $Z = 10^{-4}$7 higher-mass ones. two-pop-py, tracking only the maximum particle size, produces planets several times more massive than either full-distribution code. Parameter trends show non-monotonic behavior with $Z = 10^{-4}$8: growth accelerates up to $Z = 10^{-4}$9 but reverses beyond that as faster-drifting large grains deplete the disc. An embedded $10^{-1}$0 planet at 10 au reduces the final mass of a 5 au embryo by only ~30%, indicating modest gap filtering even for giant companions.

2D radial–vertical comparison

The 2D comparison targets sedimentation-driven coagulation ($10^{-1}$1) versus turbulence-dominated coagulation ($10^{-1}$2), run for $10^{-1}$3 yr on identical column-density structures. Despite completely different numerical methods — GPU Smoluchowski solving, Monte Carlo superparticles, and semi-analytic fluids — the three codes agree well on vertically integrated size distributions, mass-averaged grain sizes, and vertical density profiles within ~2 pressure scale heights. All reproduce the characteristic bimodal, sedimentation-driven pattern in the weak-turbulence outer disc, where grains grow while settling and rain out into the midplane. Differences concentrate in the upper atmosphere: mcdust suffers from finite particle-count resolution there (partially alleviable with more superparticles), and TriPoD cannot capture the tapered, non-power-law distributions that emerge at altitude. DustPy's assumption of vertically constant Stokes numbers fails in the dynamical growth phase, underestimating grain sizes when settling-driven coagulation operates. These results validate approximate fluid models like TriPoD for hydrodynamic coupling, while identifying cuDisc as the reference tool for hydrostatic structure calculations and Monte Carlo methods for chemistry requiring particle histories.

Limitations and open questions

The paper is candid about scope. All simulations neglect bouncing (which Dominik & Dullemond showed shrinks maximum sizes and depletes micron grains), porosity evolution and compaction, ice composition and icelines, and any disc dispersal mechanism — photoevaporation could trigger a late generation of planetesimals as gas densities drop. Gas evolution is purely viscous despite magnetized disc winds being a viable alternative, and the effect of wind-driven evolution on dust depletion is left unquantified. The SI criteria are applied outside their derivation regime (monodisperse, fixed pressure gradient), and no validated polydisperse clumping threshold exists. Numerically, high fragmentation velocities ($10^{-1}$4) produce boundary-driven artifacts or instability in all three 1D codes, deep gaps ($10^{-1}$5) become computationally intractable, and extreme minimum grain sizes cause runtime or accuracy problems. Whether the recalibrated TriPoD parameters generalize beyond the setups tested here remains open, as does the question of how the mono-disperse-derived SI thresholds translate to genuinely growing, polydisperse populations.

Conclusion

This work establishes that the choice of dust evolution prescription is largely inconsequential for integrated disc observables but decisive for planet formation predictions. DustPy and TriPoD form a mutually consistent pair across most of parameter space; two-pop-py is reliable only for smooth, gap-free discs and fails qualitatively in the presence of pressure bumps. The demonstration that pebble accretion amplifies sub-factor-of-two differences in dust evolution into factor-of-several differences in planetary mass, and that the representative-Stokes-number convention can dominate SI-based planetesimal formation maps, provides concrete guidance for interpreting existing population synthesis results and for designing future coupled dust–planet formation models.