A Well-Posed Kelvin-Helmholtz Instability Test and Comparison

Recently, there has been a significant level of discussion of the correct treatment of Kelvin-Helmholtz instability in the astrophysical community. This discussion relies largely on how the KHI test is posed and analyzed. We pose a stringent test of the initial growth of the instability. The goal is to provide a rigorous methodology for verifying a code on two dimensional Kelvin-Helmholtz instability. We ran the problem in the Pencil Code, Athena, Enzo, NDSPHMHD, and Phurbas. A strict comparison, judgment, or ranking, between codes is beyond the scope of this work, though this work provides the mathematical framework needed for such a study. Nonetheless, how the test is posed circumvents the issues raised by tests starting from a sharp contact discontinuity yet it still shows the poor performance of Smoothed Particle Hydrodynamics. We then comment on the connection between this behavior to the underlying lack of zeroth-order consistency in Smoothed Particle Hydrodynamics interpolation. We comment on the tendency of some methods, particularly those with very low numerical diffusion, to produce secondary Kelvin-Helmholtz billows on similar tests. Though the lack of a fixed, physical diffusive scale in the Euler equations lies at the root of the issue, we suggest that in some methods an extra diffusion operator should be used to damp the growth of instabilities arising from grid noise. This statement applies particularly to moving-mesh tessellation codes, but also to fixed-grid Godunov schemes.

💡 Research Summary

The paper presents a rigorously defined two‑dimensional Kelvin‑Helmholtz instability (KHI) test that focuses on the initial linear growth phase, and uses it to benchmark five widely used astrophysical fluid codes: Pencil Code, Athena, Enzo, NDSPHMHD, and Phurbas. The authors argue that many previous KHI tests suffer from an ill‑posed initial condition—a sharp contact discontinuity—that introduces infinite gradients and prevents convergence of numerical derivatives. To avoid this, they construct smooth initial profiles for density and x‑velocity using exponential transitions over a small length scale L = 0.025, with densities ρ₁ = 1.0, ρ₂ = 2.0 and velocities U₁ = 0.5, U₂ = ‑0.5. A small sinusoidal perturbation Vy = 0.01 sin(4πx) seeds the instability. The pressure is uniform (p = 2.5) and an ideal‑gas equation of state with γ = 5/3 is employed. The computational domain is a unit square with periodic boundaries; simulations are run at three resolutions (128², 256², 512²).

All codes solve the inviscid compressible Euler equations, allowing a direct comparison of discretization strategies (finite‑difference, finite‑volume, unstructured mesh, and mesh‑free Lagrangian). The authors introduce two diagnostic quantities that can be evaluated consistently across all methods. The first is the global amplitude M of the y‑velocity mode, defined through a weighted convolution with sin(4πx) and an exponential envelope that isolates the shear layer. For grid‑based codes the sums are taken over cell centers; for SPH the particle contributions are weighted by the smoothing length raised to the dimensionality (h_i^p); for unstructured meshes a generic quadrature weight w_i is used. The second diagnostic is the maximum y‑direction kinetic energy density, ½ ρ Vy², taken over all computational elements. The amplitude M is relatively insensitive to noise, while the kinetic‑energy maximum is highly sensitive and therefore useful for exposing numerical artifacts.

Results show that the fixed‑grid Godunov codes (Athena, Enzo) and the high‑order Pencil Code reproduce the analytic linear growth rate (≈ exp(4.384 t) for the incompressible limit) with excellent convergence as resolution increases; at 512² the deviation is below a few percent. In contrast, traditional Smoothed Particle Hydrodynamics (SPH) exhibits a markedly slower growth due to its inherent zeroth‑order inconsistency, which manifests as an artificial surface‑tension‑like force at the shear interface. Even when higher‑order kernels (e.g., quintic) or Godunov‑SPH formulations are employed, the convergence remains poorer than that of grid‑based methods.

The study also highlights that codes with very low numerical diffusion—particularly moving‑mesh Voronoi schemes and low‑diffusion fixed‑mesh Godunov solvers—can generate secondary Kelvin‑Helmholtz billows from grid‑scale noise. Because the Euler equations lack a physical diffusion scale, such secondary structures are numerical artifacts rather than physical phenomena. The authors therefore recommend adding a modest, controlled diffusion operator (e.g., a low‑order hyper‑diffusion or explicit filter) to damp spurious secondary modes, especially in moving‑mesh implementations where mesh reconstruction can amplify asymmetries.

Overall, the paper makes four key contributions: (1) a well‑posed KHI test with smooth initial conditions that enables convergence studies; (2) a unified diagnostic framework (global mode amplitude and maximum kinetic energy) applicable to grid, unstructured, and mesh‑free codes; (3) a quantitative comparison of five representative astrophysical fluid solvers, demonstrating the superior performance of high‑order Godunov schemes and the limitations of traditional SPH; and (4) practical guidance on mitigating secondary numerical instabilities via controlled diffusion. These results provide a robust benchmark for future code development and validation in astrophysical fluid dynamics, where accurate modeling of shear‑driven instabilities is essential for studying mixing, turbulence, and a wide range of astrophysical phenomena.