Periodic three-body orbits with vanishing angular momentum

Periodic three-body orbits with vanishing angular momentum in the Jacobi-Poincaré “strong” potential

V. Dmitrašinović, Luka V. Petrović and Milovan Šuvakov Institute of Physic Belgrade, Belgrade University, Pregrevica 118, Zemun, P.O.Box 57, 11080 Beograd, Serbia
Faculty of Physics, University of Vienna, Boltzmanngasse 5, 1090 Vienna, Austria
present address: ETH Zürich, Otto-Stern-Weg 1, 8093 Zürich, Switzerland,,

Moore [9] and Montgomery [10] have argued that planar periodic orbits of three bodies moving in the Jacobi-Poincaré, or the “strong” pairwise potential , can have all possible topologies. Here we search systematically for such orbits with vanishing angular momentum and find 24 topologically distinct orbits, 22 of which are new, in a small section of the allowed phase space, with a tendency to overcrowd, due to overlapping initial conditions. The topologies of these 24 orbits belong to three algebraic sequences defined as functions of integer . Each sequence extends to , but the separation of initial conditions for orbits with becomes practically impossible with a numerical precision of 16 decimal places. Nevertheless, even with a precision of 16 decimals, it is clear that in each sequence both the orbit’s initial angle and its period approach finite values in the asymptotic limit (). Two of three sequences are overlapping in the sense that their initial angles occupy the same segment on the circle and their asymptotic values are (very) close to each other. The actions of these orbits rise linearly with the index that describes the orbit’s topology, which is in agreement with the Newtonian case. We show that this behaviour is consistent with the assumption of analyticity of the action as a function of period.

45.50.Jf, 05.45.-a, 95.10.Ce
: J. Phys. A: Math. Gen.
  • August 2017

Keywords: three-body problem, periodic orbits,strong potential, nonlinear dynamics

1 Introduction

One singular homogeneous potential stands out for its properties: the attractive pairwise , also known as the Jacobi-Poincaré, or the “strong” potential has long been known to lead to a simplification of the classical N-body problem due to its dilation/conformal symmetry. In the 19th century Jacobi discovered its dilational symmetry [1] and Poincaré [2] discussed some properties of its action integral. More recently [3, 4, 5] it was shown that in one spatial dimension, this potential leads to integrable dynamics of the celebrated Calogero-Moser model.

In the meantime, this potential has acquired a physical significance: 1) it is part and parcel of post-Newtonian gravity, responsible for the gravitational collapse [6]; and 2) it is central to Efimov three-body bound states in molecular physics, see Refs. [7, 8] and references therein. Thus, it appears important that the classical dynamics in this potential be studied in more dimensions than one.

In two spatial dimensions the classical dynamics of this potential has been studied only sporadically, mostly emphasizing mathematical aspects, see Refs. [10, 11, 13], with only one significant numerical study [9]. It is known that this potential leads to a reduction of the number of degrees of freedom by unity, but not to integrability of the system, at least for unequal masses, [13]. That reduction implies a significant simplification, at least for periodic orbits.

The authors of Refs. [9, 10] have argued, along different lines, that periodic three-body orbits with all possible topologies, which is a countably infinite set, must exist in this potential; that is not the case for the Newtonian potential. Here we note that before the present paper, Moore [9] had found two numerical periodic three-body orbits in the strong potential, one of which (the figure-eight) exists in the Newtonian potential and was rediscovered in Ref. [11]. Moore’s second periodic orbit (figure-eight on the shape sphere, see Sect. 3) in the strong potential does not exist in Newtonian gravity, however. Thus, a plethora of previously unknown types of periodic orbits appears to be beckoning. Recent numerical searches for periodic orbits in the Newtonian potential have produced more than 200 new topologically distinct periodic orbits, see Refs. [14, 15, 16, 17, 18, 19, 20, 21]. The next natural step is to apply the same methods to the three-body problem in the strong potential.

Now, mathematical existence theorems do not necessarily tell one where to search and find the respective solutions, merely that they exist. Fortunately, in the case of the strong potential the phase space of initial conditions has fewer dimensions than the one in the case of Newtonian potential. Moreover, it has been shown that every periodic orbit (with one exception, that has non-zero angular momentum, however) must pass through a collinear configuration at least once during its period, [12, 11], which further reduces the space of initial configurations. These facts suggest that the numerical search for periodic orbits in the strong potential should be simpler than in the Newtonian one.

Here we report the results of our search: we found 24 periodic solutions, 22 of which are new,111The two orbits previously discovered by Moore, Ref. [9], are the figure-8 orbit in configuration space, with topology and the figure-8 in the shape space, with topology . falling into three algebraically defined (infinite) sequences. Many other candidates for periodic orbits have been detected, but could not be clearly separated within our self-imposed numerical precision requirements, for reasons that will be explained shortly.

