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

TL;DR

This work addresses momentum non-conservation and parallelization challenges in self-gravity calculations for SPH with adaptive softening. It implements a parallel, momentum-conserving symmetric fast multipole method by duplicating leaf-node interactions in a dual-tree traversal, preserving Newton's third law while maintaining comparable accuracy to existing methods. The approach yields machine-precision linear momentum conservation, improved angular momentum behavior, and robust orbitalPhase fidelity in binary polytrope and common-envelope tests, with performance scales similar to the previous solver. The method is now the default in Phantom, enabling accurate, scalable gravity calculations for large SPH simulations with adaptive softening.

Abstract

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.

Figures (9)

• Figure 1: Multipole methods for computing self-gravity on particles. Left columns shows computation of a single force applied to a target particle, shown in blue. Right column shows forces applied to all particles. The first row shows the direct sum method where every pairwise interaction is computed. The multipole method (second row) uses multipole expansions directly applied on particles of the system to reduce the number of interaction. The fast multipole method (FMM; third row) applies multipole expansion on nodes that hold particles instead to further reduce the number of interactions. The last row corresponds to the symmetric fast multipole method (SFMM) where only pairwise node-node interactions satisfy Newton's third law.
• Figure 2: Mean relative error (compared to a direct sum) versus the maximum opening angle computed in a Plummer sphere (top) and a homogeneous sphere (bottom) of $10^5$ equal-mass particles, comparing FMM (red) to SFMM (black). The range boundaries around each curve represent the maximum error (upper) and the tenth percentile (lower). The close overlap between both methods suggests a similar overall accuracy.
• Figure 3: Wall time (in seconds) spent computing gravitational forces in a Plummer sphere of $10^5$ particles for a range of maximum opening angles. The two multipole methods are shown in red (FMM) and black (SFMM). The reference time of the direct algorithm is plotted in green and is unrelated to the critical opening angle.
• Figure 4: Wall time (in seconds) spent computing gravitational forces in a Plummer sphere with $10^4$ to $10^6$ particles with $\theta_\mathrm{crit}=0.5$. The two multipole methods, FMM in red and SFMM in black, are compared to the direct sum algorithm in green. Both multipole method follows their expected $\mathcal{O}(N\log N)$ scaling while the direct sum is several order of magnitude more computationally expensive following a $\mathcal{O}(N^2)$.
• Figure 5: Twenty orbits of a binary system of polytropic stars. The render of the logarithmic column density of each timestep is stacked to construct the path of the binary components. With the non-symmetric tree (left), a small but accumulative non-conservation error causes the orbit to drift with time, while this drift is removed with the new implementation (right).