GENERAL ARTICLE

An Introduction to the Classical Three-Body Problem

From Periodic Solutions to Instabilities and Chaos

Govind S Krishnaswami and Himalaya Senapati

The classical three-body problem arose in an attempt to un-

derstand the effect of the Sun on the Moon’s Keplerian orbit
around the Earth. It has attracted the attention of some of the
best physicists and mathematicians and led to the discovery
of ‘chaos’. We survey the three-body problem in its historical
context and use it to introduce several ideas and techniques
that have been developed to understand classical mechanical
(Left) Govind Krishnaswami
systems. is on the faculty of the
Chennai Mathematical
Institute. His research
concerns various problems in
1. Introduction theoretical and mathematical
physics.
The three-body problem is one of the oldest problems in classical
dynamics that continues to throw up surprises. It has challenged (Right) Himalaya Senapati is
scientists from Newton’s time to the present. It arose in an at- a PhD student at the Chennai
Mathematical Institute. He
tempt to understand the Sun’s effect on the motion of the Moon
works on dynamical systems
around the Earth. This was of much practical importance in ma- and chaos.
rine navigation, where lunar tables were necessary to accurately
determine longitude at sea (see Box 1).
The study of the three-body problem led to the discovery of the
planet Neptune (see Box 2), it explains the location and stability
of the Trojan asteroids and has furthered our understanding of the
stability of the solar system [1]. Quantum mechanical variants of Keywords
the three-body problem are relevant to the helium atom and water Kepler problem, three-body prob-
lem, celestial mechanics, classi-
molecule [2].
cal dynamics, chaos, instabilities.
The three-body problem admits many ‘regular’ solutions such as
the collinear and equilateral periodic solutions of Euler and La-
grange as well as the more recently discovered figure-8 solution.
On the other hand, it can also display chaos as serendipitously dis- Euler had gone blind
covered by Poincaré. Though a general solution in closed form when he developed much
of his lunar theory!

RESONANCE | January 2019 https://doi.org/10.1007/s12045-019-0760-1 87

 GENERAL ARTICLE

is not known, Sundman while studying binary collisions, discov-

ered an exceptionally slowly converging series representation of
solutions in fractional powers of time.
The three-body problem
The importance of the three-body problem goes beyond its appli-
provides a context in
which to study the cation to the motion of celestial bodies. As we will see, attempts
development of classical to understand its dynamics have led to the discovery of many phe-
dynamics as well as a nomena (e.g., the abundance of periodic motions, resonances (see
window into several
Box 3), homoclinic points, collisional and non-collisional singu-
areas of mathematics
(geometry, calculus and larities, chaos and KAM tori) and techniques (e.g., Fourier series,
dynamical systems). perturbation theory, canonical transformations and regularization
of singularities) with applications across the sciences.

Box 1. Longitude Act

The Longitude Act (1714) of the British Parliament offered £20,000 for a method to determine the longitude
at sea to an accuracy of half a degree. This was important for marine navigation at a time of exploration
of the continents. In the absence of accurate clocks that could function at sea, a lunar table along with the
observed position of the Moon was the principal method of estimating the longitude. Leonhard Euler, Alexis
Clairaut and Jean-Baptiste d’Alembert competed to develop a theory accounting for solar perturbations to
the motion of the Moon around the Earth. For a delightful account of this chapter in the history of the
three-body problem, including Clairaut’s explanation of the annual 40◦ rotation of the lunar perigee (which
had eluded Newton), see [3]. Interestingly, Clairaut’s use of Fourier series in the three-body problem (1754)
predates their use by Joseph Fourier in the analysis of heat conduction!

Box 2. Discovery of Neptune

The French mathematical astronomer Urbain Le Verrier (1846) was intrigued by the discrepancies between
the observed and Keplerian orbits of Mercury and Uranus. He predicted the existence of Neptune (as was
widely suspected) and calculated its expected position based on its effects on the motion of Uranus around
the Sun (the existence and location of Neptune was independently inferred by John Adams in Britain). The
German astronomer Johann Galle (working with his graduate student Heinrich d’Arrest) discovered Nep-
tune within a degree of Le Verrier’s predicted position on the very night that he received the latter’s letter.
It turned out that both Adams’ and Le Verrier’s heroic calculations were based on incorrect assumptions
about Neptune, they were extremely lucky to stumble upon the correct location!

88 RESONANCE | January 2019

 GENERAL ARTICLE

Box 3. Orbital Resonances

The simplest example of an orbital resonance occurs when the periods of two orbiting bodies (e.g., Jupiter
and Saturn around the Sun) are in a ratio of small whole numbers (T S /T J ≈ 5/2). Resonances can enhance
their gravitational interaction and have both stabilizing and destabilizing effects. For instance, the moons
Ganymede, Europa and Io are in a stable 1 : 2 : 4 orbital resonance around Jupiter. The Kirkwood gaps
in the asteroid belt are probably due to the destabilizing resonances with Jupiter. Resonances among the
natural frequencies of a system (e.g., Keplerian orbits of a pair of moons of a planet) often lead to difficulties
in naive estimates of the effect of a perturbation (say of the moons on each other).

2. Review of the Kepler Problem

If we ignore the non-zero size of celestial bodies, Newton’s sec- As preparation for the
three-body problem, we
ond law for the motion of two gravitating masses states that begin by reviewing some
key features of the
(r2 − r1 ) (r1 − r2 )
m1 r̈1 = α and m2 r̈2 = α . (1) two-body problem.
|r1 − r2 |3 |r1 − r2 |3

Here, α = Gm1 m2 measures the strength of the gravitational at-

traction and ‘dots’ denote time derivatives. This system has six
degrees of freedom, say the three Cartesian coordinates of each
mass r1 = (x1 , y1 , z1 ) and r2 = (x2 , y2 , z2 ). Thus, we have a sys-
tem of 6 nonlinear (due to division by |r1 − r2 |3 ), second-order or-
dinary differential equations (ODEs) for the positions of the two
masses. It is convenient to switch from r1 and r2 to the center of
mass (CM) and relative coordinates An advantage of the CM
and relative coordinates
m1 r1 + m2 r2 is that in the absence of
R= and r = r2 − r1 . (2)
m1 + m2 external forces, the CM
moves at constant
In terms of these, the equations of motion (EOM) become velocity, which can be
chosen to vanish by
M R̈ = 0 and mr̈ = −αr/|r|3 . (3) going to a frame moving
with the CM.
Here, M = m1 + m2 is the total mass and m = m1 m2 /M the ‘re-
duced’ mass. The motion of the relative coordinate r decouples
from that of R and describes a system with three degrees of free-
dom r = (x, y, z). Expressing the conservative gravitational force
in terms of the gravitational potential V = −α/|r|, the equation

RESONANCE | January 2019 89

 GENERAL ARTICLE

for the relative coordinate r becomes

 
∂V ∂V ∂V
ṗ ≡ mr̈ = −∇r V = − , , , (4)
∂x ∂y ∂z
where p = mṙ is the relative momentum. Taking the dot product
with the ‘integrating factor’ ṙ = ( ẋ, ẏ, ż), we get
   
d 1 2 ∂V ∂V ∂V dV
mṙ · r̈ = mṙ = − ẋ + ẏ + ż = − , (5)
dt 2 ∂x ∂y ∂z dt

which implies that the energy E ≡ 1

