# KAM theory for the three-body problem

This article has not yet been published; it may contain inaccuracies, unapproved changes, or be unfinished.

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.

## 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.

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.

• 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.
• 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.