The obtained results are interesting for the following reasons: (i) the distribution of periodic orbits in the (phase) sub-space of initial conditions is far from uniform: there are vast regions that are entirely devoid of periodic orbits, whereas several small regions are severely overcrowded with (infinitely) many closely packed - indeed overlapping - orbits; (ii) all periodic orbits discovered thus far belong to three sequences, indexed by integer ; each sequence having an accumulation point of its own, both as a function of the initial angle and as a function of the period ; (iii) the sensitivity to “errors”, i.e. to differences between the fiducial and the actual values of the initial conditions, rapidly increases with the sequential index ; (iv) the (almost) linear dependence of the scale-invariant period on the sequential index , or equivalently on topology, that was observed in the Newtonian case [22], breaks down here spectacularly - the periods do not increase linearly, but approach finite asymptotic accumulation points (a different one for each sequence); (v) a linear relation between the initial conditions - angles - and the periods holds in each sequence separately in the asymptotic limit ; (vi) a linear relation between the (scale-invariant) action and the sequential index , or equivalently the topology, holds in the asymptotic limit . 222In other homogeneous potentials one linear dependence would imply the other, but not here, see Sect. 4.1.

We are not aware of another system of ordinary differential equations with such, or similar distribution of periodic solutions. At least one (vi) property of these solutions, the asymptotically linear relation, as well as the existence of small non-linear corrections that disappear in the limit, can be explained in terms of analytic properties of the action.

This paper consists of six sections: after some preliminaries in Sect. 2.1, wherein we remind the readers about the basic facts of periodic orbits in the Jacobi-Poincaré “strong” potential must have zero total energy and constant hyper-radius, we present our search method in Sect. 2. Then, in Sect. 3 we show our results - the 22 new orbits that together with two previously known orbits form three sequences. In Sect. 4 we display the action “quantization” regularities within the three sequences. In Sect. 5 we show some mathematical argument based on analyticity of the obtained solutions. Finally, in Sect. 6 we summarize and draw our conclusions. A contains a description of the return proximity function’s dependence on the numerical accuracy of initial conditions. B contains detailed calculation of the action of the periodic three-body orbit in the (1D) Calogero-Moser model. In C we discuss the analytic properties of the action expressed in terms of complex variables.

2 Search for periodic solutions

Two-dimensional motions in the strong potential can be divided into (at least) three types: 1) positive energy () leads either to infinite expansion (for positive values of the time derivative of the hyper-radius, ), or to collapse of the system (for negative values of the time of the hyper-radius, ); 2) negative energy () always leads to the collapse of the system; and 3) zero energy () motions with: (a) non-zero values of the time derivative of the hyper-radius, ) lead to either collapse or infinite expansion, whereas (b) only vanishing time derivative of the hyper-radius, ) leads to dilation-invariant motion(s), which includes, but is not limited to, periodic ones. Thus periodic orbits occupy a small subset of all possible motions. Further, vanishing angular momentum condition ensures that the corresponding three-body orbit must be in a plane.

Thus, in case (3.b), the motions of the system can be described by two dimensionless degrees of freedom, such as two angles parametrizing a unit-radius sphere in three-dimensional “shape space”. There are three points on this sphere that correspond to two-body collisions, and thus to infinite kinetic and potential energies; these three points determine the topology of a periodic solution. The “volume” of such a configuration space (i.e. the area of the sphere) is finite, and the corresponding phase space (excluding the collision points) is also compact, which makes the subspace of periodic orbits finite. Moreover, due to the at least one syzygy theorem [12], the initial configuration space may be further reduced to just the equator on the shape sphere, excluding the three collision points.

2.1 Preliminaries

The planar, or two-dimensional (2D) three-body motions in the potential display several peculiarities that we list below.

1) All periodic solutions in the strong potential three-body problem must have exactly zero energy and the hyper-radius (the “scalar moment of inertia”) must stay constant at all times. These facts follow from the Lagrange-Jacobi identity, or the virial theorem.

2) We use the fact, shown in Ref. [12] for the Newtonian potential, and in Ref. [11] for the strong potential, that all but one periodic orbits must cross the equator on the shape sphere at least once during their periods (“the syzygy theorem”).

These two facts reduce the dimensionality of the initial state phase space and allow us to use a point on the equator of the shape sphere as an initial configuration. This facilitates the search for periodic solutions. Moreover, the virial theorem also implies that the action of a periodic orbit is not proportional to the (vanishing) energy of the periodic solution.

“Normalizing” orbital periods by spatially scaling solutions to the same energy level is not possible in this potential, due to the vanishing energy. Consequently, the constant appearing in Kepler’s third law - the “scale-invariant period” is trivial here, as it always equals zero. Instead, one keeps the scale, or the hyperradius of all orbits invariant at one fixed value. The periods depend on this scale, as , i.e., the product is scale invariant. The minimized action of periodic orbits is scale invariant, however, and can be used instead of the above “scale-invariant period” , to “measure” the topological dependence of orbits in the case of strong potential, in the sense of Refs. [20, 22].

2.2 Initial conditions

