Orbits in the problem of two fixed centers on the sphere

TL;DR

This paper establishes a mathematical isomorphism between the two fixed center problem on a sphere and its planar counterpart, enabling analysis of spherical orbits through elliptic integrals and functions.

Contribution

It introduces a novel isomorphism via gnomonic projections that links spherical and planar two fixed center problems, facilitating bifurcation and dynamical analysis.

Findings

Bifurcation diagrams of the spherical problem are derived from planar systems.

Spherical orbits are expressed using Jacobi elliptic functions.

Quadratures are transformed into elliptic integrals for analysis.

Abstract

A trajectory isomorphism between the two Newtonian fixed center problem in the sphere and two associated planar two fixed center problems is constructed by performing two simultaneous gnomonic projections in $S^2$. This isomorphism converts the original quadratures into elliptic integrals and allows the bifurcation diagram of the spherical problem to be analyzed in terms of the corresponding ones of the planar systems. The dynamics along the orbits in the different regimes for the problem in $S^2$ is expressed in terms of Jacobi elliptic functions.

Equations

$$X^2+Y^2+Z^2=R^2$$

$$\mathcal{U}({\theta }_1,{\theta }_2)=-\frac{{\gamma }_1}{R}\mathrm{cotan}{\theta }_1-\frac{{\gamma }_2}{R}\mathrm{cotan}{\theta }_2,$$

$$\mathcal{H}=\frac{1}{2R^2}\left(L_X^2+L_Y^2+L_Z^2\right)-\frac{1}{R}\left(\frac{{\gamma }_1(\sigma Z+\bar{\sigma }X)}{\sqrt{R^2-(\sigma Z+\bar{\sigma }X{)}^2}}+\frac{{\gamma }_2(\sigma Z-\bar{\sigma }X)}{\sqrt{R^2-(\sigma Z-\bar{\sigma }X{)}^2}}\right),$$

$$\Omega =\frac{1}{2R^2}\left(L_X^2+{\sigma }^2{L}_Y^2\right)-\frac{\sigma }{R}\left(\frac{{\gamma }_{1}Z}{\sqrt{R^2-(\sigma Z+\bar{\sigma }X{)}^2}}+\frac{{\gamma }_{2}Z}{\sqrt{R^2-(\sigma Z-\bar{\sigma }X{)}^2}}\right).$$

$$\mathcal{U}({\theta }_1,{\theta }_2)=-\frac{({\gamma }_1+{\gamma }_2)\sin \frac{{\theta }_1+{\theta }_2}{2}\cos \frac{{\theta }_1+{\theta }_2}{2}+({\gamma }_1-{\gamma }_2)\sin \frac{{\theta }_2-{\theta }_2}{2}\cos \frac{{\theta }_2-{\theta }_1}{2}}{R\left({\sin }^2\frac{{\theta }_1+{\theta }_2}{2}-{\sin }^2\frac{{\theta }_2-{\theta }_1}{2}\right)}$$

$$U=\sin \frac{{\theta }_1+{\theta }_2}{2},V=\sin \frac{{\theta }_2-{\theta }_1}{2};-\bar{\sigma }<V<\bar{\sigma },\bar{\sigma }<U<1$$

$$X=\frac{R}{\bar{\sigma }}UV,Y^2=\frac{R^2}{{\sigma }^2{\bar{\sigma }}^2}(U^2-{\bar{\sigma }}^2)({\bar{\sigma }}^2-V^2),Z^2=\frac{R^2}{{\sigma }^2}(1-U^2)(1-V^2)$$

$${\mathcal{U}}_{+}(U,V)=-\frac{1}{R(U^2-V^2)}\left(({\gamma }_1+{\gamma }_2)U\sqrt{1-U^2}+({\gamma }_1-{\gamma }_2)V\sqrt{1-V^2}\right),$$

$${\mathcal{U}}_{-}(U,V)=-\frac{1}{R(U^2-V^2)}\left(-({\gamma }_1+{\gamma }_2)U\sqrt{1-U^2}+({\gamma }_1-{\gamma }_2)V\sqrt{1-V^2}\right).$$

$${\mathcal{H}}_{\pm }=\frac{1}{2R^2(U^2-V^2)}\left((U^2-{\bar{\sigma }}^2)(1-U^2)p_U^2+({\bar{\sigma }}^2-V^2)(1-V^2)p_V^2\right)+{\mathcal{U}}_{\pm }(U,V).$$

$${\mathcal{H}}_{\pm }\left(\frac{\partial S}{\partial U},\frac{\partial S}{\partial V},U,V\right)+\frac{\partial S}{\partial t}=0$$

$${\mathcal{H}}_{\pm }\to \frac{{\gamma }_1+{\gamma }_2}{R}{\mathcal{H}}_{\pm },t\to \frac{\sqrt{R^3}}{\sqrt{{\gamma }_1+{\gamma }_2}}t,p_{U,V}\to \sqrt{R({\gamma }_1+{\gamma }_2)}{p}_{U,V}$$

$$\gamma =\frac{{\gamma }_2}{{\gamma }_1+{\gamma }_2},$$

$$S_{\pm }(t;U,V)$$

$$\mathcal{G}=\mathcal{H}-\Omega .$$

$$\frac{dU}{d\varsigma }$$

$$\frac{dV}{d\varsigma }$$

$$\frac{dU}{d\varsigma }$$

$$\frac{dV}{d\varsigma }$$

$$\mathcal{U}(X,Y,Z)=-\frac{1}{R}\left(\frac{{\gamma }_1(\sigma Z+\bar{\sigma }X)}{\sqrt{R^2-(\sigma Z+\bar{\sigma }X{)}^2}}+\frac{{\gamma }_2(\sigma Z-\bar{\sigma }X)}{\sqrt{R^2-(\sigma Z-\bar{\sigma }X{)}^2}}\right).$$

$$\ddot{X}=-\frac{\partial \mathcal{U}}{\partial X}+\lambda X,\ddot{Y}=-\frac{\partial \mathcal{U}}{\partial Y}+\lambda Y,\ddot{Z}=-\frac{\partial \mathcal{U}}{\partial Z}+\lambda Z,$$

