Computational advances and challenges in simulations of turbulence and star formation

TL;DR

The paper surveys computational advances in simulating turbulence and star formation, detailing numerical codes, parallel optimization, and the integration of MHD, gravity, radiation transfer, and cosmic-ray transport. It highlights the challenges of achieving ISM-like Reynolds numbers, the necessity of high Jeans-resolution and robust sink-particle treatments, and the crucial role of protostellar jets and radiative feedback in shaping star formation and the IMF. The review compares RT methods (FLD, M1, VET, and Monte Carlo), discusses subgrid turbulence and dissipation diagnostics, and outlines gravity solvers (FFT, multi-grid, tree) and their coupling to gas and stars. It argues that future progress will come from multi-physics, multi-scale approaches, including hybrid fluid-particle methods and GPU-accelerated exascale computing, to bridge small-scale SN/shock physics with galaxy evolution and to improve observational diagnostics for validation.

Abstract

We review recent advances in the numerical modeling of turbulent flows and star formation. An overview of the most widely used simulation codes and their core capabilities is provided. We then examine methods for achieving the highest-resolution magnetohydrodynamical turbulence simulations to date, highlighting challenges related to numerical viscosity and resistivity. State-of-the-art approaches to modeling gravity and star formation are discussed in detail, including implementations of star particles and feedback from jets, winds, heating, ionization, and supernovae. We review the latest techniques for radiation hydrodynamics, including ray tracing, Monte Carlo, and moment methods, with comparisons between the flux-limited diffusion, moment-1, and variable Eddington tensor methods. The final chapter summarizes advances in cosmic-ray transport schemes, emphasizing their growing importance for connecting small-scale star formation physics with galaxy-scale evolution.

Figures (12)

• Figure 1: Weak-scaling tests for three codes: arepo (diamonds), comparing CPU (turquoise) and GPU (blue) implementations of radiation transport (RT) in simulations including hydrodynamics and gravity, run on the 'Raven' system; athena (circles), comparing athena 4.2 (orange) with the modernized athena$^{++}$ (red) for pure MHD on a Cray XC50; and flash (stars), comparing the public version (black) with an optimized hybrid-precision version (magenta; Sect. \ref{['sec:hybrid-prec']}) for MHD turbulence on 'SuperMUC-NG'. For each code, the two modes shown are directly comparable (but not across different codes). Substantial performance gains are evident: GPU acceleration ( arepo), code modernization ( athena), and hybrid precision ( flash).
• Figure 2: Zoom into the magnetic current structures (top panels) and into the gas density (bottom left), together with the power spectrum of the magnetic-to-kinetic energy ratio $E_\mathrm{mag}/E_\mathrm{kin}$ (bottom right) in a supersonic MHD turbulence simulation with $10,\!080^3$ grid cells from BeattieEtAl2025. Supersonic turbulence produces density contrasts spanning several orders of magnitude and intricate current structures across scales, ultimately forming strongly magnetized plasma structures on small scales. Accurately modeling these requires powerful numerical schemes and extreme resolution. A key finding is that while the magnetic field has only a modest impact on large scales (low $k$), magnetic energy always approaches equipartition with kinetic energy on small scales ($k>k_\mathrm{eq}$). Numerical convergence of the energy content below $k_\mathrm{eq}$ requires resolutions $\gtrsim5,\!000^3$. Adapted from Figs. 1 and 2, and Suppl. Fig. 3 in BeattieEtAl2025.
• Figure 3: Kinetic Reynolds number ($\mathrm{Re}$, top) and magnetic Reynolds number ($\mathrm{Rm}$, bottom) as functions of the number of resolution elements ($N$) along one side of a cubic domain, shown for the subsonic (left) and supersonic (right) turbulence regimes. Data points represent simulations with explicit dissipation, i.e., $\mathrm{Re}$ and $\mathrm{Rm}$ values set by the authors of the respective works (see legend): GFN02 GotohFukayamaNakano2002, SCT+04 SchekochihinEtAl2004, HB04 HaugenBrandenburg2004dyn, HBD04 HaugenBrandenburgDobler2004, HBM04 HaugenBrandenburgMee2004, MB06 MeeBrandenburg2006, SIC+07 SchekochihinEtAl2007, FCS+11 FederrathEtAl2011, BR19 BrandenburgRempel2019, AFT+21 AchikanathEtAl2021, SF21 SetaFederrath2021, KBS+22 KrielEtAl2022, GKW+22 GalishnikovaEtAl2022, BFK+23 BeattieEtAl2023, and KBF+25 KrielEtAl2025 for the subsonic regime, and HBM04, FCS+11, FSB+14 FederrathSchoberBovinoSchleicher2014, SF21, and KBF+25 for the supersonic regime. The lines show $\mathrm{Re}$-$N$ and $\mathrm{Rm}$–$N$ relations derived by ShivakumarFederrath2025, with parameters listed in each panel. These relations provide estimates of the maximum $\mathrm{Re}$ and $\mathrm{Rm}$ achievable for a given $N$ due solely to numerical dissipation. Shaded regions indicate the associated uncertainties, accounting for methodological variations as well as differences across 14 numerical schemes, including different grid-based methods and SPH.
• Figure 4: Comparison of double-precision, single-precision, and hybrid-precision schemes for supersonic MHD turbulence. The hybrid scheme (blue line) matches the accuracy of double precision (red dotted line), conserving gas mass (left) and momentum (right), while single precision (green line) exhibits significant errors. At the same time, the hybrid scheme provides a computing speed-up and reduces memory use by a factor of $\sim2$ relative to double precision.
• Figure 5: Jeans resolution study of accretion disk and outflow formation, $1000\,\mathrm{yr}$ after the formation of the protostar, with edge-on slices through the disk. Panels show identical simulations but with Jeans resolutions of $\mathrm{J_{res}}=2$, 4, 8, 16, 32, and 64 cells per $\lambda_\mathrm{J}$. For $\mathrm{J_{res}}=2$, artificial fragmentation produces four sink particles (annotated as $N_\mathrm{sink}$), while only a single star forms for $\mathrm{J_{res}}\geq4$, confirming the TrueloveEtAl1997 criterion. However, disk and outflow structure, as well as the accretion rate (measured by the fraction of gas accreted, i.e., the star formation efficiency after $1000\,\mathrm{yr}$), only converge for $\mathrm{J_{res}}\gtrsim30$. Figure adapted from FederrathEtAl2014.