As stated above, the two main distinctions of the Jacobi-Poincaré “strong” potential are that: 1) the energy of all periodic orbits must be zero ; 2) the hyper-radius, or “overall size”, of all periodic orbits must remain constant at all times . These two conditions reduce the number of degrees of freedom of the planar equal-mass three-body problem from three to two. The conditions of 1) vanishing angular momentum, and 2) equal masses of all three bodies, are our own choice, in this paper.

There are independent kinematic variables that define the state vector (in phase space) of a three-body system (in two dimensions): for each of three bodies there are two coordinates ( and ) for each body, and two components ( and ) of its velocity vector. We set , which does not reduce the generality of our results, as certain scaling rules hold, see Ref. [23]. The scaling rules allow one to obtain solutions for any (real, positive) value of and/or of the common mass from the solutions presented here. Of course, distinct-mass orbits cannot be obtained from the equal-mass limit ones. The choice of center-of-mass system () cuts down this number to eight independent variables, as there are (only) two independent relative coordinate vectors and two corresponding velocities. The choice of vanishing angular momentum reduces this number down to seven.

We choose the so-called Euler initial configuration, the three bodies being collinear, say on the x-axis, with the distance between the bodies number “one” and number “two” equaling two units (2), and with the body number “three” at the mid-point between bodies number “one” and number “two”. That sets the (initial state) hyper-radius at . As the hyper-radius must remain constant during periodic motions, we are left with six independent variables. The conditions , , and put together imply

thus inferring that only one velocity two-vector is independent in this choice of initial conditions.

Finally, demanding , leaves the system with two independent variables: the angle between the -axis and the velocity 2-vector , and the overall size . Due to the zero-energy condition the size of the system, which has already been set at by our choice of initial positions, determines the value of the initial kinetic energy as , thus leaving the angle as the only free variable.

This means that, in order to find periodic orbits passing through the Euler point, the only variable that can be varied in this sub-space of initial conditions is the angle between the two components of the vector : .

2.3 Search Method

The return proximity function in phase space is defined as the absolute minimum of the distance from the initial condition by , where


is the distance (Euclidean norm) between two 12-vectors in phase space (the Cartesian coordinates and velocities of all three bodies without removing the center-of-mass motion). We define the return time as the time for which this minimum is reached. Numerical minimization of this function has been used in several successful searches for periodic orbits in the Newtonian potential [14, 15, 16, 17, 18, 19, 20, 21]. Those searches were conducted in a two-dimensional subspace, because the hyper-radius is variable in the Newtonian three-body periodic orbits. The choice of the initial conditions has been explained in Sect. 2.2.

In order to look for periodic solutions numerically, we have discretized the search window in the one-dimensional subspace and calculated the return proximity function for each grid point up to some pre-defined upper limit on the integration time , which was set at . We shall see that the period does not grow with the length of the orbit’s “word”, but rather approaches an asymptotic limit that depends on the algebraic structure of the “word”.

2.4 Solving the equations of motion

We did the calculation in two stages: at first we used a fourth-order Runge-Kutta-Fehlberg integrator, for the full set of three-body equations in Cartesian coordinates, which, of course, does not conserve the energy, or the hyper-radius. When that integrator started showing its limitations we wrote a new integrator for the two true dynamical variables (two hyperangles on the shape sphere), thus eliminating two constants of motion (the hyper-radius and the angular momentum) from the start. The third constant of motion is the energy , which must vanish () for all periodic orbits; that constraint can be “hard-wired” into the code using the Hamiltonian formalism that is manifestly symplectic and leads to a (much) higher accuracy even with a Runge-Kutta algorithm. That allowed us to determine the accuracy of the previous RKF calculation, as shown in detail in A.

It is well known that the planar three-body dynamics can be expressed in terms of hyper-spherical coordinates [24] as a function of the hyper-radius and the shape-sphere unit vector:


As stated already in Sect. 2.2, the hyper-radius is a constant of motion for all periodic solutions in this potential. We can assume without loss of generality that . Therefore, the configuration space is two dimensional in this case. We choose the polar angle and the azimuthal angle as the dynamical variables on the shape-sphere. The phase space is four dimensional with these two angles as the generalized coordinates and their two conjugate generalized momenta. Therefore, the equations of motion are four first-order Hamilton’s equations.

One can implement the zero-energy constraint into these four equations and thus eliminate one variable and one equation. The initial-state variables defined in Sect. 2.2 can be written as functions of a single parameter which was defined as the angle between the two components of the velocity vector : , or equivalently as the angle of the angular velocity, measured with respect to the equator, on the shape sphere:




is the hyper-angular potential as a function of and . One can easily check that this change of variables identically satisfies , by calculating


where the kinetic energy is given by the formula:


Three variables , , and completely define the state of the system. Using this set of variables we now have three first-order differential equations of motion, instead of four; they are:


where and are:


The initial conditions defined previously in Sect. 2.2 correspond to , , with taken as a free parameter. Using the variables , , and , the return proximity function can be redefined as follows


We solved the equations of motion using an explicit fifth-order Runge-Kutta algorithm [25] implemented in python library SciPy [26] with best relative tolerance of . We note that the Hamiltonian (first-order) nature of the equations of motion ensures exact conservation of energy at each step of the calculation, regardless of the particular algorithm used for the numerics.

