# KAM theory for the three-body problem

Prof. Jacques Féjoz accepted the invitation on 16 March 2010

**KAM theory for the three-body problem** refers to applications of the theory of Kolmogorov-Arnold-Moser to the three-body problem in celestial mechanics. It aims at proving the existence of quasi-periodic motions for three particles undergoing Newtonian attraction. If periodic orbits are conjecturally dense, it is also likely that their union has Lebesgue measure zero. In contrast, quasi-periodic solutions have a positive measure in various regions of the phase and parameter spaces where the three-body problem is nearly integrable. The main known examples of such regions are:

- The
**planetary problem**, which is the asymptotic regime of the three-body problem when two bodies, the**planets**, have small masses compared to the third body, the**Sun**. It is also assumed that the planets have negative energies and do not come close to each other. In this regime, each planet primarily undergoes the attraction of the Sun, and thus the system is the perturbation of two uncoupled Kepler problems (as it will turn out below, this point is not completely obvious, since the Sun is attracted by both planets, and there is a subtelty if the inner planet gets close to the Sun, thus making the Sun move very fast close to the collision). This problem also makes sense for any number of planets around the Sun. It is important from the astronomical viewpoint, since it is a good model of the Sun-Jupiter-Saturn system, the mass of Jupiter, the heaviest planet of our Solar system, being approximately 1/1000th that of our Sun. V. Arnold gave the first proof of existence of such motions, in the neighborhood of coplanar circular Keplerian ellipses.