$${\Pi }_{+}:x=\frac{R}{Z}X,y=\frac{R}{Z}Y;{\Pi }_{-}:x=\frac{R}{-Z}X,y=\frac{R}{-Z}Y.$$

$$x_1\equiv x,x_2\equiv \frac{y}{\sigma }.$$

$$x_1^{\prime \prime }(\tau )=-\frac{\partial {\mathcal{V}}_{+}}{\partial x_1},x_2^{\prime \prime }(\tau )=-\frac{\partial {\mathcal{V}}_{+}}{\partial x_2}$$

$${\mathcal{V}}_{+}(x_1,x_2)=-\frac{{\alpha }_1}{\sqrt{{\left(x_1-a\right)}^2+x_2^2}}-\frac{{\alpha }_2}{\sqrt{{\left(x_1+a\right)}^2+x_2^2}}$$

$$d\tau =\frac{R^2}{Z^2}dt$$

$$x_1^{\prime \prime }(\tau )=-\frac{\partial {\mathcal{V}}_{-}}{\partial x_1},x_2^{\prime \prime }(\tau )=-\frac{\partial {\mathcal{V}}_{-}}{\partial x_2}$$

$${\mathcal{V}}_{-}(x_1,x_2)=\frac{{\alpha }_2}{\sqrt{{\left(x_1-a\right)}^2+x_2^2}}+\frac{{\alpha }_1}{\sqrt{{\left(x_1+a\right)}^2+x_2^2}}$$

$${\mathrm{h}}_{\pm }=\frac{1}{2}\left(p_1^2+p_2^2\right)+{\mathcal{V}}_{\pm }(x_1,x_2).$$

$$x_i\to a{x}_i,p_i\to \frac{\sqrt{{\alpha }_1+{\alpha }_2}}{\sqrt{a}}{p}_i,\tau \to \frac{\sqrt{a^3}}{\sqrt{{\alpha }_1+{\alpha }_2}}\tau ,{\mathrm{h}}_{\pm }\to \frac{{\alpha }_1+{\alpha }_2}{a}{\mathrm{h}}_{\pm };\alpha =\frac{{\alpha }_2}{{\alpha }_1+{\alpha }_2}=\gamma$$

$$u=\frac{\sqrt{(x_1+1{)}^2+x_2^2}+\sqrt{(x_1-1{)}^2+x_2^2}}{2},v=\frac{\sqrt{(x_1+1{)}^2+x_2^2}-\sqrt{(x_1-1{)}^2+x_2^2}}{2}$$

Full text

MSC2010 numbers: 70F15, 70H06.

Keywords: spherical two-center problem, separation of variables, spheroconical coordinates, elliptic coordinates.

1 Introduction

The two fixed center problem on the two-dimensional sphere goes back to Killing [21], and in modern times to Kozlov and Harin [22], who proved the separability of the problem, thus its integrability, in sphero-conical coordinates. These coordinates on $S^{2}$ were introduced by Liouville [25], and independently by Neumann [27] in one of the first examples of dynamics in spaces of constant curvature, and they are closely related to elliptic coordinates, in fact, the first system of coordinates plays a role in the dynamics on the sphere completely similar to the second system with respect to the Euclidian case. Integrability and Hamilton-Jacobi separability in sphero-conical coordinates has been constructed for different physical systems defined on the sphere, see, for instance, [11]. In particular, the authors analyzed in this context the Neumann problem and the Garnier system on $S^{2}$ in order to study solitary waves in one-dimensional nonlinear $S^{2}$-sigma models, see [7] and [8]. A detailed historical review of several systems defined in spaces of constant curvature, including open problems, has been recently performed in [14], where a precise bibliography is contained.

The two fixed center problem on the sphere is the superposition of two Kepler problems on $S^{2}$. An explicit expression for the second constant of motion for this problem and also for some generalizations was given in [26, 12]. In [13] Borisov and Mamaev, inspired by a previous work of Albouy and Stuchi [1, 2], established a trajectory isomorphism between the orbits lying in the half-sphere that contains the two attractive centers and the bounded orbits of an associated planar system of two attractive centers.

This isomorphism is constructed using a gnomonic projection from $S^{2}$ to the tangent plane at the middle point between the centers. The idea of relating planar and spherical problems in the framework of general force fields by considering the gnomonic projection goes back to Appell [9, 10], who also explained in this context the previous results of Serret [30] about the one fixed center problem on the sphere. Almost one century later Higgs [19] rediscovered these techniques. More recently, Albouy has developed these ideas and extended their scope to a general projective dynamics [1, 2, 3, 4, 5].

In this work we extend the results of [13] to the whole sphere, i.e., we establish a trajectory isomorphism between the complete set of orbits of the original problem and the corresponding one to two associated planar problems. The underlying idea is to identify each trajectory crossing the equator with the conjunction of two planar unbounded orbits, one of the two attractive center problem and another one for the system of two repulsive centers.

The extended trajectory isomorphism allows us to understand the bifurcation diagram of the spherical problem, previously analyzed in [31, 32, 33], as the superposition of the diagrams corresponding to the planar problem of two attractive centers [34] plus those associated to two repulsive centers [29]. This point of view permits the identification of the domains of allowable motions and the different cases for orbits, already described in [33], in terms of their partners for the planar diagrams. Finally, these identifications lead in a natural way to the determination of the quadratures (elliptic integrals) for the parametric equations of the orbits in terms of a local time. Using adequately the properties of elliptic integrals, the quadratures are inverted to obtain explicit formulas in terms of Jacobi elliptic functions for the different types of orbits of the spherical problem.

The existence of closed orbits in the sphere is guaranteed for the case of commensurability between the involved periods of the Jacobi functions.

The structure of the paper is as follows: The problem is presented in Section 2 using sphero-conical coordinates on $S^{2}$. In Section 3 the extended trajectory isomorphism is defined, and thus the quadratures are converted into elliptic integrals. The bifurcation diagram for the spherical problem is constructed from the diagrams of the two associated planar problems in Section 4. Finally, in Section 5, the process of inversion of elliptic integrals in $S^{2}$ is detailed, and general expressions for the solutions are shown.

The complete list of analytical expressions for the different types of orbits in $S^{2}$, in terms of a local time, is included in the Appendix, together with a gallery of figures for all the significative cases.