The return proximity function has been calculated using linear interpolation between numerically obtained values of , , and . This interpolation method creates additional errors that cannot be explicitly calculated, but can be numerically investigated as in A.

2.5 Topological classification of orbits in the shape space

Any newly found periodic three-body orbit must be identified and classified so as to be distinguished from previously discovered orbits. For that purpose we use Montgomery’s topological classification [27]: he noticed the connection between the “fundamental group of a two-sphere with three punctures”, i.e. the “free group on two letters” , and the conjugacy classes of the “projective coloured/pure braid group” of three strands . The utility of this classification becomes apparent when we identify the “two-sphere with three punctures” with the shape-space sphere and the three two-body collision points with the punctures.

Graphically, this method amounts to classifying closed curves according to their “topologies” on a sphere with three punctures. A stereographic projection of this sphere onto a plane, using one of the punctures as the “north pole” effectively removes that puncture to infinity, and reduces the problem to one of classifying closed curves in a plane with two punctures. That leads to the aforementioned free group on two letters , where (for definiteness) denotes a clockwise “full turn” around the right-hand-side puncture, and denotes the counter-clockwise full turn around the other puncture, see Ref. [14]. For better legibility we denote their inverses by capitalized letters , .

Of course, there need not be only one solution with a particular topology: indeed orbits with different values of the angular momentum and identical topology define one (continuous) topological family of orbits. In the Newtonian potential there are sometimes multiple orbits333And sometimes none. with identical topologies and angular momenta, but, as we shall see, that does not happen here.

A specific sequence of letters, or a word, is not the only possible description of a periodic orbit, because there is no preferred initial point on a periodic orbit; therefore any other word that can be obtained by a cyclic permutation of the letters in the original word is an equally good description of such an orbit. For example the conjugacy class of the free group element contains also the element . The set of all cyclically permuted words is the aforementioned conjugacy class of a free group element (word). Thus, each family of orbits is associated with the conjugacy class of a free group element.

Moreover, the time-reversed orbits correspond to physically identical solutions, but their free group elements and their conjugacy classes are generally different. So, for example, families of orbits described by a and A are equivalent, but families ab and AB are not because the inverse of ab is BA, not AB.

There is finally one last ambiguity concerning all non-cyclically permutation symmetric orbits [22]. One may apply the stereographic projection of the shape sphere onto a plane using any one of the three two-body collision points as the North Pole, and in that way one may obtain three (at least in principle) different words. When the orbit is cyclically permutation symmetric, this ambiguity disappears, for all other orbits one may switch to the three-symbol labeling scheme [20] to resolve this issue rigorously. But, for all practical purposes one need not go to such lengths, as it is sufficient to make sure that the same two-body collision point (“puncture”) has been used as the North Pole in the stereographic projection for all the orbits treated. Then, the same algebraic relations must hold among the orbits’ words, irrespective of the choice of the puncture.

For a working algorithm that “reads” an orbit’s word, see the Appendix to Ref. [19].

3 Results

In Fig.1 we show the dependence on the initial angle in the range . There one can see many overlapping peaking structures some of them broad (“hills”), others quite narrow (“trees”, or “spikes”), in roughly four distinct regions (“groves”). Several properties of the newly discovered orbits shine through immediately:

  1. periodic orbits exist only in several finite, small segments of the (quarter-circular) perimeter [in the angular range (in radians)], rather than on the whole quarter-circle ;

  2. We found around 100 peaking structures (at the present level of precision) in this region, only 24 of which kept increasing their values as the precision was improved, see Appendix A. These 24 initial conditions appear in three distinct regions (“groves”) in Fig. 1, and are tabulated in Tables 1, 2, 3.

  3. The remaining (roughly 80) peaking structures are interspersed through all distinct regions (“groves”). These structures (probably) correspond to (quasi-periodic) orbits that do not pass exactly through the Euler point, but come close to it. Therefore, their values cannot be increased beyond some limiting value. Note that one “grove”, , consists entirely of such “quasi-periodic” orbits. The corresponding periodic orbits must be searched for in a two-parameter space (angle and the “displacement” from the Euler point), which is beyond our present scope.

  4. The periodic orbits found thus far can be classified into three sequences with well-defined algebraic structures describing the topologies of the orbits, [27], very much as in the Newtonian case, Ref. [20]. The topologies of solutions belonging to three sequence are defined as , , , with integer .

  5. The periods and initial angles of the -th periodic orbit converge to one of three accumulation points, depending on the sequence, see Fig. 2. It follows immediately that there is no linear relation between the ordering number and the period of the -th orbit in any of the three sequences, in contrast to the Newtonian case.

  6. This accumulation of orbits’ initial conditions makes the search for higher-, i.e., beyond some (fairly small number, such as 10) , orbits in the same sequence very difficult, if not completely impossible. This is because for some value of the index the difference between two sequential initial conditions becomes comparable with, or smaller than the numerical accuracy, see A. It should therefore be clear that finding further sequences with longer periods will be severly limited by the available numerical accuracy. Thus, at the present time we can only conjecture that initial conditions of orbits in sequences with topologies described by , and , are “hidden behind” the two already observed ones with .

  7. Not one topological-power (with topology of the form , ) orbit has been found, thus far, in this potential, in contrast with the Newtonian potential. If that were a general rule, it would explain the absence of the term from the , sequence in Table 3, as this term would be the second power of .

  8. Only one orbit has been found for each topology, i.e., there are no multiple periodic orbits with the same topology, again in contrast to the Newtonian three-body problem. Uniqueness of the figure-eight orbit has been proven in Ref. [10], so the uniqueness of the here observed orbits leads us to conjecture that all periodic orbits in this potential might be unique.

