# Great Inequality of Jupiter and Saturn I: The Planetary Three Body Problem,
Heliocentric development by Lagrange multipliers, Perturbation Theory
Formulation

Jonathan Tot
Department of Mathematics and Statistics, Dalhousie University,
Halifax, Nova Scotia, Canada B3H 4R2
<EMAIL_ADDRESS>

S.R. Valluri
Department of Physics and Astronomy, University of Western Ontario
and Mathematics, King’s University College
London, Ontario, Canada N6A 3K7
<EMAIL_ADDRESS>

P.C. Deshmukh
CAMOST, IIT Tirupati and IISER Tirupati,
Tirupati, Andhra Pradesh 517619, India
<EMAIL_ADDRESS>
###### Abstract

In this paper, we undertake to present a self-contained and thorough analysis
of the gravitational three body problem, with anticipated application to the
Great Inequality of Jupiter and Saturn. The analysis of the three body
Lagrangian is very convenient in heliocentric coordinates with Lagrange
multipliers, the coordinates being the vector-sides $\vec{r}_{i},\,i=1,2,3$ of
the triangle that the bodies form. In two dimensions to begin with, the
equations of motion are formulated into a dynamical system for the polar
angles $\theta_{i}$, angular momenta $\ell_{i}$ and eccentricity vectors
$\vec{e}_{i}$. The dynamical system is simplified considerably by change of
variables to certain auxiliary vector $\vec{f}_{i}=\hat{r}_{i}+\vec{e}_{i}$.
We then begin to formulate the Hamiltonian perturbation theory of the problem,
now in three dimensions. We first give the geometric definitions for the
Delaunay action-angle variables of the two body problem. We express the three
body Hamiltonian in terms of Delaunay variables in each sector $i=1,2,3$,
revealing that it is a nearly integrable Hamiltonian. We then present the KAM
theory perturbative approach that will be followed in future work, including
the modification that will be required because the Hamiltonian is degenerate.

##

###  I Introduction

The three body-problem and the Great Inequality of Jupiter and Saturn are much
studied from the 18th and 19th centuries. They are part of the history of the
development of Celestial Mechanics, inaugurated by the publication of Newton’s
Principia in 1687 [1], and many of the great names that come down to us from
that time were heavily involved, most chiefly the mathematical techniques
proposed and analyses completed by Leonhard Euler, Joseph-Louis Lagrange, and
finally Pierre-Simon Laplace, in his Théorie de Jupiter et de Saturne, written
in 1785 [2]. These advances were enabled by increasing observational accuracy
during the time of John Flamsteed, the first Astronomer Royal, and
astronomical forecasts of such figures as Cassini and Halley [3, p.16].

Following the grand realization, due to Laplace, that the observed
discrepancies of Jupiter and Saturn’s motions from the essentially Keplerian
predictions could actually be accounted for by the mutual Newtonian
gravitation of the planets, other researchers into the 19th century continued
to make improved theories and contributions to the methods, perhaps most
notably the work of G.W. Hill [4] and the series manipulations of Charles
Delaunay [5].

The study of the three body problem, as a mathematical problem in its own
right, also began with the widespread acknowledgement of Newtonian universal
gravitation. Euler (1767) and Lagrange (1772) both initially found particular
periodic solutions, which today are understood to be associated with central
configurations of the system [6]. Much of the early focus was given to the
restricted three body problem, in which one mass is negligible relative to the
other two, so that the larger masses form a Keplerian system, and the third
body orbits in the gravitational field of the first two. It was Euler who
first set the (circular) restricted three body problem in rotating
coordinates, and Lagrange who demonstrated the existence of equilibrium points
within the rotating frame, today known as the five Lagrange points, which
correspond to five orbits for the third mass with the same orbital frequency
as the two-body system.

In the late 1800s, Poincaré studied the general three body problem. Famously,
his work won the prize competition for the 60th birthday of Oscar II, the King
of Sweden and Norway in 1885. Poincaré’s work on the three body problem led
him to consider what are now known as Poincaré sections and first return maps,
and these insights ultimately led to the development of Kolmogorov-Arnold-
Moser theory in the mid-20th century [7]. In the large, three-volume Les
Méthodes Nouvelles de la Mécanique Céleste (1892-99) [8], Poincaré saw, in the
unpredictable nature of three-body problem solutions, the first glimpses of
chaotic dynamics, which dominates much of dynamical systems analysis today
[9].

In this work we present a self-contained analysis of the planetary three body
problem, in which two lighter masses orbit a heavier central body, with
particular application to the Sun-Jupiter-Saturn system. Working with modern
mathematical notation and physical terminology, one aim of this work is
pedagogical in nature, that the subject matter would be more accessible to
modern audiences. For readers who are new to the subject, we put classical
terminology in italics upon their first occurrence and definition in the text.
Jupiter and Saturn’s conjunction in December 2020 gained much attention in the
media [11, 11], so popular presentation of this work should serve to foster
public engagement in mathematics and the sciences.

Crucial to our analysis is the work of Brouke and Lass [10], who formulate the
Lagrangian problem in terms of the three vector-sides of the triangle that the
bodies form, using Lagrange multipliers. In Part I, we present this treatment
of the three body problem in heliocentric coordinates. In §II we show how
total energy and total angular momentum are conserved in this scheme. Staying
in two dimensions in §III, we transform the equations of motion into a system
of first order ODEs for the polar angles, angular momenta and eccentricity
vectors $\vec{e}_{i}$ of the three sectors $i=1,2,3$ of the model. We also
demonstrate how the constraint can be employed to remove the third sector
entirely, leaving equations for only $i=1,2$, corresponding to the planets.
Auxiliary vectors $\vec{f}_{i}=\hat{r}_{i}+\vec{e}_{i}$ for each sector are
introduced in §IV, which make a considerable simplification to the algebraic
form of the dynamical equations. The geometrical properties of the auxiliary
vector is explored. We also present alternate forms of these equations in §V,
in terms of the polar representations $(e_{i},\beta_{i})$, $(f_{i},\psi_{i})$
of eccentricity and auxiliary vectors, respectively.

In Part II we return to three dimensions, and move toward the perturbational
analysis of this problem. In §VI we begin with the geometric definition and
construction of the Delaunay action-angle variables for the two-body problem,
with particular focus on the mean anomaly. In §VIII we present Hamilton’s
equations for the problem, in terms of the the perturbing function
$\mathbf{R}=-\vec{\lambda}\cdot(\vec{r}_{1}+\vec{r}_{2}+\vec{r}_{3})$, where
$\vec{\lambda}$ is the Lagrange multiplier,
$\vec{r}_{1}+\vec{r}_{2}+\vec{r}_{3}=0$ being the constraint. Finally, in §IX
we present the basic approach and the setup of Kolmogorov-Arnold-Moser (KAM)
theory for this problem.

##  Part I The Planetary Three Body Problem

Let the masses of light bodies be $m_{J},m_{S}$ respectively, and $M$ the mass
of a heavier body. Let $\mu=m_{S}/m_{J}$, expected to be $\mathcal{O}(1)$,
while $\epsilon=m_{J}/M$ is small $\epsilon\ll 1$. For Jupiter, Saturn and the
sun these are

$\begin{array}[]{cc}\begin{array}[]{l}M=1.989\times 10^{30}\text{ kg}\\\
m_{J}=1.898\times 10^{27}\text{ kg}\\\ m_{S}=5.683\times 10^{26}\text{
kg}\end{array}&\begin{array}[]{l}\epsilon=9.54..\times 10^{-4}\approx
10^{-3}\\\ \mu=0.2994..\approx 0.3\end{array}\end{array}$ (1)

Let the positions of the bodies in an inertial reference frame be $\vec{X}$
for mass $M$, and $\vec{x}_{J},\vec{x}_{S}$ for $m_{J},m_{S}$. Then the
Lagrangian is

$\mathcal{L}_{\text{in}}=\frac{1}{2}\left(M\left\|\dot{\vec{X}}\right\|^{2}+m_{J}\norm{\dot{\vec{x}}_{J}}^{2}+m_{S}\norm{\dot{\vec{x}}_{S}}^{2}\right)+G\left(\frac{Mm_{J}}{\|\vec{x}_{J}-\vec{X}\|}+\frac{Mm_{S}}{\norm{\vec{x}_{S}-\vec{X}}}+\frac{m_{J}m_{S}}{\norm{\vec{x}_{S}-\vec{x}_{J}}}\right)$
(2)

We will make a change of variables to: the center of mass and heliocentric
coordinates

$\displaystyle\vec{R}$
$\displaystyle=\frac{M\vec{X}+m_{J}\vec{x}_{J}+m_{S}\vec{x}_{S}}{M_{\Sigma}}$
(3) $\displaystyle\vec{r}_{J}$ $\displaystyle=\vec{x}_{J}-\vec{X}$ (4)
$\displaystyle\vec{r}_{S}$ $\displaystyle=\vec{X}-\vec{x}_{S}$ (5)

where $M_{\Sigma}=M+m_{J}+m_{S}$ is the total mass. With these the Lagrangian
becomes

$\displaystyle\mathcal{L}_{\text{in}}=$
$\displaystyle\frac{1}{2}M_{\Sigma}\norm{\dot{\vec{R}}}^{2}+\frac{1}{2}\frac{M}{M_{\Sigma}}\left\\{m_{J}\norm{\dot{\vec{r}}_{J}}^{2}+m_{S}\norm{\dot{\vec{r}}_{S}}^{2}+\frac{m_{J}m_{S}}{M}\norm{\dot{\vec{r}}_{J}+\dot{\vec{r}}_{S}}^{2}\right\\}$
$\displaystyle+\frac{m_{J}M}{M_{\Sigma}}\cdot
GM_{\Sigma}\left(r_{J}^{-1}+\frac{m_{S}/m_{J}}{r_{S}}+\frac{m_{S}/M}{\norm{\vec{r}_{J}+\vec{r}_{S}}}\right)$
(6)

We see the center-of-mass coordinate $\vec{R}$ is cyclic; it’s dynamics
decouple from the other variables, and so can be ignored. We will consequently
drop the first term of (6). Then taking a factor $m_{J}$ out of the kinetic
terms, we can see both remaining terms are proportional to
$\tilde{m}_{J}=\frac{M}{M_{\Sigma}}m_{J}$, which we may regard as a reduced
mass of ‘Jupiter’. Then a reduced and normalized Lagrangian is

$\displaystyle\mathbf{L}_{\text{in}}=\frac{\mathcal{L}_{\text{in}}}{\tilde{m}_{J}}=\frac{1}{2}\left(\norm{\dot{\vec{r}}_{J}}^{2}+\mu\norm{\dot{\vec{r}}_{S}}^{2}+\epsilon\mu\norm{\dot{\vec{r}}_{J}+\dot{\vec{r}}_{S}}^{2}\right)+\alpha\left(r_{J}^{-1}+\mu
r_{S}^{-1}+\frac{\epsilon\mu}{\norm{\vec{r}_{J}+\vec{r}_{S}}}\right)$ (7)

where $\alpha=GM_{\Sigma}$ is the gravitational parameter. Often $\alpha$ is
taken to be 1. However, we will work in unit’s of Jupiter’s average distance
and Jupiter’s year, so that we should take $\alpha=R^{3}\omega^{2}=4\pi^{2}$.

At this stage we recognize 1) that the first and second terms of each
parentheses, what can be called the ‘Jupiter’ and ‘Saturn’ terms, are just
like the terms for the displacement vector of a two body problem, with large
central mass at the origin, and 2) that the third terms, proportional to
$\epsilon\mu$, also look like two-body problem terms, but with displacement
vector $\vec{r}_{SJ}=\pm(\vec{r}_{J}+\vec{r}_{S})$. Taking the minus-sign
option, which gives $\vec{r}_{SJ}=\vec{x}_{S}-\vec{x}_{J}$, we find we are
dealing with a Lagrangian for three independent two-body problems, with
vectors $\vec{r}_{J},\vec{r}_{S}$ and $\vec{r}_{SJ}$ that satisfy the
condition $\vec{r}_{J}+\vec{r}_{S}+\vec{r}_{SJ}=0$. Looking back at our
definition of these vectors in terms of inertial coordinates
$\vec{X},\vec{x}_{J},\vec{x}_{S}$, the constraint is satisfied identically:

$(\vec{x}_{J}-\vec{X})+(\vec{X}-\vec{x}_{S})+(\vec{x}_{S}-\vec{x}_{J})\equiv
0$ (8)

We see that the three vectors are the sides of the triangle that the three
bodies form. The system is analogous to three bodies, of relative masses
$1,\mu$ and $\epsilon\mu$ with respect to the first mass, all orbiting a
central, stationary mass located at the origin, such that the gravitational
parameter for each two-body problem is $\alpha=GM_{\Sigma}$. Of particular
note is that these three supposed bodies orbiting a central mass do not
gravitate to each other. The constraint is maintained by forces that are
sourced by Lagrange multipliers. We thus consider the modified Lagrangian

$\mathbf{L}_{\lambda}=\sum_{i}\left\\{\mu_{i}\left(\frac{1}{2}\norm{\vec{r}_{i}}^{2}+\alpha/r_{i}\right)+\vec{\lambda}\bm{\cdot}\vec{r}_{i}\right\\}$
(9)

where the sum is on $i=J,S,SJ$ or simply $i=1,2,3$, and
$\mu_{i}=1,\mu,\epsilon\mu$. The additional terms are
$\vec{\lambda}\bm{\cdot}\sum_{i=1}^{3}\vec{r}_{i}$, so that the constraint
equation is

$\sum_{i=1}^{3}\vec{r}_{i}\,=0$ (10)

The Euler-Lagrange equations are

$\mu_{i}\ddot{\vec{r}}_{i}=-\mu_{i}\frac{\alpha\,\vec{r}_{i}}{r_{i}^{3}}+\vec{\lambda}\,,\quad
i=1,2,3.$ (11)

The dynamical version of the constraint equation is