2 The two Newtonian centers problem in $S^{2}$

We consider the problem of a unit mass lying on the sphere $S^{2}$ of radius $R$, viewed as immersed in the Euclidean space ${\mathbb{R}}^{3}$ with Cartesian coordinates $(X,Y,Z)$:

[TABLE]

under the influence of the superposition of two Kepler potentials on $S^{2}$, i.e., the potential:

[TABLE]

where $\theta_{1}$ and $\theta_{2}$ denote the great circle angles between the location of the centers $F_{1}$ and $F_{2}$, see Fig. 1, and a given point $P$ on $S^{2}$, in such a way that $R\theta_{1}$ and $R\theta_{2}$ are the orthodromic distances from $F_{1}$ and $F_{2}$ to $P$, respectively. $\gamma_{1}$ and $\gamma_{2}$ are the strengths of the centers, where we have considered $0<\gamma_{2}\leq\gamma_{1}$, i.e., the test mass feels the presence of two attractive centers in $F_{1}$ and $F_{2}$, and correspondingly two repulsive centers at their antipodal points $\bar{F}_{1}$ and $\bar{F}_{2}$. Without loss of generality, the chosen points, notation and orientation are shown in Fig. 1. Thus, the Cartesian coordinates of $F_{1}$ and $F_{2}$ are $(R\sin\theta_{f},0,R\cos\theta_{f})=(R\bar{\sigma},0,R\sigma)$ and $(-R\sin\theta_{f},0,R\cos\theta_{f})=(-R\bar{\sigma},0,R\sigma)$, respectively. Parameters $\sigma=\cos\theta_{f}$ and $\bar{\sigma}=\sin\theta_{f}$ have been introduced in order to alleviate the notation.

This problem is completely integrable, see, e.g., [21, 22], there exist two constants of motion, the Hamiltonian:

[TABLE]

where $\vec{L}=\vec{X}\times\vec{P}$, $\vec{P}=(P_{X},P_{Y},P_{Z}),\vec{X}=(X,Y,Z)$, and the second invariant:

[TABLE]

This constant of motion is slightly different, but equivalent to the invariant obtained by Borisov and Mamaev in [26, 12]. Potential (1) can be rewritten as

[TABLE]

in such a way that it is natural to introduce an á la Euler version of sphero-conical coordinates on $S^{2}$, i.e.,

[TABLE]

Coordinate lines with fixed $U$ or $V$ resemble “spherical ellipses” or “spherical hyperbolas”, respectively, with foci $F_{1}$ and $F_{2}$, with the understanding that “spherical hyperbolas” are nothing more than “spherical ellipses” with respect to the pair of foci $\bar{F}_{1}$ and $F_{2}$ or $F_{1}$ and $\bar{F}_{2}$.

The change of coordinates

[TABLE]

is a four-to-one map because of the ambiguities in the signs of $Y$ and $Z$. Obviously, coordinates $U$ and $V$ are dimensionless.

Potential (4) is written in these sphero-conical coordinates with two different expressions depending on the hemisphere that it is considered. For $S^{2}_{+}=\{(X,Y,Z)\in S^{2},\ Z\geq 0\}$, we have

[TABLE]

whereas in $S^{2}_{-}=\{(X,Y,Z)\in S^{2},\ Z\leq 0\}$ the potential reads:

[TABLE]

Note that both expressions (6) and (7) coincide on the Equator $Z=0$, or $U=1$. Thus, Hamiltonian (2) has also to be split into two different expressions:

[TABLE]

The Hamilton-Jacobi equations coming from (8)

[TABLE]

are separable into two ordinary differential equations if we look for solutions of the form: $S_{\pm}(t;U,V)=S_{t}(t)+S_{U\pm}(U)+S_{V}(V)$. Introducing nondimensional variables

[TABLE]

and defining the parameter

[TABLE]

the complete solution of (9) is

[TABLE]

where $H$ and $G$ are the values of the constants of motion: ${\cal H}=H$, ${\cal G}=G$; ${\cal G}$ is the separation constant related to $\Omega$ and ${\cal H}$, (3) and (2), by the expression

[TABLE]

Given the local time $\varsigma$ by $d\varsigma=\frac{dt}{U^{2}-V^{2}}$, the standard separation procedure leads us to the first-order equations

[TABLE]

for the problem in the Northern hemisphere $S_{+}^{2}$, and

[TABLE]

for the Southern $S_{-}^{2}$ one.

A direct attack to the quadratures involved looks apparently cumbersome and, as far as we know, they are not solved in the literature. Nevertheless, some of the qualitative and topological properties of these orbits have been analyzed in [31, 32, 33].

3 Trajectory isomorphism between the spherical and two different planar problems

Following Borisov & Mamaev [13], we go back to Cartesian coordinates $(X,Y,Z)$ where the potential (1) can be written as

[TABLE]

The corresponding Newton equations for this problem are

[TABLE]

where the dots represent derivatives with respect to the physical (dimensional) time $t$ and $\lambda$ is the Lagrange multiplier. In [13] it was proved that the gnomonic projection from $S_{+}^{2}$ to the tangent plane $\Pi_{+}$ at the point $(0,0,R)$, together with a linear transformation in $\Pi_{+}$, maps Newton equations (15) to the Newton equations of an associated problem of two attractive centers in ${\mathbb{R}}^{2}$.

Here, we shall also consider simultaneously another gnomonic projection, from $S_{-}^{2}$ to the tangent plane $\Pi_{-}$, at $(0,0,-R)$. The projected coordinates $(x,y)$ are given in the two planes by

[TABLE]

for $Z\neq 0$. The Equator is mapped to the infinity in both $\Pi_{+}$ and $\Pi_{-}$ planes.

We will use throughout the paper the following criteria: uppercase letters describe magnitudes and variables specifically defined in the sphere, whereas lowercase letters will be associated to the planar cases. Following [13], we perform in $\Pi_{+}$ the linear transformation:

[TABLE]

The Newton equations (15) for potential (14) are rewritten in transformed projected coordinates $(x_{1},x_{2})$ on $\Pi_{+}$ as

[TABLE]

[TABLE]

where the primes denote the derivative with respect to a new time $\tau$ defined by

[TABLE]