Figure 1: The negative logarithm of the return proximity function resulting from a scan over the initial angle range .
0 0.444445
1.02824692798800732 1.21269962942932730
0.96097510395224572 1.41012664240185859
0.94454479748942055 1.46028397633386153
0.93814716192692216 1.48000532746712521
0.93500390397660027 1.48973256473400473
0.93322897503736779 1.49523614367717128
0.93212904360729043 1.49865060080713763
0.93140041145929431 1.50091405578946335
0.93089288086015565 1.50249147597823884
0.93052521736517435 1.50363454265253282

Table 1: Initial condition (angle ), and period of three-body orbits in the sequence , , in the “strong force”, with coupling constant is taken as and equal masses as . The total energy is zero (), so the initial velocities are determined by a single angle . We also show the lowest return proximity value , the period and the free group element (f.g.e.). The orbit denoted by 0 is the Calogero-Schubart colliding one.
1.02824692798800732 1.21269962942932730
1.17049635743816727 2.51487201239182445
1.16542247378727315 2.98924613111245119
1.16451391284537076 3.07222781280329249
1.16420804317607640 3.10000877826202093
1.16406923772482696 3.11258622399824114
1.16399470076300249 3.11933231322032567
1.16395009046894526 3.12336590267499581

Table 2: Initial condition (angle ), and period of three-body orbits in the sequence , in the “strong force” potential, with the same values of parameters as in Table 1. We also show the lowest return proximity value , the period and the free group element (f.g.e.).
1.02824692798800732 1.21269962942932730
0.96053567047115329 2.61984543973343476
0.94436142254818101 2.66897947579246475
0.93805794284507649 2.68839038409076103
0.93495448623723931 2.69799863788414607
0.93319891411979261 2.70345025077037748
0.93210944940321705 2.70684099103309794
0.93138694948558265 2.70908665447520480
Table 3: Initial condition (angle ), and period of three-body orbits in the sequence , , in the “strong force” potential, with the same values of parameters as in Table 1. We also show the lowest return proximity value , the period and the free group element (f.g.e.). Note that the orbit is missing.
Figure 2: (online colour) The initial angle as a function of the period , for orbits in the three sequences, Tables 1,2,3 together with quadratic fits (colored continuous lines). Three accumulation points are clearly visible, as well as the fact that one sequence obstructs the “visibility” of another.

In Figs. 3 and 4 we show two solutions on the shape sphere and in configuration space. These are typical solutions, insofar as they do not approach any of the three collision points. Indeed, all the orbits found so far fall between an inner and an outer envelope, both in the real and in the shape space: the inner envelope precludes the orbit(s) from getting too close to a(ny) collision point.


Figure 3: (color on line) Orbits on the shape-space sphere: (left hand side) “figure-8 on the shape-sphere” - ; (right hand side) orbit .


Figure 4: (color on line) Orbits in configuration space: (left-hand side) “figure-8 on the shape-sphere” - ; (right-hand side) orbit .

4 Topological dependence of the action

4.1 Some properties of the action

For periodic orbits in a homogeneous potential with (arbitrary) degree of homogeneity , the (minimized) action can be related to the energy and period of the orbit as follows

Note that the case is singular: both the numerator and the denominator on the right-hand-side of this identity are zero. Consequently, the action is not linearly proportional to the period , and the two must be separately determined. Therefore, the periods are not constrained by any “quantization” of the action.

The following approximate linear relation for the minimized actions , of two orbits described by topologies, or free group “words”, and ,

has been found numerically in the () Newtonian potential, Refs. [19, 22, 20].

In an potential holds, and the (minimized) action is scale invariant. Then

ought to hold and should be tested; the trouble with this specific example is that only one orbit in each topological sector of the strong potential three-body problem has been found (thus far) - there are no known examples of topological powers of any orbit in the potential. Therefore, we cannot test this relation directly.

A more general relation, that the actions of all orbits in the same sequence are proportional to the “length” of the word , i.e., to the sum of , was conjectured in Ref. [22]. Thus, one ought to see a correlation between the actions and lengths of orbits’ words

This conjecture will be tested for the sequence in Sect. 4.2.1, for sequence in Sect. 4.2.2 and for the sequence in Sect. 4.2.3.

The action S can be directly evaluated as