$\sum_{i=1}^{3}\ddot{\vec{r}}_{i}\,=0$ (12)

If, in addition to this, we have initial conditions that satisfy both

$\sum_{i=1}^{3}\vec{r}_{i}\,=0\text{ and }\sum_{i=1}^{3}\dot{\vec{r}}_{i}\,=0$
(13)

then the condition (10) will be satisfied for all time. This allows us to
solve for the Lagrange multipliers by taking linear combination of the
equations (11). We have

$0=\sum_{i=1}^{3}\ddot{\vec{r}}_{i}=-\alpha\left(\sum_{i=1}^{3}\frac{\hat{r}_{i}}{r_{i}^{2}}\right)+\left(\sum_{i=1}^{3}\frac{1}{\mu_{i}}\right)\vec{\lambda}$
(14)

and thus

$\vec{\lambda}=\alpha\delta\sum_{i=1}^{3}\frac{\hat{r}_{i}}{r_{i}^{2}}$ (15)

where the coefficient $\delta$ is the reciprocal of the sum of reciprocal
masses

$\delta=\left(\sum_{i=1}^{3}\frac{1}{\mu_{i}}\right)^{-1}=\frac{\epsilon\mu}{1+\epsilon+\epsilon\mu}=\frac{m_{S}}{M_{\Sigma}}=\frac{\tilde{m}_{S}}{M}=2.8536..\times
10^{-4}$ (16)

The equations (11) are

$\ddot{\vec{r}}_{i}=A_{ij}\frac{\hat{r}_{j}}{r_{j}^{2}}\quad\text{(summation
on {\hbox{j}})}$ (17)

where the matrix $A$ of coefficients is

$\displaystyle A$
$\displaystyle=\alpha\left(\begin{array}[]{ccc}-\frac{1+\epsilon}{1+\epsilon+\epsilon\mu}&\delta&\delta\\\
\frac{\epsilon}{1+\epsilon+\epsilon\mu}&-\frac{1+\epsilon\mu}{1+\epsilon+\epsilon\mu}&\frac{\epsilon}{1+\epsilon+\epsilon\mu}\\\
(1+\epsilon+\epsilon\mu)^{-1}&(1+\epsilon+\epsilon\mu)^{-1}&-\frac{\epsilon(1+\mu)}{1+\epsilon+\epsilon\mu}\end{array}\right)$
(21) $\displaystyle=G\left(\begin{array}[]{ccc}-(M+m_{J})&m_{S}&m_{S}\\\
m_{J}&-(M+m_{S})&m_{J}\\\ M&M&-(m_{J}+m_{S})\end{array}\right).$ (25)

That each column sums to $0$ corresponds to
$\ddot{\vec{r}}_{1}+\ddot{\vec{r}}_{2}+\ddot{\vec{r}}_{3}=0$.

###  II Constants of Motion

####  II.1 Energy

The conjugate linear momenta are

$\vec{p}_{i}=\nabla_{\dot{\vec{r}}_{i}}\mathbf{L}_{\lambda}=\mu_{i}\dot{\vec{r}}_{i}$
(26)

and by Legendre transform, the Hamiltonian is

$\displaystyle\mathbf{H}_{\lambda}$
$\displaystyle=\sum_{i}\vec{p}_{i}\bm{\cdot}\dot{\vec{r}}_{i}\,-\mathbf{L}_{\lambda}$
(27)
$\displaystyle=\sum_{i}\left\\{\frac{p_{i}^{2}}{\mu_{i}}-\left(\frac{1}{2}\frac{p_{i}^{2}}{\mu_{i}}+\mu_{i}\frac{\alpha}{r_{i}}+\vec{\lambda}\bm{\cdot}\vec{r}_{i}\right)\right\\}$
(28)
$\displaystyle=\sum_{i}\left\\{\frac{p_{i}^{2}}{2\mu_{i}}-\mu_{i}\frac{\alpha}{r_{i}}-\vec{\lambda}\bm{\cdot}\vec{r}_{i}\right\\}$
(29)

thus when the constraint (10) is satisfied, the following (reduced) energy is
conserved

$\xi=\frac{E}{\tilde{m}_{J}}=\sum_{i}\mu_{i}\left(\frac{1}{2}\norm{\dot{\vec{r}}_{i}}^{2}-\alpha/r_{i}\right)=\sum_{i}\mu_{i}\xi_{i}$
(30)

where $\xi_{i}=\norm{\dot{\vec{r}}_{i}}^{2}/2-\alpha/r_{i}$ is the specific
energy for each sector of the model.

####  II.2 Angular Momentum

Now we put the dynamical vectors $\vec{r}_{i}$ in polar coordinates. Up to
this point, we could have been in three dimensions, but for now we work in two
dimensions. Each vector is $\vec{r}_{i}(t)=r_{i}(t)\hat{\theta}_{i}(t)$, where
$\hat{r}(\theta)=(\cos\theta,\,\sin\theta)^{\text{T}}$

$\displaystyle\norm{\dot{\vec{r}}_{i}}^{2}=\dot{r}_{i}^{2}+r_{i}^{2}\dot{\theta}_{i}^{2}$
(31) $\displaystyle\implies\quad\text{momenta conjugate to }\theta_{i}\text{
are }$ $\displaystyle
h_{i}=\frac{\partial\,\mathbf{L}_{\lambda}}{\partial\dot{\theta}_{i}}=\mu_{i}r_{i}^{2}\dot{\theta}_{i}$
(32) and the Euler-Lagrange equations are
$\displaystyle\dot{h}_{i}=\frac{\partial\,\mathbf{L}_{\lambda}}{\partial\theta_{i}}=\frac{\partial}{\partial\theta_{i}}\left(\vec{\lambda}\bm{\cdot}\vec{r}_{i}\right)$
(33)
$\displaystyle\hphantom{\dot{h}_{i}}=\vec{\lambda}\bm{\cdot}r_{i}\hat{\theta}_{i}$
(34)
$\displaystyle\hphantom{\dot{h}_{i}}=\vec{\lambda}\bm{\cdot}r_{i}T\hat{r}_{i}$
(35) $\displaystyle\hphantom{\dot{h}_{i}}=\vec{\lambda}\bm{\cdot}T\vec{r}_{i}$
(36)

where $\hat{\theta}$ is the vector function
$\hat{\theta}(\varphi)=(-\sin\varphi,\,\cos\varphi)^{\text{T}}$,
$\hat{\theta}_{i}$ is the evaluation $\hat{\theta}(\theta_{i}(t))$, and $T$ is
the $2\times 2$ matrix $T=\left[\begin{array}[]{lr}0&-1\\\
1&0\end{array}\right]$, which is ccw-rotation by $\pi/2$, so that
$\hat{\theta}=T\hat{r}$. This shows that the total angular momentum satisfies

$\dot{h}=\sum_{i=1}^{3}\dot{h}_{i}=\vec{\lambda}\bm{\cdot}T\sum_{i}\vec{r}_{i}=0\text{
by the constraint}.$ (37)

The presence of the rotation matrix $T$ in $\dot{h}$ shows that conservation
of angular momentum is due to the fact that the system as a whole is invariant
under global rotation.

Going forward, it will be best to use specific angular momenta

$\ell_{i}=h_{i}/\mu_{i}=r_{i}^{2}\dot{\theta}_{i}$ (38)

which gives

$\dot{\ell}_{i}=\frac{1}{\mu_{i}}\vec{\lambda}\bm{\cdot}r_{i}\hat{\theta}$
(39)

Using the Lagrange multiplier (15), this is

$\displaystyle\dot{\ell}_{i}$
$\displaystyle=\frac{\alpha\delta}{\mu_{i}}\left(\sum_{j}\frac{\hat{r}_{j}}{r_{j}^{2}}\right)\bm{\cdot}r_{i}\hat{\theta}_{i}$
(40) $\displaystyle\therefore\dot{\ell}_{i}$
$\displaystyle=\frac{\alpha\delta}{\mu_{i}}r_{i}\sum_{j\neq
i}r_{j}^{-2}\sin(\theta_{j}-\theta_{i})$ (41)

we shall label these functions
$\dot{\ell}_{i}=\tau_{i}=\tau_{i}(r_{j},\theta_{j},\ell_{j})$; $\tau$ for
torque. Notice that these equations are perturbations for $i=1,2$, since
$\mu_{1,2}\sim\mathcal{O}(1)$, but this is not the case for $i=3$;
$\mu_{3}=\epsilon\mu\sim\mathcal{O}(\delta)$, so the coefficient in the
equation for $\dot{\ell}_{3}$ is $\mathcal{O}(1)$.

#####  II.2.1 Angular Momentum in Three Dimensions

In three dimensions, we have the vector (total) angular momentum

$\vec{h}=\sum_{i}\vec{r}_{i}\times\vec{p}_{i}=\sum_{i}\mu_{i}\vec{r}_{i}\times\dot{\vec{r}}_{i}$
(42)

the time-derivative of which is

$\displaystyle\dot{\vec{h}}$
$\displaystyle=\sum_{i}\vec{r}_{i}\times\dot{\vec{p}}_{i}$ (43)
$\displaystyle=\sum_{i}\vec{r}_{i}\times\left(-\mu_{i}\alpha\frac{\hat{r}_{i}}{r_{i}^{2}}+\vec{\lambda}\right)$
(44) $\displaystyle=\left(\sum_{i}\vec{r}_{i}\right)\times\vec{\lambda}$ (45)

so angular momentum is conserved by the constraint.

###  III Returning to the Equations of Motion

$\ddot{\vec{r}}_{i}=A_{ij}\frac{\hat{r}_{j}}{r_{j}^{2}}$ (46)

With $\vec{r}_{i}(t)=r_{i}(t)\hat{\theta_{i}(t)}$, we develop these equations
in polar coordinates. We will work out the $\hat{r}_{i}$\- and
$\hat{\theta}_{i}$-components for the $i^{\text{th}}$ equation of (89).

$\displaystyle\dot{\vec{r}}_{i}=\dot{r}_{i}\hat{r}_{i}+r_{i}\dot{\theta}_{i}\hat{\theta}_{i}$
(47) $\displaystyle\text{so
}\norm{\dot{\vec{r}}_{i}}^{2}=\dot{r}_{i}^{2}+r_{i}^{2}\dot{\theta}_{i}^{2}\,\,,\,\text{
and }$
$\displaystyle\ddot{\vec{r}}_{i}=(\ddot{r}_{i}-r_{i}\dot{\theta}_{i}^{2})\hat{r}_{i}+(r_{i}\ddot{\theta}_{i}+2\dot{r}_{i}\dot{\theta}_{i})\hat{\theta}_{i}\,.$
(48)

The $\hat{\theta}_{i}$-component of $\ddot{\vec{r}}_{i}$ is nothing other than
$\dot{\ell}_{i}/r_{i}$; indeed the $\hat{\theta}_{i}$-components of (46) are
just the torque equations (41) derived above. That leaves us with the
$\hat{r}_{i}$-components

$\ddot{r}_{i}-r_{i}\dot{\theta}^{2}_{i}=\left(\sum_{j=1}^{3}A_{ij}\frac{\hat{r}_{j}}{r_{j}^{2}}\right)\bm{\cdot}\hat{r}_{i}$
(49)

with $\dot{\theta}_{i}=\ell_{i}/r_{i}^{2}$ and
$A_{ii}=\alpha(-1+\delta/\mu_{i})$, while $A_{ij}=\alpha\delta/\mu_{i}$ for
$j\neq i$, this is

$\displaystyle\ddot{r}_{i}=\frac{\ell_{i}^{2}}{r_{i}^{3}}+\frac{A_{ii}}{r_{i}^{2}}+\frac{\alpha\delta}{\mu_{i}}\sum_{j\neq
i}\frac{\cos(\theta_{j}-\theta_{i})}{r_{j}^{2}}$ (50)
$\displaystyle\therefore\quad$
$\displaystyle\ddot{r}_{i}=\frac{\ell_{i}^{2}}{r_{i}^{3}}-\alpha
r_{i}^{-2}+\frac{\alpha\delta}{\mu_{i}}\sum_{j=1}^{3}\frac{\cos(\theta_{j}-\theta_{i})}{r_{j}^{2}}$
(51)

Here, in analogy to the definition we have for
$\tau_{i}=\dot{\ell}_{i}=(\alpha\delta/\mu_{i})\,r_{i}\sum_{j}\sin(\theta_{j}-\theta_{i})/r_{j}^{2}$,
we define three functions $\sigma_{i}$ as

$\displaystyle\sigma_{i}$
$\displaystyle=\frac{1}{\mu_{i}}\vec{\lambda}\bm{\cdot}\vec{r}_{i}$
$\displaystyle=\frac{\alpha\delta}{\mu_{i}}\,r_{i}\sum_{j=1}^{3}\frac{\cos(\theta_{j}-\theta_{i})}{r_{j}^{2}}$
(52)

then what we have is

$\ddot{r}_{i}=\frac{\ell_{i}^{2}}{r_{i}^{3}}-\alpha
r_{i}^{-2}+\frac{\sigma_{i}}{r_{i}}$ (53)

and define the right-hand sides as functions $a_{i}=a_{i}(\theta,\ell,r)$

####  III.1 Eccentricity or Laplace-Runge-Lenz Vectors

At this stage, we may formulate our differential equations as a system of 12
first order ODEs

$\displaystyle\dot{\theta_{i}}$ $\displaystyle=\Omega_{i}=\ell_{i}/r_{i}^{2}$
(54) $\displaystyle\dot{\ell}_{i}$
$\displaystyle=\tau_{i}=\frac{\alpha\delta}{\mu_{i}}r_{i}\sum_{j\neq
i}r_{j}^{-2}\sin(\theta_{j}-\theta_{i})$ (55) $\displaystyle\dot{r}_{i}$
$\displaystyle=v_{i}$ (56) $\displaystyle\dot{v}_{i}$
$\displaystyle=a_{i}=\frac{\ell_{i}^{2}}{r_{i}^{3}}-\alpha
r_{i}^{-2}+\frac{\sigma_{i}}{r_{i}}$ (57)

