The figure eight is one of the simplest three-body gravitational choreographies known. The purpose of this presentation is to display details of the figure eight as explicit functions of time as a step towards explicit parametrization of the orbit.

The Hamiltonian in what will be called direct variables is

$H=\frac{m}{2}({\stackrel{\xb7}{x}}_{1}^{2}+{\stackrel{\xb7}{y}}_{1}^{2}+{\stackrel{\xb7}{x}}_{2}^{2}+{\stackrel{\xb7}{y}}_{2}^{2}+{\stackrel{\xb7}{x}}_{3}^{2}+{\stackrel{\xb7}{y}}_{3}^{2})-\frac{k}{{r}_{12}}-\frac{k}{{r}_{23}}-\frac{k}{{r}_{31}}=E$

with equal masses and the radial differences

$\begin{array}{l}{r}_{12}=\sqrt{({x}_{1}-{x}_{2}{)}^{2}+({y}_{1}-{y}_{2}{)}^{2}}\\ {r}_{23}=\sqrt{({x}_{2}-{x}_{3}{)}^{2}+({y}_{2}-{y}_{3}{)}^{2}}\\ {r}_{31}=\sqrt{({x}_{3}-{x}_{1}{)}^{2}+({y}_{3}-{y}_{1}{)}^{2}}\end{array}$

The equations of motion in direct variables are

$\begin{array}{l}m{\stackrel{\xb7\xb7}{x}}_{1}=-k\frac{{x}_{1}-{x}_{2}}{{r}_{12}^{3}}+k\frac{{x}_{3}-{x}_{1}}{{r}_{31}^{3}}\\ m{\stackrel{\xb7\xb7}{x}}_{2}=-k\frac{{x}_{2}-{x}_{3}}{{r}_{23}^{3}}+k\frac{{x}_{1}-{x}_{2}}{{r}_{12}^{3}}\\ m{\stackrel{\xb7\xb7}{x}}_{3}=-k\frac{{x}_{3}-{x}_{1}}{{r}_{31}^{3}}+k\frac{{x}_{2}-{x}_{3}}{{r}_{23}^{3}}\\ \\ m{\stackrel{\xb7\xb7}{y}}_{1}=-k\frac{{y}_{1}-{y}_{2}}{{r}_{12}^{3}}+k\frac{{y}_{3}-{y}_{1}}{{r}_{31}^{3}}\\ m{\stackrel{\xb7\xb7}{y}}_{2}=-k\frac{{y}_{2}-{y}_{3}}{{r}_{23}^{3}}+k\frac{{y}_{1}-{y}_{2}}{{r}_{12}^{3}}\\ m{\stackrel{\xb7\xb7}{y}}_{3}=-k\frac{{y}_{3}-{y}_{1}}{{r}_{31}^{3}}+k\frac{{y}_{2}-{y}_{3}}{{r}_{23}^{3}}\end{array}$

Adding triples of equations gives the constraints

${\stackrel{\xb7\xb7}{x}}_{1}+{\stackrel{\xb7\xb7}{x}}_{2}+{\stackrel{\xb7\xb7}{x}}_{3}=0\phantom{\rule{7em}{0ex}}{\stackrel{\xb7\xb7}{y}}_{1}+{\stackrel{\xb7\xb7}{y}}_{2}+{\stackrel{\xb7\xb7}{y}}_{3}=0$

which can be translated into constraints on the center of mass

$\begin{array}{r}{x}_{1}+{x}_{2}+{x}_{3}=0\\ {\stackrel{\xb7}{x}}_{1}+{\stackrel{\xb7}{x}}_{2}+{\stackrel{\xb7}{x}}_{3}=0\end{array}\phantom{\rule{7em}{0ex}}\begin{array}{r}{y}_{1}+{y}_{2}+{y}_{3}=0\\ {\stackrel{\xb7}{y}}_{1}+{\stackrel{\xb7}{y}}_{2}+{\stackrel{\xb7}{y}}_{3}=0\end{array}$

by suitable choices of the initial conditions.

The appearance of only differences in the potential energy and hence right-hands sides of equations of motion motivates the introduction of difference variables:

