method is that we only solve Kepler’s equation once for each particle/shadow-particle pair (one of the most time-consuming steps in a Wisdom-Holman integrator for small particle numbers).

The MEGNO is then straightforwardly computed from the variations as [@Cincotta2003] $$\begin{aligned}
Y(t) = \frac{2}{t} \int_0^t  t' \frac{\sum_{i=0}^{N-1} \dot{\pmb{\delta}_i}(t') \cdot \pmb{\delta}_i(t')}{\sum_{i=0}^{N-1} \pmb{\delta}_i^2(t')} dt'.\end{aligned}$$ If $Y(t)\rightarrow \infty$, then the system is chaotic. For quasi-periodic orbits, the MEGNO converges to a finite value, $Y(t)\rightarrow 2$ [e.g. @Hinse2010].

One can obtain the Lyapunov characteristic number (LCN), the inverse of the Lyapunov timescale, from the time evolution of the MEGNO via a linear least square fit to $Y(t)$.

Kepler Problem with Variations {#sec:kepler}
------------------------------

In this section, we summarize how to solve the two-body Kepler problem numerically, including the variational equations. Although the solution has been known since the days of Newton, the transcendental nature of Kepler’s equation does not admit a closed-form mathematical expression.

We closely follow the work of [@MikkolaInnanen1999] where the reader can find additional information that we have left out. The equivalent one-body Hamiltonian for the Kepler problem is $$\begin{aligned}
H_{\rm Kepler} &=& \frac12 \mathbf{v}^2 - \frac{M}{|\mathbf{r}|},\end{aligned}$$ where $M$ is the total mass of the two bodies. For consistency with [@MikkolaInnanen1999], we have 