Now we introduce a change of variables from $(r_{i},v_{i})$ to eccentricity
vectors, in each sector

$\vec{e}=\dot{\vec{r}}\times\vec{\ell}/\alpha-\hat{r}.$ (58)

This is a normalization of the Laplace-Runge-Lenz vector

$\vec{A}=\vec{p}\times\vec{L}-m^{2}\alpha\,\hat{r}=m^{2}\alpha\,\vec{e}$ (59)

where $\vec{p}$ is linear momentum and $\vec{L}$ is the (dimensionfull)
angular momentum. In the two-body problem, these vectors are constants of
motion. The eccentricity vector has magnitude equal to that of the
eccentricity of the Keplerian orbit and points in the direction of periapsis.
Working in two dimensions, our specific angular momenta are out of the plane,
$\vec{\ell}=\ell\hat{z}$, so that in polar coordinates

$\displaystyle\vec{e}$
$\displaystyle=(v\,\hat{r}+r\,\dot{\theta}\,\hat{\theta})\times(\ell\,\hat{z}/\alpha)-\hat{r}$
(60) $\displaystyle=\left(\frac{\ell^{2}/\alpha}{r}-1\right)\hat{r}-\ell
v/\alpha\,\hat{\theta}$ (61)

(where $\ell/r^{2}$ has been substituted for $\dot{\theta}$). This defines the
eccentricity vector in terms of it’s polar components

$e^{r}=\frac{\ell^{2}/\alpha}{r}-1,\quad e^{\theta}=-\ell\,v/\alpha$ (62)

the reverse change of variables being

$\displaystyle r$ $\displaystyle=\frac{\ell^{2}/\alpha}{1+e^{r}}$ (63)
$\displaystyle v$ $\displaystyle=-\alpha\frac{e^{\theta}}{\ell}.$ (64)

Equation (63) is precisely the form of a Keplerian elliptic orbit if $\ell$
and $\vec{e}$ are constant, as in the two-body problem, since

$\displaystyle r=\frac{\ell^{2}/\alpha}{1+e^{r}}$
$\displaystyle=\frac{\ell^{2}/\alpha}{1+\vec{e}\bm{\cdot}\hat{r}}$ (65)
$\displaystyle=\frac{\ell^{2}/\alpha}{1+e^{x}\cos\theta+e^{y}\sin\theta}$ (66)
$\displaystyle=\frac{\ell^{2}/\alpha}{1+e\cos(\theta-\beta)}$ (67)

where $\beta$ is the angle of $\vec{e}\,$ from the positive $x$-axis, called
the longitude of periapsis. Equation $(\ref{osculating})$ thus describes the
osculating orbit to the trajectory, which is the elliptic orbit that a body
would follow given it’s instantaneous angular momentum and eccentricity. As a
further note, if we take the time-derivative of (63), writing $\tau$ for
$\dot{\ell}$, we find

$v=\dot{r}=\frac{1}{\alpha}\left(\frac{2\ell\tau}{1+e^{r}}-\frac{\ell^{2}}{(1+e^{r})^{2}}\left(\dot{\vec{e}}\bm{\cdot}\hat{r}+\frac{\ell}{r^{2}}\vec{e}\bm{\cdot}\hat{\theta}\right)\right)=\frac{1}{\alpha}\left(\frac{2\ell\tau}{1+e^{r}}-\frac{r^{2}}{\ell^{2}}\dot{\vec{e}}\bm{\cdot}\hat{r}\right)-\alpha\frac{e^{\theta}}{\ell}$
(68)

which reduces to precisely (64) if $\tau$ and $\dot{\vec{e}}$ are 0.

####  III.2 System of Equations using Eccentricity vectors

We can differentiate the definitions (58) and use the system (54-57) to derive
first order differential equations for the eccentricity vectors
$\dot{\vec{e}}_{i}=\dots(\theta,\ell,\vec{e}\,)$. First working in one sector
(no indices $i$), writing $\tau$ for $\dot{\ell}$, we have

$\displaystyle\dot{\vec{e}}$
$\displaystyle=\frac{1}{\alpha}\left(\frac{2\ell\tau}{r}-\cancel{\frac{\ell^{2}v}{r}}\,\right)\hat{r}+\left(\frac{\ell^{2}/\alpha}{r}-1\right)\frac{\ell}{r^{2}}\hat{\theta}-(\tau\,v+\ell\,\ddot{r})\,\hat{\theta}/\alpha+\cancel{\frac{\ell
v}{\alpha}\frac{\ell}{r^{2}}\hat{r}}$ (69)
$\displaystyle=\frac{2\ell\tau}{\alpha\,r}\hat{r}-\frac{\ell}{\alpha}\left(\ddot{r}-\frac{\ell^{2}}{r^{3}}+\alpha
r^{-2}+v\frac{\tau}{\ell}\right)\hat{\theta}.$ (70) Substituting $v=-\alpha
e^{\theta}/\ell$ and $\ddot{r}=a=\ell^{2}/r^{3}-\alpha r^{-2}+\sigma/r$ gives
$\displaystyle\dot{\vec{e}}$
$\displaystyle=\frac{2\ell\tau}{\alpha\,r}\hat{r}+\left(\frac{\tau}{\ell}\,e^{\theta}-\frac{\ell\sigma}{\alpha\,r}\right)\hat{\theta}.$
(71) Pair half of the first term with the second $\hat{\theta}$-term
$\displaystyle\dot{\vec{e}}$
$\displaystyle=\frac{\ell\tau}{\alpha\,r}\hat{r}+\frac{\tau}{\ell}e^{\theta}\hat{\theta}+\frac{\ell}{\alpha\,r}\left(\tau\hat{r}-\sigma\hat{\theta}\right).$
(72) In the first term we substitute $\ell/\alpha\,r=(1+e^{r})/\ell$
$\displaystyle\dot{\vec{e}}$
$\displaystyle=\frac{\tau}{\ell}(1+e^{r})\hat{r}+\frac{\tau}{\ell}e^{\theta}\hat{\theta}+\frac{\ell}{\alpha\,r}\left(\tau\hat{r}-\sigma\hat{\theta}\right)$
(73)
$\displaystyle=\frac{\tau}{\ell}\left(\hat{r}+e^{r}\hat{r}+e^{\theta}\hat{\theta}\right)+\frac{\ell}{\alpha\,r}\left(\tau\hat{r}-\sigma\hat{\theta}\right).$
(74) Thus we find $\displaystyle\dot{\vec{e}}$
$\displaystyle=\frac{\tau}{\ell}\left(\vec{e}+\hat{r}\right)+\frac{\ell}{\alpha\,r}\left(\tau\,\hat{r}-\sigma\,\hat{\theta}\right).$
(75)

The combination $(\tau\,\hat{r}-\sigma\,\hat{\theta})/r$ is, in each sector

$\displaystyle(\tau_{i}\,\hat{r}_{i}-\sigma_{i}\,\hat{\theta}_{i})/r_{i}=$
$\displaystyle\left(\frac{\alpha\delta}{\mu_{i}}\sum_{j=1}^{3}\frac{\sin(\theta_{j}-\theta_{i})}{r_{j}^{2}}\right)\left[\begin{array}[]{c}\cos\theta_{i}\\\
\sin\theta_{i}\end{array}\right]$ (78)
$\displaystyle-\left(\frac{\alpha\delta}{\mu_{i}}r_{i}\sum_{j=1}^{3}\frac{\cos(\theta_{j}-\theta_{i})}{r_{j}^{2}}\right)\left[\begin{array}[]{c}-\sin\theta_{i}\\\
\cos\theta_{i}\end{array}\right]$ (81) $\displaystyle=$
$\displaystyle\frac{\alpha\delta}{\mu_{i}}\,\sum_{j=1}^{3}r_{j}^{-2}\left[\begin{array}[]{c}\sin\theta_{j}\\\
-\cos\theta_{j}\end{array}\right]$ (84) $\displaystyle=$
$\displaystyle-\frac{\alpha\delta}{\mu_{i}}\,\sum_{j=1}^{3}r_{j}^{-2}\hat{\theta}_{j}=-\frac{\alpha\delta}{\mu_{i}}\,\sum_{j=1}^{3}r_{j}^{-2}\,T\hat{r}_{j}$
(85)
$\displaystyle\therefore\,(\tau_{i}\,\hat{r}_{i}-\sigma_{i}\,\hat{\theta}_{i})/r_{i}=$
$\displaystyle-T\vec{\lambda}/\mu_{i}$ (86)

So finally, our differential equations for $\dot{\vec{e}}_{i}$ are

$\dot{\vec{e}}_{i}=\frac{\tau_{i}}{\ell_{i}}\left(\vec{e}_{i}+\hat{r}_{i}\right)-\frac{\ell_{i}}{\alpha\mu_{i}}\,T\vec{\lambda}.$
(87)

Observe that these equations are $\propto\,\delta/\mu_{i}$, so that
$\dot{\vec{e}}_{i}$ are $\mathcal{O}(\delta)$ for $i=1,2$, i.e. Jupiter and
Saturn, while $\dot{\vec{e}}_{3}$ is $\mathcal{O}(1)$. The system of twelve
first order ODEs for the variables $(\theta_{i},\ell_{i},\vec{e}_{i})$ is

$\displaystyle\dot{\theta_{i}}$ $\displaystyle=\Omega_{i}=\ell_{i}/r_{i}^{2}$
(88) $\displaystyle\dot{\ell}_{i}$
$\displaystyle=\tau_{i}=\frac{\alpha\delta}{\mu_{i}}r_{i}\sum_{j\neq
i}r_{j}^{-2}\sin(\theta_{j}-\theta_{i})$ (89) $\displaystyle\dot{\vec{e}}_{i}$
$\displaystyle=\frac{\tau_{i}}{\ell_{i}}\left(\vec{e}_{i}+\hat{r}_{i}\right)-\frac{\ell_{i}}{\alpha\mu_{i}}T\vec{\lambda}\quad\text{for
}i=1,2,3$ (90) $\displaystyle\text{where
}r_{i}=\ell_{i}^{2}/\left[\alpha(1+e^{r}_{i})\right]$

Mixing between the sectors enters through the torques $\tau_{i}$ and
$\vec{\lambda}$.

####  III.3 Reduction by the Constraint to Two Sectors

Given that solutions will satisfy the constraint $0=\sum_{i}\vec{r}_{i}$, we
can write down algebraic/trigonometric expressions for the variables of the
3rd sector $\vec{r}_{3}=-\vec{r}_{1}-\vec{r}_{2}$ in terms of the 1st and 2nd
sectors. Substituting these relations into the $i=1,2$ equations would leave
the system (88-90) for $i=1,2$ only, and the equations for
$\dot{\ell}_{i},\dot{\vec{e}}_{i}$ would all be $\mathcal{O}(\delta)$.

In particular, if $\vec{r}_{3}=-\vec{r}_{1}-\vec{r}_{2}$, then we can write
expressions for the radius and trig ratios of the argument of $\vec{r}_{3}$ in
terms of those of $\vec{r}_{1},\vec{r}_{2}$. First of all, by the cosine-law

$r_{3}^{2}=r_{1}^{2}+r_{2}^{2}+2r_{1}r_{2}\cos(\theta_{2}-\theta_{1})$ (91)

and then from basic trigonometry, we can express
$\cos\theta_{3},\sin\theta_{3}$ as follows

$\displaystyle\cos\theta_{3}=-\frac{r_{1}\cos\theta_{1}+r_{2}\cos\theta_{2}}{r_{3}}$
(92)
$\displaystyle\sin\theta_{3}=-\frac{r_{1}\sin\theta_{1}+r_{2}\sin\theta_{2}}{r_{3}}$
(93)

These may then be worked into the equations (89,90), still using
$r_{i}=\alpha^{-1}\ell_{i}^{2}/(1+e^{r}_{i})$ but now only for $i=1,2$. This
results in a system of first order ODEs
$\dot{\theta}_{1}=\Omega_{1},\,\dot{\theta}_{2}=\Omega_{2},\,\dot{\ell}_{1}=\tau_{1},\,\dot{\ell}_{2}=\tau_{2},\dot{\vec{e}}_{1}=\vec{v}_{1},\dot{\vec{e}}_{2}=\vec{v}_{2}$
where the angular velocities $\Omega_{i}$ are $\mathcal{O}(1)$, but the
torques $\tau_{i}$ and eccentricity-velocities $\vec{v}_{i}$ are
$\mathcal{O}(\delta)$, facilitating a multiple-scales analysis.

###  IV Auxiliary vectors $\vec{f}=\hat{r}+\vec{e}$

The expressions for the right-hand-sides of (88-90), especially the torques
and eccentricity-velocities, in terms of the variables
$(\theta,\ell,e^{r},e^{\theta})$, are very large and cumbersome. Written as
rational functions, the numerators and denominators expanded out, are
respectively 16 and 9 terms for the torques and 222 and 9 terms for the
eccentricity-velocities, not including terms within $3/2$-roots of the
denominators.

The prevalence of the combination $\hat{r}_{i}+\vec{e}_{i}$ presents an
opportunity for simplification. Observe that the denominator of the expression
for the osculating orbit (63) is the $\hat{r}_{i}$-component of this vector

$r=\frac{\ell^{2}/\alpha}{1+e^{r}}=\frac{\ell^{2}/\alpha}{(\hat{r}+\vec{e}\,)\cdot\hat{r}}.$
(94)

The combination $\hat{r}_{i}+\vec{e}_{i}$ also occurs in the DEs for
eccentricity vectors (87). If we change variables from eccentricity vectors
$\vec{e}_{i}$ to these auxiliary vectors
$\vec{f}_{i}=\hat{r}_{i}+\vec{e}_{i}$, the modified equations become