$\begin{array}{l}{x}_{12}={x}_{1}-{x}_{2}\\ {x}_{23}={x}_{2}-{x}_{3}\\ {x}_{31}={x}_{3}-{x}_{1}\end{array}\phantom{\rule{4em}{0ex}}\begin{array}{l}{y}_{12}={y}_{1}-{y}_{2}\\ {y}_{23}={y}_{2}-{y}_{3}\\ {y}_{31}={y}_{3}-{y}_{1}\end{array}\phantom{\rule{4em}{0ex}}\begin{array}{l}{r}_{12}^{2}={x}_{12}^{2}+{y}_{12}^{2}\\ {r}_{23}^{2}={x}_{23}^{2}+{y}_{23}^{2}\\ {r}_{31}^{2}={x}_{31}^{2}+{y}_{31}^{2}\end{array}$

These differences will also be called dual variables for a reason to be made patent. The equations of motion in differences are somewhat simpler,

$\begin{array}{l}m{\stackrel{\xb7\xb7}{x}}_{12}=-3k\frac{{x}_{12}}{{r}_{12}^{3}}+X\\ m{\stackrel{\xb7\xb7}{x}}_{23}=-3k\frac{{x}_{23}}{{r}_{23}^{3}}+X\\ m{\stackrel{\xb7\xb7}{x}}_{31}=-3k\frac{{x}_{31}}{{r}_{31}^{3}}+X\end{array}\phantom{\rule{5em}{0ex}}\begin{array}{l}m{\stackrel{\xb7\xb7}{y}}_{12}=-3k\frac{{y}_{12}}{{r}_{12}^{3}}+Y\\ m{\stackrel{\xb7\xb7}{y}}_{23}=-3k\frac{{y}_{23}}{{r}_{23}^{3}}+Y\\ m{\stackrel{\xb7\xb7}{y}}_{31}=-3k\frac{{y}_{31}}{{r}_{31}^{3}}+Y\end{array}$

where the additional terms, identical for triples of equations, are

$X=k(\frac{{x}_{12}}{{r}_{12}^{3}}+\frac{{x}_{23}}{{r}_{23}^{3}}+\frac{{x}_{31}}{{r}_{31}^{3}})\phantom{\rule{7em}{0ex}}Y=k(\frac{{y}_{12}}{{r}_{12}^{3}}+\frac{{y}_{23}}{{r}_{23}^{3}}+\frac{{y}_{31}}{{r}_{31}^{3}})$

For a choreography these functions are periodic with a period equal to one-third of that of the choreography.

The constraints on the center of mass are carry over to the dual variables by subtracting each constraint from itself:

$\begin{array}{r}{x}_{12}+{x}_{23}+{x}_{31}=0\\ {\stackrel{\xb7}{x}}_{12}+{\stackrel{\xb7}{x}}_{22}+{\stackrel{\xb7}{x}}_{31}=0\end{array}\phantom{\rule{7em}{0ex}}\begin{array}{r}{y}_{12}+{y}_{23}+{y}_{31}=0\\ {\stackrel{\xb7}{y}}_{12}+{\stackrel{\xb7}{y}}_{23}+{\stackrel{\xb7}{y}}_{31}=0\end{array}$

To write the energy entirely in terms of differences, first note that with existing constraints

$({\stackrel{\xb7}{x}}_{1}+{\stackrel{\xb7}{x}}_{2}+{\stackrel{\xb7}{x}}_{3}{)}^{2}=0={\stackrel{\xb7}{x}}_{1}^{2}+{\stackrel{\xb7}{x}}_{2}^{2}+{\stackrel{\xb7}{x}}_{3}^{2}+2({\stackrel{\xb7}{x}}_{1}{\stackrel{\xb7}{x}}_{2}+{\stackrel{\xb7}{x}}_{1}{\stackrel{\xb7}{x}}_{3}+{\stackrel{\xb7}{x}}_{2}{\stackrel{\xb7}{x}}_{3})$

followed by

