stable and unstable periodic orbits.

There are two ways to calculate the LCN or the MEGNO. Conceptually the simplest is to integrate an additional shadow particle for each body in the simulation, i.e. a particle with slightly perturbed initial conditions. One can then directly measure the divergence of each particle’s path from its shadow. The second approach is to consider each body’s six-dimensional displacement vector $\pmb{\delta_i}$ from its shadow (in both position and velocity) as a dynamical variable. One can then obtain differential equations for each $\pmb{\delta_i}$ vector by applying a variational principle to the trajectories of the original bodies. We choose to follow the latter approach, as it is both faster and numerically more robust [@Tancredi2001]. In this scheme, one can imagine shadow particles with phase-space coordinates $\pmb \xi^s_i = \pmb\xi_i + \pmb\delta_i$, where $\pmb\xi_i=(\mathbf{r_i}, \mathbf{v_i})$ is the phase-space coordinate of the $i$-th original particle. Initially, we set each component of $\pmb \delta_i$ to a small value.

We follow the work of [@MikkolaInnanen1999] who describe how to efficiently couple the variational equations to the original equations of motion. This allows us to construct a symplectic integrator for the variational equations (a symplectic tangent map). An important advantage of this 