2 mṙ
2 + V or Hamiltonian
p2
2m + V is conserved. The relative angular momentum L = r ×
The conservation of mṙ = r × p is another constant of motion as the force is central :
angular momentum in a L̇ = ṙ × p + r × ṗ = 0 + 0. The constancy of the direction of L
central force is a implies planar motion in the CM frame: r and p always lie in the
consequence of rotation
invariance: V = V(|r|) is ‘ecliptic plane’ perpendicular to L, which we take to be the x-y
independent of polar and plane with origin at the CM (see Figure 1). The Kepler problem
azimuthal angles. More is most easily analyzed in plane-polar coordinates r = (r, θ) in
generally, Noether’s which the energy E = 12 mṙ2 + Veff (r) is the sum of a radial kinetic
theorem relates
continuous symmetries energy and an effective potential energy Veff = L2z /(2mr2 ) + V(r).
to conserved quantities. Here, Lz = mr2 θ̇ is the vertical component of angular momentum
and the first term in Veff is the centrifugal ‘angular momentum
barrier’. Since L (and therefore Lz ) is conserved, Veff depends
only on r. Thus, θ does not appear in the Hamiltonian: it is a
‘cyclic’ coordinate. Conservation of energy constrains r to lie
between ‘turning points’, i.e., zeros of E − Veff (r) where the radial
velocity ṙ momentarily vanishes. One finds that the orbits are
Keplerian ellipses for E < 0 along with parabolae and hyperbolae
for E ≥ 0: r(θ) = ρ(1 +  cos θ)−1 [4, 5]. Here, ρ = L2z /mα is the
radius of the circular orbit corresponding to angular momentum
α
Lz ,  the eccentricity and E = − 2ρ (1 −  2 ) the energy.
In addition to E and L, the Laplace–Runge–Lenz (LRL) vector
The LRL vector A points A = p × L − mα r̂ is another constant of motion. Thus, we have 7
along the semi-major conserved quantities: energy and three components each of L and
axis from the CM to the
A. However, a system with three degrees of freedom has a six-
perihelion and its
magnitude determines dimensional phase space (space of coordinates and momenta, also
the eccentricity of the called the state space) and if it is to admit continuous time evolu-
orbit. tion, it cannot have more than 5 independent conserved quantities.

90 RESONANCE | January 2019

 GENERAL ARTICLE

Figure 1. Keplerian el-

lipse in the ecliptic plane of
motion showing the constant
LRL vector A. The constant
angular momentum L points
out of the ecliptic plane.

The apparent paradox is resolved once we notice that E, L and A

are not all independent; they satisfy two relations:

A 2 − m 2 α2
L · A = 0 and E= . (6)
2mL2

Newton used the solution of the two-body problem to understand Wolfgang Pauli (1926)
the orbits of planets and comets. He then turned his attention derived the quantum
to the motion of the Moon around the Earth. Lunar motion is mechanical spectrum of
the Hydrogen atom
significantly affected by the Sun. For instance, A is not conserved using the relation
and the lunar perigee rotates by 40◦ per year. Thus, he was led to between E, L2 and A2
study the Moon–Earth–Sun three-body problem. before the development
of the Schrödinger
equation. Indeed, if we
3. The Three-Body Problem postulate circular Bohr
orbits which have zero
eccentricity (A = 0) and
We consider the problem of three point masses (ma with position
quantized angular
vectors ra for a = 1, 2, 3) moving under their mutual gravitational momentum L2 = n2 2 ,
2
attraction. This system has 9 degrees of freedom, whose dynam- then En = − 2mα2 n2 where
ics is determined by 9 coupled second order nonlinear ODEs: α = e2 /4π0 is the
electromagnetic
d2 ra  rb − ra analogue of Gm1 m2 .
ma = Gma mb for a = 1, 2 and 3. (7)
dt 2
ba
|rb − ra |3

As before, the three components of momentum P = a ma ṙa ,

three components of angular momentum L = a ra × pa and
energy
1  Gm m
3
a b
E= ma ṙ2a − ≡T +V, (8)
2 a=1 |r
a<b a
− rb |

RESONANCE | January 2019 91

 GENERAL ARTICLE

Figure 2. Jacobi vectors

J1 , J2 and J3 for the three-
body problem. O is the ori-
gin of the coordinate system
while CM12 is the center of
mass of particles 1 and 2.

furnish 7 independent conserved quantities. Lagrange used these

conserved quantities to reduce the above EOM to 7 first order
ODEs (see Box 4).

Box 4. Lagrange’s Reduction From 18 to 7 Equations

The 18 phase space variables of the three-body problem (components of r1 , r2 , r3 , p1 , p2 , p3 ) satisfy 18 first
order ODEs ṙa = pa , ṗa = −∇ra V. Lagrange (1772) used the conservation laws to reduce these ODEs to a
system of 7 first order ODEs. Conservation of momentum determines 6 phase space variables comprising

the location RCM and momentum P of the CM. Conservation of angular momentum L = ra × pa and
energy E lead to 4 additional constraints. By using one of the coordinates as a parameter along the orbit (in
place of time), Lagrange reduced the three-body problem to a system of 7 first order nonlinear ODEs.

Jacobi vectors (see Figure 2) generalize the notion of CM and

relative coordinates to the three-body problem [6]. They are de-
A quadratic form fined as J1 = r2 − r1 ,

a,b a Qab rb is diagonal
r
if Qab = 0 for a  b. m1 r1 + m2 r2 m1 r1 + m2 r2 + m3 r3
J2 = r3 − and J3 = . (9)
Here, M1−1 = m−1 1 + m2
−1 m1 + m2 m1 + m2 + m3
is the reduced mass of
the first pair, J3 is the coordinate of the CM, J1 the position vector of m2 rel-
M2−1 = (m1 + m2 )−1 + m−13
ative to m1 and J2 that of m3 relative to the CM of m1 and m2 .
is the reduced mass of A nice feature of Jacobi vectors is that the kinetic energy T =
m3 and the (m1 , m2 ) 1 
2 a=1,2,3 ma ṙa and moment of inertia I = a=1,2,3 ma ra , regarded
2 2
system and
M3 = m1 + m2 + m3 the as quadratic forms, remain diagonal:
total mass. 1  
T = Ma J̇2a and I= Ma J2a . (10)
2 1≤a≤3 1≤a≤3

92 RESONANCE | January 2019

 GENERAL ARTICLE

Figure 3. (a, b) Euler

collinear periodic solutions
with three masses traversing
Keplerian ellipses/circles
with one focus at the CM.
The constant ratios of
separations are functions of
the mass ratios alone. (c)
Lagrange’s periodic solution
Moreover, just as the potential energy −α/|r| in the 2-body prob- with three bodies at vertices
lem is a function only of the relative coordinate r, here the poten- of equilateral triangles.
tial energy V may be expressed entirely in terms of J1 and J2 :
Gm1 m2 Gm2 m3 Gm3 m1 m1,2
V=− − − where μ1,2 = .
|J1 | |J2 − μ1 J1 | |J2 + μ2 J1 | m1 + m2
(11)
Thus, the components of the CM vector J3 are cyclic coordinates
in the Hamiltonian H = T + V. In other words, the center of mass
motion (J̈3 = 0) decouples from that of J1 and J2 .
An instantaneous configuration of the three bodies defines a tri-
angle with masses at its vertices. The moment of inertia about
the center of mass ICM = M1 J21 + M2 J22 determines the size of
the triangle. For instance, particles suffer a triple collision when
ICM → 0 while ICM → ∞ when one of the bodies flies off to
infinity.

4. Euler and Lagrange Periodic Solutions

The planar three-body problem is the special case where the masses In Lagrange’s solutions,
always lie on a fixed plane. For instance, this happens when the bodies lie at vertices of
CM is at rest (J̇3 = 0) and the angular momentum about the CM equilateral triangles
while they are collinear
vanishes (LCM = M1 J1 × J̇1 + M2 J2 × J̇2 = 0). In 1767, the in Euler’s solutions. In
Swiss scientist Leonhard Euler discovered simple periodic solu- both cases, the force on
tions to the planar three-body problem where the masses are al- each body is always
ways collinear, with each body traversing a Keplerian orbit about toward the common
center of mass and
their common CM. The line through the masses rotates about the proportional to the
CM with the ratio of separations remaining constant (see Fig- distance from it.
ure 3). The Italian/French mathematician Joseph-Louis Lagrange
rediscovered Euler’s solution in 1772 and also found new periodic

RESONANCE | January 2019 93

 GENERAL ARTICLE

solutions where the masses are always at the vertices of equilat-

eral triangles (see Figure 3c).

In Lagrange’s equilateral It is convenient to identify the plane of motion with the complex
solution, the size and plane C and let the three complex numbers za=1,2,3 (t) denote the
angular orientation of positions of the three masses at time t. For e.g., the real and imag-
the triangle may change
with time. In the limiting inary parts of z1 denote the Cartesian components of the position
case of zero angular vector r1 of the first mass. In Lagrange’s solutions, za (t) lie at
momentum, the three the vertices of an equilateral triangle while they are collinear in
bodies move toward or Euler’s solutions. In both cases, the force on each body is always
away from their CM
along straight lines. toward the common CM and proportional to the distance from it.
These implosion or For instance, the force on m1 in a Lagrange solution is:
explosion solutions are
called Lagrange r2 − r1 r3 − r1
F1 = Gm1 m2 + Gm1 m3
homotheties. |r2 − r1 |3 |r3 − r1 |3
Gm1
= (m1 r1 + m2 r2 + m3 r3 − M3 r1 ) , (12)
d3

