Sundman’s Power Series Solution to the Three-Body Problem (Part I)

July 12, 2024 4 minute read Advanced

Summer research 2024: the three-body problem · Part 7 of 11

A weekly research log from my summer project at the University of Edinburgh, supervised by Prof. Jenni Smillie.

← Problem of Two Fixed Centres (Part II) · Sundman's Power Series Solution to the Three-Body Problem (Part II) →

In 1912 Karl F. Sundman found a convergent power-series solution to the three-body problem, valid whenever the total angular momentum is non-zero. In principle this is a general solution, but its convergence is so slow that it is of no practical use. This post outlines how it works.

Key points of Sundman’s paper

At a double (binary) collision, the distance to the third, non-colliding body stays bounded.
A triple collision can only occur if the total angular momentum (the “area integral”) is zero.
At a double collision, the velocities and accelerations of the colliding bodies tend to infinity.
The unit vector along the line joining the two colliding bodies tends to a limit.
The singularity of a collision between a given pair of bodies can be removed by introducing a new independent variable \(u\).
A different independent variable \(w\) removes all the double-collision singularities at once, allowing any number and combination of double collisions.
The coordinates of the three bodies can then be expanded as power series in \(w\).
These series converge for all time, although extremely slowly.

The regularising variable \(w\)

Where the singularities come from

The variable \(u\) regularises a collision between one particular pair of bodies, say \(m_0\) and \(m_1\). The equations of motion contain terms such as \(1/r_0\) and \(1/r_1\), where the \(r_i\) are mutual distances, so a collision involving a different pair still produces a singularity. We need a variable that treats all three pairs symmetrically.

Definition

Sundman’s variable \(w\) is defined by

\[dt= \Gamma\, dw, \qquad w=0 \text{ at } t=0,\]

where

\[\Gamma = \left(1-e^{-r_0/l}\right)\left(1-e^{-r_1/l}\right)\left(1-e^{-r_2/l}\right)\]

for a suitable constant \(l\), and \(r_0, r_1, r_2\) are the three mutual distances. The definition is symmetric in the three distances, so every double collision is regularised in the same way.

Convergence in a strip

After a long analysis (Yeomans gives a complete account), one finds that the coordinates are holomorphic functions of \(w\) in a strip around the real axis,

\[\lvert \operatorname{Im} w\rvert < \Omega,\]

for some \(\Omega > 0\): a strip of width \(2\Omega\).

Diagram of the strip of half-width Omega around the real w-axis in which the solution is holomorphic.

Mapping the strip to a disc

A power series converges in a disc, not a strip, so we map the strip conformally onto the unit disc.

Write the strip as \(S = \{s : \lvert\operatorname{Im} s\rvert < \beta\}\). The exponential map

\[g(s) = e^{\pi s / (2\beta)}\]

sends \(S\) onto the right half-plane \(\{\operatorname{Re} z > 0\}\). Points of the right half-plane are closer to \(1\) than to \(-1\), i.e. \(\lvert z-1\rvert < \lvert z+1\rvert\), so the Möbius transformation

\[h(z)=\frac{z-1}{z+1}\]

maps the right half-plane onto the unit disc. The composition \(h \circ g\) therefore maps \(S\) onto the unit disc. With \(\beta = \Omega\), the new variable \(\tau = h(g(w))\) satisfies

\[w = \frac{2 \Omega}{\pi} \log\frac{1+ \tau}{1- \tau},\]

and the coordinates become power series in \(\tau\) that converge for \(\lvert\tau\rvert < 1\). The whole real time axis, \(-\infty < w < \infty\), corresponds to the real interval \(-1 < \tau < 1\) (a good exercise to check).

This is also the source of the slow convergence: late times are squeezed towards the boundary \(\lvert\tau\rvert = 1\) of the disc of convergence, where power series converge very slowly.

When does it work?

The series exist whenever the total angular momentum is non-zero,

\[\mathbf{c} = m_0 (\mathbf{r}_0 \times \mathbf{v}_0) + m_1 (\mathbf{r}_1\times \mathbf{v}_1) + m_2 ( \mathbf{r}_2 \times \mathbf{v}_2) \neq \mathbf{0} .\]

That covers essentially every case, including double collisions, but not triple collisions, which require \(\mathbf{c}=\mathbf{0}\).

The catch is the convergence rate: Belorizky estimated that something like \(10^{8\,000\,000}\) terms would be needed for useful accuracy in a simple case. Part 8 goes through the estimate.

Summary

Sundman found an independent variable \(w\) in which every double collision is regular. The solution is analytic in a strip of width \(2\Omega\) around the real \(w\)-axis, which a conformal map sends to the unit disc in a new variable \(\tau\). The resulting power series in \(\tau\) solve the three-body problem for all time whenever the angular momentum is non-zero, but they converge far too slowly to be useful.

Progress

Alongside the reading, I have been simulating the Broucke–Hadjidemetriou–Hénon (BHH) orbits. Here is one, with its curve on the shape sphere seen from two angles:

Trajectories of a Broucke–Hadjidemetriou–Hénon periodic orbit in the plane.

The BHH orbit's shape curve on the shape sphere, first viewing angle.

The BHH orbit's shape curve on the shape sphere, second viewing angle.

Credits

This post is based on D. K. Yeomans’s exposition of Sundman’s regularisation of the three-body problem, NASA technical memorandum X-640-66-481 (1966), and on notes by Peter Senchyna.

Changes on 25 October 2026

What was wrong.

The map from the right half-plane to the unit disc was written \(h(z) = \lvert z-1\rvert / \lvert z+1\rvert\), which is a number, not a map. It is the Möbius transformation \(h(z) = (z-1)/(z+1)\).
The strip was written \(\operatorname{Im} s \le \beta\) instead of \(\lvert\operatorname{Im} s\rvert < \beta\), and the condition on \(w\) used \(\lvert w - \bar w\rvert\), which is twice \(\lvert\operatorname{Im} w\rvert\).
A reference link no longer worked.

What changed. All are fixed, and the reference now points to the copy on NASA’s Technical Reports Server.