$\displaystyle r_{i}$
$\displaystyle=\frac{\ell_{i}^{2}}{\alpha\,f_{i}^{r}}=\frac{\ell_{i}^{2}}{\alpha\,\vec{f}_{i}\cdot\hat{r}_{i}}=\frac{\ell_{i}^{2}}{\alpha\,f_{i}\cos(\psi_{i}-\theta_{i})}$
(95) $\displaystyle\dot{\vec{f}}_{i}$
$\displaystyle=\frac{\tau_{i}}{\ell_{i}}\vec{f}_{i}-\frac{\ell_{i}}{\alpha\mu_{i}}(T\vec{\lambda})+\frac{d}{dt}\hat{r}_{i}$
(96)
$\displaystyle=\frac{\tau_{i}}{\ell_{i}}\vec{f}_{i}-\frac{\ell_{i}}{\alpha\mu_{i}}(T\vec{\lambda})+\Omega_{i}\,\hat{\theta}_{i}$
(97)

where $f_{i}=\norm{\vec{f}_{i}}$ and $\psi_{i}$ is the argument of
$\vec{f}_{i}$ (such that $\vec{f}_{i}/f_{i}=\hat{r}(\psi_{i})$). Let the point
be laboured, that given this definition of $\vec{f}_{i}$, the differential
equation is

$\dot{\vec{f}}_{i}=\frac{d}{dt}\hat{r}_{i}+\mathcal{O}(\delta).$ (98)

So to leading order, the vector $\vec{f}_{i}$ will be very nearly equal to the
unit vector $\hat{r}_{i}(t)=\hat{r}(\theta_{i}(t))$. Indeed, the solution is
exactly $\vec{f}_{i}=\hat{r}_{i}+\vec{e}_{i}$. So solutions will have
$f_{i}\sim 1$ and $\psi_{i}\sim\theta_{i}$ for small initial eccentricity, at
least for a finite duration after initial conditions.

####  IV.1 Geometric relationship of the Auxiliary vectors

(a)

(b)

Figure 1: These figures shows the elliptic geometry of the auxiliary vector.
In a) $\vec{f}=\hat{r}+\vec{e}$, where $\vec{e}$ is the eccentricity vector,
shown to the left of $A$. This ellipse has an eccentricity of 0.44 and semi-
major axis equal to 1.02. The points $O,O^{\prime}$ are the foci, and $C$ is
the center. The point $P$ represents the location of a planet, so the vector
from $O$ to $P$ is $\vec{r}$. The point $E$ on the semi-major axis is the
location of $\vec{r}-r\vec{f}=-r\vec{e}$. The point $F$ on the blue curve is
the coordinate of $\vec{r}-\vec{f}$, and the blue curve is the locus of such
points, for all possible points $P$ on the ellipse. The black interval along
the $x$-axis is the locus of points $E$. As a subset of the $x$-axis, it is
the interval $[ae(1-e),ae(1+e)]$. In b) the ellipse has eccentricity $e=0.6$
and $a=1$; the equation $a\hat{r}=-a\vec{e}+a\vec{f}$ is demonstrated, where
$a\hat{r}=\vec{OQ}$ is the vector position of the orbit $\vec{OP}$ continued
to circle about $O$, and $\vec{CQ}$ is the auxiliary vector $a\vec{f}$.

For ellipses with even moderate eccentricity, up to $\sim 0.5$, the vector
$-\vec{f}=-\hat{r}-\vec{e}$ points, to leading order, from the position of the
planet towards the center of the ellipse. Indeed, the vector
$r\vec{f}=\vec{r}+r\vec{e}$ is such that

$\vec{r}-r\vec{f}=-r\vec{e},$ (99)

while the coordinate of the center of the ellipse is $-a\vec{e}$, where $a$ is
the semi-major axis, as demonstrated in Fig. 1(a). So the degree to which
these coincide is the degree to which $r$ and $a$ agree. In terms of
eccentricity, semi-major axis and true anomaly $\nu=\theta-\beta$, the
relationship is

$r=\frac{a(1-e^{2})}{1+e\cos\nu}.$ (100)

At minimum $re=ae(1-e)$, while the maximum value is $ae(1+e)$. Thus the
location $\vec{r}-r\vec{f}$ lies within a segment of the major-axis which is a
length $a\,e^{2}$ on either side of the ellipse center, as shown in Fig. 1(a).
The error of $-\vec{f}$ pointing from the location of the planet to the center
of the orbit is $\mathcal{O}(e^{2})$. Specifically, the difference in angle of
$-r\vec{e}$ vs. $-a\vec{e}$ as seen from the position $\vec{r}$—in other
words, the angle $\angle CPE$—is $e^{2}\sin\nu\cos\nu$ to leading order in
$e$.

Moreover, if we consider the equation

$a\hat{r}=-a\vec{e}+a\vec{f},$ (101)

then we see the following geometry: construct an circle around the focus of
the ellipse, with radius $a$. If we continue the vector $\vec{r}$ from the
focus out to radius $a$, we reach the point $a\hat{r}$ on the circle. The
position of the ellipse centre from the focus is the first term $-a\vec{e}$.
Thus we see from (101) that the vector $a\vec{f}$ points from the centre of
the ellipse to the position of the orbit projected radially to radius $a$, as
shown in Fig. 1(b).

####  IV.2 The Dynamical System in terms of the Auxiliary vectors

What might simplify the equations the most is to write the equations for
$\dot{\vec{f}}_{i}$ in polar form, for the components
$\dot{f}_{i}^{r}=\frac{d}{dt}(\hat{r}_{i}\cdot\vec{f}_{i})$ and
$\dot{f}_{i}^{\theta}=\frac{d}{dt}(\hat{\theta}_{i}\cdot\dot{\vec{f}}_{i})$.
That is, we have

$\displaystyle\vec{f}_{i}$
$\displaystyle=f_{i}^{r}\hat{r}_{i}+f_{i}^{\theta}\hat{\theta}_{i}$ (102)
$\displaystyle\text{and }\dot{\vec{f}}_{i}$
$\displaystyle=\dot{f}_{i}^{r}\hat{r}_{i}+f_{i}^{r}\left(\Omega_{i}\hat{\theta}_{i}\right)+\dot{f}_{i}^{\theta}\hat{\theta}_{i}+f_{i}^{\theta}\left(-\Omega_{i}\hat{r}_{i}\right)$
(103)
$\displaystyle=\left(\dot{f}_{i}^{r}-\Omega_{i}f_{i}^{\theta}\right)\hat{r}_{i}+\left(\dot{f}_{i}^{\theta}+\Omega_{i}f_{i}^{r}\right)\hat{\theta}_{i}$
(104)

The equation (97) for $\dot{\vec{f}}_{i}$ becomes

$\displaystyle\dot{f}_{i}^{r}$
$\displaystyle=\frac{\tau_{i}}{\ell_{i}}f_{i}^{r}+\Omega_{i}f_{i}^{\theta}-\frac{\ell_{i}}{\alpha\mu_{i}}\left(\hat{r}_{i}\cdot
T\vec{\lambda}\right)$ (105) $\displaystyle\dot{f}_{i}^{\theta}$
$\displaystyle=\frac{\tau_{i}}{\ell_{i}}f_{i}^{\theta}-\Omega_{i}f_{i}^{r}-\frac{\ell_{i}}{\alpha\mu_{i}}\left(\hat{\theta}_{i}\cdot
T\vec{\lambda}\right)+\Omega_{i}.$ (106)

We know the components of $T\vec{\lambda}$ from (86)111Multiplying (86) by
$-T$, we can also find
$\vec{\lambda}=\frac{\mu_{i}}{r_{i}}\left(\sigma_{i}\hat{r}_{i}+\tau_{i}\hat{\theta}_{i}\right)$.

$\displaystyle\hat{r}_{i}\cdot(-T\vec{\lambda})=\mu_{i}\,\tau_{i}/r_{i}$ (107)
$\displaystyle\hat{\theta}_{i}\cdot(-T\vec{\lambda})=-\mu_{i}\,\sigma_{i}/r_{i}$
(108)

which gives

$\displaystyle\dot{f}_{i}^{r}$
$\displaystyle=\frac{\tau_{i}}{\ell_{i}}f_{i}^{r}+\frac{\ell_{i}\tau_{i}}{\alpha\,r_{i}}+\Omega_{i}f_{i}^{\theta}$
(109) $\displaystyle\dot{f}_{i}^{\theta}$
$\displaystyle=\frac{\tau_{i}}{\ell_{i}}f_{i}^{\theta}-\frac{\ell_{i}\sigma_{i}}{\alpha\,r_{i}}+\Omega_{i}(1-f_{i}^{r}).$
(110)

####  IV.3 Final substitutions

Things may be made yet more concise. The osculating orbits are given by
$r_{i}=\ell_{i}^{2}/\alpha f_{i}^{r}$, and the prevalent coefficients in
(109,110) become

$\displaystyle\Omega_{i}$
$\displaystyle=\frac{\ell_{i}}{r_{i}^{2}}=\alpha^{2}\frac{{f_{i}^{r}}^{2}}{\ell_{i}^{3}}$
(111) $\displaystyle\frac{\ell_{i}}{\alpha\,r_{i}}$
$\displaystyle=\frac{f_{i}^{r}}{\ell_{i}}$ (112)

so that we find

$\displaystyle\dot{f}_{i}^{r}$
$\displaystyle=2\frac{\tau_{i}}{\ell_{i}}f_{i}^{r}+\alpha^{2}\frac{{f_{i}^{r}}^{2}}{\ell_{i}^{3}}f_{i}^{\theta}$
(113) $\displaystyle\dot{f}_{i}^{\theta}$
$\displaystyle=\left(-\frac{\sigma_{i}}{\ell_{i}}+\alpha^{2}\frac{{f_{i}^{r}}(1-f_{i}^{r})}{\ell_{i}^{3}}\right)f_{i}^{r}+\frac{\tau_{i}}{\ell_{i}}f_{i}^{\theta}.$
(114)

For completeness, we give the expressions for $\sigma$ and $\tau$ in these
terms

$\displaystyle\tau_{i}$
$\displaystyle=\frac{\alpha^{2}\delta}{\mu_{i}}\frac{\ell_{i}^{2}}{f_{i}^{r}}\sum_{j\neq
i}\frac{{f_{j}^{r}}^{2}}{\ell_{j}^{4}}\sin(\theta_{j}-\theta_{i})$ (115)
$\displaystyle\sigma_{i}$
$\displaystyle=\frac{\alpha^{2}\delta}{\mu_{i}}\left[\frac{f_{i}^{r}}{\ell_{i}^{2}}+\frac{\ell_{i}^{2}}{f_{i}^{r}}\sum_{j\neq
i}\frac{{f_{j}^{r}}^{2}}{\ell_{j}^{4}}\cos(\theta_{j}-\theta_{i})\right].$
(116)

Equations (115,116) are best read in the 3-sector version of this problem
(that is, if one does not eliminate $r_{3}$ and $\theta_{3}$). If one wishes
to work in only two sectors $i=1,2$, the substitutions (91,92,93) should be
made into (41) and (52).

We make the substitutions (111,112) into (87) and (97) as well, which become

$\displaystyle\dot{\vec{e}}_{i}$
$\displaystyle=\frac{\tau_{i}}{\ell_{i}}(\vec{e}_{i}+\hat{r}_{i})+\frac{f_{i}^{r}}{\ell_{i}}\left(\tau_{i}\hat{r}_{i}-\sigma_{i}\hat{\theta}_{i}\right)$
(117) $\displaystyle\dot{\vec{f}}_{i}$
$\displaystyle=\dot{\vec{e}}_{i}+\Omega_{i}\hat{\theta}_{i}$
$\displaystyle=\frac{\tau_{i}}{\ell_{i}}\vec{f}_{i}+\frac{f_{i}^{r}}{\ell_{i}}\left(\tau_{i}\hat{r}_{i}-\sigma_{i}\hat{\theta}_{i}\right)+\alpha^{2}\frac{{f_{i}^{r}}^{2}}{\ell_{i}^{3}}\hat{\theta}_{i}$
(118)

In full then, the system of equations, in terms of the polar components of the
auxiliary vectors $f_{i}^{r},f_{i}^{\theta}$ are

$\displaystyle\dot{\theta_{i}}$
$\displaystyle=\Omega_{i}=\frac{\alpha^{2}{f_{i}^{r}}^{2}}{l_{i}^{3}}$ (119)
$\displaystyle\dot{\ell}_{i}$ $\displaystyle=\tau_{i}=(\ref{taus})$ (120)
$\displaystyle\dot{f}_{i}^{r}$
$\displaystyle=2\frac{\tau_{i}}{\ell_{i}}f_{i}^{r}+\alpha^{2}\frac{{f_{i}^{r}}^{2}}{\ell_{i}^{3}}f_{i}^{\theta}$
(121) $\displaystyle\dot{f}_{i}^{\theta}$
$\displaystyle=\left(-\frac{\sigma_{i}}{\ell_{i}}+\alpha^{2}\frac{{f_{i}^{r}}(1-f_{i}^{r})}{\ell_{i}^{3}}\right)f_{i}^{r}+\frac{\tau_{i}}{\ell_{i}}f_{i}^{\theta}.$
(122)

Numerical solutions of (119-122) for Jupiter’s eccentricity vector are shown
in Fig.

(a)

(b)

(c)

(d)

Figure 2: Numerical solutions of the system (119-122) for Jupiter’s
eccentricity vector
$\vec{e}_{1}=(f_{1}^{r}(t)-1)\hat{r}(\theta_{1}(t))+f_{1}^{\theta}(t)\hat{\theta}(\theta_{i}(t))$,
over different time-scales, with initial conditions based on the J2000
Epoch[15]. 2(a) shows the solution for 5 Jupiter-years (Jyrs), which is the
cycle of the Jupiter-Saturn conjunctions, and reveals a repeating triangular
pattern. 2(b) is for 240 Jyrs, showing that the small ‘scallops’ of the final
plot happen on the scale of every 100 Jyrs. 2(c) plots from $t=0$ to $t=5000$
and shows that the large ‘scallops’ are on the order of 4500 Jyrs. 2(d) plots
for 27 kJyrs, and reveals a hexagonal pattern; indeed, 27,000/4500=6.