where d = |r2 − r1 | = |r3 − r1 | is the side-length of the equilateral

triangle and M3 = m1 + m2 + m3 . Recalling that rCM = (m1 r1 +
m2 r2 + m3 r3 )/M3 , we get

Gm1 rCM − r1
F1 = M3 (rCM − r1 ) ≡ Gm1 δ1 , (13)
d 3 |rCM − r1 |3

1 Indeed, M3 (rCM − r1 ) = where δ1 = M3 |rCM − r1 |3 /d3 is a function of the masses alone1 .

m2 (r2 − r1 ) + m3 (r3 − r1 ) ≡ Thus, the equation of motion for m1 ,
m2 b + m3 c where b and c
are two of the sides of the
equilateral triangle of length d.
m1 r̈1 = Gm1 δ1 (rCM − r1 )/(|rCM − r1 |3 ) , (14)
 to |(rCM − r1 )/d| =
This leads
(1/M3 ) m22 + m23 + m2 m3 takes the same form as in the two-body Kepler problem [see (1)].
which is a function of masses The same applies to m2 and m3 . So if za (0) denote the initial
alone.
positions, the curves za (t) = z(t)za (0) are solutions of Newton’s
equations for three bodies provided z(t) is a Keplerian orbit for
an appropriate two-body problem. In other words, each mass
traverses a rescaled Keplerian orbit about the common centre of
mass. A similar analysis applies to the Euler collinear solutions
as well: locations of the masses is determined by the requirement
that the force on each one is toward the CM and proportional to
the distance from it (see Box 5 on central configurations).

94 RESONANCE | January 2019

 GENERAL ARTICLE

Box 5. Central Configurations

Three-body configurations in which the acceleration of each particle points towards the CM and is propor-
tional to its distance from the CM (ab = ω2 (RCM − rb ) for b = 1, 2, 3) are called ‘central configurations’.
A central configuration rotating at angular speed ω about the CM automatically satisfies the equations of
motion (7). Euler collinear and Lagrange equilateral configurations are the only central configurations in
the three-body problem. In 1912, Karl Sundmann showed that triple collisions are asymptotically central
configurations.

5. Restricted Three-Body Problem

The restricted three-body problem is a simplified version of the

three-body problem where one of the masses m3 is assumed much
smaller than the primaries m1 and m2 . Thus, m1 and m2 move The Sun–Earth–Moon
in Keplerian orbits which are not affected by m3 . In the pla- system provides an
example of the restricted
nar circular restricted three-body problem, the primaries move in
three-body problem
fixed circular orbits around their common CM with angular speed where we further have
Ω = (G(m1 +m2 )/d3 )1/2 given by Kepler’s third law and m3 moves m2 = m E
m 1 = m S .
in the same plane as m1 and m2 . Here, d is the separation be-
tween the primaries. This system has 2 degrees of freedom asso-
ciated to the planar motion of m3 , and therefore, a 4-dimensional
phase space just like the planar Kepler problem for the reduced
mass. However, unlike the latter which has three conserved quan-
tities (energy, z-component of angular momentum and direction
of LRL vector) and is exactly solvable, the planar restricted three-
body problem has only one known conserved quantity, the ‘Ja-
cobi integral’, which is the energy of m3 in the co-rotating (non-
inertial) frame of the primaries: The ‘Roche’ effective
    potential Veff , named
1 1 1 m1 m2 after the French
E = m3 ṙ + m3 r φ̇ − m3 Ω r −Gm3
2 2 2 2 2
+ ≡ T +Veff .
2 2 2 r1 r2 astronomer Édouard
(15) Albert Roche, is a sum
Here, (r, φ) are the plane polar coordinates of m3 in the co-rotating of centrifugal and
gravitational energies
frame of the primaries with origin located at their CM while r1
due to m1 and m2 .
and r2 are the distances of m3 from m1 and m2 (see Figure 4).
A system with n degrees of freedom needs at least n constants
of motion to be exactly solvable. For the restricted three-body

RESONANCE | January 2019 95

 GENERAL ARTICLE

Figure 4. The secondary

m3 in the co-rotating frame
of primaries m1 and m2
in the restricted three-body
problem. The origin is lo-
cated at the CM of m1 and
m2 which coincides with the
CM of the system since
m3
m1,2 .
problem, Henri Poincaré (1889) proved the nonexistence of any
conserved quantity (other than E) that is analytic in small mass
ratios (m3 /m2 and (m3 + m2 )/m1 ) and orbital elements (J1 , M1 J̇1 ,
A Hamiltonian system J2 and M2 J̇2 ) [7, 8, 9]. This was an extension of a result of
with n degrees of Heinrich Bruns who had proved in 1887 the nonexistence of any
freedom is exactly new conserved quantity algebraic in Cartesian coordinates and
solvable in the sense of
Liouville if it possesses momenta for the general three-body problem [10]. Thus, roughly
n independent conserved speaking, Poincaré showed that the restricted three-body problem
quantities in involution, is not exactly solvable. In fact, as we outline in §7., he discovered
i.e., with vanishing
that it displays chaotic behavior.
pairwise Poisson
brackets (see Boxes 6 Euler and Lagrange points (denoted L1−5 ) of the restricted three-
and 10). body problem are the locations of a third mass (m3
m1 , m2 ) in
the co-rotating frame of the primaries m1 and m2 in the Euler and
Lagrange points L1−5 are Lagrange solutions (see Figure 5). Their stability would allow
also called libration an asteroid or satellite to occupy a Lagrange point. Euler points
(literally, balance)
L1,2,3 are saddle points of the Roche potential while L4,5 are max-
points.
ima (see Figure 6). This suggests that they are all unstable. How-
ever, Veff does not include the effect of the Coriolis force since
it does no work. A more careful analysis shows that the Cori-
olis force stabilises L4,5 . It is a bit like a magnetic force which
does no work but can stabilise a particle in a Penning trap. Euler
points are always unstable while the Lagrange points L4,5 are sta-
ble to small perturbations iff (m1 + m2 )2 ≥ 27m1 m2 [11]. More
generally, in the unrestricted three-body problem, the Lagrange
Stable ‘Halo’ orbits equilateral solutions are stable iff:
around Euler points have
been found numerically. (m1 + m2 + m3 )2 ≥ 27(m1 m2 + m2 m3 + m3 m1 ). (16)

The above criterion due to Edward Routh (1877) is satisfied if

96 RESONANCE | January 2019

 GENERAL ARTICLE

Figure 5. The positions of

Euler (L1,2,3 ) and Lagrange
(L4,5 ) points when m1 
m2  m3 . m2 is in an
approximately circular orbit
around m1 . L3 is almost di-
ametrically opposite to m2
and a bit closer to m1 than
m2 is. L1 and L2 are sym-
metrically located on either
side of m2 . L4 and L5 are
equidistant from m1 and m2
and lie on the circular orbit
of m2 .

Figure 6. Level curves

of the Roche effective po-
tential energy Veff of m3
in the co-rotating frame of
the primaries m1 and m2 in
the circular restricted three-
body problem for G = 1,
m1 = 15, m2 = 10 and
m3 = .1. Lagrange points
L1−5 are at extrema of Veff .
The trajectory of m3 for a
given energy E must lie in
the Hill region defined by
Veff (x, y) ≤ E. E.g., for
E = −6, the Hill region is
one of the masses dominates the other two. For instance, L4,5 the union of two neighbour-
for the Sun–Jupiter system are stable and occupied by the Trojan hoods of the primaries and
asteroids. a neighborhood of the point
at infinity. The lobes of the
∞-shaped level curve pass-
ing through L1 are called
6. Planar Euler Three-Body Problem Roche’s lobes. The sad-
dle point L1 is like a moun-
Given the complexity of the restricted three-body problem, Euler tain pass through which ma-
(1760) proposed the even simpler problem of a mass m moving in terial could pass between the
lobes.

RESONANCE | January 2019 97

 GENERAL ARTICLE

