Finally On An Ethereal Orrery


Over the course of the year, my fellow students and I have been experimenting with an ethereal orrery which models the motion of heavenly bodies using nought but Sir N-----'s laws of gravitation and motion[1][2][3]. Whilst the consequences of those laws are not generally subject to solution by mathematical reckoning, we were able to approximate them with a scheme that admitted errors of the order of the sixth power of the steps in time by which we advanced the positions of those bodies[4][5].
We have thus far employed it to model the solar system itself, uniformly distributed bodies of matter and the accretion of bodies that are close to Earth's orbit about the Sun. Whilst we were most satisfied by its behaviour, I should now like to report upon an altogether more surprising consequence of its engine's action.

Sir N-----'s laws imply that a particle \(P_i\), having a mass of \(m_i\) and a vector position of \(\mathbf{x}_i\), will be subjected by a particle \(P_j\), with a mass of \(m_j\) and a position of \(\mathbf{x}_j\), to a force of
\[ \mathbf{F}_{i,j} = \frac{G m_i m_j}{\left|\mathbf{x}_j-\mathbf{x}_i\right|^3} \times \left(\mathbf{x}_j-\mathbf{x}_i\right)\\ \]
for some constant value \(G\) and where the vertical bars stand for the magnitude, or length, of the vector lying between them, which in this case will be the distance between \(P_i\) and \(P_j\). Furthermore, \(P_i\) will have a vector acceleration toward \(P_j\) of
\[ \mathbf{a}_{i,j} = \frac{\alpha \mathbf{F}_{i,j}}{m_i} = \frac{\alpha \,G m_j}{\left|\mathbf{x}_j-\mathbf{x}_i\right|^3} \times \left(\mathbf{x}_j-\mathbf{x}_i\right) \]
for some constant \(\alpha\), and, since such accelerations accumulate additively, its total acceleration should it be one of several particles would be
\[ \mathbf{a}_i = \sum_{j \neq i} \mathbf{a}_{i,j} = \alpha \, G \sum_{j \neq i} \frac{m_j}{\left|\mathbf{x}_j - \mathbf{x}_i\right|^3} \times \left(\mathbf{x}_j - \mathbf{x}_i\right) \]
where \(\sum\) is the summation sign.

We fashioned our orrery with the units of terameter, hellagram, megasecond and yottanewton given by
\[ \begin{align*} 1 \, \mathrm{Tm} &= 10^{12} \, \mathrm{m}\\ 1 \, \mathrm{Hg} &= 10^{27} \, \mathrm{g} = 10^{24} \, \mathrm{kg}\\ 1 \, \mathrm{Ms} &= 10^6 \, \mathrm{s}\\ 1 \, \mathrm{YN} &= 10^{24} \, \mathrm{N} \end{align*} \]
for distance, mass, time and force respectively, under which the constants \(G\) and \(\alpha\) are given by
\[ \begin{align*} G &\approx 6.67408 \times 10^{-11}\\ \alpha &= 1 \end{align*} \]
where \(\approx\) is read as approximately equal to.

Listing 1 shows the implementation of our ethereal orrery's mathematical engine, with accelerations figuring the accelerations of the particles according to Sir N-----'s laws, and substep and step applying the frankly incomprehensible scheme to reckon their consequent motions.

Listing 1: Our Orrery's Engine
var G = 6.67408e-11;

