Periodic threebody orbits with vanishing angular momentum in the JacobiPoincaré “strong” potential
Abstract
Moore [9] and Montgomery [10] have argued that planar periodic orbits of three bodies moving in the JacobiPoincaré, 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.
pacs:
45.50.Jf, 05.45.a, 95.10.Ce
August 2017
Keywords: threebody 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 JacobiPoincaré, or the “strong” potential has long been known to lead to a simplification of the classical Nbody 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 CalogeroMoser model.
In the meantime, this potential has acquired a physical significance: 1) it is part and parcel of postNewtonian gravity, responsible for the gravitational collapse [6]; and 2) it is central to Efimov threebody 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 threebody 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 threebody orbits in the strong potential, one of which (the figureeight) exists in the Newtonian potential and was rediscovered in Ref. [11]. Moore’s second periodic orbit (figureeight 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 threebody 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 nonzero 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,^{1}^{1}1The two orbits previously discovered by Moore, Ref. [9], are the figure8 orbit in configuration space, with topology and the figure8 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 selfimposed 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) subspace 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 scaleinvariant 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 (scaleinvariant) action and the sequential index , or equivalently the topology, holds in the asymptotic limit . ^{2}^{2}2In 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 nonlinear 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 JacobiPoincaré “strong” potential must have zero total energy and constant hyperradius, 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 threebody orbit in the (1D) CalogeroMoser model. In C we discuss the analytic properties of the action expressed in terms of complex variables.
2 Search for periodic solutions
Twodimensional 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 hyperradius, ), or to collapse of the system (for negative values of the time of the hyperradius, ); 2) negative energy () always leads to the collapse of the system; and 3) zero energy () motions with: (a) nonzero values of the time derivative of the hyperradius, ) lead to either collapse or infinite expansion, whereas (b) only vanishing time derivative of the hyperradius, ) leads to dilationinvariant 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 threebody 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 unitradius sphere in threedimensional “shape space”. There are three points on this sphere that correspond to twobody 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 twodimensional (2D) threebody motions in the potential display several peculiarities that we list below.
1) All periodic solutions in the strong potential threebody problem must have exactly zero energy and the hyperradius (the “scalar moment of inertia”) must stay constant at all times. These facts follow from the LagrangeJacobi 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 “scaleinvariant 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 “scaleinvariant 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 JacobiPoincaré “strong” potential are that: 1) the energy of all periodic orbits must be zero ; 2) the hyperradius, 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 equalmass threebody 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 threebody 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, distinctmass orbits cannot be obtained from the equalmass limit ones. The choice of centerofmass 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 socalled Euler initial configuration, the three bodies being collinear, say on the xaxis, with the distance between the bodies number “one” and number “two” equaling two units (2), and with the body number “three” at the midpoint between bodies number “one” and number “two”. That sets the (initial state) hyperradius at . As the hyperradius 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 twovector 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 2vector , and the overall size . Due to the zeroenergy 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 subspace 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
(1) 
is the distance (Euclidean norm) between two 12vectors in phase space (the Cartesian coordinates and velocities of all three bodies without removing the centerofmass 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 twodimensional subspace, because the hyperradius is variable in the Newtonian threebody 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 onedimensional subspace and calculated the return proximity function for each grid point up to some predefined 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 fourthorder RungeKuttaFehlberg integrator, for the full set of threebody equations in Cartesian coordinates, which, of course, does not conserve the energy, or the hyperradius. 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 hyperradius 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 “hardwired” into the code using the Hamiltonian formalism that is manifestly symplectic and leads to a (much) higher accuracy even with a RungeKutta 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 threebody dynamics can be expressed in terms of hyperspherical coordinates [24] as a function of the hyperradius and the shapesphere unit vector:
(2) 
As stated already in Sect. 2.2, the hyperradius 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 shapesphere. 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 firstorder Hamilton’s equations.
One can implement the zeroenergy constraint into these four equations and thus eliminate one variable and one equation. The initialstate 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:
(3) 
(4) 
where
(5) 
is the hyperangular potential as a function of and . One can easily check that this change of variables identically satisfies , by calculating
(6) 
where the kinetic energy is given by the formula:
(7) 
Three variables , , and completely define the state of the system. Using this set of variables we now have three firstorder differential equations of motion, instead of four; they are:
(8)  
(9)  
(10) 
where and are:
(11)  
(12) 
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
(13) 
We solved the equations of motion using an explicit fifthorder RungeKutta algorithm [25] implemented in python library SciPy [26] with best relative tolerance of . We note that the Hamiltonian (firstorder) 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 threebody 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 twosphere 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 “twosphere with three punctures” with the shapespace sphere and the three twobody 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 righthandside puncture, and denotes the counterclockwise 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 orbits^{3}^{3}3And 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 timereversed 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 noncyclically permutation symmetric orbits [22]. One may apply the stereographic projection of the shape sphere onto a plane using any one of the three twobody 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 threesymbol 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 twobody 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:

periodic orbits exist only in several finite, small segments of the (quartercircular) perimeter [in the angular range (in radians)], rather than on the whole quartercircle ;

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.

The remaining (roughly 80) peaking structures are interspersed through all distinct regions (“groves”). These structures (probably) correspond to (quasiperiodic) 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 “quasiperiodic” orbits. The corresponding periodic orbits must be searched for in a twoparameter space (angle and the “displacement” from the Euler point), which is beyond our present scope.

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.

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 .

Not one topologicalpower (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 .

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 threebody problem. Uniqueness of the figureeight 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.
f.g.e.  

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 
f.g.e.  

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 
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 
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.
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 righthandside 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 threebody 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,
(14) 
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 actiontopology 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 ^{4}^{4}4The cubic fit is and the quartic one . Then we use these fits to predict: 1) the single lowerlying () orbit’s action; 2) the higherlying ( [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.
f.g.e.  

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 
Thence we see that: 1) all four fits agree with/correctly predict the order of magnitude of the (the CalogeroSchubart 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 
4.2.2 Sequence
f.g.e.  

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 
The linear fit, excluding the point, is the quadratic one .
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
f.g.e.  

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 
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?
5 Analytical arguments
Calculus with one complex variable has been used to evaluate the action of the twobody problem, see Refs. [23, 28]. Application of few complex variables methods to the planar fewbody problem(s) in celestial mechanics can be traced back (at least) to Siegel and Moser [29]. That means two complex variables for the planar threebody problem. In the case of the strong potential, conservation of the hyperradius, together with the zero angular momentum condition, eliminate one of two complex variables, see C. ^{5}^{5}5Analyticity of solutions to the threebody equations of motion in the Newtonian potential, has been proven by Sundman, Ref. [30, 31], subject to the condition that there are no threebody collisions, albeit allowing for twobody 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 shapesphere variables can be recast as a Cauchy closedcontour integral over a single complex variable, see pp. 2230 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:
(15) 
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.
(16) 
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 nonvanishing 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 freegroup 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 twobody 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 zeroangular momentum CalogeroSchubart 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 nonvanishing 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 threebody orbits in the strong threebody problem, and found 24 orbits, 22 of which are new. These solutions neatly fill in the lowest entries in Montgomery’s topological classification of threebody orbits and contain several solutions with topologies that do not exist in the Newtonian potential. These 24 orbits fall into three algebraically welldefined 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 zeroangularmomentum periodic threebody 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 scaleinvariant 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’ configurationspace 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 threebody potential: (1) the postNewtonian 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.
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.
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 .
Appendix B Action of the periodic threebody orbit in the onedimensional “strong” potential
In the onedimensional (1D) “strong” potential CalogeroMoser model, Refs. [3, 35, 5], the periodic threebody problem has only one variable  the meridional angle (“longitude”) of a point on the equator of the shapesphere, 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
(17) 
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:
(18) 
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 halfperiod is
b.2 Action of the SchubartCalogeroMoser 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 halfperiod
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