Unlike in the restricted the gravitational potential of two fixed masses m1 and m2 . Initial
three-body problem, in conditions can be chosen so that m always moves on a fixed plane
the Euler three-body, the containing m1 and m2 . Thus, we arrive at a one-body problem
rest-frame of the
primaries is an inertial with two degrees of freedom and energy:
frame, so there are no
centrifugal or Coriolis
forces. This 1
μ1 μ2
E = m ẋ2 + ẏ2 − − . (17)
simplification allows the 2 r1 r2
Euler three-body
problem to be exactly
solved. Here, (x, y) are the Cartesian coordinates of m, ra the distances of
m from ma and μa = Gma m for a = 1, 2 (see Figure 7).
Just as the Kepler problem simplifies in plane-polar coordinates
(r, θ) centered at the CM, the Euler three-body problem simplifies
in an elliptical coordinate system (ξ, η). The level curves of ξ and
η are mutually orthogonal confocal ellipses and hyperbolae (see
Figure 7) with the two fixed masses at the foci 2 f apart:

x = f cosh ξ cos η and y = f sinh ξ sin η. (18)

Here, ξ and η are like the radial distance r and angle θ, whose
level curves are mutually orthogonal concentric circles and radial
rays. The distances of m from m1,2 are r1,2 = f (cosh ξ ∓ cos η).

Figure 7. Elliptical coor-

dinate system for the Euler
three-body problem. Two
masses are at the foci (± f, 0)
of an elliptical coordinate
system with f = 2 on the x-y
plane. The level curves of ξ
and η (confocal ellipses and
hyperbolae) are indicated.

98 RESONANCE | January 2019

 GENERAL ARTICLE

Box 6. Canonical Transformations

We have seen that the Kepler problem is more easily solved in polar coordinates and momenta (r, θ, pr , pθ )
than in Cartesian phase space variables (x, y, px , py ). This change is an example of a canonical trans-
formation (CT). More generally, a CT is a change of canonical phase space variables (q, p) →
(Q(p, q, t), P(p, q, t)) that preserves the form of Hamilton’s equations. For one degree of freedom, Hamil-
ton’s equations q̇ = ∂H ∂p
and ṗ = − ∂H∂q
become Q̇ = ∂K∂P
and Ṗ = − ∂Q∂K
where K(Q, P, t) is the new
Hamiltonian (for a time independent CT, the old and new Hamiltonians are related by substitution:
H(q, p) = K(Q(q, p), P(q, p))). The form of Hamilton’s equations is preserved provided the basic Pois-
son brackets (PB) do not change i.e.,

{q, p} = 1, {q, q} = {p, p} = 0 ⇒ {Q, P} = 1, {Q, Q} = {P, P} = 0. (19)

Here, the Poisson bracket of two functions on phase space f (q, p) and g(q, p) is defined as

∂ f ∂g ∂ f ∂g
{ f (q, p), g(q, p)} = − . (20)
∂q ∂p ∂p ∂q

For one degree of freedom, a CT is simply an area and orientation preserving transformation of the
q-p phase plane. Indeed, the condition {Q, P} = 1 simply states that the Jacobian determinant J =

det ∂Q , ∂Q | ∂P
∂q ∂p
, ∂P = 1 so that the new area element dQ dP = J dq d p is equal to the old one. A CT
∂q ∂p
can be obtained from a suitable generating function, say of the form S (q, P, t), in the sense that the equa-
tions of transformation are given by partial derivatives of S :

∂S ∂S ∂S
p= , Q= and K=H+ . (21)
∂q ∂P ∂t

For example, S = qP generates the identity transformation (Q = q and P = p) while S = −qP generates a
rotation of the phase plane by π (Q = −q and P = −p).

The above confocal ellipses and hyperbolae are Keplerian orbits

when a single fixed mass (m1 or m2 ) is present at one of the foci
(± f, 0). Remarkably, these Keplerian orbits survive as orbits of
the Euler three-body problem. This is a consequence of Bonnet’s
theorem, which states that if a curve is a trajectory in two separate
force fields, it remains a trajectory in the presence of both. If v1
and v2 are the speeds of the Keplerian
 trajectories when only m1
or m2 was present, then v = v21 + v22 is the speed when both are
present.

RESONANCE | January 2019 99

 GENERAL ARTICLE

Box 7. Hamilton–Jacobi Equation

The Hamilton–Jacobi (HJ) equation is an alternative formulation of Newtonian dynamics. Let i = 1, . . . , n

label the degrees of freedom of a mechanical system. Cyclic coordinates qi (i.e., those that do not appear
in the Hamiltonian H(q, p, t) so that ∂H/∂qi = 0) help to understand Newtonian trajectories, since their
conjugate momenta pi are conserved ( ṗi = ∂H ∂qi
= 0). If all coordinates are cyclic, then each of them evolves
∂H ∂H
linearly in time: qi (t) = qi (0) + ∂p i
t. Now time-evolution is even simpler if ∂p i
= 0 for all i as well, i.e.,
if H is independent of both coordinates and momenta! In the HJ approach, we find a CT from old phase
space variables (q, p) to such a coordinate system (Q, P) in which the new Hamiltonian K is a constant
(which can be taken to vanish by shifting the zero of energy). The HJ equation is a nonlinear, first-order
partial differential equation (PDE) for Hamilton’s principal function S (q, P, t) which generates the canonical
∂S ∂S ∂S
transformation from (q, p) to (Q, P). As explained in Box 6, this means pi = ∂q i , Q = ∂P and K = H + ∂t .
j
j
Thus, the HJ equation  
∂S ∂S
H q, ,t + =0 (22)
∂q ∂t
is simply the condition for the new Hamiltonian K to vanish. If H is time-independent, we may ‘separate’
the time-dependence of S by writing S (q, P, t) = W(q, P) − Et where the ‘separation constant’ E may be
interpreted as energy. Thus, the time independent HJ-equation for Hamilton’s characteristic function W is
 
∂W
H q, = E. (23)
∂q

2
∂W
1
E.g., for a particle in a potential V(q), it is the equation 2m ∂q
+ V(q) = E. By solving (23) for W, we find
the desired canonical transformation to the new conserved coordinates Q and momenta P. By inverting the
relation (q, p) → (Q, P) we find (qi (t), p j (t)) given their initial values. W is said to be a complete integral
of the HJ equation if it depends on n constants of integration, which may be taken to be the new momenta
P1 , . . . , Pn . When this is the case, the system is said to be integrable via the HJ equation. However, it is
seldom possible to find such a complete integral. In favourable cases, separation of variables can help to
solve the HJ equation (see Box 8).

When the primaries Bonnet’s theorem however does not give us all the trajectories of
coalesce at the origin the Euler three-body problem. More generally, we may integrate
( f → 0), Whittaker’s the equations of motion by the method of separation of variables
constant reduces to the
conserved quantity L2 of in the Hamilton–Jacobi equation (see [12] and Boxes 6, 7 and
the planar 2-body 8). The system possesses two independent conserved quantities:
problem. energy and Whittaker’s constant [2, 10]

w = L1 · L2 + 2m f (−μ1 cos θ1 + μ2 cos θ2 )

= m2 r12 r22 θ̇1 θ̇2 + 2 f m (−μ1 cos θ1 + μ2 cos θ2 ) . (24)

100 RESONANCE | January 2019

 GENERAL ARTICLE

Here, θa are the angles between the position vectors ra and the
positive x-axis and L1,2 = mr1,22 θ̇ ẑ are the angular momenta
1,2
about the two force centers (Figure 7). Since w is conserved,
it Poisson commutes with the Hamiltonian H. Thus, the planar
Euler three-body problem has two degrees of freedom and two
conserved quantities in involution. Consequently, the system is
integrable in the sense of Liouville.

Box 8. Separation of Variables

In the planar Euler three-body problem, Hamilton’s characteristic function W depends on the two ‘old’
elliptical coordinates ξ and η. The virtue of elliptical coordinates is that the time-independent HJ equation
can be solved by separating the dependence of W on ξ and η: W(ξ, η) = W1 (ξ) + W2 (η). Writing the energy
(17) in elliptical coordinates (18) and using pξ = W1 (ξ) and pη = W2 (η), the time-independent HJ equation
(23) becomes
W  (ξ)2 + W2 (η)2 − 2m f (μ1 + μ2 ) cosh ξ − 2m f (μ1 − μ2 ) cos η
E= 1 . (25)
2m f 2 (cosh2 ξ − cos2 η)
Rearranging,

W12 − 2Em f 2 cosh2 ξ − 2m f (μ1 + μ2 ) cosh ξ = −W22 − 2Em f 2 cos2 η + 2m f (μ1 − μ2 ) cos η. (26)