###  V Other formulations

For completeness, we present the following forms of the equations for both
$\dot{\vec{e}}_{i}$ and $\dot{\vec{f}}_{i}$. First of all, for
$\dot{\vec{e}}_{i}$ in the polar coordinates as above

$\displaystyle\dot{e}_{i}^{r}$
$\displaystyle=2\frac{\tau_{i}}{\ell_{i}}(1+e^{r}_{i})+\alpha^{2}\frac{(1+e_{i}^{r})^{2}}{\ell_{i}^{3}}e_{i}^{\theta}$
(123) $\displaystyle\dot{e}_{i}^{\theta}$
$\displaystyle=-\left(\frac{\sigma_{i}}{\ell_{i}}+\alpha^{2}\frac{e_{i}^{r}(1+e_{i}^{r})}{\ell_{i}^{3}}\right)(1+e_{i}^{r})+\frac{\tau_{i}}{\ell_{i}}e_{i}^{\theta}\,.$
(124)

Of course, these are related to (113,114) as the components of
$\vec{e}_{i},\vec{f}_{i}$ are related by

$f_{i}^{r}=1+e_{i}^{r},\quad f_{i}^{\theta}=e_{i}^{\theta}.$ (125)

Another option, using either the eccentricity or auxiliary vectors, is to
express them in their own polar coordinates; that is,

$\displaystyle\vec{e}_{i}(t)$ $\displaystyle=e_{i}(t)\hat{r}(\beta_{i}(t))$
(126) $\displaystyle\vec{f}_{i}(t)$
$\displaystyle=f_{i}(t)\hat{r}(\psi_{i}(t)).$ (127)

Here, $e_{i}$ are the eccentricities (of the osculating orbits) and
$\beta_{i}$ are the longitudes of periapsis, whereas $f_{i}$ are the lengths
of the auxiliary vectors and $\psi_{i}$ are the angles which $\vec{f}_{i}$
make with the fixed reference (positive $x$-axis, for instance). These
definitions give

$\displaystyle\dot{\vec{e}}_{i}=\dot{e}_{i}\,\hat{r}_{\beta_{i}}+e_{i}\,\dot{\beta}_{i}\,\hat{\theta}_{\beta_{i}}$
(128)
$\displaystyle\dot{\vec{f}}_{i}=\dot{f}_{i}\,\hat{r}_{\psi_{i}}+f_{i}\,\dot{\psi_{i}}\,\hat{\theta}_{\psi_{i}}$
(129)

where here the notation $\hat{r}_{\times}$ has been used for
$\hat{r}(\times(t))$. From these, we derive differential equations for norms
($e_{i},f_{i}$) and arguments ($\beta_{i},\psi_{i}$) by taking
$\hat{r}_{\beta_{i}},\hat{\theta}_{\beta_{i}}$-components222These are derived
with the following helpful properties of the $\hat{r},\hat{\theta}$ functions:
$\hat{r}(a)\cdot\hat{r}(b)=\hat{\theta}(a)\cdot\hat{\theta}(b)=\cos(a-b)$,
$\hat{r}(a)\cdot\hat{\theta}(b)=\sin(a-b)$. of (117) and
$\hat{r}_{\psi_{i}},\hat{\theta}_{\psi_{i}}$-components of (118)

$\displaystyle\dot{e}_{i}$
$\displaystyle=\hat{r}_{\beta_{i}}\cdot\dot{\vec{e}}_{i}=\frac{\tau_{i}}{\ell_{i}}\left[e_{i}+(1+f_{i}^{r})\cos(\theta_{i}-\beta_{i})\right]+\frac{\sigma_{i}}{\ell_{i}}f_{i}^{r}\sin(\theta_{i}-\beta_{i})$
(130) $\displaystyle e_{i}\,\dot{\beta}_{i}$
$\displaystyle=\hat{\theta}_{\beta_{i}}\cdot\dot{\vec{e}}_{i}=\frac{\tau_{i}}{\ell_{i}}(1+f_{i}^{r})\sin(\theta_{i}-\beta_{i})-\frac{\sigma_{i}}{\ell_{i}}f_{i}^{r}\cos(\theta_{i}-\beta_{i})$
(131)

while the equations for $f_{i},\psi_{i}$ are

$\displaystyle\dot{f}_{i}=\hat{r}_{\psi_{i}}\cdot\dot{\vec{f}}_{i}$
$\displaystyle=\frac{\tau_{i}}{\ell_{i}}\left(f_{i}+f_{i}^{r}\cos(\psi_{i}-\theta_{i})\right)+\frac{f_{i}^{r}}{\ell_{i}}\left(\alpha^{2}\frac{f_{i}^{r}}{\ell_{i}^{2}}-\sigma_{i}\right)\sin(\psi_{i}-\theta_{i})$
(132) $\displaystyle
f_{i}\,\dot{\psi}_{i}=\hat{\theta}_{\psi_{i}}\cdot\dot{\vec{f}}_{i}$
$\displaystyle=-\frac{\tau_{i}f_{i}^{r}}{\ell_{i}}\sin(\psi_{i}-\theta_{i})+\frac{f_{i}^{r}}{\ell_{i}}\left(\alpha^{2}\frac{f_{i}^{r}}{\ell_{i}^{2}}-\sigma_{i}\right)\cos(\psi_{i}-\theta_{i}).$
(133)

The systems
$\\{\dot{\theta}_{i},\dot{\ell}_{i},\dot{e}_{i},\dot{\beta}_{i}\\}$ and
$\\{\dot{\theta}_{i},\dot{\ell}_{i},\dot{f}_{i},\dot{\psi}_{i}\\}$ can be
closed by the relations

$f_{i}^{r}=f_{i}\cos(\psi_{i}-\theta_{i})=1+e_{i}^{r}=1+e_{i}\cos(\theta_{i}-\beta_{i}).$
(134)

##  Part II Perturbation Theory for the Three Body Problem

###  VI Delaunay Variables - an interpretation of the Mean Anomaly

We will not here reproduce the whole definition and derivation of action-angle
variables for the two-body problem, nor the canonical transformation to what
are known as the Delaunay variables. Suffice it now to say, that the Delaunay
variables are a set of action-angle variables for the two-body problem, and
that using these variables the Hamiltonian depends on none of the conjugate
coordinates, revealing that the Hamiltonian is integrable, and the dynamics
incredibly simple. Excellent references for the material can be found in [13,
14].

The two-body problem, of masses $m$ and $M$ orbiting each other subject to an
attractive inverse-square-of-distance force $\vec{F}=-k/r^{3}\vec{r}$, where
$\vec{r}$ is the relative position of one body to the other, reduces to a
“1-body” problem of a mass $\mu=mM/(m+M)$, called the reduced mass, moving
about a fixed centre by the same inverse-square force, and the displacement of
the reduced mass from the centre is equal to the relative separation of the
original bodies. This problem has the equation of motion
$\mu\ddot{\vec{r}}=-k/\absolutevalue{\vec{r}}^{3}\vec{r}$, and this may be
given by the Lagrangian
$\mathcal{L}=\tfrac{1}{2}\,\mu|\dot{\vec{r}}|^{2}+k/\absolutevalue{\vec{r}}$;
equivalently the Hamiltonian
$\mathcal{H}=|\vec{p}|^{2}/(2\mu)-k/\absolutevalue{\vec{r}}$, where the linear
momentum is $\vec{p}=\nabla_{\dot{\vec{r}}}\mathcal{L}=\mu\dot{\vec{r}}$. The
only parameter essential to the problem is the ratio $\alpha=k/\mu$, and in
the case of the gravitational two-body problem this is $GmM/\mu=G(M+m)$,
called the standard gravitation parameter (sometimes ‘standard’ is dropped).

Many things are well known of this problem and it’s solutions. It is straight
forward to confirm that the angular momentum
$\vec{L}=\mu\,\vec{r}\times\vec{v}$ is conserved, and that consequently the
motion is planar. It is also well known that the orbits $\\{\vec{r}(t)|\,t\in
T\\}$ are conic sections, including ellipses for negative energies
$E=\mu\,v^{2}/2-k/r$ (with the fixed ‘centre’ at one of the foci of the
ellipse), which will be our focus. The physical size of elliptical orbits is
characterized by the semi-major axis $a$, and it is also well known that the
period of elliptical trajectories is given vy Kepler’s Third Law: the square
of an orbital period is proportional to the cube of the length of the semi-
major axis

$T^{-2}a^{3}=\alpha/4\pi^{2}$ (135)

or equivalently, $a^{3}\omega^{2}=\alpha$, where
$\omega=2\pi/T=\sqrt{\alpha/a^{3}}$, the angular frequency associated with
period $T$, is called the mean motion.

Before describing the Delaunay variables, we introduce various ‘specific’
quantities—quantities that are made ‘massless’, dividing by the chosen mass
scale. Thus the specific energy $\xi=v^{2}/2-\alpha/r$, and specific angular
momentum $\vec{\ell}=\vec{r}\times\vec{v}$, $\ell=|\vec{\ell}|$. These are
both constants of the motion, given by $\xi=-\alpha/(2a)$ and
$\ell=\sqrt{\alpha a(1-e^{2})}$.

Figure 3: The construction of the eccentric and mean anomalies $E,M$.
$E=\angle ACQ$, while $M$ is defined such the area $\absolutevalue{AOQ}$
bounded by the circle is equal to the sector $\absolutevalue{ACX}$.

With respect to a fixed orthonormal frame $\\{\hat{x},\hat{y},\hat{z}\\}$, the
Delaunay momenta are: the $z$-component of angular momentum $\ell^{z}$, the
total angular momentum $\ell$, and finally a third action $j$, related to both
the energy and semi-major axis, given as $j=\sqrt{\alpha
a}=\alpha/\sqrt{-2\xi}$. For lack of another name, I will refer to this as the
impulse of the orbit. The coordinates are as follows333It should be noted that
the Delaunay variables completely describe the 3D orientation of the
elliptical orbit, as well as the progression of motion within the orbit. The
coordinates $\Phi$ and $\eta$ describe the axis of rotation of the orbital
plane relative to the reference $xy$-plane and the position of periapsis
within the orbital plane. From the momenta $\ell^{z},\ell,j$ can be determined
$a,e$ and $\alpha$—in other words, the size, shape and speed of motion within
the orbit. The actual rotation of the orbital plane seems to be missing; that
is, the inclination $\zeta$. But since the angular momentum is perpendicular
to the orbital plane, we can determine the inclination by
$\cos(\zeta)=\ell^{z}/\ell$.. Conjugate to $\ell^{z}$ is the longitude of
ascending node, indicated with $\Phi$. This is the angle from the positive
$x$-axis to the ray along which the plane of the orbit intersects the
$xy$-plane and where $\dot{z}>0$—this direction is called the ascending node.
Conjugate to $\ell$ is the argument of periapsis, $\eta$: the angle from the
ascending node to the ray connecting the central body to the position of
closest approach along the orbit, or periapsis, measured in the orbital plane.
The final angle, conjugate to the impulse $j$, is the so-called mean anomaly
$M$. Unlike the previous two coordinates, $M$ is not so straightforward to
define geometrically. We must first construct the eccentric anomaly, denoted
$E$, by considering a circle of radius $a$ which circumscribes the ellipse,
tangent at periapsis and apoapsis ($A$ and $B$ respectively). The construction
is shown in Fig. 3. The centre $C$ of this circle coincides with the centre of
the ellipse. We describe the position $P$ of the orbiting body around the
ellipse by the angle subtended at the focus $O$ to periapsis, called the true
anomaly and denoted $\nu$

$\displaystyle\vec{OP}=\vec{r}=r(\nu)$ $\displaystyle\hat{r}(\nu)\equiv
r(\nu)\left(\cos\nu,\sin\nu\right)$ (136) $\displaystyle r(\nu)$
$\displaystyle=\frac{a(1-e^{2})}{1+e\cos\nu}.$ (137)

From $P$ we project to the auxiliary circle perpendicular to the major axis,
arriving at a point $Q$. The eccentric anomaly is the angle about $C$ from
periapsis to $Q$, $E=\angle ACQ$. Now, the ellipse can be recovered from the
circle by a contraction of factor $(1-e^{2})^{1/2}$ parallel to the minor
axis. This reveals that, working in the orbital plane and relative to the
centre $C$, the position $P$ is

$\vec{CP}=\left(a\cos E,b\sin E\right),$ (138)

the axes being aligned with the major and minor axes, and $b=a\sqrt{1-e^{2}}$
is the semi-minor-axis.

Since the separation of the centre and focus is $CO=ae$, we have the relations

$\displaystyle\cos E$
$\displaystyle=\frac{(1-e^{2})\cos\nu}{1+e\cos\nu}+e=\frac{\cos\nu-e}{1+e\cos\nu}$
(139) $\displaystyle\sin E$
$\displaystyle=\frac{\sqrt{1-e^{2}}\sin\nu}{1+e\cos\nu}$ (140)
$\displaystyle\tan\tfrac{E}{2}$
$\displaystyle=\sqrt{\frac{1-e}{1+e}}\tan\tfrac{\nu}{2}$ (141)

and the inverse relations

$\displaystyle\cos\nu$ $\displaystyle=\frac{\cos E-e}{1-e\cos E}$ (142)
$\displaystyle\sin\nu$ $\displaystyle=\frac{\sqrt{1-e^{2}}\sin E}{1-e\cos E}.$
(143)

These relations give the separation of the orbit as

$OP=r(E)=a(1-e\cos E).$ (144)

We move to constructing the mean anomaly by employing Kepler’s Second Law: we
know that equal area is swept out by the orbit in equal times. This can be
expressed as the fact that the time-rate-of-change (t.r.o.c) of the area
bounded by the ellipse $\mathcal{A}=\absolutevalue{AOP}$ is constant; in
particular

