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

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 symmetrization 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 vertical bar Delta E/E vertical bar similar to 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-GPU(max) similar to 4 x 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 regularized integrator with post-Newtonian corrections for close encounters and binaries in the code is straightforward.

Automated summary

The study introduces a novel symplectic integrator suitable for simulating systems with a large dynamic range. By computing an additional gradient term, strict positive integrator sub-steps are achieved in the fourth-order symplectic integrator. The integration algorithm demonstrates momentum conservation through synchronous force calculations and kick operations.