Since the LHS and RHS are functions only of ξ and η respectively, they must both be equal to a ‘separation
constant’ α. Thus, the HJ PDE separates into a pair of decoupled ODEs for W1 (ξ) and W2 (η). The latter
may be integrated using elliptic functions. Note that Whittaker’s constant w (24) may be expressed as
w = −2m f 2 E − α.

More generally, in the three-dimensional Euler three-body prob-

lem, the mass m can revolve (non-uniformly) about the line join-
ing the force centers (x-axis) so that its motion is no longer con-
fined to a plane. Nevertheless, the problem is exactly solvable
as the equations admit three independent constants of motion in
involution: energy, Whittaker’s constant and the x component of
angular momentum [2].

RESONANCE | January 2019 101

 GENERAL ARTICLE

7. Some Landmarks in the History of the Three-Body Prob-

lem

Developments that arose The importance of the three-body problem lies in part in the de-
from attempts to solve velopments that arose from attempts to solve it [7, 8]. These have
the three-body problem
had an impact all over astronomy, physics, and mathematics.
have had an impact
across the sciences. Can planets collide, be ejected from the solar system or suffer
significant deviations from their Keplerian orbits? This is the
question of the stability of the solar system. In the 18th cen-
tury, Pierre-Simon Laplace and J L Lagrange obtained the first
significant results on stability. They showed that to first order in
the ratio of planetary to solar masses (M p /MS ), there is no un-
bounded variation in the semi-major axes of the orbits, indicating
stability of the solar system. Siméon Denis Poisson extended this
Haretu’s results result to second order in M p /MS . However, in what came as a
essentially put an end to surprise, the Romanian Spiru Haretu (1878) overcame significant
the hope of proving the technical challenges to find secular terms (growing linearly and
stability or instability of
the solar system using a quadratically in time) in the semi-major axes at third order! This
perturbative approach. was an example of a perturbative expansion, where one expands a
physical quantity in powers of a small parameter (here the semi-
major axis was expanded in powers of M p /MS
1). Haretu’s
result however did not prove instability as the effects of his secu-
lar terms could cancel out (see Box 9 for a simple example). But
it effectively put an end to the hope of proving the stability or
instability of the solar system using such a perturbative approach.

The scale of Delaunay’s The development of Hamilton’s mechanics and its refinement in
hand calculations is the hands of Carl Jacobi was still fresh when the French dynami-
staggering: he applied a cal astronomer Charles Delaunay (1846) began the first extensive
succession of 505
canonical use of canonical transformations (see Box 6) in perturbation the-
transformations to a 7th ory [13]. He arrived at the EOM for m3 in Hamiltonian form
order perturbative using 3 pairs of canonically conjugate orbital variables (3 angular
treatment of the momentum components, the true anomaly, longitude of the as-
three-dimensional
elliptical restricted cending node and distance of the ascending node from perigee).
three-body problem. He obtained the latitude and longitude of the Moon in trigono-
metric series of about 450 terms with secular terms eliminated
(see Box 9). It wasn’t till 1970–71 that Delaunay’s heroic calcu-

102 RESONANCE | January 2019

 GENERAL ARTICLE

lations were checked and extended using computers at the Boeing

Scientific Laboratories [13]!

Box 9. Poincaré–Lindstedt Method

The Poincaré-Lindstedt method is an approach to finding series solutions to a system such as the anharmonic
oscillator ẍ + x + gx3 = 0, which for small g, is a perturbation of the harmonic oscillator m ẍ + kx = 0 with
mass m = 1 and spring constant k = 1. The latter admits the periodic solution x0 (t) = cos t with initial
conditions x(0) = 1, ẋ(0) = 0. For a small perturbation 0 < g
1, expanding x(t) = x0 (t) + gx1 (t) + · · · in
powers of g leads to a linearized equation for x1 (t)

ẍ1 + x1 + cos3 t = 0. (27)

However, the perturbative solution

 
1 3
x(t) = x0 + gx1 + O(g2 ) = cos t + g (cos 3t − cos t) − t sin t + O(g2 ) (28)
32 8

is unbounded due to the linearly growing secular term (−3/8)t sin t. This is unacceptable as the energy
E = 12 ẋ2 + 12 x2 + 14 gx4 must be conserved and the particle must oscillate between turning points of the
potential V = 12 x2 + g4 x4 . The Poincaré–Lindstedt method avoids this problem by looking for a series
solution of the form
x(t) = x0 (τ) + g x̃1 (τ) + · · · , (29)

where τ = ωt with ω = 1 + gω1 + · · · . The constants ω1 , ω2 , · · · are chosen to ensure that the coefficients of
the secular terms at order g, g2 , · · · vanish. In the case at hand, we have

x(t) = cos(t + gω1 t) + g x̃1 (t) + O(g2 ) = cos t + g x̃˜1 (t) + O(g2 ) where x̃˜1 (t) = x̃1 (t) − ω1 t sin t. (30)

x̃˜1 satisfies the same equation (27) as x1 did, leading to

 
1 3
x̃1 (t) = (cos 3t − cos t) + ω1 − t sin t. (31)
32 8

The choice ω1 = 3/8 ensures cancellation of the secular term at order g, leading to the approximate bounded
solution  

3 g
x(t) = cos t + gt + (cos 3t − cos t) + O g2 . (32)
8 32

The Swede Anders Lindstedt (1883) developed a systematic method

to approximate solutions to nonlinear ODEs when naive pertur-
bation series fail due to secular terms (see Box 9). The technique

RESONANCE | January 2019 103

 GENERAL ARTICLE

was further developed by Poincaré. Lindstedt assumed the series

to be generally convergent, but Poincaré soon showed that they
are divergent in most cases. Remarkably, nearly 70 years later,
Kolmogorov, Arnold and Moser showed that in many of the cases
where Poincaré’s arguments were inconclusive, the series are in
fact convergent, leading to the celebrated KAM theory of inte-
grable systems subject to small perturbations (see Box 10).

Box 10. Action-Angle Variables and Invariant Tori

Time evolution is particularly simple if all the generalized coordinates θ j are cyclic so that their conjugate
∂H
momenta I j are conserved: I˙j = − ∂θ j = 0. A Hamiltonian system with n degrees of freedom is integrable in

the sense of Liouville if it admits n canonically conjugate ({θ j , Ik } = δkj ) pairs of phase space variables (θ j , I j )
with all the θ j cyclic, so that its Hamiltonian depends only on the momenta, H = H(I). Here, the Kronecker
symbol δkj is equal to one for j = k and zero otherwise. It follows that the ‘angle’ variables θ j evolve linearly
in time (θ j (t) = θ j (0) + ω j t) while the momentum or ‘action’ variables I j are conserved. Here, ω j = θ̇ j = ∂H ∂I j
are n constant frequencies. Typically, the angle variables are periodic, so that the θ j parametrize circles.
The common level sets of the action variables I j = c j are, therefore, a family of tori that foliate the phase
space. Recall that a torus is a Cartesian product of circles. For instance, for one degree of freedom, θ1 labels
points on a circle S 1 while for 2 degrees of freedom, θ1 and θ2 label points on a 2-torus S 1 × S 1 which looks
like a vada or doughnut. Trajectories remain on a fixed torus determined by the initial conditions. Under
a sufficiently small and smooth perturbation H(I) + gH  (I, θ), Andrei Kolmogorov, Vladimir Arnold and
Jürgen Moser showed that some of these ‘invariant’ tori survive provided the frequencies ωi are sufficiently
‘non-resonant’ or ‘incommensurate’ (i.e., their integral linear combinations do not get ‘too small’).

George William Hill was motivated by discrepancies in the lunar

perigee calculations. His celebrated paper on this topic was pub-
lished in 1877 while working with Simon Newcomb at the Amer-
Simon Newcomb’s ican Ephemeris and Nautical Almanac . He found a new family of
project of revising all the periodic orbits in the circular restricted (Sun–Earth–Moon) three-
orbital data in the solar body problem by using a frame rotating with the Sun’s angular ve-
system established the
missing 42 in the 566 locity instead of that of the Moon. The solar perturbation to lunar
centennial precession of motion around the Earth results in differential equations with pe-
Mercury’s perihelion. riodic coefficients. He used Fourier series to convert these ODEs
This played an important to an infinite system of linear algebraic equations and developed a
role in validating
Einstein’s general theory theory of infinite determinants to solve them and obtain a rapidly
of relativity. converging series solution for lunar motion. He also discovered