where represents the periodic solution to the e.o.m. with given energy . We may factor out the constant term (due to ) in front of the integral,


and evaluate the remaining integral with our solutions. In spite of the explicit scale () dependence, this integral is scale invariant, due to the time ’s compensating scale dependence.

4.2 Action “quantization”

4.2.1 Sequence

We test the linear action-topology regularity for the sequence, Fig. 5, with data from Table 4. At first, we fit the five points ( [1,5]) with linear, quadratic, cubic and quartic polynomials: The linear fit to orbits with [1,5] is the quadratic one 444The cubic fit is and the quartic one . Then we use these fits to predict: 1) the single lower-lying () orbit’s action; 2) the higher-lying ( [6,10]) orbits’ actions S and compare them with the actual values in Table 5, shown in Fig. 5. One can immediately see that the cubic and the quartic fits deviate significantly from the data at higher values of , Fig. 5, thus suggesting that all terms with powers higher than the second are inappropriate.

0.88889 1.09861 - -
2.4254 14.4191 12.2256 -
2.82025 27.4681 13.0453 0.8197
2.92057 40.231 12.7238 -0.3215
2.96001 52.9039 12.6713 -0.0525
2.97947 65.5371 12.6356 -0.0357
2.99047 78.1493 12.547 -0.0886
2.9973 90.7493 12.7089 0.1619
3.00183 103.341 12.5883 -0.1206
3.00498 115.928 12.587 -0.0013
3.00727 128.51 12.582 -0.005
Table 4: The action , the period , the difference of action , the difference of differences and the free group element (f.g.e.) of three-body orbits in the sequence. The entry denoted by 0 is the Calogero-Schubart colliding orbit whose period and action were calculated analytically in B.

Thence we see that: 1) all four fits agree with/correctly predict the order of magnitude of the (the Calogero-Schubart orbit’s) action, where the calculated value , see B, satisfies the inequality ; 2) the largest deviations of and from their mean values are for , the calculation of which involves .; 3) the actual values for [6,10] generally fall between the linear and the quadratic fits: , see Table 4, and Fig. 5.

Note that the constraint is in (rough) agreement with the in Table 4, and that the fitted value of generally agrees with the sign and the order of magnitude of . As increases, the values of and approach constants. All of this suggests that the value of the action asymptotically approaches a linear function of , as .

The linear fit to all 11 orbits in sequence I, and the complete quadratic one, are in slightly better agreement with the data than the ones without the first term.

S f.g.e.
0.88889 1.09861 1.82486 1.38896 0.988914
2.99047 78.1493 78.3531 77.9172 78.3172
2.9973 90.7493 91.1078 90.236 91.4361
3.00183 103.341 103.862 102.43 105.002
3.00498 115.928 116.617 114.5 119.158
3.00727 128.51 129.372 126.445 134.046
Table 5: Predictions (linear , quadratic , cubic ), of the action S of three-body orbits in the sequence. The first orbit () is the Calogero-Schubart colliding one.

