$m'_i$ is the reduced mass of $m_i$ and $M_{i-1}$. Explicit expressions for the momenta can be found by evaluating the time derivative r[Eq. \[eq:p\].]{}

The Jacobi coordinates above are relative coordinates for $i=1,\ldots, N-1$. For the $0$-th coordinate, a different convention is used, $$\begin{aligned}
\mathbf{r}'_0 = \mathbf{R}_{N-1},\quad\quad\quad m'_0 = M_{N-1},\quad\quad\quad \mathbf{p}'_0 = \sum_{j=0}^{N-1} \mathbf{p}_j.\end{aligned}$$ Thus, $\mathbf{r}'_0$ points towards the centre-of-mass of the entire system, $\mathbf{p}'_0$ is the total momentum and $m'_0$ is the total mass.

Hamiltonian Splitting {#sec:splitting}
---------------------

After some algebra, we can rewrite the Hamiltonian in Eq. \[eq:H\] in terms of the conjugate momenta of the Jacobi coordinates [e.g. @solarsystemdynamics; @SussmanWisdom2001]. We only rewrite the kinetic term and keep the potential term expressed as a function of the 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|}.\end{aligned}$$ Note that the kinetic term is still diagonal, i.e. there are no cross terms involving $\mathbf{p}_i\mathbf{p}_j$ with $i\neq j$. Next, we add and subtract the term $$\begin{aligned}
\operatorname{\mathcal{H}}_\pm = \sum_{i=1}^{N-1}\frac{G m'_i M_{i} }{|\mathbf{r}'_i|}. \label{eq:hpm}\end{aligned}$$ After grouping terms in the Hamiltonian, we arrive at $$\begin{aligned}
\mathcal{H} &=& 
\underbrace{\frac{{\mathbf{p}'_0}^2}{2m'_0}\vphantom{\sum_{i=1}^{N-1} \frac{{\mathbf{p}'_i}^2}{2m'_i} }}_{\mathcal{H}_{0}} 
+\underbrace{\sum_{i=1}^{N-1} \frac{{\mathbf{p}'_i}^2}{2m'_i} -\sum_{i=1}^{N-1}\frac{G m'_i M_{i} }{|\mathbf{r}'_i|}}_{\mathcal{H}_{\rm Kepler}} \nonumber \\
&&  + \underbrace{\sum_{i=1}^{N-1}\frac{G m'_i M_{i} }{|\mathbf{r}'_i|}- \sum_{i=0}^{N-1} \sum_{j=i+1}^{N-1} \frac{Gm_im_j}{|\mathbf{r}_i-\mathbf{r}_j|}}_{\mathcal{H}_{\rm Interaction}}.\label{eq:hamsplit}\end{aligned}$$ The 