time to the power of one half.

We have also sped up the integrator through various improvements to the integrator’s Kepler-solver. Our implementation allows for the evolution of variational equations (to determine whether orbits are chaotic) at almost no additional cost. Additionally, we implement so-called symplectic correctors up to order eleven to increase the accurary [@Wisdom1996], allow for arbitrary unit choices, and do not tie the integration to a particular frame of reference.

We make our integrator, which we call [[WHFast]{}]{}, publicly available in its native C99 implementation and as an easy-to-use python module.

The remainder of this paper is structured as follows. We first summarize the concepts and algorithms used in this paper, including Jacobi Coordinates, our choice of Hamiltonian splitting, the symplectic Wisdom-Holman map, symplectic correctors and the variational equations in Sect. \[sec:background\]. We then go into detail discussing the improvements we have made to these algorithms in Sect. \[sec:improvements\]. Numerical tests are presented in Sect. \[sec:numericalresults\] before we conclude in Sect. \[sec:conclusions\].

Background {#sec:background}
==========

The Hamiltonian $\mathcal{H}$ of the gravitational $N$-body system can be written as the sum of kinetic and potential terms in Cartesian coordinates $$\begin{aligned}
\mathcal{H} &=& \sum_{i=0}^{N-1} \frac{\mathbf{p}_i^2}{2m_i} - \sum_{i=0}^{N-1} \sum_{j=i+1}^{N-1} \frac{Gm_im_j}{|\mathbf{r}_i-\mathbf{r}_j|}.\label{eq:H}\end{aligned}$$ One 