The three-body problem

In the summer of 2024 I spent eight weeks on the gravitational three-body problem in the School of Physics and Astronomy at the University of Edinburgh, supervised by Prof. Jenni Smillie and funded by a Career Development Scholarship. The project had two aims: to reproduce and visualise important orbits, and to understand two landmark results, the proof that the figure-eight orbit exists and the convergence of Sundman’s series.

This page goes through the main ideas, with small simulations running in your browser. The full write-up is in the final report (PDF), and I kept a research log while the project ran. One part became a short student journal paper on Euler’s problem of two fixed centres.

1The problem

Three point masses attract each other by Newtonian gravity:

For two bodies this reduces to one body moving in a central force, and the orbits are conic sections: Kepler’s ellipses, parabolas and hyperbolas. With three bodies there is no such reduction. The conserved quantities (energy, momentum and angular momentum) are not enough to pin down the motion, and many motions are chaotic: two starting points that differ by a tiny amount end up far apart.

Throughout this page the three masses are equal and I use units with \(G = m = 1\).

2A short history

• 1687. Newton’s Principia treats the motion of the Moon under the pull of both the Earth and the Sun, the first three-body problem.
• 1767. Euler finds solutions in which the three bodies stay on a rotating straight line.
• 1772. Lagrange finds solutions in which they stay at the corners of a rotating equilateral triangle.
• 1887. Bruns shows that there are no further conserved quantities of a simple algebraic kind.
• 1890. Poincaré’s prize memoir for King Oscar II uncovers the extreme sensitivity to initial conditions that we now call chaos.
• 1912. Sundman finds a power series that solves the problem for all time (section 6).
• 1993. Moore finds the figure-eight orbit numerically.
• 2000. Chenciner and Montgomery prove that it exists (section 5).
• 2013. Šuvakov and Dmitrašinović find thirteen new families of periodic orbits, with names such as butterfly, moth and yin-yang.
• 2017. Li and Liao find more than six hundred more.

3Integrating without losing energy

With no formula, the orbits have to be computed numerically. The simplest method, Euler’s, moves the positions and velocities forward along their current rates of change. It is easy to write and it fails: the energy drifts steadily, so orbits slowly spiral outwards or inwards.

Gravity is a Hamiltonian system, with energy \(H = T(\mathbf p) + V(\mathbf q)\), and its time evolution preserves areas in phase space (it is a symplectic map). A symplectic integrator is a time step that preserves them too. Splitting each step into kicks, in which only the velocities change, and drifts, in which only the positions change, gives such steps automatically, because each piece is an exact Hamiltonian flow on its own:

One kick and one drift give symplectic Euler. With \(c = (\tfrac12, \tfrac12)\) and \(d = (1, 0)\) this is the velocity Verlet method. Ruth’s third-order scheme and Forest and Ruth’s fourth-order scheme use three and four stages with carefully chosen coefficients. A symplectic integrator does not conserve the energy exactly, but it exactly conserves a slightly modified Hamiltonian, so the energy error stays bounded instead of drifting. The classical fourth-order Runge–Kutta method is just as accurate per step but not symplectic, and its error keeps growing.

4Periodic orbits

Most three-body motions are chaotic, but special starting points give orbits that repeat exactly. A family studied by Šuvakov and Dmitrašinović starts with the three bodies on a line, at \((-1, 0)\), \((1, 0)\) and \((0, 0)\), with velocities \((a, b)\), \((a, b)\) and \((-2a, -2b)\), so that the total momentum and the angular momentum are both zero. Different values of \((a, b)\) give the figure-eight, butterflies, moths and yin-yangs. Šuvakov and Dmitrašinović found them by scanning the \((a, b)\) plane for starting points whose motion comes back close to where it began.

The last part of my project was a smaller search of the same kind. I varied \(a\) and \(b\) from 0 to 1, integrated each orbit with a fixed time step, and kept the ones that stayed bounded with an energy error below 1. Looking back, that test said more about the integrator than about the orbits: a fixed step fails during close encounters, so a large energy error mostly meant that the step could not follow one.

The map below redoes the search. Each of 401 × 401 starting velocities is integrated to \(t = 100\) with the Forest–Ruth method and a step that shrinks during close encounters,

where \(r_{\min}\) is the smallest distance between two of the bodies. For each orbit the map records whether a body escapes, which I count as moving more than 20 from the centre of mass. If none does, it records the closest return: the smallest distance, over time, between the state and the starting state, measured over all positions and velocities. Bright points come back almost exactly, so they are candidates for periodic orbits. Orbits whose energy error exceeded \(10^{-3}\) are marked as not resolved. These make up 4 % of the grid, most of them at small \(b\), where the bodies start almost in a line and pass very close to each other.

Outside the dashed arc the total energy, \(3(a^2+b^2) - \tfrac52\), is positive, so the three bodies cannot all stay together, and 99.8 % of the orbits there lose a body within 100 time units. Inside the arc a third of them do. Fifteen of the sixteen known orbits fall on grid points that stay bound and come back to within 0.06 of their starting state. The goggles orbit lies in such a thin strip that its nearest grid point escapes. The figure-eight sits in a broad bright patch: it is stable, so starting points near it stay near it. Most of the others are unstable, and in the live view a small error grows until the orbit breaks apart. The goggles orbit loses a body after about 125 time units, and butterfly I after about 140.