Figure 5: S(n) “curve” for the sequence of solutions with =0,…,10 from Table 4 together with four different polynomial fits (linear, quadratic, cubic and quartic (blue, green, magenta and red, respectively) to the =1,…,5 orbits, the values at being predictions. One can see that all four fits predict the point correctly, but only the linear and the quadratic fits (blue and green, respectively) bracket the higher- values.

4.2.2 Sequence

We repeat this procedure for the sequence and find similar results in Table 6 and Fig. 6.

2.4254 14.4191 - -
5.02974 22.5645 8.1454 -
5.97849 36.0148 13.4503 5.3049
6.14445 48.8613 12.8465 -0.6038
6.18929 61.4806 12.6193 -0.2272
6.22509 74.2112 12.7306 0.1113
6.23852 86.8308 12.6196 -0.111
6.2463 99.433 12.6022 -0.0174
Table 6: The action , the period , the difference of action , the difference of differences and the free group element (f.g.e.) of three-body orbits in the sequence (II). The entry denoted by () is also the first orbit from sequence I.

The linear fit, excluding the point, is the quadratic one .

Figure 6: The action S as a function of for the sequence of solutions with =0,…,7 from Table 6 together with two lowest-order polynomial fits (the same color code as in Fig. 5).

Including the point, the fits become and . Note that the magnitudes of in these fits generally agree with those in the sequence , with being an exception, for obvious reasons.

4.2.3 Sequence

We repeat this procedure for the sequence and find similar results in Table 7 and Fig. 7.

2.4254 14.4191 - -
5.23969 41.8824 13.7316 -
5.33796 54.6432 12.7608 -0.97085
5.37678 67.3197 12.6765 -0.0843
5.39600 79.9482 12.6285 -0.048
5.40845 92.6465 12.6983 0.0698
5.41368 105.229 12.5825 -0.1158
5.41817 117.864 12.6350 0.0525
Table 7: The action , the period , the difference of action , the difference of differences and the free group element (f.g.e.) of three-body orbits in the sequence. The entry denoted by 0 is the Calogero-Schubart colliding orbit whose period and action were calculated analytically in Appendix B.

The linear fit, excluding the point, is the quadratic one . Including the point, the fits become and . Once again the same remarks apply: the values of in these fits are generally close to those in the previous two sequences. The manifest question is: why?

Figure 7: The action S as a function of for the sequence of solutions with =0,2,…,8 from Table 7 together with two lowest-order polynomial fits (the same color code as in Fig. 5).

5 Analytical arguments

Calculus with one complex variable has been used to evaluate the action of the two-body problem, see Refs. [23, 28]. Application of few complex variables methods to the planar few-body problem(s) in celestial mechanics can be traced back (at least) to Siegel and Moser [29]. That means two complex variables for the planar three-body problem. In the case of the strong potential, conservation of the hyper-radius, together with the zero angular momentum condition, eliminate one of two complex variables, see C. 555Analyticity of solutions to the three-body equations of motion in the Newtonian potential, has been proven by Sundman, Ref. [30, 31], subject to the condition that there are no three-body collisions, albeit allowing for two-body collisions. Analogous proofs do not exist in the case of the strong potential, to our knowledge, but we shall nevertheless proceed, subject to the assumption that the periodic solutions to the strong potential e.o.m. are also analytic, at least within some region of the complex plane, that is, perhaps, determined by collisions. As we deal only with collisionless orbits here, we believe to be within the assumed region of analyticity.

The action of a periodic orbit can be written as a line integral over a closed path on a shape sphere parametrized by two angles in Eq. (14). Such an integral over the shape-sphere variables can be recast as a Cauchy closed-contour integral over a single complex variable, see pp. 22-30 in Ref. [32], or C, with three poles, because the sphere can be stereographically projected into a (real) plane, and the (real) plane can be replaced by a single complex variable [32]. Thus we have a closed line integral in a complex plane with three poles (one of which is at infinity) - its value is determined by Cauchy’s residue theorem and by the topology of the integration path as follows:


Then, wee see that the value of the integral is naturally proportional to the number of times the path goes around either of two poles, which is the basic hallmark of our numerical results.


where the term(s) arise (only) if the difference between the residues at the -th and the st turn does not equal zero.

Note, however, that: (a) there is an additive constant in the fits to the sequence (I), above - what is its source in terms of complex integration? (b) the non-vanishing value of the quadratic term in the fit(s) indicates that the value of the residue changes with , yielding the term - but, why should that happen?

(a) The contribution in sequence (I) simply cannot originate from a pole, as its contribution would have to be multiplied by , and consequently would vanish at . This is not necessarily true for sequences (II) and (III), as their free-group elements, and , respectively, indicate that their terms still involve some contribution from the two poles in the complex plane (without the point at infinity).

In the case of sequence (I), can be due to the pole at infinity, which is circumscribed only once, or from crossing a branch cut. A branch cut would also explain the different changes of the residue with the changing , but that cut would have to be logarithmic in order to change the residues at infinitely many values of . Now, no branching points are manifest in Eq. (14), but two have been explicitly found, associated with the two-body collison poles, and the logarithmic nature of the branch cut shown, by a calculation of the action of one particular solution, see B.

In B we have evaluated the action of the collisional zero-angular momentum Calogero-Schubart orbit, Refs. [3, 4, 35, 33, 34], and show that it provides: 1) a logarithmic branch cut; 2) the correct size of the constant term in the power series, that is related to the discontinuity of the integral at the branch cut.

(b) The size of the change of the residue due to the “crossing of the cut” can not be estimated with our present knowledge, but if it is small compared with the actual value of the residue, , then it would explain both (1) the approximately equal slopes of the graphs, for different sequences of orbits; and (2) the branch cut implies a non-vanishing quadratic term that varies from one value of to another. The above arguments suggest that the analyticity of solutions could be sufficient to explain the remarkable, numerically observed behaviour of the action.

6 Summary and Conclusions

We have searched for periodic planar three-body orbits in the strong three-body problem, and found 24 orbits, 22 of which are new. These solutions neatly fill in the lowest entries in Montgomery’s topological classification of three-body orbits and contain several solutions with topologies that do not exist in the Newtonian potential. These 24 orbits fall into three algebraically well-defined sequences.

Two properties of the discovered orbits shine through immediately: 1) periodic orbits exist only in (several) small regions of the allowed phase space; 2) the orbits discovered thus far form three sequences with respect to the algebraic description of their topologies; 3) the orbits’ periods and initial conditions converge to, in this case three (thus far) accumulation points, which makes them difficult to disentangle beyond some fairly small number (about 10 when working to a precision of 16 decimal places) of orbits in each sequence. We are not aware of another system of ordinary differential equations that has a similar structure of solutions.

Thus we have shown that the problem of finding zero-angular-momentum periodic three-body orbits in the strong potential is amenable to a numerical study, albeit with some perhaps unexpectedly severe limitations imposed by the structure of the solutions themselves. The task of finding an efficient method of disentangling further sequences of orbits and their accumulation points remains as a challenge for future work.