- The
**lunar problem**(also sometimes called the**hierarchical problem**, or**Hill's problem**), which is the asymptotic regime of the three-body problem when one of the bodies, the "Sun", is far from the other two, the "Earth" and the "Moon". It is a good model for the real Sun-Earth-Moon system. In the first approximation, the Sun only feels the attraction of one fictitious body located at the center of mass of the Earth and the Moon, whereas the Earth and the Moon primarily undergo their own mutual attraction. See below.

- The neighborhood of relative equilibria (the three-body problem failing to have any equilibrium).

- The neighborhood of periodic orbits (which Poincaré worshiped as the main breach through which the three-body problem can be attacked).

- The neighborhood of other quasiperiodic motions.

In all these special, nearly integrable regimes of the three-body problem, KAM theory asserts that there exist quasiperiodic motions, provided some transversality condition is satisfied. Without this transversality condition, all the quasiperiodic motions of the integrable approximation may vanish according to the perturbation. Below we give a simple example of computation showing that this transversality condition holds in the lunar problem in the neighborhood of circular coplanar Keplerian motions.

## Contents |

## Statement for the lunar problem

The Hamiltonian of the three-body problem is
\[H = \sum_{0 \leq j \leq 2} \frac{\|p\|^2}{2m_j} - \sum_{0 \leq j < k
\leq 2} \frac{m_jm_k}{\|q_j-q_k\|},\] where \(q_j \in \R^3\) is the
position of the \(j\)-th body and \(p_j\in \R^3\) is its impulsion. We
will prove the existence of quasiperiodic motions in the
**hierarchical** regime, where two bodies
(say, \(q_0\) and \(q_1\)) revolve around each other while the third body
revolves, far away, around the center of mass of the two
primaries.

**Theorem** *There exist a set of initial conditions of positive
Lebesgue measure leading to quasiperiodic solutions, arbitrarily close
to Keplerian, coplanar, circular motions, with semi major axis ratio
arbitrarily small.*

The hurried reader may simplify the following discussion by focusing on the plane invariant subproblem.

## Reduction by the symmetry of translations

Let \((Q_0,Q_1,Q_2,P_0,P_1,P_2)\) be the Jacobi coordinates, defined by: \[ \begin{cases} Q_0 = q_0\\ Q_1 = q_1 - q_0\\ Q_2 = q_2 - \sigma_0 q_0 - \sigma_1 q_1, \end{cases} \begin{cases} P_0 = p_0 + p_1 + p_2\\ P_1 = p_1 + \sigma_1 p_2\\ P_2 = p_2, \end{cases} \] where \(1/\sigma_0 = 1 + m_1/m_0\) and \(1/\sigma_1 = 1 + m_0/m_1\). \(P_0\) is the total linear momentum, which can be assumed equal to \(0\) without loss of generality. Besides, \(H\) does not depend on \(Q_0\). So, \((Q_1,Q_2,P_1,P_2)\) is a symplectic coordinate system on the phase space reduced by the symmetry of translation, and the equations read \[ \begin{cases} \dot Q_i = \partial_{P_i} H \quad (i=1,2)\\ \dot P_i =-\partial_{Q_i} H. \end{cases} \] A direct computation shows that \[H = \sum_{1 \leq i \leq 2} \frac{\|P_i\|^2}{2\mu_i} - \sum_{0 \leq i < j \leq 2} \frac{m_im_j}{\|q_i-q_j\|},\] with \[M_0 = m_0, \quad M_1= m_0 + m_1 \quad \mbox{and} \quad \frac{1}{\mu_i} = \frac{1}{M_{i-1}} + \frac{1}{m_i}.\]

## Keplerian approximation

One can split \(H\) into two parts \[H = \mathrm{Kep}\, + \mathrm{Rem}\,\] where \[\mathrm{Kep}\, = \sum_{1 \leq i \leq 2} \left(\frac{\|P_i\|^2}{2\mu_i} - \frac{\mu_iM_i}{\|Q_i\|}\right)\] is a sum of two uncoupled Kepler problems, and \[\mathrm{Rem}\, = \frac{\mu_2M_2}{\|Q_2\|} - \frac{m_0m_2}{\|q_2 - q_0\|} - \frac{m_1m_2}{\|q_2 - q_1\|}\] is the remainder.

Let us assume that the two terms of \(\mathrm{Kep}\,\) are negative so
that each body \(Q_i\) under the flow of \(\mathrm{Kep}\,\) describes
a Keplerian ellipse. Let
\((\ell_i,L_i,g_i,G_i,\theta_i,\Theta_i)_{i=1,2}\) be the associated
Delaunay coordinates. These coordinates are symplectic and analytic
over the open set where motions are non-circular and
non-horizontal (see Les méthodes nouvelles or, as a
modern reference, (Féjoz, 2013)); since we will precisely
be interested in a neighborhood of circular coplanar motions, these
variables are only intermediate coordinates for computations. One
shows that
\[\mathrm{Kep}\, = - \sum_{i=1,2} \frac{\mu_i^3M_i^2}{2L_i^2}.\]
The Keplerian frequencies are
\[\kappa_i = \frac{\partial \mathrm{Kep}\,}{\partial L_i} =
\frac{\mu_i^3M_i^2}{L_i^3} = \frac{\sqrt{M_i}}{a_i^{3/2}},\]
so that the *Keplerian frequency map*
\[\kappa : (L_1,L_2) \mapsto (\kappa_1, \kappa_2)\]
is a diffeomorphim \((\R_+^\times)^2 \to (\R_+^\times)^2\). Due to the
fact that the Keplerian part depends only on \(2\) of the action
variables, solutions of the Keplerian approximation are quasiperiodic
with atmost \(2\) independant frequencies. This degeneracy has been
interpreted as a hidden \(SO(4)\)-symmetry for each planet, whose
momentum map is given partly by the eccentricity vector. How the
Keplerian ellipses slowly rotate and deform will be determined by
mutual attractions. This degeneracy is specific to the Newtonian and
elastic potentials, as Bertrand's theorem asserts.

## Quadrupolar approximation

In the hierarchical regime (\(a_1\ll a_2\)) the dominating term of the remainder is \[ \mathrm{Main}\, := - \mu_1 m_2 P_2(\cos \theta) \frac{\|X_1\|^2}{\|X_1\|^3}, \] with \(P_2(c) = \frac{1}{2}(3c^2-1)\) (second Legendre polynomial) and \(\theta = \widehat{Q_1,Q_2}\).

Since the Keplerian frequencies satisfy \(\kappa_1 \gg \kappa_2\), we
may average out the fast, Keplerian angles \(\ell_1\) and \(\ell_2\)
successively, thus without small
denominators (Féjoz, 2002), (Jefferys, 1966).
The *quadrupolar Hamiltonian* is
\[
\mathrm{Quad}\, = \int_{\mathbb{T}^2} \mathrm{Main}\, \, \frac{d\ell_1
\, d\ell_2}{4\pi^2};
\]
It is the dominating interaction term which rules the slow
deformations of the Keplerian ellipses. It naturally defines a
Hamiltonian on the space of pairs of Keplerian ellipses with fixed
semi major axes, which is called the *secular space*.

After reduction by the symmetry of rotations (e.g. with Jacobi's reduction of the nodes, which consists in fixing the angular momentum vector, say, vertically, and quotienting the so-obtained codimension-3 Poisson submanifold by rotations around the angular momentum), the secular space has \(4\) dimensions, with coordinates \((g_1,G_1,g_2,G_2)\) outside coplanar or circular motions.

**Lemma** * The quadrupolar system \(\mathrm{Quad}\,\) is
integrable.*

Indeed, it happens that \(\mathrm{Quad}\,\) does not depend on the argument
\(g_2\) of the pericenter of the outer ellipse (but the next higher
order term, the ``octupolar term*, does), thus proving its*
integrability:
\[
\mathrm{Quad}\, = - \frac{\mu_{1} m_2 a_{1}^2}{8 a_2^3
\left(1-e_2^2\right)^{3/2}}
\left[
\begin{array}[c]{l}
\left(15 e_1^2\,\cos ^2 g_1 -12 e_1^2 -3 \right)\,\sin
^2(i_2-i_1)\\
+3 e_1^2+2
\end{array}
\right],
\]
where \(i_j\) is the inclination of the ellipse of \(Q_j\) with respect to
the Laplace plane (Lidov-Ziglin, 1976); the Hamiltonian in the
plane problem is simply obtained by letting \(i_1 = i_2 =0\).

We now need to estimate the frequencies and the torsion of the quadrupolar system, somewhere in the secular space. Lidov-Ziglin (Lidov-Ziglin, 1976) have established the bifurcation diagram of the system, and proved the existence of \(5\) regimes in the parameter space, according to the number of equilibrium points of the reduced quadrupolar system. Here, for the sake of simplicity we will localize our study in some regular region (i.e. a region with a uniform action-angle coordinate system). We will focus on the neighborhood of the origin of the secular space, i.e. circular horizontal Keplerian ellipses. See (Laskar-Robutel, 1995) for more details on the computations.

**Lemma (Lagrange, Laplace)** * The first quadrupolar system has
a degenerate elliptic singularity at the origin of the secular space,
whose normal frequency vector is \[\alpha_{\mathrm{Quad}\,} (0) = -
\frac{3a_1^2}{4a_2^3}
\begin{pmatrix}
\Lambda_1^{-1}\\
\Lambda_2^{-1}\\
-\Lambda_1^{-1}-\Lambda_2^{-1}\\
0
\end{pmatrix}.
\]
*

Indeed, the following steps lead to the wanted expansion of \(\mathrm{Quad}\,\):

- Using elementary geometry, express \(\cos \theta_{12}\) in terms of

the elliptic elements and the true anomalies. Then substitute the variable \(u_1\) for \(v_1\), using the relations \[\cos(v_1) = \frac{a_1}{\|X_1\|} (\cos u_1 - e_1) \quad \mbox{and} \quad \sin v_1 = \frac{a_1}{\|X_1\|} \sqrt{1-e_1^2} \sin u_1.\]

- Multiply \(\mathrm{Main}\,\) by the Jacobian of the change of angles

\[\frac{d\ell_1 \, d\ell_2}{du_1\, dv_2} = \frac{\|X_1\|}{a_1} \, \frac{\|X_2\|^2}{a_2\sqrt{1-e_2^2}}.\]

- Express the distances to the Sun in terms of the inner eccentric

anomaly \(u_1\) and outer true anomaly \(v_2\): \[\|X_1\| = a_1(1-e_1\cos u_1) \quad \mbox{and} \quad \|X_2\| = \frac{a_2(1-e_2^2)}{1+e_2\cos v_2},\] and expand at the second order with respect to eccentricitites and inclinations (even powers vanish; the fourth order will be needed later).

- The obtained expression is trigonometric polynomial in \(u_1\)

and \(v_2\). Average it.

- Switch to the Poincaré coordinates \((\xi_j,\eta_j,p_j,q_j)\), which

are symplectic and analytic in the neighborhood of circular horizontal ellipses, and are defined by the relations \[ \begin{cases} \xi_j + i \eta_j = \sqrt{2L_j}\sqrt{1-\sqrt{1-e_j^2}}\, e^{-i(g_j+\theta_j)}\\ p_j + i q_j = \sqrt{2L_j}\sqrt{\sqrt{1-e_j^2}(1-\cos i_j)} \, e^{-i \theta_j}; \end{cases} \] here we use the notation \(\Lambda_j = L_j = \mu_j \sqrt{M_ja_j}\).

The invariance of \(\mathrm{Quad}\,\) by horizontal rotations entails
that, as proved by Lagrange and Laplace, there exist two quadratic
forms \(Q_h\) and \(Q_v\) (indices \(h\) and \(v\) here stand for
``horizontal* and ``vertical*) on \(\R^2\) such that
\[\mathrm{Quad}\, = - \frac{3a_1^2}{8a_2^3} \left( Q_h(\xi) +
Q_h(\eta) + Q_v(p) + Q_v(q) + O_4(\xi,\eta,p,q) \right).\]
The computation shows that
\[
\begin{cases}
Q_h(\xi) = \frac{\xi_1^2}{L_1} + \frac{\xi_2^2}{L_2}\\
Q_v(p) = - \frac{p_1^2}{L_1} - \frac{p_2^2}{L_2} +
\frac{2p_1p_2}{\sqrt{L_1L_2}}
\end{cases}
\]
The horizontal part is already in diagonal form. The vertical part
\(Q_v\) is diagonalized by the orthogonal operator of \(\R^2\)
\[\rho = \frac{1}{\sqrt{\Lambda_1 + \Lambda_2}}
\begin{pmatrix}
\phantom{-} \sqrt{\Lambda_2} &\sqrt{\Lambda_1}\\
-\sqrt{\Lambda_1} &\sqrt{\Lambda_2}\\
\end{pmatrix}.
\]
This operator of \(\R^2\) lifts to a symplectic operator
\begin{align*}
\tilde \rho \; : \; &(x_1,y_1,x_2,y_2,x_3,y_3,x_4,y_4)\\
&\mapsto (\xi_1,\eta_1,\xi_2,\eta_2,p_1,q_1,p_2,q_2) =
(x_1,y_1,x_2,y_2, p_1,q_1,p_2,q_2)
\end{align*}
with \[(p_1,p_2) = \rho \cdot (x_3,x_4) \quad \mbox{and} \quad
(q_1,q_2)= \rho \cdot (y_3,y_4).\]
and
\[
\mathrm{Quad}\, = - \frac{a_1^2}{4a_2^3} \left(
\begin{array}[c]{l}
1 +
\frac{3}{2\Lambda_1}(x_1^2 + y_1^2) + \frac{3}{2\Lambda_2}(x_2^2 +
y_2^2) - \\
\frac{3}{2} \left( \frac{1}{\Lambda_1} +
\frac{1}{\Lambda_2} \right) \left( x_3^2 + y_3^2 \right) +
O_4(x,y)
\end{array}\right),
\]
thus showing that the origin is elliptic (as mentioned above), and
degenerate (since there is no term in \(x_4^2+y_4^2\)).

Switching (outside the origin) to symplectic polar coordinates \((\tilde\varphi_j,\tilde r_j)_{j=1,...,4}\) defined by \[x_j + i y_j = \sqrt{2\tilde r_j} \, e^{-i\tilde\varphi_j},\] one gets the wanted expression of \(\alpha_{\mathrm{Quad}\,}(0) = \frac{\partial \mathrm{Quad}\,}{\partial \tilde r} (0)\). \(\Box\)

It is an exercise (e.g. using generating functions) to check that all the changes of coordinates we have made on the secular space lift to changes of coordinates in the full phase space, up to adequately modifying the mean longitude. This does not change the Keplerian frequencies.

The quadrupolar frequency vector \(\alpha_{\mathrm{Quad}\,}(0)\) calls for some comments:

- Due to the \(SO(3)\)-symmetry, rotations of the two inner ellipses

around a horizontal axis leave \(\mathrm{Quad}\,\) invariant. Hence the infinitesimal generators of such rotations (last two columns of the matrix \(\tilde\rho\) above) span an eigenplane of the quadratic part, with eigenvalue \(0\). This explains for the vanishing last component of the normal frequency vector (for all \(r\)'s for that matter): \[ \alpha_{\mathrm{Quad}\,}(0)_4 = 0 \qquad (\forall \Lambda_1, \Lambda_2).\]

- Unexpectedly, the sum of the frequencies vanishes:

\[ \sum_{1 \leq j \leq 4} \alpha_{\mathrm{Quad}\,}(0)_j = 0 \qquad (\forall \Lambda_1, \Lambda_2). \]

- The local image of the map \((\Lambda_1, \Lambda_2) \mapsto \alpha_{\mathrm{Quad}\,}(0)\) thus lies in a \(2\)-plane of \(\R^4\)

but in no line, since the map \[(\Lambda_1, \Lambda_2) \mapsto - \frac{3a_1^2}{4a_2^3} \begin{pmatrix} \Lambda_1^{-1}\\ \Lambda_2^{-1} \end{pmatrix} = - \frac{3}{4} \frac{M_2^3}{M_1^2} \frac{\mu_2^6}{\mu_1^4} \begin{pmatrix} \Lambda_1^3 \Lambda_2^{-6}\\ \Lambda_1^4 \Lambda_2^{-7} \end{pmatrix}\] is a diffeomorphism. Hence, additional resonances may always be removed by slightly shifting \(\Lambda_1\) and \(\Lambda_2\).

## Torsion

Resonances a priori prevent from eliminating resonant terms in the Lindstedt (or Birkhoff) normal form of \(\mathrm{Quad}\,\). And they will not disappear by adjusting the \(\Lambda_j\)'s. But, as the following lemma shows, there are no resonant terms at the second order.

Let \[\mathcal{L}(2) = \left\{(\Lambda_1,\Lambda_2) \in \R^2, \; \forall |k|\leq 4, \; k_1 \Lambda_1^{-1} + k_2 \Lambda_2^{-1} \neq 0 \right\}\] be the open set of values of \((\Lambda_1,\Lambda_2)\) for which the horizontal first quadrupolar frequency vector satisfies no resonance of order \(\leq 4\).

**Lemma** *If the parameters \((\Lambda_1,\Lambda_2)\) belong to
\(\mathcal{L}(2)\), \(\mathrm{Quad}\,\) has a non-resonant Lindstedt
normal form at order \(2\) i.e., there exist coordinates
\((\varphi_j,r_j)_{j=1,...,4}\), tangent to \((\tilde\varphi_j,\tilde
r_i)_{j=1,...,4}\), such that \[\mathrm{Quad}\, = \mathrm{cst} +
\alpha_{\mathrm{Quad}\,}(0) \cdot r + \frac{1}{2}
\tau_{\mathrm{Quad}\,} \cdot r^2 + O(r^3).\] Besides, the torsion is
\[\tau_{\mathrm{Quad}\,} =
\left(
\begin{array}[c]{ccc|c}
&&&0\\
&\bar\tau_{\mathrm{Quad}\,}&&0\\
&&&0\\
\hline
0 &0 &0 &0
\end{array}
\right)
\]
with \(\bar\tau_{\mathrm{Quad}\,} =\)
\[- \frac{a_1^2}{8a_2^3}\frac{1}{\Lambda_1^2 \Lambda_2^2}
\begin{pmatrix}
-6\,\Lambda_2^2&18\,\Lambda_1\,\Lambda_2&-24\,\Lambda_2^2-
18\,\Lambda_1\,\Lambda_2\\
18\,\Lambda_1\,\Lambda_2&24\,\Lambda_1^2
&-18\,\Lambda_1\,\Lambda_2-24\,\Lambda_1^2\\
-24\,\Lambda_2^2- 18\,\Lambda_1\,\Lambda_2&-18\,
\Lambda_1\,\Lambda_2-24\,\Lambda_1^2&6
\,\Lambda_2^2+18\,\Lambda_1\,\Lambda_2+6\,\Lambda_1^2
\end{pmatrix}.
\]
*

We carry out the same computation as before, now up to the order \(2\) in the \(\tilde r_j\)'s. The truncated expression is a trigonometric polynomial in the angles \(\tilde\varphi_j\), of degree \(\leq 4\). Eliminating non-resonant monomials, i.e. functions of \(k \cdot \varphi\) with \(k \cdot \alpha_{\mathrm{Quad}\,}(0)\), is a classical matter. Two kinds of terms cannot be eliminated by averaging:

- Monomials in the angle \(4\tilde\varphi_4\)

- Monomials in \(\tilde\varphi_1 + \tilde\varphi_2 + \tilde\varphi_3 + \tilde\varphi_4\)

Such monomials actually cannot occur in the expansion, due to the invariance by rotations (they would not satisfy the d'Alembert relation (Chierchia-Pinzari, 2011 or Malige-Robutel-Laskar, 2002). <ref>Thank you to Gabriella Pinzari for clarifying this assertion to me.</ref>

A direct computation leads to the given expression of the torsion \(\tau_{\mathrm{Quad}\,}\). \(\Box\)

Note that the torsion \(\tau_{\mathrm{Quad}\,}\), as a function of
\(\Lambda_1\) and \(\Lambda_2\), extends analytically outside
\(\mathcal{L}(2)\) (as often do first order normal forms). This allows
us to define the *quadrupolar frequency map*
\[\alpha_{\mathrm{Quad}\,} : r \mapsto \alpha_{\mathrm{Quad}\,} +
\frac{1}{2} \tau_{\mathrm{Quad}\,} \cdot r,\]
a first order approximation of the normal frequencies.

**Proposition** * The first quadrupolar frequency map has
constant rank \(3\) and, in restriction to the symplectic submanifold
obtained by fixing vertically the direction of the angular momentum,
is a local diffeomorphism. *

To see this, we denote by \(C = (C_x,C_y,C_z) \in \R^3\) the angular momentum of the first two planets. The submanifold \(\mathcal{V}\) of vertical angular momentum, has equation \(C_x = C_y = 0\). It is a symplectic, codimension-\(2\) submanifold, transverse to the Hamiltonian vector fields \(X_{C_x}\) and \(X_{C_y}\) of \(C_x\) and \(C_y\). Since it is invariant by the flow of \(\mathrm{Quad}\,\), its tangent space has equations, in the coordinates \((x_j,y_j)_{j=1,...,4}\), \(x_4 = y_4 = 0\). So, the upper left \(3 \times 3\) submatrix \(\tilde\tau_{\mathrm{Quad}\,}\) of \(\tau_{\mathrm{Quad}\,}\) is the Hessian of the restriction of \(\mathrm{Quad}\,\) to \(\mathcal{V}\). In order to conclude, one merely needs to notice that the determinant of the torsion \(\tilde\tau_{\mathrm{Quad}\,}\): \[\det \tilde\tau_{\mathrm{Quad}\,} = - \left(\frac{a_1^2}{a_2^3}\right)^3 \frac{27}{64 \Lambda_1^2 \Lambda_2^2} \left( 39 \Lambda_1^2 + 39 \Lambda_1 \Lambda_2 + 4 \Lambda_2^2 \right)\] is non-zero. \(\Box\)

So, \(\mathrm{Quad}\,\) (adequately truncated) has a non-degenerate quasiperiodic dynamics in the three degrees of freedom corresponding to coordinates \((\psi_j,s_i)_{j=1,2,3}\).

Therefore, there exist invariant tori which have frequencies of the following order (assuming \(a_1 = O(1)\) and \(a_2 \to \infty\)): \[\alpha = O(1, a_2^{-3/2}, a_2^{-3}, a_2^{-7/2}, a_2^{-3}, 0).\]

## References

- Arnold, V. I. (1963).
*Russian Math. Survey*18: 13-40.

- Arnold, V. I. (1963).
*Russian Math. Survey*18: 85-191.

- Bertrand, J. (1873). Théorème relatif au mouvement d'un point attiré vers un centre fixe.
*C. R. Math. Acad. Sci. Paris*73: 849-853.

- Chierchia, A. (2007). KAM stability and Celestial Mechanics.
*Memoirs of the Amerc. Math. Soc.*187: viii+134.

- Chierchia, L. and Pinzari, G. (2011). Planetary Birkhoff normal forms.
*Journal of Modern Dynamics*5: 623-664.

- Féjoz, J. (2002). Global Secular Dynamics in the Planar Three-Body Problem.
*Celestial Mechanics and Dynamical Astronomy*84: 159-195.

- Féjoz, J. (2013). On action-angle coordinates and the Poincaré coordinates.
*Regular and Chaotic Dynamics*18: 703-718.

- Jefferys, W. H. and Moser, J. (1966). Quasi-periodic solutions for the three-body problem.
*The Astronomical Journal*71: 568-578.

- Laskar, J. and Robutel, Ph. (1995). Stability of the planetary three-body problem. {I}. Expansion of the planetary Hamiltonian.
*Celestial Mech. Dynam. Astronom.*62: 193-217.

- Lidov, M. L. and Ziglin, S. L. (1976). Non-restricted double-averaged three body problem in {H}ill's case.
*Celestial Mech.*13: 471-489.

- Malige, F.; Robutel, Ph. and Laskar, J. (2002). Partial Reduction in the N-Body Planetary Problem using the Angular Momentum Integral.
*Celestial Mechanics and Dynamical Astronomy*84: 283-316.

- Poincaré, H. (1892-1899).
*Les méthodes nouvelles de la mécanique céleste*. Gauthier-Villars.

## Further reading

- Arnold, V. I. (1989).
*Mathematical methods of classical mechanics*. Springer-Verlag. ISBN 978-0-387-96890-2.

- Arnold, V. I.; Kozlov, V. V. and Neishadt, A. I. (2006).
*Mathematical aspects of classical and celestial mechanics*. Springer-Verlag. ISBN 978-3-540-28246-4.

- Meyer, K. R. and Hall, G. R. (1992).
*Introduction to Hamiltonian Dynamical Systems and the N-Body Problem*. Springer-Verlag. ISBN 978-0-387-09723-7.

- Knauf, A. (1989).
*Mathematische Physik: Klassische Mechanik*. Springer-Verlag. ISBN 978-3-642-20977-2.

- Moser, J. (1973).
*Stable and random motions in dynamical systems. With special emphasis on celestial mechanics*. Annals of Mathematics Studies, Princeton. ISBN 9780691089102.

- Siegel, C. L. and Moser, J. (1971).
*Lectures on celestial mechanics*. Springer-Verlag. ISBN 9780387054193.

- Sternberg, S. (1969).
*Celestial Mechanics*. W.A. Benjamin. ISBN 9780805391039.

## External links

## See also

Averaging, Aubry-Mather theory, Chaos, Computational celestial mechanics, Dynamical Systems, Hamiltonian Dynamics, Hamiltonian Normal Forms, Kolmogorov-Arnold-Moser theory, Normal Forms, Symplectic maps, Three body problem,