function accelerations(p) {
 var n = p.length;
 var i, j, dx, dy, d2, g, ax, ay;

 for(i=0;i<n;++i) {
  p[i].ax = 0;
  p[i].ay = 0;

 for(i=1;i<n;++i) {
  for(j=0;j<i;++j) {
   dx = p[j].x-p[i].x;     dy = p[j].y-p[i].y;
   d2 = dx*dx+dy*dy;       g  = G/(d2*Math.sqrt(d2));
   ax = g*(p[j].x-p[i].x); ay = g*(p[j].y-p[i].y);

   p[i].ax += p[j].m*ax;   p[i].ay += p[j].m*ay;
   p[j].ax -= p[i].m*ax;   p[j].ay -= p[i].m*ay;

function substep(p, dt) {
 var n = p.length;
 var i;

 for(i=0;i<n;++i) {
  p[i].x += p[i].vx*dt + 0.5*p[i].ax*dt*dt;
  p[i].y += p[i].vy*dt + 0.5*p[i].ay*dt*dt;

  p[i].vx += 0.5*p[i].ax*dt;
  p[i].vy += 0.5*p[i].ay*dt;


 for(i=0;i<n;++i) {
  p[i].vx += 0.5*p[i].ax*dt;
  p[i].vy += 0.5*p[i].ay*dt;

var C = [0.784513610477560, 0.235573213359357, -1.17767998417887, 1.31518632068391];

function step(p, dt) {
 var i;
 for(i=0;i<7;++i) substep(p, dt*C[Math.min(i, 6-i)]);

Deck 1 applies it to model the solar system, plotting the planets upon a quadratic scale of distance from the Sun, Sun in orange, Mercury in black, Venus in grey, Earth in green, Mars in red, Jupiter in light brown, Saturn in reddish brown, Uranus in cyan and Neptune in blue, the operation of which you may bring to a conclusion with the halt button.

Deck 1: Our Ethereal Orrery

You will no doubt recall that we had added a facility for determining when bodies collide so that they shouldn't approach infinite accelerations as the distances between them approach zero. During our final investigation we were not especially concerned with the outcomes of such events, but rather whether they might have happened at all, and so we removed it for the sake of economy of calculation.


Specifically, our interest had been drawn to periodic motions of point masses, being particles with positive mass but zero volume, that are quite unlike those of the solar system which are largely governed by the tremendous mass of the Sun.

Figure 1: A Regular Polygon Of Particles
For example, let us suppose that \(n\) point masses are located at the vertices of a regular polygon centred upon the origin, such as those illustrated by figure 1, and rotating about it at a constant rate so that angle between the line from the origin to the \(j^\mathrm{th}\) particle and the positive \(x\) axis and its first and second derivatives with respect to time are given by
\[ \begin{align*} \theta_j &= \theta_0 + \frac{2\pi j}{n}\\ \frac{\mathrm{d}}{\mathrm{d}t} \theta_j &= \theta^\prime\\ \frac{\mathrm{d}^2}{\mathrm{d}t^2} \theta_j &= \theta^{\prime\prime} = 0 \end{align*} \]
If that polygon has a radius of \(r\) then the Cartesian coordinates of the \(j^\mathrm{th}\) particle should be
\[ \mathbf{x}_j = \begin{pmatrix} r \cos \theta_j\\ r \sin \theta_j \end{pmatrix} \]
its velocity should be
\[ \mathbf{v}_j = \frac{\mathrm{d}}{\mathrm{d}t} \mathbf{x}_j = \begin{pmatrix} \frac{\mathrm{d}}{\mathrm{d}t} r \cos \theta_j\\ \frac{\mathrm{d}}{\mathrm{d}t} r \sin \theta_j \end{pmatrix} = \begin{pmatrix} -r \theta^\prime \times \sin \theta_j\\ \phantom{-}r \theta^\prime \times \cos \theta_j \end{pmatrix} \]
and its acceleration should be
\[ \mathbf{a}_j = \frac{\mathrm{d}}{\mathrm{d}t} \mathbf{v}_j = \begin{pmatrix} -r \theta^{\prime\prime} \times \sin \theta_j -r {\theta^\prime}^2 \times \cos \theta_j \\ \phantom{-}r \theta^{\prime\prime} \times \cos \theta_j -r {\theta^\prime}^2 \times \sin \theta_j \end{pmatrix} = \begin{pmatrix} -r {\theta^\prime}^2 \times \cos \theta_j\\ -r {\theta^\prime}^2 \times \sin \theta_j \end{pmatrix} \]
Let us further suppose that the first particle lies upon the positive \(x\) axis so that \(\theta_0\) equals zero and hence
\[ \begin{align*} \mathbf{x}_0 &= \begin{pmatrix} r \cos 0\\ r \sin 0 \end{pmatrix} = \begin{pmatrix} r\\ 0 \end{pmatrix}\\ \mathbf{v}_0 &= \begin{pmatrix} -r \theta^\prime \times \sin 0\\ \phantom{-}r \theta^\prime \times \cos 0 \end{pmatrix} = \begin{pmatrix} 0\\ r \theta^\prime \end{pmatrix}\\ \mathbf{a}_0 &= \begin{pmatrix} -r {\theta^\prime}^2 \times \cos 0\\ -r {\theta^\prime}^2 \times \sin 0 \end{pmatrix} = \begin{pmatrix} -r {\theta^\prime}^2\\ 0 \end{pmatrix} \end{align*} \]
Now the vector from that particle to the \(j^\mathrm{th}\) and its magnitude are given by
\[ \begin{align*} \mathbf{x}_j - \mathbf{x}_0 &= \begin{pmatrix} r \cos \theta_j - r\\ r \sin \theta_j \end{pmatrix} = r \begin{pmatrix} \cos \theta_j - 1\\ \sin \theta_j \end{pmatrix} \\ \left|\mathbf{x}_j - \mathbf{x}_0\right| &= r \left(\left(\cos \theta_j - 1\right)^2 + \sin^2 \theta_j\right)^\frac12 = r \left(\cos^2 \theta_j - 2 \cos \theta_j + 1 + \sin^2 \theta_j\right)^\frac12 = r \left(2 - 2 \cos \theta_j\right)^\frac12 \end{align*} \]
since the square of the cosine plus the square of the sign of the same angle always equals one, and Sir N-----'s laws consequently imply that its acceleration toward the \(j^\mathrm{th}\) equals
\[ \begin{align*} \mathbf{a}_{0,j} &= \frac{Gm}{\left|\mathbf{x}_j - \mathbf{x}_0\right|^3} \times \left(\mathbf{x}_j - \mathbf{x}_0\right) = \frac{Gm}{r^3 \left(2 - 2 \cos \theta_j\right)^\frac32} \times r \begin{pmatrix} \cos \theta_j - 1\\ \sin \theta_j \end{pmatrix}\\ &= -\frac{Gm}{r^2} \begin{pmatrix} \dfrac{1-\cos \theta_j}{\left(2 - 2 \cos \theta_j\right)^\frac32}\\ -\dfrac{\sin \theta_j}{\left(2 - 2 \cos \theta_j\right)^\frac32} \end{pmatrix} = -\frac{Gm}{r^2} \begin{pmatrix} \dfrac{1}{2\left(2 - 2 \cos \theta_j\right)^\frac12}\\ -\dfrac{\sin \theta_j}{\left(2 - 2 \cos \theta_j\right)^\frac32} \end{pmatrix} \end{align*} \]
Since the other particles lie equally above and below the \(x\) axis, their forces upon the first perpendicular to it, which is to say along the \(y\) axis, necessarily cancel each other out and we therefore need only consider its acceleration along the former
\[ a_{0,j} = -\frac{Gm}{2^\frac32 r^2} \times \frac{1}{\left(1 - \cos \theta_j\right)^\frac12} \]
Rearranging the trigonometrical half-angle formula
\[ \cos \theta = \frac{1 - \tan^2 \frac12\theta}{1 + \tan^2 \frac12\theta} \]
\[ \begin{align*} \cos \theta &= \frac{1 - \frac{\sin^2 \frac12\theta}{\cos^2 \frac12\theta}}{1 + \frac{\sin^2 \frac12\theta}{\cos^2 \frac12\theta}} = \frac{\cos^2 \frac12\theta - \sin^2 \frac12\theta}{\cos^2 \frac12\theta + \sin^2 \frac12\theta} = \left(1 - \sin^2 \tfrac12\theta\right) - \sin^2 \tfrac12\theta = 1 - 2\sin^2 \tfrac12\theta\\ 1 - \cos \theta &= 2\sin^2 \tfrac12\theta \end{align*} \]
implying that
\[ a_{0,j} = -\frac{Gm}{2^\frac32 r^2} \times \frac{1}{\left(2 \sin^2 \frac{\theta_j}{2}\right)^\frac12} = -\frac{Gm}{4 r^2 \sin \frac{\theta_j}{2}} = -\frac{Gm}{4 r^2 \sin \frac{\pi j}{n}} \]
The first particle's acceleration along the \(x\) axis is trivially given by
\[ a_0 = \sum_{j=1}^{n-1} a_{0,j} \]
from which we may deduce that
\[ \begin{align*} -r {\theta^\prime}^2 &= -\frac{Gm}{4 r^2}\sum_{j=1}^{n-1} \frac{1}{\sin \frac{\pi j}{n}}\\ \theta^\prime &= \left(\frac{Gm}{4 r^3}\sum_{j=1}^{n-1} \frac{1}{\sin \frac{\pi j}{n}}\right)^\frac12 \end{align*} \]
Finally, defining
\[ \Sigma = \sum_{j=1}^{n-1} \frac{1}{\sin \frac{\pi j}{n}} \]
the first particle's velocity must equal
\[ \mathbf{v}_0 = \begin{pmatrix} 0\\ r \theta^\prime \end{pmatrix} = \begin{pmatrix} 0\\ \left(\dfrac{Gm\Sigma}{4r}\right)^\frac12 \end{pmatrix} \]
and we can figure those of the remaining particles with the matrix-vector product
\[ \begin{align*} \mathbf{v}_j &= \begin{pmatrix} \cos \theta_j & -\sin \theta_j\\ \sin \theta_j & \cos \theta_j \end{pmatrix} \times \mathbf{v}_0 = \begin{pmatrix} \cos \frac{2 \pi j}{n} & -\sin \frac{2 \pi j}{n}\\ \sin \frac{2 \pi j}{n} & \cos \frac{2 \pi j}{n} \end{pmatrix} \times \begin{pmatrix} 0\\ \left(\dfrac{Gm\Sigma}{4r}\right)^\frac12 \end{pmatrix}\\ &= \begin{pmatrix} -\left(\dfrac{Gm\Sigma}{4r}\right)^\frac12 \sin \frac{2 \pi j}{n}\\ \phantom{-}\left(\dfrac{Gm\Sigma}{4r}\right)^\frac12 \cos \frac{2 \pi j}{n} \end{pmatrix} \end{align*} \]
We may conclude that, if we give the first particle an initial velocity of \(v_0\) in the direction of the positive \(y\) axis, the polygonal system of particles shall exhibit periodic motion if
\[ \begin{align*} v_0 &= \left(\frac{Gm\Sigma}{4r}\right)^\frac12\\ v_0^2 &= \frac{Gm\Sigma}{4r}\\ m &= \frac{4 v_0^2 r}{G \Sigma} \end{align*} \]
To choose values for \(r\) and \(v_0\) my fellow students and I took inspiration from the initial conditions that our ethereal orrery had for the Earth
\[ \begin{align*} m &= 5.97 \, \mathrm{Hg}\\ x &= 0.152 \, \mathrm{Tm}\\ y &= 0 \, \mathrm{Tm}\\ v_x &= 0 \, \mathrm{Tm}/\mathrm{Ms}\\ v_y &= 0.0293 \, \mathrm{Tm}/\mathrm{Ms} \end{align*} \]
which we simplified to
\[ \begin{align*} r &= 0.15\\ v_0 &= 0.03 \end{align*} \]
Now, the rate at which the particles rotate about the origin is given by
\[ \theta^\prime = \frac{v_0}{r} \]
and so they shall return to their initial positions after a time of
\[ T = \frac{2\pi}{\theta^\prime} = \frac{2 \pi r}{v_0} \]
which we have used to divide a full cycle of a pentagon of particles into five hundred steps in deck 2.

Deck 2: A Single Cycle

Whilst our orrery certainly confirms that our reckoning has been sound, it also demonstrates that such systems are chronically unstable, with the slightest deviation from their proper positions and velocities throwing them into complete disarray. For example, deck 3 allows you to run it for several cycles, with the scheme that we are using to approximate their motion introducing small errors at each step; small errors that will accumulate until they are sufficient to disrupt the periodic motion entirely.

Deck 3: Several Cycles

The truly surprising aspect of this breakdown is how suddenly it occurs. The deviations of the particles from their correct positions and velocities are barely perceptible until moments before the whole system falls apart!

Fortunately, we can erase those errors by noting that after one \(n^\mathrm{th}\) part of a cycle the first particle must be in the position and have the velocity that the second occupied and had at the outset, that the second must likewise have the position and velocity that had been had by the third and so on and so forth. This gives us the means to rectify the states of the particles before they are overwhelmed by approximation error, as demonstrated by deck 4.

Deck 4: Correcting The Errors

To show just how effective this strategy is, you might increase the number of particles in both deck 3 and deck 4 to several dozens.

Pieces Of Eight

The surprising behaviour to which I alluded was the discovery by one of my fellow students[6] of a stable motion in which three particles traverse a figure of eight, as demonstrated by deck 5.

Deck 5: A Figure Of Eight

Note that approximation errors in the solution of the equations of motion causes this trio of particles to slowly rotate about the origin, as demonstrated by deck 6.

Deck 6: Precession

Having found that such a curious motion was possible, my fellow students and I commenced a systematic search for others of a similar nature, with another of us discovering an unstable motion involving four particles[7], as illustrated by deck 7.

Deck 7: Four Bodies

Another of my fellows set to cataloguing the initial conditions that give rise to such motions[8], of which deck 8 is an example involving five bodies.

Deck 8: Five Bodies

These are the only such \(n\)-body choreographies for which we have constructed decks, although we are aware of other student bodies that have done so, albeit by employing Fourier transforms to explicitly define the curves that they follow rather than by implicitly defining them through the application of Sir N-----'s laws, so as to better cope with their general chronic instability.


[1] On An Ethereal Orrery,, 2019.

[2] Further On An Ethereal Orrery,, 2019.

[3] Further Still On An Ethereal Orrery,, 2019.

[4] Yoshida, H., Construction of higher order symplectic integrators, Physics Letters A, Vol. 150, No. 5, 6, 7, Elsevier, 1990.

[5] Hut, P. & Makino, J., The 2-Body Problem: Higher-Order Integrators,, 2007.

[6] Moore, C., Braids in Classical Gravity, Physical Review Letters, Vol. 70, 1993.

[7] Simó, C., New Families of Solutions in N-Body Problems, Proceedings of the ECM , 2000.

[8] Dyck, J., Periodic solutions to the n-Body Problem, Masters thesis, University of Manitoba, 2015.

Leave a comment

This site requires HTML5, CSS 2.1 and JavaScript 5 and has been tested with

Chrome Chrome 26+
Firefox Firefox 20+
Internet Explorer Internet Explorer 9+