$\frac{d}{dt}(2\mathcal{A})=\ell.$ (145)

Owing to the aforementioned contraction, the area bounded by the auxiliary
circle $\tilde{\mathcal{A}}=\absolutevalue{AOQ}$ also grows uniformly with
time

$\frac{d}{dt}(2\tilde{\mathcal{A}})=(1-e^{2})^{-1/2}\frac{d}{dt}(2\mathcal{A})=\sqrt{\alpha
a}=j.$ (146)

Thus, we may see the significance of $j$ as (twice) the t.r.o.c. of orbit-area
projected to the auxiliary circle, just as angular momentum $\ell$ is twice
the t.r.o.c. of area swept in the orbit itself.

The mean anomaly $M$ is defined as an angle in the auxiliary circle, measured
at the center, say to a point $X$ on the circle, $M=\angle ACX$, such that the
resulting sector has the same area as
$\tilde{\mathcal{A}}=\absolutevalue{AOQ}$. At a constant radius $a$, we see
that the uniform growth of area $\tilde{\mathcal{A}}$ implies a constant rate-
of-change for $M$

$\displaystyle 2\absolutevalue{ACX}=a^{2}M\text{, so }a^{2}\dot{M}$
$\displaystyle=\frac{d}{dt}(2\tilde{\mathcal{A}})=\sqrt{\alpha a}$ (147)
$\displaystyle\implies\quad\dot{M}$
$\displaystyle=\sqrt{\frac{\alpha}{a^{3}}}=\omega.$ (148)

The mean and eccentric anomalies can be link by a straightforward calculation

$\displaystyle a^{2}M=2\absolutevalue{ACX}=2\absolutevalue{AOQ}$
$\displaystyle=2(\absolutevalue{ACQ}-\absolutevalue{OCQ})$ (149)
$\displaystyle=a^{2}E-2(ae)(a\sin E)/2$ (150) $\displaystyle\therefore\quad M$
$\displaystyle=E-e\sin E.$ (151)

We thus have the relation

$E-e\sin E=\omega(t-\tau)$ (152)

where $\tau$ is a time when the orbit is at periapsis. A tantalizing equation,
but unfortunately this cannot be inverted for the eccentric anomaly as
elementary functions of time. If this could be done, we could write the true
anomaly explicitly as a function of the mean anomaly, and thus of time. As it
is, the relationship is implicit, although of course $\nu(M;e)$ is formally a
well-defined function

$\nu(M;e)=2\arctan\left(\sqrt{\frac{1+e}{1-e}}\tan\frac{E(M;e)}{2}\right);\quad
0\leq e<1$ (153)

where $E(M;e)$ is the inverse444It can be seen that this inverse is a well-
defined function as follows: for $\absolutevalue{e}\leq 1$ the function
$M(E;e)=E-e\sin E$ is increasing except for isolated points. Where $M(E;e)$ is
increasing, the inverse $E(M;e)$ is a differentiable function. Even for $e=1$,
on any domain $dM/de\geq 0$, with equality only at
$E=2k\pi,\,\,k\in\mathbb{Z}$. This gives $E(M;e)$ a continuous and increasing
function, with $dE/dM\rightarrow+\infty$ for $M\in 2\pi\mathbb{Z}$. of (151)
for given eccentricity $0\leq e\leq 1$. We can see that (151) pairs $E=M=k\pi$
for $k\in\mathbb{Z}$. If we consider $E(M;e)=M+\vartheta(M;e)$, this is
$\vartheta=0$ for $M=k\pi$. Furthermore, $\vartheta(M;e)$ is $2\pi$-periodic
in $M$, and an odd function. We can thus express $\vartheta$ as a
$\sin$-Fourier series, with coefficients that are functions of $e$.

$E(M;e)=M+\vartheta(M;e)=M+\sum_{k\geq 1}c_{k}^{E}(e)\sin(kM)$ (154)

These Fourier coefficients are $c_{k}^{E}(e)=2J_{k}(ke)/k$, where $J_{k}(z)$
are the Bessel functions of the first kind. We can see that the relationship
$\nu(E)$ is just the same kind of relationship

$\nu(E;e)=E+\sum_{k\geq 1}c_{k}^{\nu}(e)\sin(kE).$ (155)

Finally, it is straightforward to confirm that composition of such functions
is closed—these functions being the sum of the identity function and some odd
$2\pi$-periodic function. We conclude that the expression for the true anomaly
in terms of the mean anomaly has the same form

$\nu(M;e)=M+\sum_{k\geq 1}C_{k}(e)\sin(kM)$ (156)

this Fourier sine-series is known as the “equation of centre”. The
coefficients have a remarkable expression in terms of Bessel functions

$C_{k}(e)=\frac{2}{k}\left\\{J_{k}(ke)+\sum_{m\geq
1}\beta^{m}\left[J_{k-m}(ke)+J_{k+m}(ke)\right]\right\\},\quad\beta=\frac{1-\sqrt{1-e^{2}}}{e}\sim
e/2+\cdots$ (157)

These coefficients are of order $C_{k}(e)\sim\mathcal{O}(e^{k})$; the Taylor
series for the first of these functions begin

$\displaystyle C_{1}(e)$
$\displaystyle=2e-\frac{1}{4}e^{3}+\frac{5}{96}e^{5}+\frac{107}{4608}e^{7}+\cdots$
(158) $\displaystyle C_{2}(e)$
$\displaystyle=\frac{5}{4}e^{2}-\frac{11}{24}e^{4}+\frac{17}{192}e^{6}+\cdots$
(159) $\displaystyle C_{3}(e)$
$\displaystyle=\frac{13}{12}e^{3}-\frac{43}{64}e^{5}+\frac{95}{512}e^{7}+\cdots$
(160) $\displaystyle C_{4}(e)$
$\displaystyle=\frac{103}{96}e^{4}-\frac{451}{480}e^{6}+\cdots$ (161)
$\displaystyle C_{5}(e)$
$\displaystyle=\frac{1097}{960}e^{5}-\frac{5957}{4608}e^{7}+\cdots$ (162)
$\displaystyle C_{6}(e)$ $\displaystyle=\frac{1223}{960}e^{6}+\cdots$ (163)
$\displaystyle C_{7}(e)$ $\displaystyle=\frac{47273}{32256}e^{7}+\cdots$ (164)

####  VI.1 Auxiliary Circular Orbit

We now discuss two different proposals for what could be considered the
“auxiliary circular orbit”, auxiliary to a given elliptical orbit. The
pedagogical idea is to be able to compare/connect the motion of a body in a
elliptic Keplerian orbit, to a corresponding position in a related circular
orbit. This has the flavour of a “method of images”: the (uniformly orbiting)
position in the circular orbit is an image of the position in the elliptic
orbit, the motion of which is non-uniform, and perhaps more difficult to
develop an intuition for. The relationship between the image-position and the
elliptical position will ultimately be the geometrical one, based on the
definition of the mean anomaly. Indeed, we will take the mean anomaly, as
determined in the elliptical orbit, to be the angular position of a body in a
circular orbit about the same centre (the focus of the ellipse).

Before proceeding, it should be pointed out that we will think primarily in
the orbital plane, or equivalently in two dimensions. In two dimensions $\Phi$
and $\eta$ themselves are undefined (better to say, undefinable), but their
sum $\beta=\Phi+\eta$ is the well defined longitude of periapsis. Delaunay
variables in two dimensions are the conjugate pairs $\ell,\beta$ and $j,M$.

Certainly, the obvious option that first presents itself is to take the orbit
with eccentricity $e^{\prime}=0$ and radius $a^{\prime}=a$. With
$j=\sqrt{\alpha a},\,\ell=\sqrt{\alpha a(1-e^{2})}$,
$\omega=\sqrt{\alpha/a^{3}}$ and $\xi=-\alpha/2a$, this choice gives
$\ell^{\prime}=j^{\prime}=j$, $\omega^{\prime}=\omega$ and $\xi^{\prime}=\xi$.
It should be noted than in this scenario we are also taking the same
gravitational parameter $\alpha^{\prime}=\alpha$, ie. the same total mass in
the ‘phantom’ circular two body set-up as in the original. This certainly has
simplicity to it’s advantage, and that the phantom orbit has the same period
as the original.

However, since the definition of $M$ is so linked to Kepler’s Second Law and
the uniform t.r.o.c. of area-sweep, it would seem desirable to consider an
auxiliary circular orbit which has the same angular momentum
$\ell^{\prime}=\ell$ as the original, as well as having the same period,
requiring $\omega^{\prime}=\omega$. With $e^{\prime}=0$, these conditions
require not only taking a circular radius different than the semi-major axis
of the elliptical orbit, but also considering a change to the gravitational
parameter $\alpha^{\prime}\neq\alpha$. These conditions are the following

$\displaystyle\sqrt{\alpha^{\prime}\,a^{\prime}}$ $\displaystyle=\sqrt{\alpha
a(1-e^{2})},\quad\sqrt{\alpha^{\prime}/{a^{\prime}}^{3}}=\sqrt{\alpha/a^{3}}$
(165) the solution to which is $\displaystyle a^{\prime}$
$\displaystyle=a(1-e^{2})^{1/4},\quad\alpha^{\prime}=\alpha(1-e^{2})^{3/4}.$
(166)

With equal t.r.o.c. of area as well as equal orbital duration, it follows that
the orbits enclose the same area: the elliptical area $\pi ab=\pi
a^{2}\sqrt{1-e^{2}}=\pi{a^{\prime}}^{2}$ is equal to the area of circle with
radius $a^{\prime}$. We thus see the sense of this auxiliary circular orbit,
as the Keplerian orbit with the same period and bounding the same area. In
this sense we really can say the orbits are of the same size.

That we have to consider a modified gravitational parameter is to say that we
must consider the auxiliary circular orbit as having a different total mass.
But this should not be problematic to us. We do well to remember that the
frequency of a Keplerian orbit depends principally on the semi-major axis and
the total mass of the bodies in the orbit (we will consider Newton’s constant
$G$ as universal). So to take a circular orbit with a radius different from a
given semi-major axis, if we want an orbit with the same frequency, we must
take a different total mass. All this is to suggest that we should consider
two orbital scenarios—in which one is both double the size in linear dimension
and has 8 times the total mass than the other—as more similar to each other
than orbits that have merely equal total mass or equal semi-major axis.

This consideration would also seem to elevate the stance of the angular
momentum $\ell$ as somehow principal over the impulse $j$: we have
$j^{\prime}=\ell^{\prime}=\ell\neq j$, $\xi^{\prime}=\xi(1-e^{2})^{1/2}$.
Indeed, we are seeing that the angular momentum is more characteristic of the
orbit, as it corresponds to the actual area-sweeping rate in the physical
orbit, and we preserve that in our ‘phantom’ circular orbit, whereas $j$ is
the projected area rate in the $a$-radius circle that we construct around the
ellipse. It is seen that this circle has far less to do with the actual
physics than the above proposed auxiliary circular orbit.

###  VII Alternate Mass parametrization

We have freedom to take a different parametrization of the masses and the
coefficients $\mu_{1,2,3}$. First note, from the mass distribution of the
three masses $m_{J},m_{S},M$, that any dimensionless ratios of masses depends
on two independent mass-ratios. We can form the following

$\displaystyle\beta_{1}=m_{J}$ $\displaystyle/M_{\Sigma}=9.54\times
10^{-4},\quad\beta_{2}=m_{S}/M_{\Sigma}=2.85\times 10^{-4}$
$\displaystyle\beta_{3}=M/M_{\Sigma}=0.99876=1-1.2384\times 10^{-3}$
$\displaystyle\text{which are constrained by
}\beta_{1}+\beta_{2}+\beta_{3}=1.$

Our freedom is in the mass scale $\tilde{m}$ we take to divide the Lagrangian
(6), giving the coefficients

$\displaystyle\mu_{1}=\frac{Mm_{J}}{M_{\Sigma}\tilde{m}},\quad\mu_{2}=\frac{Mm_{S}}{M_{\Sigma}\tilde{m}},\quad\mu_{3}=\frac{m_{J}m_{S}}{M_{\Sigma}\tilde{m}}$
(167)

If we take $\tilde{m}=\lambda m_{J}$ for some $\lambda>0$ (originally we had
$\lambda=M/M_{\Sigma}=\beta_{3}$), then these are

$\displaystyle\mu_{1}=\frac{\beta_{3}}{\lambda},\quad\mu_{2}=\frac{\beta_{2}\beta_{3}}{\beta_{1}\lambda},\quad\mu_{3}=\frac{\beta_{2}}{\lambda},\quad\delta=\frac{\beta_{2}\beta_{3}}{\lambda}$
(168)

We can choose $\lambda$ such that
$\delta=1-\beta_{3}=\beta_{1}+\beta_{2}=(m_{J}+m_{S})/M_{\Sigma}$, in which
case $\mu_{3}=\delta/\beta_{3}=\delta/(1-\delta)$. This is

$\displaystyle\lambda=\beta_{2}\beta_{3}\frac{M_{\Sigma}}{m_{J}+m_{S}}=\beta_{3}\frac{m_{S}}{m_{J}+m_{S}}$
(169)

and this gives

$\mu_{1}=\frac{m_{J}+m_{S}}{m_{S}},\quad\mu_{2}=\frac{m_{J}+m_{S}}{m_{J}}$
(170)

which are related by $1/\mu_{1}+1/\mu_{2}=1$. In terms of the original
parameter $\mu=m_{S}/m_{J}\sim 0.3$, these are $\mu_{1}=1+\mu^{-1}$ and
$\mu_{2}=1+\mu$. For the masses of Jupiter and Saturn, these are
$\mu_{1}=4.340=4\tfrac{1}{3}+6.452\times 10^{-3}$ and
$\mu_{2}=1.2994=1.3-5.796\times 10^{-4}$. This new mass scale is
$\tilde{m}=\frac{m_{S}}{m_{J}+m_{S}}\frac{M}{M_{\Sigma}}m_{J}=\frac{Mm_{J}m_{S}}{M_{\Sigma}(m_{J}+m_{S})}$.
We will redefine $\epsilon\equiv\mu_{3}=\delta/(1-\delta)=(m_{J}+m_{S})/M$.
Notice that $\delta$ is always $0\leq\delta\leq 1$. For the planetary regime,
certainly we would say $m_{J}+m_{s}\leq M$, so that $\delta\leq 1/2$ and
$\epsilon\leq 1$.