$\begin{array}{l}{\stackrel{\xb7}{x}}_{12}^{2}+{\stackrel{\xb7}{x}}_{23}^{2}+{\stackrel{\xb7}{x}}_{31}^{2}=2({\stackrel{\xb7}{x}}_{1}^{2}+{\stackrel{\xb7}{x}}_{2}^{2}+{\stackrel{\xb7}{x}}_{3}^{2})-2({\stackrel{\xb7}{x}}_{1}{\stackrel{\xb7}{x}}_{2}+{\stackrel{\xb7}{x}}_{1}{\stackrel{\xb7}{x}}_{3}+{\stackrel{\xb7}{x}}_{2}{\stackrel{\xb7}{x}}_{3})\\ \phantom{{\stackrel{\xb7}{x}}_{12}^{2}+{\stackrel{\xb7}{x}}_{23}^{2}+{\stackrel{\xb7}{x}}_{31}^{2}}=3({\stackrel{\xb7}{x}}_{1}^{2}+{\stackrel{\xb7}{x}}_{2}^{2}+{\stackrel{\xb7}{x}}_{3}^{2})\end{array}$

with corresponding expressions in *y*-variables. The energy can thus be written in dual variables as

$E=\frac{m}{6}({\stackrel{\xb7}{x}}_{12}^{2}+{\stackrel{\xb7}{y}}_{12}^{2}+{\stackrel{\xb7}{x}}_{23}^{2}+{\stackrel{\xb7}{y}}_{23}^{2}+{\stackrel{\xb7}{x}}_{31}^{2}+{\stackrel{\xb7}{y}}_{31}^{2})-\frac{k}{{r}_{12}}-\frac{k}{{r}_{23}}-\frac{k}{{r}_{31}}$

where the only change is a factor of one third on the kinetic energy. This concise form unfortunately cannot be used to generate accurate equations of motion.

While only the total energy is constant, it might be interesting to see the contributions from each body separately. In dual variables the definition of component energies is straightforward,

${E}_{ij}=\frac{m}{6}({\stackrel{\xb7}{x}}_{ij}^{2}+{\stackrel{\xb7}{y}}_{ij}^{2})-\frac{k}{{r}_{ij}}$

where indices are assigned cyclically. In direct variables the potential energy must be split equally over the bodies involved, so that again with cyclical index assignment one can define

${E}_{i}=\frac{m}{2}({\stackrel{\xb7}{x}}_{i}^{2}+{\stackrel{\xb7}{y}}_{i}^{2})-\frac{k}{2}(\frac{1}{{r}_{ij}}+\frac{1}{{r}_{ki}})$

The remaining constant of the motion is total angular momentum defined in direct variables by

$L=m({x}_{1}{\stackrel{\xb7}{y}}_{1}-{y}_{1}{\stackrel{\xb7}{x}}_{1}+{x}_{2}{\stackrel{\xb7}{y}}_{2}-{y}_{2}{\stackrel{\xb7}{x}}_{2}+{x}_{3}{\stackrel{\xb7}{y}}_{3}-{y}_{3}{\stackrel{\xb7}{x}}_{3})$

The can be written entirely in differences by noting

$\begin{array}{l}{x}_{12}{\stackrel{\xb7}{y}}_{12}+{x}_{23}{\stackrel{\xb7}{y}}_{23}+{x}_{31}{\stackrel{\xb7}{y}}_{31}\\ \phantom{\rule{4em}{0ex}}=3({x}_{1}{\stackrel{\xb7}{y}}_{1}+{x}_{2}{\stackrel{\xb7}{y}}_{2}+{x}_{3}{\stackrel{\xb7}{y}}_{3})-({x}_{1}+{x}_{2}+{x}_{3})({\stackrel{\xb7}{y}}_{1}+{\stackrel{\xb7}{y}}_{2}+{\stackrel{\xb7}{y}}_{3})\\ \phantom{\rule{4em}{0ex}}=3({x}_{1}{\stackrel{\xb7}{y}}_{1}+{x}_{2}{\stackrel{\xb7}{y}}_{2}+{x}_{3}{\stackrel{\xb7}{y}}_{3})\end{array}$

and a corresponding expression with dotted and undotted variables interchanged. The total angular momentum can thus be written in dual variables as

