Problem of Two Fixed Centres (Part I)
From the archive (2024). Corrected in October 2026; see the note at the end.
Summer research 2024: the three-body problem · Part 4 of 11
A weekly research log from my summer project at the University of Edinburgh, supervised by Prof. Jenni Smillie.
← Ellipses and the Two-Body Problem · Understanding and Visualising the Shape Sphere →
This post introduces the problem of two fixed centres, also known as Euler’s three-body problem: a small body moving under the gravity of two masses held fixed in place. Unlike the general three-body problem, it is integrable.
Momentum conservation fixed
In part 2 I noted that the total momentum in my simulation was drifting. The cause turned out to be a bug: an accumulator for the momentum components was not reset to zero at each time step. With the fix, the remaining error is at the level of floating-point round-off, about \(10^{-15}\):
Set-up
The two fixed masses \(m_\pm\) sit on the \(z\)-axis at \(z = \pm b\), and the third body moves in the \(x\)–\(z\) plane.
Simulation
To simulate the problem, I held two of the masses fixed and let the third move freely:
The Lagrangian
The natural coordinates are prolate spheroidal coordinates \((R, \sigma)\), which play the role of polar coordinates but with two foci:
\[x = \pm \sqrt{R^2-b^2}\, \sin \sigma, \qquad z = R \cos \sigma .\]Writing \(\mathbf{r}(R, \sigma) = (x, z)\), the partial derivatives are
\[\frac{\partial x}{\partial R} = \pm \frac{R \sin \sigma}{\sqrt{R^2 - b^2}}, \quad \frac{\partial z}{\partial R} = \cos \sigma, \qquad \frac{\partial x}{\partial \sigma} = \pm \sqrt{R^2 - b^2}\, \cos \sigma, \quad \frac{\partial z}{\partial \sigma} = -R \sin \sigma,\]so (taking the \(+\) sign)
\[\mathbf{r}_R = \frac{\partial \mathbf{r}}{\partial R} = \begin{pmatrix} \frac{R \sin \sigma}{\sqrt{R^2 - b^2}} \\ \cos \sigma \end{pmatrix}, \qquad \mathbf{r}_\sigma = \frac{\partial \mathbf{r}}{\partial \sigma} = \begin{pmatrix} \sqrt{R^2 - b^2}\, \cos \sigma \\ -R \sin \sigma \end{pmatrix} .\]The metric components \(g_{RR}= \mathbf{r}_R \cdot \mathbf{r}_R\), \(g_{R\sigma} = \mathbf{r}_R\cdot\mathbf{r}_\sigma\) and \(g_{\sigma\sigma} = \mathbf{r}_\sigma\cdot\mathbf{r}_\sigma\) come out as
\[g_{RR}=\frac{R^2-b^2 \cos ^2 \sigma}{R^2-b^2}, \qquad g_{R \sigma}=0, \qquad g_{\sigma \sigma}=R^2-b^2 \cos ^2 \sigma .\]The distances to the two centres are \(R \mp b\cos\sigma\), so the potential energy per unit mass is
\[V=-\mu\, \frac{R+\beta b \cos \sigma}{R^2-b^2 \cos ^2 \sigma}, \qquad \mu = G(m_+ + m_-), \qquad \beta=\frac{m_{+}-m_{-}}{m_{+}+m_{-}} .\]The Lagrangian per unit mass is therefore
\[\mathcal{L} = T-V =\frac{1}{2}\, \frac{R^2-b^2 \cos ^2 \sigma}{R^2-b^2}\, \dot{R}^2+\frac{1}{2}\left(R^2-b^2 \cos ^2 \sigma\right) \dot{\sigma}^2+\mu\, \frac{R+\beta b \cos \sigma}{R^2-b^2 \cos ^2 \sigma} .\]Jacobi elliptic functions
The equations of motion separate in these coordinates. Following the book, the angular motion reduces (after a change of independent variable) to an equation of the form
\[S^{\prime 2}=\left(1-S^2\right)\left(1-k^2 S^2\right),\]whose solution is
\[S=\operatorname{sn}(f+\omega \mid k).\]Here \(\operatorname{sn}\) is one of the Jacobi elliptic functions. Trigonometric functions are defined with reference to a circle, and the Jacobi elliptic functions generalise them to other conic sections, in particular the ellipse.
To verify the solution, use
\[\frac{d}{du} \operatorname{sn}(u \mid k) = \operatorname{cn}(u \mid k) \operatorname{dn}(u \mid k)\]together with the identities
\[\operatorname{sn}^2(u \mid k) + \operatorname{cn}^2(u \mid k) = 1, \qquad \operatorname{dn}^2(u \mid k) + k^2 \operatorname{sn}^2(u \mid k) = 1 .\]With \(u = f + \omega\), these give \(\operatorname{cn}^2 u = 1 - S^2\) and \(\operatorname{dn}^2 u = 1 - k^2 S^2\), so
\[S'^2 = \operatorname{cn}^2(u \mid k)\, \operatorname{dn}^2(u \mid k) = (1 - S^2)(1 - k^2 S^2) .\]Progress
At the end of the first week, the simulation code is in good shape, ahead of my plan, and I have started thinking about the structure of the final report.
Credits
This post is based on the book Integrable Systems in Celestial Mechanics by Diarmuid Ó’Mathúna. Parts of the material on Jacobi elliptic functions are from Wikipedia.
Changes on 25 October 2026
What was wrong. Nothing in the results.
What changed. The Jacobi elliptic function identities behind the last step are shown, the modulus is written \(k\) instead of \(k_{01}\), and the coordinates are introduced more carefully.