and we have introduced the parameters: $a=R\frac{\bar{\sigma}}{\sigma}$, $\alpha_{1}=\frac{\gamma_{1}}{\sigma^{2}}$ and $\alpha_{2}=\frac{\gamma_{2}}{\sigma^{2}}$.

In a similar way, the Newton equations (15) restricted to $S_{-}^{2}$ can be projected into $\Pi_{-}$ using (16) and, after applying transformation (17), the equations

[TABLE]

[TABLE]

are obtained.

Note that ${\cal V}_{-}(x_{1},x_{2})$ in $\Pi_{-}$ is nothing more than the planar potential of two repulsive centers, where the roles of the points $(\pm a,0)$, and thus the strengths of the centers in modulus, are interchanged with respect to the attractive potential ${\cal V}_{+}(x_{1},x_{2})$ in $\Pi_{+}$.

Thus, while the restriction of the Newton equations (15) to the Northern hemisphere $S_{+}^{2}$ is equivalent to the Newton equations (18) for a planar problem of two attractive centers with potential (19), the restriction to the Southern hemisphere $S_{-}^{2}$ is tantamount to a planar problem of two repulsive centers with potential (20).

Bounded orbits of the attractive planar problem are in a one-to-one correspondence with the orbits of the spherical problem that lie in $S_{+}^{2}$. However, trajectories of the spherical problem crossing the equator have to be described in this projected picture by two pieces: an unbounded orbit of the attractive planar problem (19) in $\Pi_{+}$ plus an (unbounded) orbit of the repulsive planar problem (20) in $\Pi_{-}$, corresponding to the parts of the orbit belonging to $S_{+}^{2}$ and $S_{-}^{2}$, respectively.

It is possible to describe in a compact form the two associated planar problems in $\Pi_{+}$ and $\Pi_{-}$, respectively, by the Hamiltonians

[TABLE]

It is adequate again to use nondimensional variables

[TABLE]

and to introduce “radial”, $u$, and “angular”, $v$, elliptic (Euler) coordinates in ${\mathbb{R}}^{2}$:

[TABLE]

[TABLE]

in such a way that the Hamiltonians (21) are written in terms of these coordinates as

[TABLE]

i.e., two standard Liouville separable systems in elliptic coordinates. It is straightforward to construct the associated first-order equations with respect to the local time $\zeta=\zeta(\tau)$ defined by

[TABLE]

and we finally obtain the following equations in the $\Pi_{+}$ plane:

[TABLE]

which solve the original problem in $S^{2}_{+}$. Correspondingly, for the Southern case we obtain in the $\Pi_{-}$ plane:

[TABLE]

where the constants of motion take the values ${\mathrm{h}}_{+}=h$ and ${\mathrm{g}}_{+}=g$ for the energy and the separation constant in $\Pi_{+}$, respectively, and ${\mathrm{h}}_{-}=\tilde{h}$ and ${\mathrm{g}}_{-}=\tilde{g}$ in $\Pi_{-}$. The quadratures involved in equations (22) and (23) are of elliptic type, and thus expressible in terms of the Jacobi elliptic functions.

It is possible to synthesize the chain of maps leading from the original problem in the sphere to the pair of planar two center problems (22) and (23) in a unique one-to-one transformation of coordinates in $S^{2}$, from sphero-conical $(U,V)$ to planar elliptic $(u,v)$, as follows:

[TABLE]

together with an equivalence, up to a constant factor, between the nondimensional local time $\varsigma$ of the spherical problem and the nondimensional local time $\zeta$ for the associated planar problems:

[TABLE]

The equator $Z=0$, or $U=1$, of $S^{2}$ is mapped by (24) into the point of infinity in the coordinate $u$.

Equation (25) establishes that $\zeta$ can be simultaneously seen as the local time for both problems, remembering its different meaning when regarded from $S^{2}$, local time associated with the “physical” time $t$ in the sphere, or from $\Pi_{\pm}$, local time corresponding to the projected (nonphysical) time $\tau$ in the planes.

Thus, (24) and (25) map directly the first-order equations (10, 11) in $S^{2}_{+}$ to equations (22) in $\Pi_{+}$, and (12, 13) in $S^{2}_{-}$ to (23) in $\Pi_{-}$ via the identifications

[TABLE]

[TABLE]

It is remarkable that in this projected picture the role of the planar energies $h$ and $\tilde{h}$ is played, up to a factor, by the projection of the second constant of motion $\Omega$ and not by the projection of the spherical Hamiltonian. This fact is a consequence of the behavior of constants of motion under central projections, as was explained in [4], see also [5].

Consequently, the transformation (24) establishes that fixing in $S^{2}$ a negative value of the constant of motion $\Omega$, the orbits of the problem lie in the $S^{2}_{+}$ hemisphere and are in a one-to-one correspondence with the bounded orbits, $h=\frac{\bar{\sigma}}{\sigma}\Omega<0$, of the planar attractive system in the $\Pi_{+}$ plane. Thus, coordinate $u$ is bounded for $\Omega<0$. However, if $\Omega\geq 0$, orbits cross the equator of $S^{2}$, and thus the portions of the orbits belonging to $S^{2}_{+}$ are described by equations (22) with planar energy $h\geq 0$, unbounded planar orbits in the attractive problem in $\Pi_{+}$, whereas the portions lying in the Southern hemisphere $S^{2}_{-}$ are determined by equations (23) with $\tilde{h}>0$, i.e., unbounded planar orbits of the repulsive problem in $\Pi_{-}$. In this case $u$ is unbounded.

4 The bifurcation diagrams

The isomorphic transformation (24) allows us to analyze the bifurcation diagram in $S^{2}$ starting from the bifurcation diagrams of the two associated planar problems. From this point of view a global bifurcation diagram for the spherical problem will be constructed out of the diagrams of two planar centers, see [34] and [29], respectively, attractive in $\Pi_{+}$ and repulsive in $\Pi_{-}$ and strengths interchanged. Thus, the results explained in [32, 33] will now be understood from a different perspective.

Both in $\Pi_{+}$ and $\Pi_{-}$ planes, i.e., the images of the North $S^{2}_{+}$ and South $S^{2}_{-}$ hemispheres, we rewrite (22) and (23) in terms of the ramification points:

[TABLE]

[TABLE]

In Figure 2, plotted for $\alpha=1/3$, we observe the bifurcation diagrams corresponding to the attractive and repulsive planar problems in $\Pi_{+}$ Fig. 2a) and $\Pi_{-}$ Fig. 2b), respectively, with strengths $\alpha_{1}$, $\alpha_{2}$, and $\tilde{\alpha}_{1}=-\alpha_{2}$, $\tilde{\alpha}_{2}=-\alpha_{1}$. Critical curves in both $\{h,g\}$ and $\{\tilde{h},\tilde{g}\}$ planes are determined by the existence of double roots in (26) and (27), see [34, 29], and shadowed areas in the diagrams are zones where motion is classically forbidden, i.e., velocities and/or momenta are imaginary.

The allowed motions in the $\Pi_{+}$-plane are of two types: (1) If $h<0$, orbits are bounded and are usually labeled as $\{t_{s},t_{s^{\prime}},t_{l},t_{p}\}$, for satellitary, lemniscatic and planetary ones, see [34]. (2) If $h\geq 0$, see [29], unbounded orbits occur, standardly labeled as $\{t_{1},t_{2},t_{3},t_{4},t_{5}\}$. Separatrices between bounded and unbounded motions live in the $\{h=0\}$ straight line.

In the $\Pi_{-}$ plane a similar, but simpler picture is found, see Fig. 2b). On the $\tilde{h}>0$ upper half-plane unbounded orbits exist in five different classes, labeled as $\{t^{\prime}_{1},t^{\prime}_{2},t^{\prime}_{3},t^{\prime}_{4},t^{\prime}_{5}\}$. In this case the line $\{\tilde{h}=0\}$ does not accommodate separatrices, but rather it corresponds to a limiting behavior of unbounded zero energy orbits reached from $\tilde{h}>0$.

The bifurcation diagram for the complete problem in $S^{2}$, Figure 3, can now be constructed from the planar ones using transformations (24) and (25). Two morphisms are induced: (1) Orbits in $\Pi_{+}$ are mapped to orbits in $S^{2}_{+}$ identifying the invariants as follows: $h=\frac{\bar{\sigma}}{\sigma}\Omega$ and $g=\frac{\sigma}{\bar{\sigma}}G$. (2) Orbits in $\Pi_{-}$ are mapped to orbits in $S^{2}_{-}$ if the invariants are translated to: $\tilde{h}=\frac{\bar{\sigma}}{\sigma}\Omega$ and $\tilde{g}=\frac{\sigma}{\bar{\sigma}}G$. The global bifurcation diagram in $S^{2}$ is thus displayed on the $\{\frac{\bar{\sigma}}{\sigma}\Omega,\frac{\sigma}{\bar{\sigma}}G\}$-plane.

Moreover, the lower half-plane of Figure 3, $\frac{\bar{\sigma}}{\sigma}\Omega<0$, is mapped one-to-one with the lower half-plane of the problem of two attractive centers in $\Pi_{+}$, Fig. 2a), as it was shown in [13], orbits lying only in $S_{+}^{2}$ are in a bijective correspondence with bounded orbits in $\Pi_{+}$. However, with an initial condition fixed, each point $\left(\frac{\bar{\sigma}}{\sigma}\Omega,\frac{\sigma}{\bar{\sigma}}G\right)$ in the upper half-plane of the global diagram represents an orbit that crosses the equator of $S^{2}$, and thus is mapped by (24) and (25) to the union of an unbounded orbit in $\Pi_{+}$ and another one in $\Pi_{-}$, with equal planar energies: $h=\tilde{h}=\frac{\bar{\sigma}}{\sigma}\Omega$.

Critical curves in Figure 3 are inherited from the corresponding ones in planar diagrams:

• •

Double roots in equations (26, 27) for the “radial” variable arise in the: blue straight line: ${\cal L}_{1}^{2}=\left\{\frac{\bar{\sigma}}{\sigma}\Omega-\frac{\sigma}{\bar{\sigma}}G-1=0\right\}$, red straight line: ${\cal L}_{1}^{1}=\left\{\frac{\bar{\sigma}}{\sigma}\Omega-\frac{\sigma}{\bar{\sigma}}G+1=0\right\}$, and green hyperbola: ${\cal L}_{1}^{3}=\left\{4\Omega G+1=0\right\}$.

• •

Analogously, double roots for the “angular” variable produce the: dashed blue straight line: ${\cal L}_{\gamma}^{1}=\left\{\frac{\bar{\sigma}}{\sigma}\Omega-\frac{\sigma}{\bar{\sigma}}G-(1-2\gamma)=0\right\}$, dashed red straight line: ${\cal L}_{\gamma}^{2}=\left\{\frac{\bar{\sigma}}{\sigma}\Omega-\frac{\sigma}{\bar{\sigma}}G+(1-2\gamma)=0\right\}$, and dashed green hyperbola: ${\cal L}_{\gamma}^{3}=\left\{4\Omega G+(1-2\gamma)^{2}=0\right\}$.

Orbits with $\Omega<0$ are naturally labeled with the inherited standard notation for bounded motion in the planar associated problem in $\Pi_{+}$. The branching points $u_{1},u_{2}$ and $v_{1},v_{2}$, understood as functions of $\Omega$ and $G$, allow us to specify the qualitative features of these orbits in $S_{+}^{2}$:

• •

Planetary orbits ($t_{p}$). There are two analytical possibilities that lead to the same type of orbits:

[TABLE]

In both cases the bounds $u=u_{1}$ and $u=u_{2}$ represent two caustics for these orbits, i.e., two “spherical ellipses” in the Northern hemisphere $S_{+}^{2}$, see Fig. 6a, that confine the planetary motion of these “circumbinary” orbits.

• •

Lemniscatic orbits ($t_{l}$). Analogously, there exist two possibilities:

[TABLE]

A unique caustic, $u=u_{2}$, appears in this case. The orbits describe a lemniscatic motion around the two centers in $S_{2}^{+}$. See Fig. 6b.

• •

Satellitary orbits ($t_{s}$): Each point $\left(\frac{\bar{\sigma}}{\sigma}\Omega,\frac{\sigma}{\bar{\sigma}}G\right)$ of this region in Figure 3 represents two possible orbits:

[TABLE]

around the stronger center, limited by the caustics: $u=u_{2}$ and $v=v_{2}$, and:

[TABLE]

around the weaker center, bounded by: $u=u_{2}$ and $v=v_{1}$. See Figure 6g.

• •

Satellitary orbits ($t_{s^{\prime}}$) around the stronger center:

[TABLE]

For this situation the motion is limited by the caustics: $u=u_{2}$ and $v=v_{2}$, see Figure 6c.

For $\Omega>0$, it is possible to extend the standard terminology, Planetary ($t_{p}$), Lemniscatic ($t_{l}$) and Satellitary ($t_{s^{\prime}}$), to the orbits that cross the equator, but have a behavior analogous to the corresponding cases restricted to the Northern hemisphere. However, two completely new types of orbits arise. There are two zones of admissible motion without partners between the orbits with $\Omega<0$, which we will call Dual Satellitary ($t_{ds}$) and Meridian Planetary ($t_{mp}$) orbits, taking into account its qualitative features.

Branching points are now identified by $\tilde{u}_{1}=-u_{2}$, $\tilde{u}_{2}=-u_{1}$ and $\tilde{v}_{1}=v_{1}$, $\tilde{v}_{2}=v_{2}$, because $h=\tilde{h}=\frac{\bar{\sigma}}{\sigma}\Omega$ and $g=\tilde{g}=\frac{\sigma}{\bar{\sigma}}G$ in order to glue continuously the two orbit pieces on the Northern and Southern hemispheres at the equator.

• •

Planetary orbits ($t_{p}$): The orbits in $S^{2}$ are composed by two pieces:

[TABLE]

Note that the limit $u\to\infty$ in both cases is nothing more than $U\to 1$, and thus the map (24) applies two unbounded curves to a finite one that crosses the equator of $S^{2}$. The Northern pieces present the caustic: $u=u_{2}$, whereas the Southern ones are limited by the “spherical ellipse”: $u=\tilde{u}_{2}$. The motion is confined between these curves in a planetary way and can be seen as the natural continuation of the $t_{p}$ orbits in $S_{+}^{2}$ with $\Omega<0$. See Figure 6d.

• •

Lemniscatic orbits ($t_{l}$): Analogously, there are two parts:

[TABLE]

in such a way that there are no caustics in $S_{+}^{2}$ and one in $S_{-}^{2}$: $u=\tilde{u}_{2}$. We find again a natural resemblance between these orbits and their partners in the $\Omega<0$ case. See Figure 6e.

• •

Satellitary orbits ($t_{s^{\prime}}$):

[TABLE]

The caustics are now: $u=\tilde{u}_{2}$ in $S_{-}^{2}$, and $v=v_{1}=\tilde{v}_{1}$ in the two hemispheres. See Figure 6f.

• •

Dual Satellitary orbits ($t_{ds}$):

[TABLE]

The $t_{ds}$ orbits present a behavior delimited by the two caustics: $v=v_{1}=\tilde{v}_{1}$ and $v=v_{2}=\tilde{v}_{2}$ in $S^{2}$, and: $u=\tilde{u}_{2}$ in the Southern hemisphere. Thus, the orbits pass between the two centers in $S_{+}^{2}$, but do not reach the South Pole. See Figure 6h.

• •

Meridian Planetary orbits ($t_{mp}$):

[TABLE]

The situation is similar to the $t_{ds}$ case, but now only the two “angular” caustics are allowable. Thus, the orbits complete the passing between the centers not only in $S_{+}^{2}$ but also in $S_{-}^{2}$. The $t_{mp}$ orbits resemble the planetary ones interchanging the surrounded centers. See Figure 6i.

Finally, the analysis should be completed with the case $\Omega=0$ whose orbits lie in the $S^{2}_{+}$ hemisphere. These can be easily described as the limit $\Omega\rightarrow 0$ in the $\Omega<0$ case. The caustic $u=u_{2}$ for the $t_{p}$, $t_{l}$ and $t_{s^{\prime}}$ orbits becomes $u_{2}\rightarrow\infty$, and thus $U(u_{2})\rightarrow 1$, i.e., the equator $Z=0$ of $S^{2}$. Consequently, the motions are completely similar to the corresponding ones in $S_{+}^{2}$, but now bounded by the equator.

5 Evaluation of the quadratures, inversion of the elliptic integrals

The resolution of the problem of two fixed centers has a long history that, apart from the original results of Euler [17, 18], includes the works of Lagrange [23], Legendre [24], Jacobi [20], Liouville [25] etc., see [28] and references therein. In more recent times Alexeev [6] has given a detailed qualitative analysis of the planar problem. Explicit analytical expressions determining the motion along the orbits are obtained by applying standard procedures that require the inversion of elliptic integrals, see, for instance, [15, 35]. The quadratures solving the two pairs of uncoupled ODEs (26) and (27) have been thoroughly discussed by several authors, see [28] and references therein, see also [16].

We shall briefly report here on the processes of quadrature evaluation/elliptic integral inversion in the context of the spherical problem, keeping in mind that the variables $(u,v)$, which appear in equations (26) and (27), should be regarded as coordinates in $S^{2}_{+}$ and $S^{2}_{-}$ through the map transformation (24), as it has been explained in the previous sections.

There are two distinctly different situations for the $\Omega<0$ or $\Omega>0$ ranges:

$\bullet$ $\Omega<0$. In this case the inversion of the elliptic integrals appearing in equations (26) is standard, we will detail only the planetary case as an example.

The range for the $u$-variable in (26) (left) is: $u_{1}<u<u_{2}$, and thus the curves: $u=u_{1}$ and $u=u_{2}$, $\forall v\in(-1,1)$, determine the two caustics. The quadrature solving the $u$-equation in (26) is

[TABLE]

where the initial condition $u(0)=u_{0}$ is assumed. The elliptic integral of the first kind $I(u)$ in (40) can be inverted by performing the following change of variable $z\to s$, see [15] case 256:

[TABLE]

where $\mathop{\rm sn}\nolimits s$ denotes the Jacobi sinus function: $\mathop{\rm sn}\nolimits s\equiv\mathop{\rm sn}\nolimits(s|k_{u}^{2})$, $g_{u}$ and the elliptic modulus $k_{u}$ are defined in terms of the turning points as

