A Comparative Study of the Streaming Instability: Unstratified Models with Marginally Coupled Grains

Abstract

The streaming instability is a leading mechanism for concentrating solids and initiating planetesimal formation in protoplanetary disks. Although numerous studies have explored its linear growth, nonlinear evolution, and implications for planet formation, the diversity of numerical methods and dust treatments used across the literature has made it difficult to assess which features of the instability are physically robust and which arise from code-dependent choices. We present the first systematic comparison of seven hydrodynamic codes--spanning finite-volume and finite-difference schemes and modeling dust either as Lagrangian particles or as a pressureless fluid--applied to the unstratified streaming instability with a dimensionless stopping time of unity. All codes reproduce the characteristic sequence of exponential growth, filament formation, and turbulent saturation, demonstrating broad agreement in the qualitative behavior of the instability. Quantitatively, however, the dust model remains the dominant source of variation at moderate resolution: particle-based simulations reach higher peak densities and exhibit broader high-density tails than fluid-based models at $512^2$ resolution, although increasing the number of particles brings their initial maximum density evolution into close agreement with that of dust-fluid models. At $1024^2$, these differences diminish substantially, indicating better agreement of the saturated-state statistics across dust treatments. In terms of computational performance, most particle implementations suffer from imbalanced parallelized loads, while execution on a GPU is at least two to three times more energy efficient and scales better at higher resolutions than on CPUs. Given the intrinsic stochasticity of this nonlinear system, only statistical diagnostics remain meaningful across codes.

Figures (14)

• Figure 1: Dust density snapshots for Problem BA at grid resolution $512^2$. Each row shows those from a single code (Section \ref{['sec:codes_numerical_methods']}), and each column shows the simulation time in local orbital periods $T$. In alphabetical order from top to bottom, upper and lower groups implement an average of $n_\mathrm{p} = 1$ particle per grid cell/point and a pressureless dust fluid, (Sections \ref{['sec:lagrangian_dust_particles']} and \ref{['sec:pressureless_dust_fluid']}) respectively. Radial $x$ and vertical $z$ coordinates are in units of the vertical gas scale height $H$. The color scale at the bottom right shows the dust density $\rho_\mathrm{d}$ relative to the initially uniform gas density $\rho_\mathrm{g,0}$ and applies to all snapshots.
• Figure 2: Maximum dust density $\max(\rho_\mathrm{d})$ as a function of time $t < 10T$ for Problem BA at grid resolution $512^2$. The units for density and time are the initially uniform gas density $\rho_\mathrm{g,0}$ and the local orbital period $T$, respectively. Using a similar cadence for each time series, dotted lines show particle codes (with average of $n_\mathrm{p} = 1$ particle per grid cell/point), solid lines show codes with pressureless dust fluids, and colors differentiate the codes (Sections \ref{['sec:lagrangian_dust_particles']}, \ref{['sec:pressureless_dust_fluid']}, and \ref{['sec:codes_numerical_methods']}, respectively).
• Figure 3: Similar to Figure \ref{['fig:BA-512_time_series_pre-10T']} except until $t_\mathrm{lim} = 100T$. The lower panel shows the simple moving average (SMA) with a sampling window of $10T$.
• Figure 4: Cumulative distribution functions of the dust density $\rho_\mathrm{d}$ for Problem BA at grid resolution $512^2$. Dotted and solid lines show codes that implement an average of $n_\mathrm{p} = 1$ particle per grid cell/point and a pressureless dust fluid, (Sections \ref{['sec:lagrangian_dust_particles']} and \ref{['sec:pressureless_dust_fluid']}) respectively, with all $\rho_\mathrm{d}$ relative to the initially uniform gas density $\rho_\mathrm{g,0}$. In the upper panel, curves show time averages over the saturated state $\overline{\rho_\mathrm{d}}$ (cf. Figure \ref{['fig:BA-512_time_series']}), colors different codes (Section \ref{['sec:codes_numerical_methods']}), and shaded areas the $1\sigma$ time variability. In the lower panel, curves show the mean $\overline{\overline{\rho_\mathrm{d}}}$ of all particle or fluid codes in the upper panel, and shaded areas the standard deviation of their respective $\overline{\rho_\mathrm{d}}$.
• Figure 5: Similar to Figure \ref{['fig:BA-512_snapshots']} except for particle codes with $n_\mathrm{p} = 9$ only.