It is instructive to return to the coefficient matrix of the equations (17)

$\displaystyle A$
$\displaystyle=G\left(\begin{array}[]{ccc}-(M+m_{J})&m_{S}&m_{S}\\\
m_{J}&-(M+m_{S})&m_{J}\\\ M&M&-(m_{J}+m_{S})\end{array}\right)$ (174)
$\displaystyle=\alpha\left(\begin{array}[]{ccc}\beta_{2}-1&\beta_{2}&\beta_{2}\\\
\beta_{1}&\beta_{1}-1&\beta_{1}\\\
\beta_{3}&\beta_{3}&\beta_{3}-1\end{array}\right)$ (178)
$\displaystyle=\alpha\left(\begin{array}[]{ccc}-1+\delta/\mu_{1}&\delta/\mu_{1}&\delta/\mu_{1}\\\
\delta/\mu_{2}&-1+\delta/\mu_{2}&\delta/\mu_{2}\\\
1-\delta&1-\delta&-\delta\end{array}\right)$ (182)

###  VIII Hamilton’s Equations

We may now proceed to consider the asymptotic analysis of the three body
problem, in the planetary case: when two bodies (the planets) orbit a third,
and their masses are also much smaller than the third, but not so much that
the gravitational attraction between them is negligible. Using the Delaunay
variables defined in the three sectors, we have the Hamiltonian

$\displaystyle\mathbf{H}_{\lambda}(J_{i},h_{i},h^{z}_{i},M_{i},\eta_{i},\Phi_{i},\vec{\lambda})$
$\displaystyle=\sum_{i=1}^{3}\left\\{\,-\mu_{i}^{3}\frac{\alpha^{2}}{2J_{i}^{2}}-\lambda\cdot\vec{r}_{i}\right\\}$
(183)
$\displaystyle=\sum_{i=1}^{3}\left\\{\,-\mu_{i}^{3}\frac{\alpha^{2}}{J_{i}^{2}}\right\\}-\vec{\lambda}\cdot\left(\vec{r}_{1}+\vec{r}_{2}+\vec{r}_{3}\right)$

where the Delaunay momenta are the relative variables (as opposed to specific)

$J_{i}=\mu_{i}j_{i}\,,\quad h_{i}=\mu_{i}\ell_{i}\,,\quad
h^{z}_{i}=\mu_{i}\ell^{z}_{i}\,.$ (184)

We notice that, with the exception of the impulses, the Delaunay variables
enter into this Hamiltonian via the geometry, in the vector positions
$\vec{r}_{i}$. The vector position in the $i^{\text{th}}$ sector is

$\vec{r}_{i}=\frac{J_{i}^{2}}{\alpha\mu_{i}^{2}}\left(1-e_{i}\cos{E_{i}}\right)\hat{r}_{i}$
(185)

The eccentricity is given by $h_{i}^{2}/J_{i}^{2}=1-e_{i}^{2}$, and the
direction vector $\hat{r}_{i}$ is derived by the following sequence of
rotations: we start with the unit vector in a fixed reference, say
$\hat{x}$-direction. This is then rotated ccw about the $z$-axis by both the
true anomaly (related to the eccentric) and the argument of periapsis,
$\nu+\eta$. At this stage we can imagine we have an elliptic orbit with
periapsis at argument $\eta$ and ascending node at the positive $x$-axis, but
no inclination. We need to rotate ccw about the $x$-axis by the inclination
$\zeta_{i}$. Note that the rotation by $\zeta_{i}$ is given by
$\cos\zeta_{i}=h^{z}_{i}/h_{i}$ and
$\sin\zeta_{i}=\sqrt{1-{h^{z}_{i}}^{2}/h_{i}^{2}}$. Finally, we rotate again
about the $z$-axis, bringing the ascending node to longitude $\Phi_{i}$. Thus
$\hat{r}_{i}$ is

$\hat{r}_{i}=R^{z}_{\Phi_{i}}R^{x}_{\zeta_{i}}R^{z}_{\nu_{i}+\eta_{i}}\left[1,0,0\right]^{\text{T}}\,.$
(186)

Now, the components as a result of the rotation by true anomaly, are given in
terms of the eccentric by

$R^{z}_{\nu_{i}}\left[1,0,0\right]^{\text{T}}=\left[\frac{\cos
E_{i}+e_{i}}{1-e_{i}\cos E_{i}},\frac{\sqrt{1-e_{i}^{2}}\sin
E_{i}}{1-e_{i}\cos E_{i}},0\right]^{\text{T}}\,.$ (187)

So we can see the $i^{\text{th}}$ position vector is

$\vec{r}_{i}=\frac{J_{i}^{2}}{\alpha\mu_{i}^{2}}\,R^{z}_{\Phi_{i}}R^{x}_{\zeta_{i}}R^{z}_{\eta_{i}}\left[\cos
E_{i}+e_{i},\sqrt{1-e_{i}^{2}}\sin E_{i},0\right]$ (188)

The components of the remaining rotation matrix are

$R^{z}_{\Phi}R^{x}_{\zeta}R^{z}_{\eta}=\left(\begin{array}[]{ccc}\cos\Phi\cos\eta-\sin\Phi\cos\zeta\sin\eta&-\cos\Phi\sin\eta-\sin\Phi\cos\zeta\cos\eta&\sin\Phi\sin\zeta\\\
\sin\Phi\cos\eta+\cos\Phi\cos\zeta\sin\eta&\cos\Phi\cos\zeta\cos\eta-\sin\Phi\sin\eta&-\cos\Phi\sin\zeta\\\
\sin\zeta\sin\eta&\sin\zeta\cos\eta&\cos\zeta\end{array}\right)$ (189)

From the Hamiltonian (183), Hamilton’s equation are

$\displaystyle\dot{M}_{i}$
$\displaystyle=\frac{\partial\mathbf{H}_{\lambda}}{\partial
J_{i}}=\mu_{i}^{3}\alpha^{2}J_{i}^{-3}+\frac{\partial\mathbf{R}}{\partial
J_{i}}$ $\displaystyle\dot{J}_{i}$
$\displaystyle=\mu_{i}\dot{j}_{i}=-\frac{\partial\mathbf{R}}{\partial M_{i}}$
(190) $\displaystyle\dot{\eta}_{i}$
$\displaystyle=\frac{\partial\mathbf{R}}{\partial h_{i}}$
$\displaystyle\dot{h}_{i}$
$\displaystyle=\mu_{i}\dot{h}_{i}=-\frac{\partial\mathbf{R}}{\partial\eta_{i}}$
(191) $\displaystyle\dot{\Phi}_{i}$
$\displaystyle=\frac{\partial\mathbf{R}}{\partial h^{z}_{i}}$
$\displaystyle\dot{h}^{z}_{i}$
$\displaystyle=\mu_{i}\dot{h}^{z}_{i}=-\frac{\partial\mathbf{R}}{\partial\Phi_{i}}$
(192) $\displaystyle 0$
$\displaystyle=\nabla_{\vec{\lambda}}\mathbf{H}_{\lambda}=\vec{S}$ (193)

for $i=1,2,3$, where $\mathbf{R}=-\vec{\lambda}\cdot\vec{S}$ is the function
by which the Hamiltonian is perturbed and
$\vec{S}=\vec{r}_{1}+\vec{r}_{2}+\vec{r}_{3}$ is the constraint. We know that
when the equations are satisfied subject to the constraint, $\vec{\lambda}$
takes the values (15) as determined in the first section, and thus also
determined by the geometry (188,189). Now, we must be careful with the partial
derivatives of the perturbation. At first blush, we should regard both
$\vec{\lambda}$ and $\vec{S}$ as functions of the coordinates. If $\chi$ is
one of the canonical variables, then the derivative of the perturbation with
respect to $\chi$ is

$-\frac{\partial\mathbf{R}}{\partial\chi_{i}}=\frac{\partial\vec{\lambda}}{\partial\chi}\cdot\vec{S}+\vec{\lambda}\cdot\frac{\partial\vec{S}}{\partial\chi}$
(194)

In particular, we compute the potential $\partial_{\chi}\vec{S}$ without
regard to the constraint, but by the function form of the constraint
$\vec{S}=\vec{r}_{1}+\vec{r}_{2}+\vec{r}_{3}$. If $\chi=\chi_{i}$ is a
variable of the $i^{\text{th}}$ sector, then this is
$\partial_{\chi_{i}}\vec{r}_{i}$. We might similarly determine the derivative
$\partial_{\chi_{i}}\vec{\lambda}$ by combining (15,188,189), but at this
point we may evaluate (194) subject to constraint $\vec{S}=0$, so that the
first term vanishes. Thus

$\frac{\partial\mathbf{R}}{\partial\chi_{i}}=-\vec{\lambda}\cdot\frac{\partial\vec{r}_{i}}{\partial\chi_{i}}.$
(195)

For some of the derivatives we have the following (suppressing indices $i$),
using $c_{\zeta}=\cos\zeta$ and $s_{\zeta}=\sin\zeta$,

$\displaystyle M$ $\displaystyle=(1-e\cos E)dE$
$\displaystyle\implies\partial_{M}$ $\displaystyle=(1-e\cos
E)^{-1}\partial_{E}\qquad\qquad\qquad\qquad\qquad\,\,\,\,$ (196)
$\displaystyle\partial_{J}$ $\displaystyle=$ $\displaystyle\partial_{J}|_{e}$
$\displaystyle\,\,+$ $\displaystyle\frac{\partial e}{\partial J}\,$
$\displaystyle\partial_{e}$ $\displaystyle=$ $\displaystyle\partial_{J}|_{e}$
$\displaystyle\,\,+$ $\displaystyle\,\,\frac{1-e^{2}}{eJ}\,$
$\displaystyle\partial_{e}$ (197) $\displaystyle\partial_{h}$ $\displaystyle=$
$\displaystyle\,\frac{\partial e}{\partial h}$ $\displaystyle\,\partial_{e}$
$\displaystyle+$ $\displaystyle\frac{\partial c_{\zeta}}{\partial h}\,$
$\displaystyle\partial_{c_{\zeta}}$ $\displaystyle\,\,+\,\,$
$\displaystyle\frac{\partial s_{\zeta}}{\partial h}\,$
$\displaystyle\partial_{s_{\zeta}}$ $\displaystyle=$
$\displaystyle\,-\frac{1-e^{2}}{eh}$ $\displaystyle\,\partial_{e}$
$\displaystyle-$ $\displaystyle\frac{c_{\zeta}}{h}\,$
$\displaystyle\partial_{c_{\zeta}}$ $\displaystyle\,\,+\,\,$
$\displaystyle\frac{c_{\zeta}^{2}}{s_{\zeta}h}\,$
$\displaystyle\partial_{s_{\zeta}}$ (198) $\displaystyle\partial_{h^{z}}$
$\displaystyle=$ $\displaystyle\frac{\partial c_{\zeta}}{\partial h^{z}}$
$\displaystyle\,\partial_{c_{\zeta}}$ $\displaystyle+$
$\displaystyle\frac{\partial s_{\zeta}}{\partial h^{z}}\,$
$\displaystyle\partial_{s_{\zeta}}$ $\displaystyle=$
$\displaystyle\frac{c_{\zeta}}{h^{z}}$ $\displaystyle\,\partial_{c_{\zeta}}$
$\displaystyle-$ $\displaystyle\frac{c_{\zeta}^{2}}{s_{\zeta}h^{z}}\,$
$\displaystyle\partial_{s_{\zeta}}.$ (199)

These will be instrumental to formulating the KAM theory for this
perturbational problem.

###  IX Setup for KAM Theory

KAM Theory is a perturbative approach for nearly integrable Hamiltonians,
which are a perturbation away from depending on only the momenta
$\bar{J}\in\mathbb{R}^{n}$

$\mathbf{H}(\bar{J},\bar{\theta};\delta)=\mathbf{H}_{0}(\bar{J})+\,\mathbf{\widehat{H}}(\bar{J},\bar{\theta};\delta)\sim\mathbf{H}_{0}(\bar{J})+\sum_{k\geq
1}\delta^{k}\mathbf{H}_{k}(\bar{J},\bar{\theta}).$ (200)

One seeks a nearly-identical canonical transformation
$\mathcal{C}:(\bar{J},\bar{\theta})\mapsto(\tilde{J},\tilde{\theta})$, given
by a generating function

$\displaystyle\mathbf{\Psi}(\tilde{J},\bar{\theta};\delta)$
$\displaystyle=\tilde{J}\cdot\bar{\theta}\,+\,\mathbf{\widehat{\Psi}}(\tilde{J},\bar{\theta};\delta)$
$\displaystyle\sim\tilde{J}\cdot\bar{\theta}\,+\,\sum_{k\geq
1}\delta^{k}\mathbf{\Psi}_{k}(\tilde{J},\bar{\theta}),$ (201)

which gives

$\displaystyle\bar{J}=\nabla_{\bar{\theta}}\mathbf{\Psi}$
$\displaystyle=\tilde{J}+\nabla_{\bar{\theta}}\mathbf{\widehat{\Psi}}(\tilde{J},\bar{\theta};\delta)$
(202) $\displaystyle\sim\tilde{J}\,+\,\sum_{k\geq
1}\delta^{k}\,\nabla_{\bar{\theta}}\mathbf{\Psi}_{k}$ (203)
$\displaystyle\tilde{\theta}=\nabla_{\tilde{J}}\mathbf{\Psi}$
$\displaystyle=\bar{\theta}+\nabla_{\tilde{J}}\mathbf{\widehat{\Psi}}(\tilde{J},\bar{\theta};\delta)$
(204) $\displaystyle\sim\bar{\theta}\,+\,\sum_{k\geq
1}\delta^{k}\,\nabla_{\tilde{J}}\mathbf{\Psi}_{k}.$ (205)