104 RESONANCE | January 2019

 GENERAL ARTICLE

Figure 8. Poincaré’s classi-

fication of zeros of a vector
field (equilibrium or fixed
points) on a plane. (a) Cen-
ter is always stable with os-
cillatory motion nearby, (b,
c) nodes and foci (or spirals)
can be stable or unstable and
new ‘tight binary’ solutions to the three-body problem where two
(d) saddles are unstable ex-
nearby masses are in nearly circular orbits around their center of cept in one direction.
mass CM12 , while CM12 and the far away third mass in turn orbit
each other in nearly circular trajectories.
The French mining-engineer/mathematician/physicist Henri Poincaré
began by developing a qualitative theory of differential equations
from a global geometric viewpoint of the dynamics on phase space.
This included a classification of the types of equilibria (zeros of
vector fields) on the phase plane (centers, nodes, foci and saddles,
see Figure 8). His 1890 memoir on the three-body problem was
the prize-winning entry in King Oscar II’s 60th birthday compe-
tition (for a detailed account see [9]). He proved the divergence
of series solutions for the three-body problem developed by De-
launay, Hugo Gyldén and Lindstedt (in many cases) and cover-
gence of Hill’s infinite determinants. To investigate the stability
of three-body motions, Poincaré defined his ‘surfaces of section’
and a discrete-time dynamics via the ‘return map’ (see Figure 9a).
A Poincaré surface S is a two-dimensional surface in phase space
transversal to trajectories. The first return map takes a point q1
on S to q2 , which is the next intersection of the trajectory through
q1 with S . Given a saddle point p on a surface S , he defined its
stable and unstable spaces W s and Wu as points on S that tend to
p upon repeated forward or backward applications of the return
map (see Figure 9b). He initially assumed that W s and Wu on a
surface could not intersect and used this to argue that the solar Homoclinic refers to the
system is stable. This assumption turned out to be false, as he property of being
discovered with the help of Lars Phragmén. In fact, W s and Wu ‘inclined’ both forward
and backward in time to
can intersect transversally on a surface at a homoclinic point if the same point.
the state space of the underlying continuous dynamics is at least
three-dimensional. What is more, he showed that if there is one

RESONANCE | January 2019 105

 GENERAL ARTICLE

Figure 9. (a) A Poincaré

surface S transversal to a
trajectory is shown. The
trajectory through q1 on S
intersects S again at q2 .
The map taking q1 to q2 is
called Poincaré’s first return
map. (b) The saddle point
p and its stable and unsta-
ble spaces W s and Wu are
shown on a Poincaré surface
homoclinic point, then there must be infinitely many accumulat-
through p. The points at
ing at p. Moreover, W s and Wu fold and intersect in a very com-
which W s and Wu intersect
are called homoclinic points, plicated ‘homoclinic tangle’ in the vicinity of p. This was the first
e.g., h0 , h1 and h−1 . Points example of what we now call chaos. Chaos is usually manifested
on W s (or Wu ) remain on W s via an extreme sensitivity to initial conditions (exponentially di-
(or Wu ) under forward and verging trajectories with nearby initial conditions).
backward iterations of the
When two gravitating point masses collide, their relative speed
return map. Thus, the for-
diverges and solutions to the EOM become singular at the colli-
ward and backward images
sion time tc . More generally, a singularity occurs when either a
of a homoclinic point under
position or velocity diverges in finite time. The Frenchman Paul
the return map are also ho-
moclinic points. In the fig- Painlevé (1895) showed that binary and triple collisions are the
ure h0 is a homoclinic point only possible singularities in the three-body problem. However,
whose image is h1 on the he conjectured that non-collisional singularities (e.g., where the
segment [h0 , p] of W s . Thus, separation between a pair of bodies goes to infinity in finite time)
Wu must fold back to inter- are possible for four or more bodies.
sect W s at h1 . Similarly, if
It took nearly a century for this conjecture to be proven, cul-
h−1 is the backward image
minating in the work of Donald Saari and Zhihong Xia (1992)
of h0 on Wu , then W s must
and Joseph Gerver (1991) who found explicit examples of non-
fold back to intersect Wu at
h−1 . Further iterations pro- collisional singularities in the 5-body and 3n-body problems for n
duce an infinite number of sufficiently large [14]. In Xia’s example, a particle oscillates with
homoclinic points accumu- ever growing frequency and amplitude between two pairs of tight
lating at p. The first exam- binaries. The separation between the binaries diverges in finite
ple of a homoclinic tangle time, as does the velocity of the oscillating particle.
was discovered by Poincaré
The Italian mathematician Tulio Levi-Civita (1901) attempted to
in the restricted three-body
avoid singularities and thereby ‘regularize’ collisions in the three-
problem and is a signature of
body problem by a change of variables in the differential equa-
its chaotic nature.

106 RESONANCE | January 2019

 GENERAL ARTICLE

Figure 10. Equal-mass

zero-angular momentum
figure-8 choreography
solution to the three-body
problem. A choreography
is a periodic solution where
all masses traverse the same
tions. For example, the ODE for the one-dimensional Kepler
orbit separated equally in
problem ẍ = −k/x2 is singular at the collision point x = 0. This time.
singularity can be regularized by introducing a new coordinate
x = u2 and a reparametrized time ds = dt/u2 , which satisfy the Solutions which could
nonsingular oscillator equation u (s) = Eu/2 with conserved en- be smoothly extended
ergy E = (2u̇2 − k)/u2 . Such regularizations could shed light on beyond collision time
(e.g., the bodies
near-collisional trajectories (‘near misses’) provided the differen- elastically collide) were
tial equations remain physically valid. called regularizable.
Those that could not
The Finnish mathematician Karl Sundman (1912) began by show- were said to have an
ing that binary collisional singularities in the three-body problem essential or transcendent
could be regularized by a repararmetrization of time, s = |t1 −t|1/3 singularity at the
where t1 is the the binary collision time [15]. He used this to find collision.

a convergent series representation (in powers of s) of the gen-

eral solution of the three-body problem in the absence of triple
collisions. The possibility of such a convergent series had been
anticipated by Karl Weierstrass in proposing the three-body prob-
lem for King Oscar’s 60th birthday competition. However, Sund-
man’s series converges exceptionally slowly and has not been of The point-particle
much practical or qualitative use. approximation to the
equations for celestial
The advent of computers in the 20th century allowed numeri- bodies of non-zero size
cal investigations into the three-body (and more generally the n- breaks down due to tidal
body) problem. Such numerical simulations have made possible effects when the bodies
get very close.
the accurate placement of satellites in near-Earth orbits as well
as our missions to the Moon, Mars and the outer planets. They
have also facilitated theoretical explorations of the three-body Sundman showed that
problem including chaotic behavior, the possibility for ejection for non-zero angular
momentum, there are no
of one body at high velocity (seen in hypervelocity stars [16])
triple collisions in the
and quite remarkably, the discovery of new periodic solutions. three-body problem.
For instance, in 1993, Chris Moore discovered the zero angular
momentum figure-8 ‘choreography’ solution. It is a stable peri-

RESONANCE | January 2019 107

 GENERAL ARTICLE

odic solution with bodies of equal masses chasing each other on

an ∞-shaped trajectory while separated equally in time (see Fig-
ure 10). Alain Chenciner and Richard Montgomery [17] proved
its existence using an elegant geometric reformulation of Newto-
nian dynamics that relies on the variational principle of Euler and
Maupertuis.

8. Geometrization of Mechanics

 The optical path length Fermat’s principle in  optics states that light rays extremize the
 dτ is proportional
n(r) optical path length n(r(τ)) dτ where n(r) is the (position de-
to dτ/λ, which is the pendent) refractive index and τ a parameter along the path. The
geometric length in units
variational principle of Euler and Maupertuis (1744) is a mechan-
of the local wavelength
λ(r) = c/n(r)ν. Here, c ical analogue of Fermat’s principle [18].  q2It states that the curve
is the speed of light in that extremizes the abbreviated action q p · dq holding energy
1
vacuum and ν the E and the end-points q1 and q2 fixed has the same shape as the
constant frequency.
Newtonian trajectory. By contrast, Hamilton’s principle of ex-
tremal action (1835) states that a trajectory going from q1 at time
t1 to q2 at time t2 is a curve that extremizes the t action. Here, the
action is the integral of the Lagrangian S = t 2 L(q, q̇) dt. Typ-
1
ically, L = T − V is the difference between kinetic and potential
energies.

