Smoothed Particle Hydrodynamics in pkdgrav3 for Shock Physics Simulations

Introduction:

Collisions are a fundamental driver of planetary system formation and evolution, shaping bodies from planetesimals to fully formed planets. Smoothed Particle Hydrodynamics (SPH) is a commonly used tool for modeling these highly non-linear, energetic events. At the same time, advances in high-performance computing (HPC) enable much higher resolution simulations, provided codes can efficiently exploit modern GPU-accelerated architectures, which have become ubiquitous especially due to the AI boom.

Most existing approaches specialize in either large-scale, gravity-dominated impacts or small-scale, strength-dominated collisions. Giant impact simulations typically neglect material strength, while small-body studies include increasingly sophisticated solid mechanics. Indeed, the transition between these regimes is gradual and depends on size, impact conditions, and material properties. Bridging this gap requires a single scalable framework that can model both hydrodynamic and solid behavior consistently at high resolution.

Here, we present such a framework within the massively parallel N-body code pkdgrav3 [1,2]. Rather than focusing on individual model components, we emphasize the range of problems accessible with this approach, drawing on both published [3–5] and ongoing work.

Methods:

Our implementation uses a density–energy formulation of SPH [6], ensuring conservation of mass, momentum, and energy, and compatibility with standard equations of state for planetary materials. The hydrodynamics module is designed for strongly shocked regimes while remaining stable at free surfaces and material interfaces. We include an interface correction [7] to suppress spurious surface tension, and we recover entropy conservation in adiabatic flows without explicitly evolving entropy [8].

The method is implemented within pkdgrav3 and directly leverages its domain decomposition, tree-based gravity solver, and GPU-accelerated particle interactions. The novel approach for neighbor finding uses the same tree as the FMM gravity, improving efficiency and scalability to very large particle numbers.

This hydrodynamics configuration already provides a robust tool for large-scale impacts between differentiated bodies and gas-rich systems where strength effects are minor. We extend this with a pressure-dependent strength model to represent solid materials. Instead of assuming purely fluid behavior, this allows material to sustain shear stresses and provides a continuous transition from strength- to gravity-dominated regimes. The strength model follows the same parallelization and GPU execution strategy, preserving scalability.

Results and Applications:

The hydrodynamics implementation was validated against standard test problems (e.g., Sod shock tube, Sedov-Taylor blast wave, Evrard collapse, Gresho-Chan vortex), showing excellent agreement with reference solutions and previous studies. We do not present these tests here and instead focus on the new applications the code enables.

From a computational perspective, the code shows good strong and weak scaling on modern GPU systems, maintaining high parallel efficiency from moderate problem sizes suitable for parameter studies up to simulations with billions of particles. This level of performance enables systematic studies of planetary collisions at intermediate resolutions, instead of limiting such simulations to individual cases.

At intermediate resolution, the code supports large parameter surveys that were previously limited to either low resolution or a small number of runs. We present two such studies [3,5] exploring the dependence of collision outcomes on impact angle, velocity, and material properties across broad parameter ranges. These simulations already reach what was considered high resolution in earlier work while allowing a much broader exploration of initial conditions.

At the same time, the framework enables ultra-high-resolution simulations for selected cases. We show three examples [4], (Meier et al., in prep), and (Reinhardt et al., in prep) reaching particle numbers sufficient to resolve ultra-low-mass components such as crusts, atmospheres, surface oceans and debris disks. These components, though small in mass, are often critical for the outcome of the collision and are poorly captured at lower resolution.

A key advantage of the framework is that it combines very high computational performance with material strength models, allowing different physical assumptions to be compared within the same numerical setup. Previous strength-capable codes were typically limited to only a few million particles, whereas the present framework enables substantially higher-resolution simulations without sacrificing physical realism. In particular, we compare simulations with and without material strength across several regimes. While large-scale outcomes in gravity-dominated impacts remain broadly similar, including strength systematically shifts the catastrophic disruption threshold. These effects persist to larger scales than often assumed and influence both small-body and planetary-scale collision outcomes.

Conclusions:

We present an SPH-based framework within pkdgrav3 that combines efficient GPU-accelerated hydrodynamics with pressure-dependent material strength. This enables consistent modeling of collisions across regimes typically treated separately, from strength-dominated small-body impacts to gravity-dominated planetary events.

The examples highlight the flexibility of the approach: the same code supports both large parameter studies at intermediate resolution and targeted ultra-high-resolution simulations. This makes it possible to move beyond single “representative” impacts toward systematic exploration of parameter space while retaining the ability to zoom in on specific cases.

Our results also show that material strength affects collision outcomes over a wider range of conditions than commonly assumed, particularly through its impact on the catastrophic disruption threshold. This underscores the importance of treating strength and gravity within a unified framework rather than switching between fundamentally different models.

Overall, the framework combines physical consistency with computational scalability, providing a practical tool for studying planetary collisions across a wide range of scales and conditions. Future work will extend the physical models and apply the code to a broader set of problems in planetary formation and evolution.

References:

[1] Meier, T. Potter, D., Reinhardt, C., Stadel, J., ApJ 1000, 266 (2026).

[2] Meier, T., Reinhardt, C., Jutzi, M., Potter, D., Stadel, J., https://doi.org/10.48550/arXiv.2603.19764 (2026).

[3] Matzkevich, Y., Reinhardt, C., Meier, T., Stadel, J., Helled, R., A&A 961, A184 (2024)

[4] Meier, T., Reinhardt, C., Shibata S., Müller, S., Stadel, J., Helled, R., ApJ 988, 7 (2025)

[5] Bussmann, M., Reinhardt, C., Gillmann, C., Meier, T., Stadel, J., Tackley, P., Helled, R., A&A 702, A106 (2025)

[6] Springel, V., & Hernquist, L., MNRAS 333 (2002)

[7] Ruiz-Bonilla, S., Borrow, J., Eke, V., Kegerreis, J., Massey, R., Sandnes, T., Teodoro, L., MNRAS 512.3 (2022)

[8] Reinhardt, C., & Stadel, J., MNRAS 467 (2017)