[TABLE]

Formula (40) is thus simplified to become a linear relation between $s_{u}$ and the local time $\zeta$ which is easily inverted:

[TABLE]

with: $g_{u}s_{u_{0}}=I(u_{0})$. Finally, recalling the last change of variable, the explicit inversion of (40) is achieved:

[TABLE]

where $s_{u}$ is defined as a function of the local time $\zeta$, $s_{u}(\zeta)$, in equation (41). Alternatively, using the properties of Jacobi elliptic functions, $u(\zeta)$ can be rewritten in terms of the Jacobi function $\mathop{\rm dn}\nolimits$ in the simpler form:

[TABLE]

We stress, by writing the inequalities characterizing this type of orbits, that the analytic expression for $u(\zeta)$ appearing in formula (42) is compelled to live inside the $(u_{1},u_{2})$ interval.

The companion expression for $v(\zeta)$, for instance, in the planetary case $-1\,<\,v\,<\,1$, is given, after a completely analogous procedure, by the expressions

[TABLE]

with

[TABLE]

Applying transformation (24) to these expressions of $u(\zeta)$ and $v(\zeta)$ and replacing the results in (5), a complete description in Cartesian coordinates of planetary orbits in the Northern hemisphere is obtained.

Analogously, all the integrals $I(u)$ and $I(v)$ solving equations (26) in the different ranges of $u$ and $v$ compatible with $\Omega<0$ can be inverted by similar techniques. The ensuing analytic expressions are assembled in the Appendix. The $u(\zeta)$ and $v(\zeta)$ functions which, respectively, solve the $u$- and $v$-dynamics are smooth, bounded between turning points, and periodic with periods, respectively, $T_{u}\propto K(k_{u}^{2})$ and $T_{v}\propto K(k_{v}^{2})$, where $K(k^{2})$ is the complete elliptic function of the first kind. The trajectories in all these cases are bounded between caustics in $S^{2}_{+}$ and dense, except if the $u$- and $v$-periods are commensurable.

$\bullet$ $\Omega>0$. The procedure is more delicate in this case essentially because the trajectories complying with the ODE pair (26) reach the equator, whereas there is admissible motion governed by (27) that also reaches the equator coming from the Southern hemisphere. Therefore, it is convenient to investigate the inversion of the quadratures of the $u$-equations of both (26) and (27) in a global form. However, in the “angular”$v$-integrals there are no differences with respect to the $\Omega<0$ range.

Let us focus on planetary orbits. An orbit of this type in $S^{2}$ is described by two pieces: the portion belonging to $S_{+}^{2}$ is a solution of equations (26) in the ranges

[TABLE]

whereas for the $S_{-}^{2}$ piece we have equations (27) and the ranges

[TABLE]

The first quadrature in (26) for the “radial” variable

[TABLE]

can be inverted with a change of variable like that explained before in the $\Omega<0$ case. The solution is

[TABLE]

where

[TABLE]

A plot of $u(\zeta)$, see Figure 4a, shows several relevant features of $u(\zeta)$. First, the function (44) presents infinite poles located at the points where $\mathop{\rm dn}\nolimits^{2}s_{u}=\frac{u_{2}-1}{u_{2}+1}$. This is an expected result if one sees $u(\zeta)$ as a solution of the planar problem of two attractive centers with $h>0$ reinterpreting $\zeta$ as the local time of this planar problem; the trajectory goes to infinity in a finite interval of the local time. However, in the sphere $S^{2}$ the sphero-conical variable $U(\zeta)$, given by (24), is bounded but exhibits finite discontinuities and reaches its maxima on the equator $U=1$ at the poles of $u(\zeta)$, see Fig. 4b. Second, it is remarkable, and a priori unexpected, that both $u(\zeta)$ and $U(\zeta)$ take negative values. The subtle interpretation of this fact is the understanding that, given the inversion problem posed by (43), its solution $u(\zeta)\equiv u(s_{u})$ solves also the complementary problem: $y<u_{1}<-1<1<u_{2}$, i.e., the inversion problem of the elliptic integral:

[TABLE]

in such a way that the inverse function $y(s)$ verifies: $y(s)=u(s_{u}+K)$, where $K=K(k_{u}^{2})$. Thus, $u(s_{u})$ defined in equation (44) represents simultaneously the genuine $u$-“radial”positive solution, $u\in(u_{2},\infty)$, and the negative $y(s)$-“radial”solution with $y\in(-\infty,u_{1})$. Note that, according to the plot in Figure 4, these two solutions occur in consecutive intervals of the local time $\zeta$.

A direct search for the solution of equation (27) in $S^{2}_{-}$, where $\tilde{u}_{1}<-1<1<\tilde{u}_{2}<u$, requires the inversion of the elliptic integral in the following equation:

[TABLE]

Having in mind that $\tilde{u}_{1}=-u_{2}$, $\tilde{u}_{2}=-u_{1}$, we can write

[TABLE]

where the change of variable $z=-w$ has been performed. Thus, we conclude that the inversion of $\tilde{I}(u)$, i.e., the “radial” solution in $S^{2}_{-}$, is tantamount to the inversion of $I^{\prime}(-u)$ and consequently to minus the negative part of $u(s_{u})$ given in (44). Therefore, we represent the “radial” solution simultaneously in both $S^{2}_{+}$ and $S^{2}_{-}$ by simply taking the absolute value $|u(\zeta)|$ of the solution given in (44). Moreover, with this identification the function $|U(\zeta)|$ is smooth, i.e., the gluing at the equator of the Northern and Southern branches of the orbits is continuous and differentiable, see Figure 5, with respect to the local time $\zeta$.

This argument is valid also for the “radial” quadratures of the rest of different types of orbits that cross the equator. Thus, the general expression for the orbits in Cartesian coordinates over the sphere $S^{2}$, using (24) in (5), can be written in a compact form valid for all the types of orbits described in the previous section as

[TABLE]

Here, $\mathrm{sg}$ denotes the sign function and $(u(\zeta),v(\zeta))$ are the solutions of equations (22) or (23).

Explicit expressions for (45) in all the different regimes are written in the Appendix.