A metric mi j on an It is well-known that the trajectory of a free particle (i.e., subject

n-dimensional to no forces) moving on a plane is a straight line. Similarly, tra-
configuration space M is jectories of a free particle moving on the surface of a sphere are
an n × n matrix at each
point q ∈ M that great circles. More generally, trajectories of a free particle mov-
determines the square of ing on a curved space (Riemannian manifold M) are geodesics
the distance (curves that extremize length). Precisely, for a mechanical system

(ds2 = ni, j=1 mi j dqi dq j ) with configuration space M and Lagrangian L = 12 mi j (q)q̇i q̇ j , La-
from q to a nearby point ∂L
q + dq. We often grange’s equations ddtpi = ∂q i are equivalent to the geodesic equa-

suppress the summation tions with respect to the ‘kinetic metric’ mi j on M:

symbol and follow the
1

convention that repeated mi j q̈ j (t) = − m ji,k + mki, j − m jk,i q̇ j (t) q̇k (t). (33)
indices are summed 2
from 1 to n.
Here, mi j,k = ∂mi j /∂qk and pi = ∂∂Lq̇i = mi j q̇ j is the momentum
conjugate to coordinate qi . For instance, the kinetic metric (mrr =
m, mθθ = mr2 , mrθ = mθr = 0) for a free particle moving on a

108 RESONANCE | January 2019

 GENERAL ARTICLE

plane may be read off from the Lagrangian L = 12 m(ṙ2 + r2 θ̇2 ) in

polar coordinates, and the geodesic equations shown to reduce to
Lagrange’s equations of motion r̈ = rθ̇2 and d(mr2 θ̇)/dt = 0.

Box 11. Gaussian Curvature

Given a point p on a surface S embedded in three dimensions, a normal plane through p is one that is
orthogonal to the tangent plane at p. Each normal plane intersects S along a curve whose best quadratic
approximation at p is called its osculating circle. The principal radii of curvature R1,2 at p are the maximum
and minimum radii of osculating circles through p. The Gaussian curvature K(p) is defined as 1/R1 R2 and
is taken positive if the centers of the corresponding osculating circles lie on the same side of S and negative
otherwise.

Remarkably, the correspondence between trajectories and geodesics

continues to hold even in the presence of conservative forces de-
rived from a potential V. Indeed, trajectories of the Lagrangian
L = T − V = 12 mi j (q)q̇i q̇ j − V(q) are reparametrized2 geodesics 2 The shapes of trajectories

of the Jacobi–Maupertuis (JM) metric gi j = (E − V(q))mi j (q) on and geodesics coincide but the
Newtonian time along trajec-
M where E = T + V is the energy. This geometric formulation of
tories is not the same as
the Euler–Maupertuis principle (due to Jacobi) follows from the the arc-length parameter along
observation that the square of the metric line element geodesics.

1 dqk dql
ds2 = gi j dqi dq j = (E − V)mi j dqi dq j = mkl mi j dqi dq j
2 dt dt
1
2 1
= mi j q̇i dq j = (p · dq)2 , (34)
2 2

so that the extremization  of p · dq is equivalent to the extrem- Loosely, the potential in
ization of arc length ds. Loosely, the potential V(q) on the con- particle mechanics plays
the role of a variable
figuration space plays the role of an inhomogeneous refractive
refractive index in
index. Though trajectories and geodesics are the same curves, the optics.
Newtonian time t along trajectories is in general different from
the arc-length parameter s along geodesics. They are related by
√
ds
dt = 2(E − V) [19].
This geometric reformulation of classical dynamics allows us to Nearby geodesics on a
assign a local curvature to points on the configuration space. For surface of negative
instance, the Gaussian curvature K of a surface at a point (see curvature diverge
exponentially.

RESONANCE | January 2019 109

 GENERAL ARTICLE

Box 11) measures how nearby geodesics behave (see Figure 11),
they oscillate if K > 0 (as on a sphere), diverge exponentially if
K < 0 (as on a hyperboloid) and linearly separate if K = 0 (as on
a plane). Thus, the curvature of the Jacobi–Maupertuis metric de-
fined above furnishes information on the stability of trajectories.
Negativity of curvature leads to sensitive dependence on initial
conditions and can be a source of chaos.
In the planar Kepler problem, the Hamiltonian (5) in the CM
frame is

p2x + p2y α
H= − where α = GMm > 0 and r2 = x2 + y2 .
2m r
(35)
A compact Riemann The corresponding JM metric line element in polar coordinates

surface is a closed, is ds2 = m E + αr dr2 + r2 dθ2 . Its Gaussian curvature K =

oriented and bounded
−Eα/2m(α + Er)3 has a sign opposite to that of energy every-
surface such as a sphere,
a torus or the surface of a where. This reflects the divergence of nearby hyperbolic orbits
pretzel. The genus of and oscillation of nearby elliptical orbits. Despite negativity of
such a surface is the curvature and the consequent sensitivity to initial conditions, hy-
number of handles: zero
perbolic orbits in the Kepler problem are not chaotic: particles
for a sphere, one for a
torus and two or more simply fly off to infinity and trajectories are quite regular. On
for higher handle-bodies. the other hand, negativity of curvature without any scope for es-
Riemann surfaces with cape can lead to chaos. This happens with geodesic motion on a
genus two or more admit
compact Riemann surface with constant negative curvature: most
metrics with constant
negative curvature. trajectories are very irregular [2].

9. Geometric Approach to the Planar Three-Body Problem

We now sketch how the above geometrical framework may be

usefully applied to the three-body problem. The configuration
space of the planar three-body problem is the space of triangles
on the plane with masses at the vertices. It may be identified
with six-dimensional Euclidean space (R6 ) with the three planar
Jacobi vectors J1,2,3 (see (9) and Figure 2) furnishing coordinates
on it. A simultaneous translation of the position vectors of all
three bodies r1,2,3 → r1,2,3 + r0 is a symmetry of the Hamiltonian
H = T + V of Eqs. (10, 11) and of the Jacobi–Maupertuis metric

110 RESONANCE | January 2019

 GENERAL ARTICLE

Figure 11. Local behav-