The space–time view is how Moore found the figure-eight in 1993. Followed through time, three bodies in a plane trace a braid, and he looked, for each braid, for the orbit with the least action.

The slider pushes bodies 1 and 2 out of the plane in opposite directions. This gives the system angular momentum about an axis that lies in the plane, so the plane itself slowly turns over, as the space view shows. The push changes the side lengths of the triangle only at second order, so it changes the shape much less. For the figure-eight, a push of 0.01 keeps the path on the shape sphere within 0.05 of the unpushed one for about 140 time units; an in-plane push of the same size moves it that far within two. The figure-eight survives pushes of up to 0.1 for at least 300 time units. After a push of 0.2, one body escapes after about 120.

The varying step also handles the very close encounters of the dragonfly and the yarn, which a fixed step of \(10^{-4}\) could not follow.

5The shape sphere and the figure-eight

Much of the motion of three bodies does not change the shape of the triangle they form: moving it around, rotating it, or scaling it up and down. Removing these leaves just the shape, and the possible shapes form a sphere, the shape sphere. Writing the positions as complex numbers \(x_1, x_2, x_3\) with their centre of mass at the origin, the Jacobi coordinates

map to a point on the unit sphere,

Collinear triangles make up the equator \(u_3 = 0\), with the three binary collisions on it, and the two equilateral triangles, one for each orientation, are the poles. On the sphere the figure-eight crosses the equator at the Euler configurations, where one body sits exactly midway between the other two, and passes over isosceles shapes in between.

Chenciner and Montgomery’s proof that the figure-eight exists works on this sphere:

1. Minimise the action \(A = \int_0^{T} \left(\tfrac12 K + U\right) dt\), where \(K\) is the kinetic term and \(U\) the sum of the inverse distances, over paths from an Euler configuration to an isosceles configuration. The minimiser is one twelfth of the orbit.
2. Show that the minimiser has no collisions: find a lower bound for the action of any path that contains a collision, and a collision-free test path whose action is smaller.
3. Use the symmetries of the equal-mass problem to join twelve copies of this arc smoothly into one closed curve.
4. Check that the result solves the equations of motion and traces a figure eight in the plane.

The figure-eight post in the research log goes through the steps in more detail.

6Sundman’s series

In 1912 Karl Sundman proved that as long as the total angular momentum is not zero, so that all three bodies cannot collide at once, the motion can be written as a power series that converges for all time. He removed the singularities of two-body collisions with a new time variable, then mapped the region where the solution is analytic onto the unit disc.

In principle this is a general solution. In practice it is useless. For a simple case, Lagrange’s equilateral solution, my estimate in the research log is that following the motion to one decimal place, even for a short time, needs about \(10^{8\,000\,000}\) terms. The research log goes through the construction in part 1 and the estimate in part 2.

7Euler’s problem of two fixed centres

If two of the bodies are held in place and only the third moves, the problem becomes Euler’s problem of two fixed centres. Unlike the full problem it is integrable: in prolate spheroidal coordinates the motion separates, and besides the energy there is a second conserved quantity. Its orbits fall into three types: satellite orbits around one centre, planetary orbits around both, and lemniscate orbits that pass between them.

Reproducing these orbits with the symplectic integrator, and checking them against the two conserved quantities, became a short student journal paper, an early and exploratory piece of work. The research log covers the coordinates and the orbit types in part I, part II and part III.

8The visualisation, with Allison Lau

Allison Lau (University of Toronto) built an interactive website on top of the project. It plays sixteen periodic orbits, computed with the project’s symplectic integrators at orders one to four, so you can see how the lower-order methods distort an orbit, alongside their paths on the shape sphere.

Open the visualisation Source code on GitHub

9References

1. J. Barrow-Green, Poincaré and the Three Body Problem (American Mathematical Society, 1997).
2. A. Chenciner and R. Montgomery, “A remarkable periodic solution of the three-body problem in the case of equal masses,” Ann. Math. 152, 881 (2000).
3. E. Forest and R. D. Ruth, “Fourth-order symplectic integration,” Physica D 43, 105 (1990).
4. E. Hairer, C. Lubich and G. Wanner, Geometric Numerical Integration, 2nd ed. (Springer, 2006).
5. X. Li and S. Liao, “More than six hundred new families of Newtonian periodic planar collisionless three-body orbits,” Sci. China Phys. Mech. Astron. 60, 129511 (2017).
6. C. Moore, “Braids in classical dynamics,” Phys. Rev. Lett. 70, 3675 (1993).
7. D. Ó'Mathúna, Integrable Systems in Celestial Mechanics (Birkhäuser, 2008).
8. R. D. Ruth, “A canonical integration technique,” IEEE Trans. Nucl. Sci. 30, 2669 (1983).
9. K. F. Sundman, “Mémoire sur le problème des trois corps,” Acta Math. 36, 105 (1912).
10. M. Šuvakov and V. Dmitrašinović, “Three classes of Newtonian three-body planar periodic orbits,” Phys. Rev. Lett. 110, 114301 (2013).
11. H. Yip and J. M. Smillie, “An analysis of the restricted Euler problem using symplectic integrators,” Edinburgh Student Journal of Science 1(3) (2025).