High order symplectic integrators for planetary dynamics and their implementation in REBOUND

Direct N-body simulations and symplectic integrators are effective tools to study the long-term evolution of planetary systems. The Wisdom-Holman (WH) integrator in particular has been used extensively in planetary dynamics as it allows for large timesteps at good accuracy. One can extend the WH method to achieve even higher accuracy using several different approaches. In this paper we survey integrators developed by Wisdom et al. (1996), Laskar & Robutel (2001, and Blanes et al. (2013). Since some of these methods are harder to implement and not as readily available to astronomers compared to the standard WH method, they are not used as often. This is somewhat unfortunate given that in typical simulations it is possible to improve the accuracy by up to six orders of magnitude (!) compared to the standard WH method without the need for any additional force evaluations. To change this, we implement a variety of high order symplectic methods in the freely available N-body integrator REBOUND. In this paper we catalogue these methods, discuss their differences, describe their error scalings, and benchmark their speed using our implementations.