FROST: a momentum-conserving CUDA implementation of a hierarchical fourth-order forward symplectic integrator

TL;DR

FROST introduces a GPU-accelerated hierarchical fourth-order symplectic integrator for N-body simulations, conserving momentum and energy with high accuracy, suitable for large-scale astrophysical systems with adaptive time-stepping.

Contribution

It presents a novel hierarchical formulation and GPU implementation of a fourth-order forward symplectic integrator, enabling efficient large-scale N-body simulations with high precision.

Findings

Conserves energy to |ΔE/E| ~ 10^{-10} in simulations.

Scales efficiently up to ~4×N/10^5 GPUs for million-body simulations.

Maintains negligible linear and angular momentum errors.

Abstract

We present a novel hierarchical formulation of the fourth-order forward symplectic integrator and its numerical implementation in the GPU-accelerated direct-summation N-body code FROST. The new integrator is especially suitable for simulations with a large dynamical range due to its hierarchical nature. The strictly positive integrator sub-steps in a fourth-order symplectic integrator are made possible by computing an additional gradient term in addition to the Newtonian accelerations. All force calculations and kick operations are synchronous so the integration algorithm is manifestly momentum-conserving. We also employ a time-step symmetrisation procedure to approximately restore the time-reversibility with adaptive individual time-steps. We demonstrate in a series of binary, few-body and million-body simulations that FROST conserves energy to a level of $|\Delta E / E| \sim 10^{-10}$ while errors in linear and angular momentum are practically negligible. For typical star cluster simulations, we find that FROST scales well up to $N_\mathrm{GPU}^\mathrm{max}\sim 4\times N/10^5$ GPUs, making direct summation N-body simulations beyond $N=10^6$ particles possible on systems with several hundred and more GPUs. Due to the nature of hierarchical integration the inclusion of a Kepler solver or a regularised integrator with post-Newtonian corrections for close encounters and binaries in the code is straightforward.