The quasiperiodicity properties of the functions (45) are inherited from the Jacobi elliptic functions through the functions $u(\zeta)$ and $v(\zeta)$: solutions (45) are products of periodic functions with different periods $T_{u}$ and $T_{v}$. Consequently, (45) will be periodic, and thus the orbits will be closed, only if $T_{u}$ and $T_{v}$ are commensurable, i.e., there exists $p,q\in\mathbb{N}^{*}$ such that

[TABLE]

otherwise the orbits will be dense inside the allowable region of $S^{2}$.

The periods $T_{u}$ and $T_{v}$ are proportional to $K(k_{u}^{2})$ and $K(k_{v}^{2})$, respectively, with a factor that depends on the concrete Jacobi functions involved in the respective expressions of $u(\zeta)$ and $v(\zeta)$. The search for a closed orbit, with the values of $p$, $q$ and $\Omega$ (or $G$) fixed, requires that the transcendental equation (46) be solved in the variable $G$ (alternatively $\Omega$). Explicit expressions for the periods and concrete examples of closed orbits for different values of $p$ and $q$ are collected in the Appendix.

Acknowledgments

The authors thank the Spanish Ministerio de Economía y Competitividad (MINECO) for financial support under grant MTM2014-57129-C2-1-P and the Junta de Castilla y León grant VA057U16.

Appendix. Explicit expressions for different types of orbits

The set of parameters that determines the problem is $R$, $\theta_{f}$, $\gamma_{1}$ and $\gamma_{2}$, but after defining nondimensional variables the strengths can be measured with only one relative quantity: $\gamma=\frac{\gamma_{2}}{\gamma_{1}+\gamma_{2}}$.

Our choice of integration constants to characterize the solutions (45) as functions of the nondimensional local time $\zeta$ introduced in (25) is: the two constants of motion $\Omega$ and $G$ and the two initial data $s_{u_{0}}$ and $s_{v_{0}}$. The dependence in $\Omega$ and $G$ is given implicitly through the values of the branching points:

[TABLE]

if $\Omega\neq 0$, and $u_{1}=\frac{\sigma}{\bar{\sigma}}G$, $v_{2}=\frac{\sigma}{\bar{\sigma}}\,\frac{-G}{(1-2\gamma)}$ for the $\Omega=0$ case.

Remember also that the following notation has been introduced throughout the paper:

[TABLE]

and so on for the rest of the Jacobi elliptic functions, where

[TABLE]

in such a way that the initial conditions are $s_{u_{0}}=s_{u}(0)$ and $s_{v_{0}}=s_{v}(0)$.

With all these considerations, the orbits for the two fixed centers problem in $S^{2}$ are:

$\Omega>0$: Orbits that cross the equator.

$\bullet$ Planetary orbits-$t_{p}$, see (35):

[TABLE]

where

[TABLE]

[TABLE]

$\bullet$ Lemniscatic orbits-$t_{l}$ (36):

[TABLE]

[TABLE]

[TABLE]

$\bullet$ Satellitary orbits-$t_{s^{\prime}}$ (37):

[TABLE]

[TABLE]

[TABLE]

$\bullet$ Dual Satellitary orbits-$t_{ds}$ (38):

[TABLE]

[TABLE]

[TABLE]

$\bullet$ Meridian Planetary orbits-$t_{mp}$ (39):

[TABLE]

[TABLE]

[TABLE]

$\bullet$ Taking into account the Jacobi functions involved in each type of solutions, the $u$- and $v$-periods for the different orbits with $\Omega>0$ are

[TABLE]

[TABLE]

[TABLE]

$\Omega<0$: Orbits that lie only in the Northern hemisphere.

$\bullet$ Planetary orbits-$t_{p}$ of type $1$, (28):

[TABLE]

[TABLE]

[TABLE]

$\bullet$ Planetary orbits-$t_{p}$ of type $2$, (29):

[TABLE]

[TABLE]

[TABLE]

$\bullet$ Lemniscatic orbits-$t_{l}$ of type $1$, (30):

[TABLE]

[TABLE]

[TABLE]

$\bullet$ Lemniscatic orbits-$t_{l}$ of type $2$, (31):

[TABLE]

[TABLE]

[TABLE]

$\bullet$ Satellitary orbits-$t_{s}$ in zone $1$, (32):

[TABLE]

[TABLE]

[TABLE]

$\bullet$ Satellitary orbits-$t_{s}$ in zone $2$, (33):

[TABLE]

[TABLE]

[TABLE]

$\bullet$ Satellitary orbits-$t_{s^{\prime}}$ (34):

[TABLE]

[TABLE]

[TABLE]

$\bullet$ The $u$- and $v$- periods in the case $\Omega<0$ are

[TABLE]

[TABLE]

[TABLE]

[TABLE]

$\Omega=0$: Orbits that lie in the Northern hemisphere bounded by the equator.

$\bullet$ Planetary orbits-$t_{p}$ in this case are   $-1<1<u_{1}<u\quad,\quad v_{2}<-1<v<1$.

[TABLE]

[TABLE]

[TABLE]

$\bullet$ Lemniscatic orbits-$t_{l}$:   $-1<u_{1}<1<u\quad,\quad v_{2}<-1<v<1$.

[TABLE]

[TABLE]

[TABLE]

$\bullet$ Satellitary orbits-$t_{s^{\prime}}$:   $-1<u_{1}<1<u\quad,\quad-1<v_{2}<v<1$.

[TABLE]

[TABLE]

[TABLE]

$\bullet$ Finally, the $u$- and $v$-periods for these orbits with $\Omega=0$ are

[TABLE]

[TABLE]

[TABLE]

Several orbits with $R=1$, $\gamma=\frac{1}{3}$ and $\theta_{f}=\frac{\pi}{6}$, one for each different situation, are represented in Figure 6. These orbits are dense in all the cases and they are depicted in the interval $\zeta\in[0,70]$.

We can see in Figure 7 several closed orbits of different types, with specification of the values of $p$, $q$ and initial data $s_{u_{0}}$ and $s_{v_{0}}$. In all the cases $p$, $q$ and the constant of motion $\Omega$ have been fixed, thus $G$ has been calculated by solving numerically equation (46).

[TABLE]

[TABLE]