$L=\frac{m}{3}({x}_{12}{\stackrel{\xb7}{y}}_{12}-{y}_{12}{\stackrel{\xb7}{x}}_{12}+{x}_{23}{\stackrel{\xb7}{y}}_{23}-{y}_{23}{\stackrel{\xb7}{x}}_{23}+{x}_{31}{\stackrel{\xb7}{y}}_{31}-{y}_{31}{\stackrel{\xb7}{x}}_{31})$

where the only change is again a factor of one third. The definition of component angular momenta is straightforward in both cases:

${L}_{i}=m({x}_{i}{\stackrel{\xb7}{y}}_{i}-{y}_{i}{\stackrel{\xb7}{x}}_{i})\phantom{\rule{5em}{0ex}}{L}_{ij}=\frac{m}{3}({x}_{ij}{\stackrel{\xb7}{y}}_{ij}-{y}_{ij}{\stackrel{\xb7}{x}}_{ij})$

The equations of motion are integrated in the browser with a mass of one third and gravitational constant of one ninth. Details of initial conditions and their application are given in this presentation. The code denotes $\stackrel{\xb7}{x}$ with $u$ and $\stackrel{\xb7}{y}$ with $v$ for concision, and constraints on the center of mass are used to improve the performance of the process of integration.

The results of integration are assembled in multiple ways for the graphs that follow, with quantities available globally in the browser JavaScript console. Colors of direct variables are chosen to match orbit segments of this presentation. Colors of dual variables are assigned according to what would be their third cyclical index.

All functions are period in the basic time span of 1.67611892375928 or a rational division thereof. Functions for individual bodies are shifted by a third of this period relative to one another. Total energy is minus one half and total angular momentum is zero.

Plots are presented first for direct variables, followed by dual variables.

Parametric plot of direct coordinates $x$ and $y$ , and the source of the name of the choreography:

Direct coordinate variables $x$ as a function of time:

Direct coordinate variables $y$ as a function of time:

Direct radial variables $\sqrt{{x}^{2}+{y}^{2}}$ as a function of time:

Sum of direct radial variables as a function of time:

Direct coordinate angular variables ${tan}^{-1}\left(\frac{y}{x}\right)$ as a function of time:

Sum of direct coordinate angular variables as a function of time:

Parametric plot of direct velocity variables $\stackrel{\xb7}{x}$ and $\stackrel{\xb7}{y}$ :

Direct velocity variables $\stackrel{\xb7}{x}$ as a function of time:

Direct velocity variables $\stackrel{\xb7}{y}$ as a function of time:

Direct speed variables $\sqrt{{\stackrel{\xb7}{x}}^{2}+{\stackrel{\xb7}{y}}^{2}}$ as a function of time:

Sum of direct speed variables as a function of time:

Direct velocity angular variables ${tan}^{-1}\left(\frac{\stackrel{\xb7}{y}}{\stackrel{\xb7}{x}}\right)$ as a function of time:

Sum of direct velocity angular variables as a function of time:

Component energies in direct variables as a function of time:

Component angular momenta in direct variables as a function of time:

Parametric plot of direct phase space $x$ and $\stackrel{\xb7}{x}$ :

Parametric plot of direct phase space $y$ and $\stackrel{\xb7}{y}$ :

Parametric plot of dual coordinates $x$ and $y$ :

Dual coordinate variables $x$ as a function of time:

Dual coordinate variables $y$ as a function of time:

Dual radial variables $\sqrt{{x}^{2}+{y}^{2}}$ as a function of time:

Sum of dual radial variables as a function of time:

Dual coordinate angular variables ${tan}^{-1}\left(\frac{y}{x}\right)$ as a function of time:

Sum of dual coordinate angular variables as a function of time:

Parametric plot of dual velocity variables $\stackrel{\xb7}{x}$ and $\stackrel{\xb7}{y}$ :

Dual velocity variables $\stackrel{\xb7}{x}$ as a function of time:

Dual velocity variables $\stackrel{\xb7}{y}$ as a function of time:

