our standard symplectic map. After the simulation has finished (or whenever we need an output) we convert back from mapping to real coordinates. The symplectic corrector operator that we use is a combination of several $\hat \operatorname{\mathcal{H}}_{\rm Interaction}(\mathit{dt})$ and $\hat \operatorname{\mathcal{H}}_{\rm Kepler}(\mathit{dt}) \circ \hat\operatorname{\mathcal{H}}_0(\mathit{dt})$ operators applied for different (positive and negative) intervals $\mathit{dt}$. If $\epsilon$ is the order of the perturbations, i.e. the mass ratio and therefore the relative magnitude of $\hat \operatorname{\mathcal{H}}_{\rm Interaction}$ compared to $\hat \operatorname{\mathcal{H}}_{\rm Kepler}$, then one can show that the use of symplectic correctors can lead to a scheme of order $O(\epsilon \mathit{dt}^{K} ) + O(\epsilon^2 \mathit{dt}^2)$ where $K$ is the order of the symplectic corrector [@MikkolaPalmer2000]. A second order Wisdom-Holman map without symplectic correctors has an energy error of order $O(\epsilon^2 \mathit{dt}^2)$. Because this coordinate transformation for the symplectic corrector is only performed for outputs and at the beginning and end of the simulation, its effect on the speed of the algorithm is negligible for sparse output.

A full derivation of the symplectic correctors would go beyond the scope of this paper and we refer the reader to [@Wisdom1996] and [@MikkolaPalmer2000]. The corrector coefficients are listed in a compact form in [@Wisdom2006].

We implement 