From particles to precision. Simulating subsonic turbulence with smoothed particle hydrodynamics

The direct numerical simulation of subsonic turbulence with smoothed particle hydrodynamics (SPH) has traditionally been hampered by zeroth-order (E0) errors, inaccurate gradient evaluations, and excessive numerical dissipation. We aim to investigate whether a modern SPH formulation can overcome these challenges by comparing its results to those obtained with state-of-the-art non-SPH codes such as AREPO and GIZMO. We used the SPH-EXA code, a highly scalable, natively GPU-accelerated, state-of-the-art SPH code equipped with an integral-based gradient estimation, advanced artificial viscosity switches with a slope-limited reconstruction, a flexible family of pairing-resistant interpolation kernels, generalized volume elements, and improved handling of density contrasts, which maximizes Lagrangian compatibility. Our results show that SPH-EXA accurately reproduces the Kolmogorov inertial range scaling in the subsonic regime with increasing resolution, closely matching the results of state-of-the-art non-SPH methods. We also identify accurate grad-h terms as critical: a noisy standard implementation can imprint spurious granulation in the density field once dissipation is sufficiently reduced. We demonstrate that, with appropriate methodological advances, SPH can achieve a level of fidelity in modeling subsonic turbulence that rivals the most advanced Eulerian and moving-mesh approaches.