The new orbits’ periods depend on the specific sequence of orbits, but always remain bounded from above, approaching different asymptotic limits, thus explicitly contradicting the conjectured linear relationship between the scale-invariant period and the index of the sequence. Furthermore, in each sequence there is a remarkable almost linear relation between the initial angle and the period , as . The orbits’ configuration-space trajectories fall within compact, clearly separated regions in the plane and have both an inner and an outer envelope. All of this put together might indicate a certain structure to the set of all periodic solutions, if not an outright solvability of the problem.

We have calculated the values of the (scale invariant) action, which equals twice the time integral over the period of the kinetic energy, for these periodic orbits and found that they are “topologically quantized”, i.e., approximately linearly proportional to the index of the sequence (or equivalently to the number of “letters” in the “words” describing their topologies), with a small quadratic correction that diminishes with increasing . The existence of both of these features of the action can be explained in terms of Cauchy’s theorem for a closed contour integral.

Our results ought to have interesting consequences for the two physical systems that contain the strong three-body potential: (1) the post-Newtonian gravity; and (2) semiclassical calculation of the Efimov effect. Unfortunately, our present results are insufficient to draw conclusions about these two physical systems, as yet. This calls for new and different methods to be invented and applied to this system.

The work of V.D. and M.Š was financed by the Serbian Ministry of Science and Technological Development under grant numbers OI 171037 and III 41011. The work of L.V.P. was conducted during the summer vacations of 2013 and 2014, and was accepted in lieu of the Praktikum in the “Scientific Computing” module of the third year of undergraduate curriculum at the Faculty of Physics of the University of Vienna. All numerical work was done on the Zefram cluster, Laboratory for gaseous electronics, Center for non-equilibrium processes, at the Institute of Physics, Belgrade.

Appendix A Numerical accuracy of initial conditions for periodic orbits in the “strong” potential

The negative logarithm of the deviation of the approximate initial angle from the exact value of a periodic orbit represents the highest decimal place of the “error” committed. In Figs. (8),(9) we show the dependence of the (negative logarithm of the) return proximity function as a function of the negative logarithm of the “error” for periodic orbits in the three sequences, respectively, during the process of minimization of the return proximity function using the gradient descent procedure. There one can see that generally speaking one needs an accuracy of about 14 decimal places in order to achieve uniform return proximity function values within one sequence.

Figure 8: The negative logarithm of the return proximity function as a function of the negative logarithm of the deviation (“error”) of the initial angle from the exact value for the first ten periodic orbits, in sequence I. Here one can graphically see the importance of the higher-lying decimal points for the value of the return proximity function.

Note that the “threshold” (minimal value of at which the starts growing) moves from 10 to 13, depending on the orbit () and the sequence. Moreover, note that some orbits, such as in sequence I, have fairly large values of even at fairly low values of .


Figure 9: Same as in Fig. 8, only for sequences II and III.

Appendix B Action of the periodic three-body orbit in the one-dimensional “strong” potential

In the one-dimensional (1D) “strong” potential Calogero-Moser model, Refs. [3, 35, 5], the periodic three-body problem has only one variable - the meridional angle (“longitude”) of a point on the equator of the shape-sphere, see Refs. [3, 4, 35, 33, 34]. This orbit corresponds to the collinear v. Schubart orbit in Newtonian gravity.

b.1 Solving the equation of motion

Energy conservation yields as the equation of motion

which is readily solved as

where . Integrate this equation over to find

which, after the simplification

leads to

Note the “dimensionless time” variable . This is identical to


Setting we have finally

as the solution.

b.1.1 Solution with

In the limit, the time and the e.o.m. turns into

The solution to this equation is

where is an integration constant; setting we find:


which is equivalent to the Ansatz Eq. (17). Note that the equation (18) does not have real solutions for time larger than , as the argument of exceeds unity. This means that the turns imaginary at the point in time, i.e., is a branch point of the solution.

We demand that be real, which implies an upper and a lower limit on the (real) time variable :

For other times, the solution to the Ansatz Eq. (17) becomes imaginary/complex. These two “boundaries” on time imply the following limits on the value of :

Thus, one cycle consists of an oscillation of from to and back. Therefore, the half-period is

b.2 Action of the Schubart-Calogero-Moser solution

The action S can be evaluated from either the time integral over one period of twice the kinetic energy, , or the negative of twice the potential energy , where (), see Refs. [4, 33, 34].

The action evaluated over one period with the solution Eq. (17) equals formally

but manifestly the integrand has poles at , i.e., at the values of the integration variable , and the integral is not defined for due to the turning imaginary (Eq. (17) implies a logarithmic branch cut in the complex- plane).

Therefore, we must either 1) make an analytic continuation of the solution for time values ; or 2) reformulate the integral within the confines of the real solution .

The action of a periodic orbit must be: 1) independent of the starting (and ending point), i.e.,

2) equal to twice the action over one half-period

due to the symmetry under time reversal of the integrand. The integral

is singular: it has a simple pole within the integration range. When we change the integration variable we find