Dual speed variables $\sqrt{{\stackrel{\xb7}{x}}^{2}+{\stackrel{\xb7}{y}}^{2}}$ as a function of time:

Sum of dual speed variables as a function of time:

Dual velocity angular variables ${tan}^{-1}\left(\frac{\stackrel{\xb7}{y}}{\stackrel{\xb7}{x}}\right)$ as a function of time:

Sum of dual velocity angular variables as a function of time:

Component energies in dual variables as a function of time:

Component angular momenta in dual variables as a function of time:

Parametric plot of dual phase space $x$ and $\stackrel{\xb7}{x}$ :

Parametric plot of dual phase space $y$ and $\stackrel{\xb7}{y}$ :

The function $X$ in dual variables as a function of time:

The function $Y$ in dual variables as a function of time:

The motivation for the term “dual” comes from the parametric plots above. In dual coordinate variables the closed figure eight becomes an open Gaussian-shaped curve. Conversely, the open velocity parametric plot in direct variables becomes a closed curve in dual variables. A variety of such relations are expected to appear between other corresponding variables.

Having all of the dynamic variables visualized allows one to speculate as to the exact functional forms of their temporal dependence as well as immediately dismiss some ideas. For example, the sum of dual radial variables appears almost constant when the origin is included but clearly is not exactly constant, since the integration of the equations of motion is good to ten significant digits. The actual behavior appears to be close to a cosine with small amplitude.

Several of the other functions appear to be either circular or Jacobi elliptic functions and in some case this can be graphically tested. The individual dual radial variables appear to be circular functions: plotting the right-hand sides of their equations of motion would in that case have the same appearance. First develop one such equation of motion from given quantities:

$\begin{array}{l}\hfill {r}_{12}^{2}={x}_{12}^{2}+{y}_{12}^{2}\phantom{\rule{5em}{0ex}}{\stackrel{\xb7}{r}}_{12}=\frac{{x}_{12}{\stackrel{\xb7}{x}}_{12}+{y}_{12}{\stackrel{\xb7}{y}}_{12}}{{r}_{12}}\hfill \\ \\ {\stackrel{\xb7\xb7}{r}}_{12}=\frac{{\stackrel{\xb7}{x}}_{12}^{2}+{x}_{12}{\stackrel{\xb7\xb7}{x}}_{12}+{\stackrel{\xb7}{y}}_{12}^{2}+{y}_{12}{\stackrel{\xb7\xb7}{y}}_{12}}{{r}_{12}}-\frac{({x}_{12}{\stackrel{\xb7}{x}}_{12}+{y}_{12}{\stackrel{\xb7}{y}}_{12}{)}^{2}}{{r}_{12}^{3}}\\ \phantom{{\stackrel{\xb7\xb7}{r}}_{12}}=\frac{{x}_{12}{\stackrel{\xb7\xb7}{x}}_{12}+{y}_{12}{\stackrel{\xb7\xb7}{y}}_{12}}{{r}_{12}}+\frac{({x}_{12}{\stackrel{\xb7}{y}}_{12}-{y}_{12}{\stackrel{\xb7}{x}}_{12}{)}^{2}}{{r}_{12}^{3}}\\ {\stackrel{\xb7\xb7}{r}}_{12}=-\frac{3k}{m{r}_{12}^{2}}+\frac{{x}_{12}X+{y}_{12}Y}{m{r}_{12}}+\frac{9{L}_{12}^{2}}{{m}^{2}{r}_{12}^{3}}\phantom{\rule{2em}{0ex}}\stackrel{?}{=}\phantom{\rule{2em}{0ex}}-{\omega}^{2}{r}_{12}\end{array}$

If this variable is indeed a simple circular function, then as indicated the right-hand side of the equation should be proportional to the variable itself apart from a constant displacement. Plotting the radial variable along with the scaled right-hand side in purple,

it is immediately clear that the simple looking function is not an exact circular. There remains the possibility of approximating its behavior and that of other variables using Jacobi elliptic functions in such a way that the approximation is always exact at extrema of the motion. This will be explored elsewhere.

*Uploaded 2018.10.15 — Updated 2018.10.18*
analyticphysics.com