dropped the primes, have scaled out $m'_i$ from $p'_i$, and rewritten Eq. \[eq:HKepler\] in non-dimensional form, i.e. the gravitational constant $G=1$ for the remainder of this paper. However, we have taken care to remove any dependence on the choice of units from our implementation, so $G$ can be freely set by the user in our implementation of the algorithm.

Our task is to find the final positions and velocities $\mathbf{r}$ and $\mathbf{v}$ of a particle evolving under this Hamiltonian for some time $\mathit{dt}$, given the initial conditions $\mathbf{r}_0$ and $\mathbf{v}_0$. Thus, we seek the effect of the operator $\hat H_{\rm Kepler}(\mathit{dt})$.

It is advantageous to solve the Kepler problem numerically using the Gauss f and g functions, which express the relevant quantities in terms of $\mathbf{r}_0$ and $\mathbf{v}_0$ [@WisdomHolman1991]. This avoids the computationally expensive conversion between Cartesian and classical orbital elements, and avoids coordinate singularities associated with circular orbits. We find that it is advantageous to use universal variables in this solution [@Stumpff1962]. This approach provides greater speed and numerical stability compared to a solution using elliptic elements. It also avoids the singularity associated with the transition from elliptic to hyperbolic motion.

To solve the analogue of Kepler’s equation for 