Momentum-conserving self-gravity in the phantom smoothed particle hydrodynamics code. Parallel dual tree traversal for the symmetric fast multipole method

Tree codes that approximate groups of distant particles with multipole expansions are the standard way to accelerate the computation of self-gravity on particles. While momentum-conserving fast multipole methods exist, parallelisation is non-trivial and previous implementations have been limited to self-gravity with fixed softening lengths. We aim for a practical, parallel version of Dehnen’s momentum-conserving Cartesian fast multipole method for the computation of the gravitational force in smoothed particle hydrodynamics (SPH) with adaptive gravitational force softening. We parallelise the dual tree walk by replicating the node-node interaction on the parents of each leaf node in the tree. While this duplicates work, it greatly simplifies the parallelisation and can be implemented with relatively minor changes from the previous non-conservative force algorithm in Phantom. We also adapt the tree opening criterion for adaptive softening lengths, such that all interactions within the softening kernel are handled pairwise (as in SPH) rather than with multipole expansions, also allowing the gravity calculation to be performed alongside the SPH force evaluation. We demonstrate that the new code conserves linear momentum to machine precision while giving similar force accuracy and computational performance to the previous (non-symmetric) fast multipole method in Phantom . The new method also gives better conservation of the angular momentum and orbital phase in a binary polytrope evolution. The symmetric fast multipole method is now the default for computing self-gravity in the public code.

💡 Research Summary

This paper presents a practical, parallel implementation of Dehnen’s momentum‑conserving Cartesian Fast Multipole Method (FMM) for self‑gravity in the Smoothed Particle Hydrodynamics (SPH) code Phantom, extending it to support adaptive gravitational softening. Traditional tree‑based gravity solvers (e.g., Barnes‑Hut) reduce the computational cost from O(N²) to O(N log N) but do not guarantee that the force between any two particles is exactly equal and opposite, leading to only approximate linear momentum conservation. The original symmetric FMM (SFMM) introduced by Dehnen ensures exact Newton’s third‑law compliance, but its recursive dual‑tree walk is intrinsically serial and has only been used with fixed Plummer‑type softening.

The authors address two major challenges: (1) parallelising the dual‑tree traversal, and (2) incorporating adaptive softening lengths that vary with the local SPH smoothing length. Their parallelisation strategy deliberately duplicates work: each leaf node is processed independently by a separate thread, which climbs the tree to the root, pushes the parent pair onto a stack, and then performs a local dual‑tree walk. Although this leads to multiple threads evaluating the same node‑node interaction, it eliminates the need for complex task‑scheduling or lock‑based synchronization, making the algorithm embarrassingly parallel on shared‑memory CPUs. Consequently, the theoretical O(N) scaling of a perfectly shared dual‑tree walk degrades to O(N log N), matching the cost of conventional tree codes while preserving exact linear momentum to machine precision.

To handle adaptive softening, the authors augment the standard multipole acceptance criterion (MAC) θ = (l_tgt + l_src)/d < θ_crit with a secondary distance check that incorporates the maximum smoothing lengths within the interacting nodes: d <