ior of nearby geodesics on a
surface depends on the sign
of its Gaussian curvature K.
(a) Nearby geodesics on a
plane (K = 0) separate lin-
early. (b) Distance between
neighbouring geodesics on a
sphere (K > 0) oscillates.
(c) Geodesics on a hyper-
bolic surface (K < 0) devi-
ate exponentially.
PPFROOLVLRQ
(XOHUSRLQWV
(XOH S QWV Ș
FROOLVLRQȘ
Ș 
PPFROOLVLRQ
OOLVLRQ
RQ
Q (Ș ʌȟ ʌ Figure 12. ‘Pair of pants’
(ʌ
V\]\J\JUHDWFLUFOH
]\J\JUHDWFLUFOH metric on shape sphere and
/DJUDQJHSRLQW/ʌʌ
Lagrange, Euler and colli-
6\]\J\ ȟ /DJUDQJHSRLQW/ ʌʌ
 sion points. (a) The nega-
/LQHV FROOLVLRQʌ
/DJUDQJHSRLQWV 
FROOLVLRQʌʌ tively curved ‘pair of pants’
(Ș ʌ
metric on the shape sphere
PPFROOLVLRQ
S2 . (b) Locations of La-
D E grange, Euler and collision
points on a geometrically
unfaithful depiction of the
shape sphere S2 . The neg-

3 ative curvature of S2 is in-
ds2 = (E − V(J1 , J2 )) Ma |dJa |2 . (36) dicated in (a). Syzygies are
a=1 instantaneous configurations
where the three bodies are
This is encoded in the cyclicity of J3 . Quotienting by transla- collinear (eclipses).
tions allows us to define a center of mass configuration space
R4 (the space of centered triangles on the plane with masses at
the vertices)
⎛ with its ⎞quotient JM metric. Similarly, rotations
⎜⎜⎜cos θ − sin θ⎟⎟⎟
Ja → ⎜⎝ ⎟⎠ Ja for a = 1, 2, 3 are a symmetry of the
sin θ cos θ
metric, corresponding to rigid rotations of a triangle about a ver-
tical axis through the CM. The quotient of R4 by such rotations Translations and
is the shape space R3 , which is the space of congruence classes rotations are symmetries
of any central
of centered oriented triangles on the plane. Translations and rota-
inter-particle potential.
tions are symmetries of any central inter-particle potential, so the

RESONANCE | January 2019 111

 GENERAL ARTICLE

dynamics of the three-body problem in any such potential admits

a consistent reduction to geodesic dynamics on the shape space
R3 . Interestingly, for an inverse-square potential (as opposed to
the Newtonian ‘1/r’ potential)
 Gm m Gm1 m2 Gm2 m3 Gm3 m1
a b
V=− =− − −
|r − rb |
a<b a
2 |J1 |2 |J2 − μ1 J1 |2 |J2 + μ2 J1 |2
(37)
with μ1,2 = m1,2 /(m1 +m2 ), the zero-energy JM metric (12) is also
invariant under the scale transformation Ja → λJa for a = 1, 2
and 3 (see Box 12 for more on the inverse-square potential and
for why the zero-energy case is particularly interesting). This al-
lows us to further quotient the shape space R3 by scaling to get
the shape sphere S2 , which is the space of similarity classes of
3 Though scaling is not a sym- centered oriented triangles on the plane3 . Note that collision con-
metry for the Newtonian gravi- figurations are omitted from the configuration space and its quo-
tational potential, it is still use-
tients. Thus, the shape sphere is topologically a 2-sphere with the
ful to project the motion onto
the shape sphere. three binary collision points removed. In fact, with the JM metric,
the shape sphere looks like a ‘pair of pants’ (see Figure 12a).

Box 12. Inverse Square Potential and its Behavior Under Scaling

The inverse-square potential is somewhat simpler than the Newtonian one due to the behavior of the Hamil-
 
tonian H = a p2a /2ma − a<b Gma mb /|ra − rb |2 under scale transformations ra → λra and pa → λ−1 pa :
H(λr, λ−1 p) = λ−2 H(r, p) [6]. The infinitesimal version (λ ≈ 1) of this transformation is generated by the

dilatation operator D = a ra · pa via Poisson brackets {ra , D} = ra and {pa , D} = −pa . Here, the PB
between coordinates and momenta are {rai , pb j } = δab δi j where a, b label particles and i, j label Cartesian
components. In terms of PBs, time evolution of any quantity f is given by f˙ = { f, H}. It follows that
Ḋ = {D, H} = 2H, so scaling is a symmetry of the Hamiltonian (and D is conserved) only when the energy
vanishes. To examine long-time behavior we consider the evolution of the moment of inertia in the CM

frame ICM = a ma r2a whose time derivative may be expressed as I˙ = 2D. This leads to the Lagrange–
Jacobi identity I¨ = { I,
˙ H} = {2D, H} = 4E or I = I(0) + I(0)
˙ t + 2E t2 . Hence when E > 0, I → ∞ as t → ∞
so that bodies fly apart asymptotically. Similarly, when E < 0 they suffer a triple collision. When E = 0,
˙ controls asymptotic behavior leaving open the special case when E = 0 and I(0)
the sign of I(0) ˙ = 0. By
−2/3
contrast, for the Newtonian potential, the Hamiltonian transforms as H(λ r, λ p) = λ H(r, p) leading
1/3 2/3

to the Lagrange–Jacobi identity I¨ = 4E − 2V. This is however not adequate to determine the long-time
behavior of I when E < 0.

For equal masses and E = 0, the quotient JM metric on the shape

112 RESONANCE | January 2019

 GENERAL ARTICLE

sphere may be put in the form

ds2 = Gm3 h(η, ξ2 ) dη2 + sin2 2η dξ22 . (38)

Here, 0 ≤ 2η ≤ π and 0 ≤ 2ξ2 ≤ 2π are polar and azimuthal

angles on the shape sphere S2 (see Figure 12b). The function h
is invariant under the above translations, rotations and scalings
and therefore a function on S2 . It may be written as v1 + v2 + v3
where v1 = ICM /(m|r2 − r3 |2 ) etc., are proportional to the inter-
particle potentials [19]. As shown in Figure 12a, the shape The shape sphere is
sphere has three cylindrical horns that point toward the three col- negatively curved at all
points other than
lision points which lie at an infinite geodesic distance. Moreover,
Lagrange and collision
this equal-mass, zero-energy JM metric (38) has negative Gaus- points. It has three horns
sian curvature everywhere except at the Lagrange and collision that point towards the
points where it vanishes. This negativity of curvature implies collision points which lie
at an infinite distance.
geodesic instability (nearby geodesics deviate exponentially) as
well as the uniqueness of geodesic representatives in each ‘free’
homotopy class, when they exist. The latter property was used by
Montgomery [17] to establish uniqueness of the ‘figure-8’ solu-
tion (upto translation, rotation and scaling) for the inverse-square
potential. The negativity of curvature on the shape sphere for
equal masses extends to negativity of scalar curvature on the CM Scalar curvature is an
configuration space for both the inverse-square and Newtonian average of the Gaussian
gravitational potentials [19]. This could help to explain instabili- curvatures in the various
tangent planes through a
ties and chaos in the three-body problem. point

Acknowledgement

We thank A Thyagaraja for useful comments. This work was

supported in part by the Infosys Foundation.

Suggested Reading
[1] J Laskar, Is the Solar System Stable?, Progress in Mathematical Physics, Vol.66,
pp.239–270, 2013.
[2] M C Gutzwiller, Chaos in Classical and Quantum mechanics, Springer-Verlag,
New York, 1990.
[3] S Bodenmann, The 18th-century Battle Over Lunar Motion, Physics Today,
Vol.63, No.1, p.27, 2010.

RESONANCE | January 2019 113

 GENERAL ARTICLE

[4] H Goldstein, C P Poole and J L Safko, Classical Mechanics, 3rd Ed., Pearson
Education, 2011.
[5] L N Hand and J D Finch, Analytical Mechanics, Cambridge Univ. Press, 1998.
[6] S G Rajeev, Advanced Mechanics: From Euler’s Determinism to Arnold’s Chaos,
Oxford University Press, Oxford, 2013.
[7] F Diacu and P Holmes, Celestial Encounters: The Origins of Chaos and Stabil-
ity, Princeton University Press, New Jersey, 1996.
[8] Z E Musielak and B Quarles, The Three-body Problem, Reports on Progress in
Physics, Vol.77, No.6, p.065901, 2014, arXiv:1508.02312.
[9] J Barrow-Green, Poincaré and the Three Body Problem, Amer. Math. Soc.,
Providence, Rhode Island, 1997.
[10] E T Whittaker, A Treatise on the Analytical Dynamics of Particles & Rigid Bod-
ies, 2nd Ed., Cambridge University Press, Cambridge, 1917, Chapt. XIV and
p.283.
[11] K R Symon, Mechanics, 3rd Ed., Addison Wesley, Philippines 1971.
[12] N Mukunda, Sir William Rowan Hamilton, Resonance, 21, No.6, p.493, 2016.
[13] M C Gutzwiller, Moon-Earth-Sun: The Oldest Three-body Problem, Reviews
of Modern Physics, Vol.70, p.589, 1998.
[14] D G Saari and Z Xia, Off to Infinity in Finite Time, Notices of the AMS, Vol.42,
p.538, 1993,
Address for Correspondence [15] C L Siegel and J K Moser, Lectures on Celestial Mechanics, Springer-Verlag,
Govind S Krishnaswami1 Berlin, p.31, 1971.
Himalaya Senapati2 [16] W R Brown, Hypervelocity Stars in the Milky Way, Physics Today, Vol.69, No.6,
Chennai Mathematical Institute p.52, 2016.
SIPCOT IT Park [17] R Montgomery, A New Solution to the Three-body Problem, Notices of the
Siruseri 603 103 AMS, Vol.48, No.5, p.471, 2001.
India. [18] C Lanczos, The Variational Principles of Mechanics, 4th Ed., Dover, New York,
Email: govind@cmi.ac.in1 p.139, 1970.
himalay@cmi.ac.in2 [19] G S Krishnaswami and H Senapati, Curvature and Geodesic Instabilities in
a Geometrical Approach to the Planar Three-body Problem, J. Math. Phys.,
Vol.57, p.102901, 2016, arXiv:1606.05091.

114 RESONANCE | January 2019