Three body free fall periodic orbits:
new remarkable features

By Alexander Gofen
December, 2021

The three body motion under the Newtonian gravitation has been intensively studied since Isaac Newton over 200 years, still presenting a challenge. Besides the special cases of elliptic motion by Euler and Lagrange, no other versions of regular motion of three bodies were known for a long time. Since emergence of computers, also numeric simulations were used to explore the orbits for various initial settings, most of which generated a chaotic motion.

So more surprising were discoveries of remarkable types of periodic motion of three bodies obtained in computer assisted researches. Such were the discovery of ...
In a typical random setting trajectories of three body motion are chaotic. That is why the discovery of the cases of periodic fall by Xiaoming Li and Shijun Liao [1] (2018) is so amazing. They have discovered hundreds of settings for resting three bodies of particular masses, which led to periodic motion: i.e. the bodies started at the given points of rest and returned back to their initial point after motion along sophisticated orbit during a period T. It's worth particular mentioning that ...
All such types of motion are cataloged and displayed as movies at the authors' site [1]. The first 30 simulations in the table present the cases of all masses equal to 1. They also come with the Demo installation of the software Taylor Center (where you can watch them in a resolution much higher than in the movies [1]).

This workshop is dedicated to studying new properties of the above mentioned 30 simulations: the properties discovered by chance in numeric experiments with these simulations performed and verified with the Taylor Center software. In addition to the remarkable properties in the bullet list above it appeared that...

  • In several of the 30 simulations the triangle formation at the second resting point is congruent to the initial triangle: mirrored and turned at a certain angle. Moreover...

  • In a few such cases the edges of the turned triangle are parallel to those of the initial triangle.

These properties are not mentioned in the original sources
[1-3], so that as of this moment (December 2021) it is not known whether the authors were aware about these properties.

The table below shows in which of the 30 simulations for equal masses [1] these earlier unknown (?) properties take place:

Parallel edges
Type of symmetry










Yes Yes Central
Yes Yes Central







Yes Yes Central


Some illustrations

Here are several illustrations of the characteristic cases mentioned above presented in high resolution of the Taylor Center software. All the 30 simulations may be viewed also in real time animation within this software as explained here.

In the pictures below the trajectories are color coded as A→A', B→B', C→C' where primed point is the second point of rest for the respective initial point. However the edges and angles of the triangle ∆ABC at the initial moment may not necessarily match the edges and angles of the second triangle even when both triangles are congruent: they can flip in various manner as the consequence of the motion.

Simulation #1: the second triangle and the first are dissimilar.

Simulation #18: the second triangle and the first are congruent.
AB = 1,                                          AC = 0.718590490501491,        BC = 0.309097112311281
A'B' = 1.00000000170251,          A'C' = 0.718590493477556,       B'C' = 0.309097111598536
A'B'/AB = 1.00000000170251, A'C'/AC =  1.00000000414153 , B'C'/BC=  0.999999997694106

Simulation #14: not only is the second triangle congruent to the first, but also their respective sides are parallel.
AB = 1,                                          AC = 0.921456266427119,        BC = 0.679289996641939
A'C' = 0.999999999906395,          A'B' = 0.921456266288986,       B'C' = 0.679289996759996
A'C'/AB = 0.999999999906395, A'B'/AC =  0.999999999850092 , B'C'/BC=  1.00000000017379
Moreover, AB||A'C', AC||A'B',  BC||B'C'.

The supporting data obtained via the Taylor Center software

The integration method for the 30 simulations used here was the same modern Taylor method which was used by the authors in a frame of
the so called Clean Numerical Simulation (CNS) [2] at a super-computer. As the authors wrote in [2], their implementation of the Taylor method admits arbitrary order and double precision. They, however, didn't mention what was the order used by them, and what the double precision means for their super computer.

In this Taylor Center software (also with an arbitrary order) we used the order 30 with maximum precision of float point numbers supported by processors Intel, which is the
10 byte float type called extended: 63-bit mantissa and 16-bit exponent.

For each of the 30 computed simulations the outputted data is comprised of the following elements:
If the congruence did take place, the data package contains also:
a3/ao1 = 1.00000000202946
a2/ao2 = 0.999999997920113
a1/ao3 = 1.00000000030637
                  ao1                                    ao2                              ao3
a1   138.847783298464°   179.999999953788°   104.361663992048°
a2   179.999999940341°   138.847783191836°   63.2094472309557°
a3   116.790552658371°   75.6383359097635°   5.33608528907246E-008°

If the
angles close to 0° or 180° are detected, the pairs of parallel edges are displayed.

Here is the actual data obtained with the computer for the 30 samples.

Notes on accuracy

We see that
the proportions expected to be 1, and the velocities expected to be zero, actually differ from the targeted values. The accuracy of those values depend on several factors: on the accuracy of the parameters provided by the author, and on the limits of accuracy in this Taylor integrator.

In the Taylor Center software the accuracy of integration (in ideal cases) may achieve up to 63 binary digits of the mantissa all being correct at every step, which corresponds to 18 correct decimal digits. Even with such ultimate 63-bit accuracy achievable at one step, the global error increases with growing number of steps (due to the rounding errors, or worse, due to catastrophic subtraction error in some problems). For example, in a test for simulation #1 integrated from 0 to its period T and back to 0, the accuracy of the method in terms of the absolute error was about 10-13 for the positions and
10-12 for the velocities in about 2500 integration steps.

Thinking about the reasons that the actual accuracy of the proportions and velocities obtained in this numerical experiment is not so good, first observation is that the values of the initial positions and the periods were specified by the authors only up to 11 decimal digits (instead of possible 18 in the PCs). Therefore, in order to see whether a more accurate value of the period T may make the velocities closer to zero at the second rest position, I resorted to the unique feature of this software – switching between different independent variables. In so doing, I switched from integration in independent variable t to integration in the velocity v3 as a new independent variable
setting the termination condition  v3 = 0 with the hope that  u1, v1, u2, v2, u3, would get closer to zero. They did not: perhaps it's an insufficient accuracy in the authors' given initial positions which caused the here observed inaccuracy in the discussed values.


1. Xiaoming Li and Shijun Liao, Movies of the Collisionless Periodic Orbits in the Free-fall Three-body Problem in Real Space or on Shape Sphere

2. Xiaoming Li and Shijun Liao, Collisionless periodic orbits in the free-fall three-body problem.

3. Xiaoming Li, Shijun Liao
Collisionless periodic orbits in the free-fall three-body problem