The goal of this transformation is that the Hamiltonian, in terms of the new
coordinates, depends only on the new momenta

$\mathbf{H}\circ\mathcal{C}^{-1}\,(\tilde{J},\tilde{\theta};\delta)=\tilde{\mathbf{H}}(\tilde{J};\delta)$
(206)

Practically speaking, if we only do this to so many terms, say truncating
$\mathbf{\widehat{H}}=\sum_{k=1}^{N}\delta^{k}\mathbf{H}_{k}$ and
$\mathbf{\widehat{\Psi}}=\sum_{k=1}^{N}\delta^{k}\mathbf{\Psi}_{k}$, then the
transformed Hamiltonian only depends on new coordinates $\tilde{\theta}$ at
the $N+1$ order in $\delta$

$\tilde{\mathbf{H}}(\tilde{J},\tilde{\theta};\delta)=\tilde{\mathbf{H}}_{N}(\tilde{J};\delta)+\mathcal{O}(\delta^{N+1})(\tilde{J},\tilde{\theta};\delta).$
(207)

We will need to decompose the perturbations into the average over angle-
variables $\bar{\theta}$

$\displaystyle\langle\mathbf{\widehat{H}}\rangle(\bar{J};\delta)$
$\displaystyle=(2\pi)^{-n}\oiint\mathbf{\widehat{H}}(\bar{J},\bar{\theta};\delta)\,\differential^{n}\bar{\theta}$
(208) $\displaystyle\langle\mathbf{H}\rangle_{k}(\bar{J};\delta)$
$\displaystyle=(2\pi)^{-n}\oiint\mathbf{H}_{k}(\bar{J},\bar{\theta};\delta)\,\differential^{n}\bar{\theta}$
(209)

and the remainders
$\mathbf{\widehat{F}}=\mathbf{\widehat{H}}-\langle\mathbf{\widehat{H}}\rangle$,
$\mathbf{F}_{k}=\mathbf{H}_{k}-\langle\mathbf{H}\rangle_{k}$.

Starting with just one order in $\delta$, $N=1$, writing
$\bar{\mathbf{\Omega}}_{0}(\bar{J})=\nabla_{\bar{J}}\mathbf{H}_{0}$ for the
frequency functions of the unperturbed Hamiltonian, we find

$\displaystyle\mathbf{H}(\bar{J},\bar{\theta};\delta)$
$\displaystyle\sim\mathbf{H}_{0}(\tilde{J}+\nabla_{\bar{\theta}}\mathbf{\widehat{\Psi}})+\delta\,\mathbf{H}_{1}(\tilde{J}+\nabla_{\bar{\theta}}\mathbf{\widehat{\Psi}},\bar{\theta})+\cdots$
(210)
$\displaystyle\sim\mathbf{H}_{0}(\tilde{J}+\delta\,\nabla_{\bar{\theta}}\mathbf{\Psi}_{1}+\cdots)+\delta\,\langle\mathbf{H}\rangle_{1}(\tilde{J}+\delta\,\nabla_{\bar{\theta}}\mathbf{\Psi}_{1}+\cdots)$
$\displaystyle\hphantom{=}\,\,+\delta\,\mathbf{F}_{1}(\tilde{J}+\delta\nabla_{\bar{\theta}}\mathbf{\Psi}_{1}+\cdots,\bar{\theta})+\cdots$
(211)
$\displaystyle\sim\mathbf{H}_{0}(\tilde{J})+\delta\,\bar{\mathbf{\Omega}}_{0}(\tilde{J})\cdot\nabla_{\bar{\theta}}\mathbf{\Psi}_{1}+\delta\,\langle\mathbf{H}\rangle_{1}(\tilde{J})+\delta\,\mathbf{F}_{1}(\tilde{J},\bar{\theta})+\mathcal{O}(\delta^{2})$
(212)

Thus we have

$\tilde{\mathbf{H}}_{1}(\tilde{J};\delta)=\mathbf{H}_{0}(\tilde{J})+\delta\,\langle\mathbf{H}\rangle_{1}(\tilde{J}).$
(213)

That is, the functional form of the transformed Hamiltonian, as a function of
the new momenta, is the sum of the unperturbed Hamiltonian and the average of
the perturbation over all angle variables, to leading order. This is achieved
by matching the remaining terms at first order

$\bar{\mathbf{\Omega}}_{0}(\tilde{J})\cdot\nabla_{\bar{\theta}}\mathbf{\Psi}_{1}(\tilde{J},\bar{\theta})+\mathbf{F}_{1}(\tilde{J},\bar{\theta})=0$
(214)

Expanding $\mathbf{\Psi}_{1},\,\mathbf{F}_{1}$ in Fourier multi-series in
$\bar{\theta}$, with coefficients
$\Psi_{1}^{\bar{k}}(\tilde{J}),\,F_{1}^{\bar{k}}(\tilde{J})$, (214) becomes
(suppressing dependence on $\tilde{J}$)

$\sum_{\bar{k}\in\mathbb{Z}^{n}\backslash\\{0\\}}\left\\{\left(i\bar{\mathbf{\Omega}}_{0}\cdot\bar{k}\,\Psi_{1}^{\bar{k}}+F_{1}^{\bar{k}}\right)e^{i\bar{k}\cdot\bar{\theta}}\right\\}=0$
(215)

So the solution is

$\mathbf{\Psi}_{1}(\tilde{J},\bar{\theta})=i\sum_{\bar{k}\in\mathcal{S}_{1}}\frac{F_{1}^{\bar{k}}(\tilde{J})}{\bar{k}\cdot\bar{\mathbf{\Omega}}_{0}(\tilde{J})}\,e^{i\bar{k}\cdot\bar{\theta}}$
(216)

where $\mathcal{S}_{1}\subset\mathbb{Z}^{n}$ is the set of multi-indices for
which $\mathbf{F}_{1}$ has non-zero Fourier coefficient. Here we finally see a
problem: that if ever the frequency vector is orthogonal to one of these
integer multi-indices $\bar{k}\cdot\bar{\mathbf{\Omega}}_{0}=0$, then the
solution (216) breaks down. This is the problem of resonance, and it requires
a modification called resonant perturbation theory. We have a special case of
this problem: for our Hamiltonian (183), the unperturbed terms

$-\mu_{1}^{3}\frac{\alpha^{2}}{J_{1}^{2}}-\mu_{2}^{3}\frac{\alpha^{2}}{J_{2}^{2}}$
(217)

do not depend on the third impulse $J_{3}$ nor any of the angular momenta,
which only come into the Hamiltonian through the perturbing function
$\mathbf{R}$. This results in the corresponding components of frequency vector
being identically zero. In other words, our Hamiltonian is degenerate, and the
solution is called degenerate perturbation theory. What we need to do, not
only decomposing the perturbation into it’s average-over-all-angles and
remainder, but further decomposing the remainder into it’s average over the
mean anomalies—the angles whose momenta the $0$th order Hamiltonian depends on
(it will be seen that $M_{3}$ can be included, and the integrable term for the
third sector can be included with the (217))—and remainder from that. If
$\mathbf{R}=-\vec{\lambda}\cdot\vec{S}=\langle\mathbf{R}\rangle+\mathbf{F}$

$\displaystyle\langle\mathbf{R}\rangle(J_{i},h_{i},h^{z}_{i})=(2\pi)^{-9}\oiint\,\mathbf{R}\,\,\differential^{3}M\,\differential^{3}\eta\,\differential^{3}\Phi$
(218)
$\displaystyle\mathbf{F}(J_{i},h_{i},h^{z}_{i},M_{i},\eta_{i},\Phi_{i})=\mathbf{R}-\langle\mathbf{R}\rangle.$
(219)

Then we further decompose $\mathbf{F}=\mathbf{F}_{M}+\tilde{\mathbf{F}}$

$\displaystyle\mathbf{F}_{M}(J_{i},h_{i},h^{z}_{i},\eta_{i},\Phi_{i})=(2\pi)^{-3}\oiint\,\mathbf{F}\,\,\differential^{3}M$
(220)
$\displaystyle\tilde{\mathbf{F}}(J_{i},h_{i},h^{z}_{i},M_{i},\eta_{i},\Phi_{i})=\mathbf{F}_{M}-\mathbf{F}.$
(221)

So $\mathbf{R}=\langle\mathbf{R}\rangle+\mathbf{F}_{M}+\tilde{\mathbf{F}}$.
Then we seek a canonical transformation to a new Hamiltonian that doesn’t
depend on the mean anomalies. This will be elaborated in future work.

###  X Conclusions

In this work, we have presented a self-contained analysis of the gravitational
three-body problem in the planetary scenario. This is done in heliocentric
coordinates with the use of Lagrange multipliers to elegantly handle the form
of the Lagrangian. First working in two dimensions, we develop the equations
of motion as a first order dynamical system in longitudes $\theta_{i}$,
angular momentum $\ell_{i}$ and eccentricity vectors $\vec{e}_{i}$. Auxiliary
vectors $\vec{f}_{i}=\hat{r}_{i}+\vec{e}_{i}$ greatly simplify the equations
algebraically, especially in terms of the components
$f_{i}^{r},f_{i}^{\theta}$.

We then give the definition and construction of the Delaunay action-angle
variables for the two body problem. We present a novel conceptualization of
the mean anomaly of an eccentric orbit as the true anomaly of a ‘phantom’ or
‘image’ body, also orbiting the central mass, but in a circular orbit with
radius reduced from the semi-major axis $a^{\prime}=a(1-e^{2})^{1/4}$ as well
as reduced gravitational parameter $\alpha^{\prime}=\alpha(1-e^{2})^{3/4}$.
These trajectories bound equal areas within the orbital plane, and the orbits
have the same orbital frequency and angular momentum, as opposed to energy. We
then develop the Hamiltonian for the problem, and investigate the geometry of
the orbital positions $\vec{r}_{i}$ in terms of the elements, in order to
express the functional form of the Hamiltonian perturbation
$\mathbf{R}=-\vec{\lambda}\cdot(\vec{r}_{1}+\vec{r}_{2}+\vec{r}_{3})$.
Hamilton’s equations are derived (190-193), and derivatives with respect to
Delaunay variables are written down in terms of derivatives with respect to
other elements $E,e,$ and $c_{\zeta},s_{\zeta}$. Finally, we considered the
Hamiltonian perturbation theory of the problem, and learned that we need to
proceed via a careful approach, seeing that our Hamiltonian is degenerate—not
varying with all of the Delaunay momenta when $\delta=0$.

There is much room for further work. The detailed work of the KAM theory
analysis for this specific problem can now be begun. Of particular interest is
the Great Inequality of Jupiter and Saturn: to identify the term or terms that
correspond to this perturbation that is unexpectedly large in both period and
amplitude. The equations in two dimensions, in terms of the auxiliary vectors,
are interesting in their own right, and might be amenable to a multiple-scales
asymptotic approach, as may be the Hamiltonian equations themselves, KAM
theory aside. Furthermore, the general-relativist corrections to this work are
of great interest. This would presumably be approached via the post-Newtonian
formalism.

### References

  * [1] Newton, Isaac. Philosophiae Naturalis Principia Mathematica. 1687, London.
  * [2] Laplace, Pierre-Simon. Théorie de Jupiter et de Saturne. Memoire de l’Academie des Sciences de Paris 1788, 33-160.
  * [3] Wilson, Curtis. The Great Inequality of Jupiter and Saturn: From Kepler to Laplace. 1985, Springer-Verlag.
  * [4] Hill, G.W. Notes on the Theories of Jupiter and Saturn. The Analyst 1881 8(2), 33-40.
  * [5] Hill, G.W. On the Extension of Delaunay’s Method in the Lunar Theory to the General Problem of Planetary Motion. Trans. Amer. Math. Soc. 1900 1(2), 205-242.
  * [6] Musielak, Z E and Quarles, B. The Three-Body Problem. Rep. Prog. Phys. 2014 77 065901.
  * [7] Arnold, V.I. Proof of a Theorem by A.N. Kolmogorov on the invariance of quasi-periodic motions under small perturbations of the Hamiltonian. Russ. Math. Survey 1963 18, 13-40.
  * [8] Poincaré, Henri. Méthodes Nouvelles de la Mécanique Céleste, vol 1-3. 1892-99, Paris: Gauthier-Villars.
  * [9] Feldman, David. Chaos and Dynamical Systems. 2019, Princeton Univ. Press. ISBN: 9780691161525
  * [10] Brouke, R and Lass, H. A Note on Relative Motion in the General Three-Body Problem. Celestial Mechanics 1973 8(1), 5-10.
  * [11] Hunt, Katie and Strickland, Ashley. “Jupiter and Saturn’s ’great conjunction’ captured in stunning images.” CTVNews.ca, Dec. 22, 2020. https://www.ctvnews.ca/sci-tech/jupiter-and-saturn-s-great-conjunction-captured-in-stunning-images-1.5241665
  * [12] Byrd, Deborah and McClure, Bruce. “All you need to know: 2020’s great conjunction of Jupiter and Saturn.” EarthSky.org, Dec. 21, 2020. https://earthsky.org/astronomy-essentials/great-jupiter-saturn-conjunction-dec-21-2020
  * [13] Celletti, Alessandra. Perturbation Theory in Celestial Mechanics. 2007, obtained from https://web.ma.utexas.edu/mp_arc/c/07/07-303.pdf
  * [14] Morbidelli, Alessandro. Modern Celestial Mechanics. 2011, obtained from https://www-n.oca.eu/morby/celmech.pdf
  * [15] Seidelmann, K.P., ed. Explanatory Supplement to the Astronomical Almanac. 1992 University Science Books, Mill Valley, California.