Examples of Complete Solvability of 2D Classical Superintegrable Systems

Classical (maximal) superintegrable systems in $n$ dimensions are Hamiltonian systems with $2n-1$ independent constants of the motion, globally defined, the maximum number possible. They are very special because they can be solved algebraically. In this paper we show explicitly, mostly through examples of 2nd order superintegrable systems in 2 dimensions, how the trajectories can be determined in detail using rather elementary algebraic, geometric and analytic methods applied to the closed quadratic algebra of symmetries of the system. We treat a family of 2nd order degenerate systems: oscillator analogies on Darboux, nonzero constant curvature, and flat spaces, related to one another via contractions, and obeying Kepler's laws. Then we treat two 2nd order nondegenerate systems, an analogy of a caged Coulomb problem on the 2-sphere and its contraction to a Euclidean space caged Coulomb problem. In all cases the symmetry algebra structure provides detailed information about the trajectories. An interesting example is the occurrence of ''metronome orbits'', trajectories confined to an arc rather than a loop, which are indicated clearly from the structure equations but might be overlooked using more traditional methods. We also treat the Post-Winternitz system, an example of a classical 4th order superintegrable system that cannot be solved using separation of variables. Finally we treat a superintegrable system, related to the addition theorem for elliptic functions, whose constants of the motion are only rational in the momenta, a system of special interest because its constants of the motion generate a closed polynomial algebra. This paper contains many new results but we have tried to present most of the materials in a fashion that is easily accessible to nonexperts, in order to provide entr\'ee to superintegrablity theory.


Introduction
Classical and quantum mechanical Hamiltonian systems that can be solved explicitly, both algebraically and analytically, and with adjustable parameters, are relatively rare and highly prized. Famous classical examples are the anharmonic oscillator (Lissajous patterns) and Kepler systems (planetary orbits) and the Hohmann transfer for orbital navigation [7]. Famous quantum examples are the Coulomb system (energy levels of the hydrogen atom, leading to the periodic table of the elements), and the quantum isotropic oscillator. The solvability of these systems is related to their symmetry, not necessarily group symmetry. This higher order symmetry is captured by the concept of superintegrability. A natural classical Hamiltonian system (with the Hamiltonian as kinetic energy plus potential energy) on an n-dimensional Riemannian space is said to be (maximally) superintegrable if it admits the maximally possible 2n − 1 functionally independent constants of the motion globally defined (usually required to be polynomial or at least rational in the momenta). Similarly a quantum Hamiltonian system H = ∆ n +V (where ∆ n is a Laplace-Beltrami operator on an n-dimensional Riemannian manifold and V is a potential function) is (maximally) superintegrable if it admits 2n − 1 algebraically independent partial differential operators commuting with H. There is a rapidly growing literature concerning these systems, e.g., [2,3,4,6,8,9,10,11,13,20,28,30,31,32,33,38,39,40,42]. In this paper we consider only classical systems and n = 2 so all of our systems admit 3 independent constants of the motion. By taking Poisson brackets of the classical constants of the motion we generate a symmetry algebra, not necessarily a Lie algebra, which is never abelian. (This contrasts with integrable Hamiltonian systems which admit n constants of the motion in involution, so that the symmetry algebra is always abelian.) It is this nonabelian symmetry algebra and its structure that is responsible for the solvability of superintegrable systems.
To be more explicit, along a specific classical trajectory each constant of the motion L j takes a fixed value j , so the trajectory can be characterized as the intersection of 2n−1 hypersurfaces L j = j in the 2n-dimensional phase space. Thus in principle the path of the trajectory can be determined algebraically, though not how it is traced out in time. Since the Poisson bracket of two constants of the motion is again a constant of the motion, the nonabelian symmetry algebra gives us relationships between the symmetries. In the quantum case each bound state eigenspace of the Hamiltonian is invariant under the action of the symmetry algebra, so that a knowledge of the irreducible representations of the symmetry algebra gives useful information about the dimensions of the eigenspaces and of the eigenvalues themselves.
In this paper we illustrate the value of superintegrabilty by studying and solving several families of classical superintegrable systems. Second degree superintegrable systems are those whose generating symmetries are all polynomials in the momenta of degree ≤ 2. All of these systems have been classified for 2 dimensions (as have the systems with nondegenerate potentials in 3 dimensions) [5,19,23]. They occur on constant curvature spaces (with 3 Killing vectors, admitting the most symmetries), on the 4 Darboux spaces (admitting 1 Killing vector) and 6 families of Koenigs spaces (admitting no Killing vectors) [25,33]. Thus, after the constant curvature spaces, the Darboux spaces admit the most symmetries. The Darboux metrics an be written as D1 : ds 2 = 4x dx 2 + dy 2 , D2 : D3 : ds 2 = e x + 1 e 2x dx 2 + dy 2 , D4(b) : ds 2 = 2 cos 2x + b sin 2 2x dx 2 + dy 2 ,

An example of a Koenigs space is
We first treat some analogs of the harmonic oscillator. (A treatment of the Kepler system and its analog on the 2-sphere from this point of view is contained in Chapter 3 of the article [33]. ) We start by studying a system on the Darboux space D4(b) with a 1-parameter potential. Then by taking limits (contractions [15,24]) of this system we obtain systems on the Darboux space D3, on the Poincaré upper half plane, the 2-sphere (the Higgs oscillator [14]), and finally the isotropic oscillator on Euclidean space.
The second class of examples are nondegenerate (3-parameter potential) systems. The first, S7 in our listing [23], can be regarded as a caged version of an analog to the Coulomb potential on the 2-sphere. It contracts to the caged Coulomb system E16 in Euclidean space.
The preceding examples can also be studied analytically via separation of variables. However, the Post-Winternitz system [36] cannot, see also [29]. It is 4th degree superintegrable but nonseparable. However, we show that the classical orbits can be found exactly via superintegrabilty.
The last examples are different. They are separable in elliptic coordinates and can be derived via an action-angle construction. The usual action-angle construction of a superintegrable and separable system requires the addition theorem for trigonometric or hyperbolic functions and leads to polynomial superintegrability, e.g., [20]. The construction here uses the addition formula for elliptic functions [35]. It leads to classical systems where some of the constants are nontrivially rational in the canonical momenta. We present two examples where, however, the systems have polynomial symmetry algebras, so we consider them as superintegrable and worthy of study. The new feature here is that the symmetry algebra closes polynomially, even though the system is only rationally superintegrable. These are the first examples of such behavior known to us. Again we can determine the trajectories exactly.
With the exception of the elliptic systems, all of the superintegrable systems in this paper have been derived and classified before; we take the structure equations as given and show how superintegrability alone leads to formulas for the trajectories. These formulas and their analysis are new. The elliptic systems have not been found elsewhere to our knowledge, so we demonstrate the procedure to derive them. This paper is partly pedagogical (the Introduction and the Appendix) but mostly new research (Sections 2-5). In the Appendix we give a brief review of fundamental definitions and results from classical Hamiltonian mechanics, adapted from [33], needed as background for our computations. The point here is to illustrate how, using the structure equations of superintegrable systems alone, we can derive and classify the trajectories via algebra. Hamilton's equations are used only to determine the periods of orbits for systems with degenerate potential. Separation of variable techniques are not used, except in the last section; our approach is easy to understand geometrically. Analogously the spectra of the corresponding quantum systems can be obtained algebraically from the structure relations, though we do not treat this here. Some recent papers, e.g., [26,27], adopt a related but different approach by using ladder operators constructed from the structure algebra to compute trajectories and spectra for systems with degenerate potential. We know of no prior treatments of the nondegenerate caged systems S7 and E16, or of the use of the structure algebra to call attention to special orbits, such as those for which R 2 = 0. We point out the contractions that relate our various superintegrable systems.

Examples of 2D 2nd degree degenerate systems
Our first examples are degenerate systems. These have 1-parameter potentials and always admit a symmetry that is a 1st degree polynomial in the momenta, hence a group symmetry that can be interpreted as invariance with respect to rotation or translation corresponds to a constant of the motion which leads to an analog of Kepler's 2nd Law. The only possibilities in 2 dimensions are constant curvature and Darboux spaces [22]. In [33] there is an example of an analog of the Kepler problem on the 2-sphere that satisfies Kepler's three laws and then by a limiting process (a contraction) goes to the Kepler system in Euclidean space. Here we will treat an analog of the harmonic oscillator on the Darboux space D4(b) and show that it obeys analogs of Kepler's laws. Then by taking contractions to superintegrable systems on the Darboux space D3, on the Poincaré upper half plane (equivalent to a system on a hyperboloid in Euclidean space) and, finally, to the isotropic oscillator system in Euclidean space, we will see that using ideas from superintegrability theory alone we can understand the basic properties of these systems: conservations laws, explicit trajectories, etc., and how they are related. Fig. 1 describes the contraction relationships that we will exploit. An important feature of degenerate systems is that they always admit 4 linearly independent symmetries that are 2nd degree in the momenta, whereas only 3 can be functionally independent. Thus there must be a relation between these symmetries. This relation emerges from the structure algebra obeyed by the symmetries.

The D4(b) oscillator
We consider the superintegrable system with a basis of constants of the motion, H, J = p y , and Here the spaces D4 are indexed by a parameter b > −2. (In the limit as b → −2 the space becomes a 2-sphere and this system becomes the Higgs oscillator [14].) The variable y can be interpreted as an angle, and the space and potential are periodic in y with period π. The constants of the motion generate an algebra under the Poisson bracket obeying the structure equations with Casimir Note that there are 4 linearly independent symmetries J 2 , H, Y 1 , Y 2 as polynomials in the momenta, but there can be only 3 functionally independent generators. The dependence relation is given by (2.6). From the first two equations (2.5) we see that (Y 1 , Y 2 ) transforms like a 2-vector under rotations about the 3-axis. The sum Y 2 1 + Y 2 2 is expressed in terms of constants of motion so the sum is constant. We set Y 2 1 + Y 2 2 ≡ κ 2 where κ is the length of the 2-vector. Thus we can choose a preferred coordinate system such that Y 1 = κ and Y 2 = 0.
To determine a trajectory, we need to express x and y in terms of the constants of the motion along the trajectory. We do this by eliminating the momenta from equations (2.1)-(2.6). We see that Y 1 and Y 2 can be expressed in an alternate form: Y 2 = − sin(2y)H + sin(2y) cosh(2x)p 2 y + cos(2y) sinh(2x)p x p y . (2.8) To eliminate p x , we multiply equation (2.7) by cos(2y), equation (2.8) by sin(2y) respectively, and add them together, which gives us Y 1 cos(2y) + Y 2 sin(2y) = −H + cosh(2x)p 2 y . Using the facts that Y 1 = κ, Y 2 = 0 and J = p y , we get the orbit equation, cosh(2x)J 2 − cos(2y)κ = H. (2.9) Since cosh(2x) is an even function in x, there will be two orbits for positive x and negative x that are symmetric with respect to x-axis. Therefore, without loss of generality, we restrict ourselves to positive x.
Analog of Kepler's second law of planetary motion. We exploit the fact that there is a 1st degree constant of the motion. We consider an orbital motion as taking place in the s 1 − s 2 plane with polar coordinates s 1 = r cos θ, s 2 = r sin θ where r = (2 cosh(2x) + b)/(2 sinh 2x) and θ = 2y. In these coordinates the metric is If α > 0 the potential is attractive to the origin in the plane; if α < 0 the potential is repulsive and the trajectories are unbounded.
Theorem 2.1. Trajectories of the D4(b) oscillator sweep out equal areas in equal times with respect to the origin in the s 1 − s 2 plane.
Proof . In the interval from some initial time 0 to time t the area swept out by the segment of the straight line connecting the origin and the object is Thus the rate at which the area is swept out is since dy/dt = (dy/dθ) · (dθ/dt) = (1/2) · (dθ/dt). To get (dy/dt), we take the Poisson bracket of H and y, which is then plug in this expression for (dy/dt) into equation (2.10), to get, (dA/dt) = J /2 which is a constant.
Analog of Kepler's third law of planetary motion.
Theorem 2.2. For α > 0 the period T of an orbit is Proof . The total area swept out as the trajectory goes through one period is However, from the second law we see that A(T ) = (T /2) · J and, using the fact that κ 2 = J 4 + H 2 + bJ 2 H − αJ 2 , the period can be expressed in the form (2.12).
Now we begin an analysis of the trajectories. We first assume that α > 0 so that the potential is attractive.
Restriction on H. Recall there is a restriction when we try to express r in terms of θ and other constants of the motion. That is, cosh(2x) = (H + cos θκ)/(J 2 ) should always be larger than 1 for any θ. Then since κ > 0, we have the restriction H − κ > J 2 implying that for a closed trajectory, H should always be positive. Here κ is the length of the 2-vector (Y 1 , Y 2 ). To be explicit, κ 2 = J 4 + H 2 + bJ 2 H − αJ 2 . Then by squaring both sides of an alternate form of the restriction equation, H − J 2 > κ, and plugging in the expression for κ 2 , we get (α − (b + 2)H)J 2 > 0. Noticing the fact that b + 2, H, and J 2 are all larger than 0, we reach the conclusion, H < α/(b + 2) for a bounded orbit. For larger values of H the trajectory is unbounded.
Case of bounded trajectory. To plot the trajectories on the s 1 − s 2 plane it is convenient to write s 1 and s 2 both in terms of θ, which is (2.13) (Note: it is useful to divide the r part of s 1 and s 2 by J 4 and introduce new constants H = H/J 2 and κ = κ/J 2 so we can eliminate one constant.) Since in s 1 and s 2 , J always appears in the form of J 2 , so for every positive J , there will always be a duplicate case for the corresponding −J , and without loss of generality, we can restrict our discussion to J > 0. A typical plot of the trajectory on the s 1 -s 2 plane is Fig. 2, and as one or more of H, J and b becomes larger and larger, they will dominate the r term and make r less susceptible to the change of cos θ. Thus the plots will look more and more circular. Furthermore, if κ becomes very large, the plots will tend to move towards the negative s 1 direction and appear in an elongated form as in Fig. 3. The plot will always be symmetric with respect to s 1 -axis, because s 1 (θ) is even and s 2 (θ) is odd in θ. Also, since the plot will always lean towards the negative s 1 side, larger r always occurs at smaller s 1 value.
To trace out an orbit in time, we would need a starting point to integrate Hamilton's equations. Conventionally, we choose a point that is closest to the origin, which we call perihelion. From the plots of the orbits and an easy analysis of equations (2.13), it is obvious that this point is the intersection of the trajectory with the positive s 1 axis. The point on the x-y plane that corresponds to this perihelion is (x 0 , 0) where x 0 = 1 2 arcCosh( κ+H J 2 ). Plugging y = 0 into equation (2.8) gives p x p y = 0. Realizing p y = J = 0, we know p x must be equal to 0. We see that the perihelion points on the phase space trajectory are uniquely determined by the constants of the motion as follows:    and we can see that if we know the perihelion point on the phase space as in equation (2.14), a set of constants of the motion can be uniquely determined.
Case of escape velocity. This is the case when H = J 2 + κ, a bifurcation point on the momentum map. In polar coordinates, as θ → π we have r → ∞, whereas for s 1 , s 2 in equations (2.13), s 1 → −∞ and s 2 → 0. A plot for s 2 vs. s 1 is given in Fig. 4.
In Fig. 4, the two tails will extend to infinity in the direction of negative s 1 -axis, and they will get closer and closer to the s 1 -axis for smaller and smaller s 1 value but they will never touch the s 1 -axis.
Case of unbounded trajectory. Here, there is more than one value of y, i.e., θ/2 that has no corresponding real value of x in equation (2.9). Example 2.3 (H ≤ J 2 ). In this case, since H ≤ J 2 , r will blow up before s 1 goes becomes negative, so the entire trajectory will be bounded on the positive-s 1 side of the plane. A boundary case (H = J 2 ) is plotted in Fig. 5. For Fig. 5 the constants are J = 10, H = 100, κ = 20, b = 0, so for 0 ≤ θ < π/2 and 3π/2 < θ ≤ 2π, corresponding real values of x exist and there is a trajectory. For θ = π/2 and 3π/2, r will goes to infinity and at the same time, s 2 will go to infinity and s 1 to 0 which is represented by the two tails in Fig. 5.
Example 2.4 (J 2 ≤ H < J 2 + κ). In this case, the trajectory will extend to the negative-s 1 side of the plane. A typical plot for an orbit is given in Fig. 6.

Contraction to the Poincaré upper half plane
Using the Hamiltonian H, (2.1) modified by a constant as /( ), and go to the limit as → 0. Then The structure equations become and the Casimir is If the potential is turned off this is essentially the Poincaré upper half plane model of hyperbolic geometry; for b = 4 it is exactly that. The space here consists of all real points (X, Y ) with Y > 0. In this limit there is no longer any periodicity. Using these identities to eliminate the momenta we arrive at the equation for the trajectories: Notice that the orbit equation is even in Y , so the trajectory will always be symmetric with respect to the X-axis. When restricting to the Poincaré upper half plane with Y > 0, the trajectory will be the upper half of the trajectory traced out by the orbit equation. If β > 0, so that the potential is attractive to the boundary Y = 0, this describes the portion of the ellipse (x 2 /A 2 ) + (y 2 /B 2 ) = 1 in the upper half plane, where If β ≤ 0 the potential is repulsive. If however, 4K 2 + (b + 2)β > 0 the trajectory (2.16) is again a portion of an ellipse. There is a special case that when β = 0, we will have A 2 = B 2 such that the trajectory will be upper half of a circle. A plot of a trajectory when β = 0 is given in Fig. 7. And an example of elliptic trajectory is given in Fig. 8.
If 4K 2 + (b + 2)β < 0 the trajectory is the upper sheet y > 0 of the hyperbola ( If 4K 2 + (b + 2)β = 0 equation (2.16) becomes degenerate. From (2.15) we now find Eliminating p Y from these two expressions we find the equation (K 2 Y 2 +X 1 ) 2 −H K 2 (b+2)X 2 = 0 for the unbounded trajectories. These are upper halves of parabolas with orbit equation Y 2 = −(X 1 /K) ± H (b + 2)X, and we notice that a certain subset of H , X 1 , K and b would correspond to two half-parabolas which are symmetric with respect to y-axis. Interestingly, when X 1 /K ≤ 0, the two half-parabolas will cross as in Fig. 11. A special case occurs when H = 0 so that the trajectory becomes a straight line. Since H = (4Y 2 · p 2 Y )/(b + 2), we see this implies that p Y = 0 so the trajectory must be parallel to the x-axis. Furthermore, if H = 0 and X 1 /K = 0, the trajectory will be exactly the y-axis. Moreover, for H = 0 and X 1 /K < 0, there will be no trajectory.
We mention that this Poincaré upper half plane model can be mapped to the unit disk, the Poincaré disk, as well as to the upper sheet of the 2-sheet hyperboloid is 3-space. Further, this model can be contracted to the flat space system H = p 2 x + p 2 y + αx, E5 in our listing [23]. Indeed, if we set We don't go into detail here but we give another detailed example in the next section.

Contraction to the D3 oscillator
The oscillator in the space Darboux 3 (D 3 ) has Hamiltonian with a single Killing vector J = p y and a symmetry algebra basis, {H, J 2 , X 1 , X 2 }. Here, β cos y e x + 1 , The structure relations are and there is the functional relation 4X 2 Theorem 2.5. The D3 oscillator system is a contraction of the D4 oscillator system.
Proof . We can get (2.17) as a limiting case of (2.1) by taking (2.18) First, it is obvious that J = 2J = p y . Then as → 0, since plugging these two identities into the Hamiltonian equation (2.1) and also using equation (2.18), we have because the terms with in front of them will all go to 0 as → 0. This H is exactly the Hamiltonian of the Darboux-3 oscillator. Similarly, by substituting the identities in equation (2.18) into equations (2.7) and (2.8), we have, as → 0. In the same sense, We notice they are exactly the same vectors as X 1 and X 2 . So we have proved that in this limiting case, Darboux-4 space contracts to Darboux-3.
Using this transformation, we can get the orbit equation for the D3 oscillator from the orbit

Embedding D3 in 3D Minkowski space
We can embed D3 as a surface in 3D Minkowski space with coordinates X, Y , Z in such a way as to preserve rotational symmetry. For example, let Then . Such embeddings are not unique. Discussion of trajectories. To see how the trajectory behaves on this surface, we impose the orbit equation e x J 2 − 2 cos(y )κ = H to X, Y and Z. From the orbit equation, we have where we ignore all the primes on the letters for convenience. Also, we treat a = H J 2 and b = 2κ as two new constants, so the number of constants is reduced by 1. Then we can express X, Y and Z all in terms of y. We notice that since e x > 0, the trajectory will be closed when a > b and unbounded when a < b.
Example 2.6. Plots of an"elliptical" shape trajectory and its overhead 2-D view are given in Figs. 12 and 13. (The size of the "cone" that appears on the graphics is determined by the minimum value we choose for e x in MATLAB. The smaller the minimum value, the larger the "cone". For these two plots, the minimum value for e x we choose is 1 and we adjust to different sizes of "cones" for other examples to make the plots clearer.)  Example 2.7. As a becomes smaller, the "ellipse" is prolonged, see Figs. 14 and 15. Seen from above the shape of the trajectory two merging ellipses. (The minimum value for e x we choose is 0.1 for these two plots.) Example 2.8 (escape velocity). For a = b, we encounter the boundary case between bounded trajectory and unbounded trajectory. A plot of the case is given in Fig. 16. The shape of the trajectory when looking from above is like two very "thin" parabolas meeting at the top of the "cone". (The minimum value for e x we choose is 0.001 for these two plots.)

Contraction of the D3 oscillator to the isotropic oscillator
The isotropic oscillator in Euclidean space has the Hamiltonian   and basis symmetries The structure equations are  The functional relation is L 2 We can obtain this system as a limit of the D3 oscillator as follows: In (2.17) we set x = x + ln( 1 ), y = y , H = 2 H, β = (2ω 2 )/ 2 . Then as → 0 we have K = p y and In terms of flat space Cartesian coordinates X = r cos θ, Y = r sin θ we have Figure 19. Elliptical trajectories.
. From the expressions of H and X 1 , we get Then from the expression of X 2 , we have p x p y = 2X 2 − ω 2 XY . Squaring this equation and equating it to the product of the two equations (2.21) we obtain the orbit equation, The shape of the orbit depends only on 8X 2 /(H − 4X 1 ), so we are really investigating the equation X 2 + Y 2 − aXY = 1 where a = 8X 2 /(H − 4X 1 ). We take the right side to be 1 since we now only care about the shape, not the scale.
When a = 0, i.e., X 2 = 0, the trajectory is just a circle. When 0 < a < 2, i.e., 0 < 8X 2 /(H − 4X 1 ) < 2, the trajectory is a tilted ellipse with y = x as the major axis. In Fig. 19, we show three trajectories for different values of a, and also the major axis. As a increases, the ellipse becomes more elongated in the major axis direction.
When a = 2, i.e., H − 4X 1 = 4X 2 , there is a bifurcation point on the momentum map. The trajectory splits into two parallel straight lines. For a > 2, i.e., 8X 2 /(H − 4X 1 ) > 2, the trajectories are hyperbolas with y = x as the symmetry axis. In Fig. 20, we present three hyperbolic trajectories, shown in darker colors as a increases. It is also possible that 8X 2 /(H − 4X 1 ) < 0, in which case the trajectories will be symmetric about the line y = −x.

The Higgs oscillator S3
The classical system S3 on the 2-sphere is determined by Hamiltonian [14] where J 1 = s 2 p 3 −s 3 p 2 and J 2 , J 3 are cyclic permutations of this expression. For computational convenience we have embedded the 2-sphere in Euclidean 3-space. Thus we can write computations, but at the end we restrict to the unit sphere: s 2 1 + s 2 2 + s 2 3 = 1 and s 1 p 1 + s 2 p 2 + s 3 p 3 = 0. The Hamilton equations for the trajectories s j (t), p j (t) in phase space are The classical basis for the constants of the motion is The structure relations are (2.24) and the Casimir relation is To analyze the classical trajectories it is convenient to replace the basis elements L 1 , L 2 in the algebra with the new basis set (2.25) transforms as a 2-vector with respect to rotations about the 3-axis. Indeed, a rotation through the angle β about the 3-axis rotates the vector by 2β. The remaining relations now become This verifies that the length κ of the 2-vector is unchanged under a rotation and shows that (S 1 , S 2 ) is similar to the Laplace-Runge-Lenz vector for the Kepler problem in Euclidean space. (However, a better superintegrable 2-sphere analog of the Kepler problem is the potential αs 3 / s 2 1 + s 2 2 , called S6 in our listing [23].) In analogy with the choice of periaptic coordinates to simplify the Kepler problem, we choose the preferred coordinate system such that the vector (S 1 , S 2 ) points along the 1-axis, i.e., S 2 = 0, S 1 = κ ≥ 0. Then equations (2.22) and (2.23) become (2.26) The last 3 equations can be solved to give Substituting these results into the square of the first equation (2.26) and simplifying, we get the result (assuming α = 0) This is the equation of a cone As 2 1 + Bs 2 2 + Cs 2 3 = 0. The trajectories lie on the intersection of this cone and the unit sphere s 2 1 + s 2 2 + s 2 3 = 1. Thus we get conic sections again, as in the Euclidean Kepler problem, but this time the sections are intersections with the unit sphere, rather than planes. The possible types of trajectory will depend on the signs of A, B, C.
The projection on the s 1 − s 2 plane is the curve (2.28) The identities will be important in the analysis of trajectories to follow. Analog of Kepler's second law of planetary motion. We see from equations (2.28) and (2.30) that for the nonzero angular momentum X the projection of the motion on the unit circle in the s 1 − s 2 plane is a segment of an ellipse, hyperbola or straight line and that none of these curves pass through the center (s 1 , s 2 ) = (0, 0) of the circle. Thus as the particle moves along its trajectory (s 1 (t), s 2 (t), s 3 (t)), the line segment connecting the projection (s 1 (t), s 2 (t)) to the center of the circle sweeps out an area. Introducing polar coordinates s 1 (t) = r(φ(t)) cos φ(t), s 2 (t) = r(φ(t)) sin φ(t) we see that in the interval from some initial time 0 to time t the area swept out is A(t) = 1 2 φ(t) φ(0) r 2 (φ)dφ. Thus the rate at which the area is swept out is dA dt = dA dφ dφ dt = 1 2 r 2 (φ(t)) dφ dt . Now note from Hamilton's equations that along the trajectory Thus dA dt = X 2 , a constant. Theorem 2.10. If the angular momentum is nonzero, the projection of the trajectory on the unit circle in the s 1 − s 2 plane sweeps out equal areas in equal times. Now we begin an analysis of the trajectories. There are two cases, depending on whether the potential is repulsive (α > 0) or attractive (α < 0).
Case 1: α > 0. This is the case of a repulsive potential. The equator of the sphere repels the particle. Thus motion is confined to a hemisphere. We start an analysis of the types of trajectories.
If X = 0 and using the fact that all trajectories will be periodic for a repulsive force, we see that there will necessarily be at least one point on the phase space trajectory for which p 2 = 0. Substituting into the phase space conditions (2.26) it is easy to see that this is possible only for s 1 = 0. Solving all the equations completely we find that these points on the phase space trajectory are uniquely determined by the constants of the motion as follows: This prescription gives us a point on each trajectory, comparable to the aphelion for the Kepler system, where we can start to trace out the orbit. It remains to analyze the possible orbits. To see what is the available parameter space for the case X = 0, note that first we must require κ ≥ 0. Noting the identity (2.29) we see that the quantities in brackets must have the same sign. However, if that sign is negative then the cone (2.27) will degenerate to a point and not intersect the sphere. Thus to obtain trajectories it is necessary that the constants of the motion On the other hand, if these conditions are satisfied, we see from equations (2.32) that there exist trajectories for which the corresponding constants of the motion are assumed. Thus the conditions are necessary and sufficient.
An Analog of Kepler's 3rd law of planetary motion. For nonzero angular momentum X the projection of the motion on the unit circle in the s 1 − s 2 plane is an ellipse (2.28) enclosing the center of the circle. The area of this ellipse is easily seen to be πX / √ H. Let T be the period of the trajectory. Thus T is the length of time for the projection of the trajectory to trace out the complete ellipse and A(T ) = πX / √ H. Since the area of the ellipse is swept out at the constant rate dA Equating these two expressions for A(t) and solving for T we find T = 2π/ √ H.
Theorem 2.11. For X = 0 the period of an orbit is T = 2π/ √ H.
If X = 0 then from (2.27) and (2.28) we see that the projection of the motion in the s 1 − s 2 plane is the straight line segment s 1 = 0. Thus this motion takes place in a plane and the trajectory is a portion of a great circle passing through the poles of the sphere. Here, κ = Case 2: α < 0. This is the case of an attractive potential, where the equator of the sphere attracts the particle. Again, motion is confined to a hemisphere. We first consider the case where X = 0. Without loss of generality we can assume X > 0. Then the projection of the motion is counter-clockwise. Now κ > 0 and from the identity (2.30) we see that there will be 3 classes of trajectories, depending on the value of H.
If H > 0 then the coefficients of s 2 1 and s 2 2 in the projection formula (2.28) must have the same sign, necessarily positive for a real trajectory. Thus and the projection will be a segment of an ellipse. It is straight-forward to check that the ellipse always intersects the circle at 4 points. All of these trajectories lead to annihilation at the equator. There will necessarily be a "perihelion" point on each trajectory: . If H < 0 then the coefficient of s 2 1 is positive and the coefficient of s 2 2 is negative in the projection formula (2.28). Thus and the projections will be segments of hyperbolas. Each branch of the hyperbola intersects the circle at exactly 2 points. Again, these trajectories lead to annihilation at the equator. Again there is a "perihelion" point on each trajectory: so the coefficient of s 2 1 is positive and the coefficient of s 2 2 is 0 in the projection formula (2.28). Thus and the projections will be segments of straight lines. s 2 1 = X 2 /(X 2 − α). Here p 1 = 0 and p 2 2 = X 2 −α, a constant in accordance with the analog of Kepler's Second Law. These trajectories lead to annihilation at the equator. There is a "perihelion" point on each trajectory: s 2 = 0, p 1 = p 3 = 0, s 2 1 = X 2 /(X 2 − α), p 2 2 = X 2 − α. Finally, we suppose X = 0, so that the motion takes place in a plane through the poles. If H − α > 0 then the motion takes place in the plane s 1 = 0. All trajectories annihilate at the equator. The projected trajectories each pass through the center of the circle, and we have p 2 2 = H − α − Hs 2 2 . If H − α < 0 then the motion takes place in the plane s 2 = 0. All trajectories annihilate at the equator. The projected trajectories do not pass through the center of the circle, and we have p 2 1 = H − α − Hs 2 1 . If H = α then there are possible trajectories on any plane passing through the poles. If we go to new rotated s 1 , s 2 coordinates such that the motion takes place in the plane s 1 = 0, then p 1 = 0 and p 2 2 = −αs 2 2 . All trajectories annihilate at the equator, except for unstable equilibria at the poles.
We can ignore the constant term −β and write such that H = H + β. This is the Higgs oscillator on the 2-sphere.
In terms of the Euclidean embedding of the sphere we have s 1 = cos y/cosh x, s 2 = sin y/cosh x, s 3 = sinh x/cosh x, so s 2 1 +s 2 2 +s 2 3 = 1. If we set s 1 = r cos θ , s 2 = r sin θ as in polar coordinates, then r = 1/cosh x and θ = y. However, we notice that in the analog of Kepler's 2nd law part of the Section 2.1, r = 2 cosh(2x) + b/(2 sinh 2x) and θ = 2y. As → 0, r = 1/(2 cosh x). Thus r = 2r and θ = θ/2, so this contraction space is a double covering of the original D4(b) space. We should be careful of this fact in the further computations. Expressed in the phase space (s 1 , s 2 , s 3 , p 1 , p 2 , p 3 ), we have , where J 1 = s 2 p 3 − s 3 p 2 and J 2 , J 3 are cyclic permutations of this expression. Considering the potential part of the Hamiltonian, we have a new basis set, We will show that in this limit, the orbit equation (2.9) will yield the orbit equation (2.28). First, we transform cosh(2x) and cos(2y) in equation (2.9) to express them in terms of s 1 and s 2 .
where we realize that κ = Y 1 which is 2 times the basis S 1 in equation (2.25), so κ = 2κ where κ is the one in equation (2.28). Therefore, we can see that (2.33) is the same orbit equation as (2.28).
Analog of Kepler's second law of planetary motion. Following the same procedure as in Section 2.1, and noticing that r = 2r and θ = θ/2, we have, dA dt = 1 2 4r 2 (θ(t)) dθ 2dt = r 2 dθ dt = 2 dA dt . Therefore, dA dt = J , which matches the expression for the Higgs oscillator. Period of the orbit. With the same procedure as in Section 2.1 and r = 2r and θ = θ/2, we find area swept out as a trajectory goes through one period T is A = 1 2 2π 0 r 2 (θ )dθ = 2 2π 0 r 2 (θ)dθ = 4A where A is the area in equation (2.12). Also, from the second law we see that A = 1 2 J T . Hence, where T is the same as equation (2.11).
In the limit as → 0, we take b = −2 into the above equation of T and get , which differs from the period in Theorem 2.11 by a factor of 1 2 , due to the fact that the contraction is a double covering of the original D4 space.

Contraction of the Higgs to the isotropic oscillator
This contraction has a simple geometric interpretation, the contraction of a 2-sphere to a plane. We can consider the Higgs oscillator as living in a 2-dimensional bounded "universe" of radius 1 in some set of units. Suppose an observer is situated in this universe "near" the attractive north pole. Our observer uses a system of units with unit length where 0 < 1 and we suppose that using these units the universe appears flat to the observer. Thus in the observer's units we have Here, 2 is so small that to the observer it appears that s 3 = 1. Thus, to the observer, it appears that the universe is the plane s 3 = 1 with local Cartesian coordinates (X, Y ). We compare the actual system on the 2-sphere with the system as it appears to the observer and we assume that 2 is so small that it can be neglected, unless we are dividing by it. Now we have Then substituting these results into the Higgs Hamiltonian equation Multiplying both sides of this equation by 2 we obtain the Hamiltonian equation for the Euclidean space isotropic oscillator in agreement with (2.20). Using the same procedure we find that the constants of the motion become K = X = xp y − yp x , and an alternate basis for the symmetries of the isotropic oscillator. Similarly, the structure equations for the Higgs oscillator go in the limit to the structure equations of the isotropic oscillator.
3 Examples of 2D 2nd degree nondegenerate systems 3. 1 The system S7 on the sphere The classical system S7 on the 2-sphere [23] is determined by the Hamiltonian where J 1 = s 2 p 3 −s 3 p 2 and J 2 , J 3 are cyclic permutations of this expression. We have embedded the 2-sphere in Euclidean 3-space, so and we can use the Poisson bracket but at the end we restrict to the unit sphere: s 2 1 + s 2 2 + s 2 3 = 1 and s 1 p 1 + s 2 p 2 + s 3 p 3 = 0. The Hamilton equations for the trajectories s j (t), p j (t) in phase space are The classical basis for the constants of the motion is The structure relations are R = {L 1 , L 2 } and
Suppose the trajectory touches the equatorial plane at r = 1, θ = θ 0 . Then taking a limit in (3.1) as the trajectory goes to the boundary we find that if the trajectory touches the equatorial plane it does so at an angle cos θ which is a solution of the quadratic equation Thus, a necessary condition for the trajectory to reach the equatorial plane is that H ≥ L 1 .
If the trajectory just touches the equatorial plane but doesn't go into the lower hemisphere then we must have H = L 1 and − a 3 2 > L 2 > a 3 2 . An example of touching the equatorial plane is Fig. 21. The different colors in the accompanying graphs correspond to the different curves (3.2), (3.3) that make up the trajectories. If the trajectory passes through the equatorial plane then we must have H > L 1 and 4L 2 1 − 4a − 2L 1 + a 2 1 ≥ 0. The angle of crossing lies in the interval and for each θ 0 in that interval the possible values of L 2 are .   If R 2 = 0 then S(r) ≡ 0. The previous analysis is correct, except the trajectory is a single arc, rather than a loop. We call this a metronome orbit. The particle moves back and forth along the arc with a fixed period. We give no more details here because we will study the analogous systems in our treatment of E16.   Case 2: a 1 > a 2 . Then the portion −1 < s 1 < −a 2 /a 1 of the s 1 -axis becomes attractive, the trajectories are not closed; each end of the trajectory impacts this interval, assuming R 2 > 0. If H > L 1 then one impact occurs in the northern hemisphere and one impact in the southern hemisphere. If H ≤ L 1 both impacts occur in the northern hemisphere.
An example for H > L 1 is Fig. 26, and for case H < L 1 is Fig. 27. If R 2 = 0 then the trajectory is a single arc that impacts the interval −1 < s 1 < −a 2 /a 1 of the s 1 -axis in the northern hemisphere.

System E16 in Euclidean space as a contraction of S7
E16 can be described as a Kepler-Coulomb system with barrier in two dimensions [23,30]. The system has Hamiltonian, In terms of the rotation generator J = xp y − yp x the defining constants of the motion are H and The time derivative of a function F of the phase space parameters along a trajectory is given byḞ ≡ ∂F ∂t = {H, F }, so in Cartesian coordinates we have the equations of motioṅ x = 2p x ,ṗ x = 1 x 3 (x 2 + y 2 ) 3/2 a 1 x 4 + a 2 y 3x 2 + 2y 2 + 2a 3 x 2 + y 2 3/2 , y = 2p y ,ṗ y = 1 (x 2 + y 2 ) 3/2 (a 1 y − a 2 ).
In polar coordinates x = r cos θ, y = r sin θ we have Setting R = {L 1 , L 2 } we can verify that Note also that H = p 2 r + a 1 r + L 1 +a 3 r 2 . In polar coordinates the equations of motion arė From these equations we see that points on a trajectory such thatθ = 0 are characterized by p θ = 0. If such a point has coordinates (θ 0 , r 0 , p r 0 , 0) we find that , Similarly, ifṙ = 0 then p r = 0 and if such a point has coordinates (θ 1 , r 1 , 0, p θ 1 ) we find in particular that r 1 = (a 1 ± a 2 1 + 4H(L 1 + a 3 ))/(2H).
Possible values taken by the constants of the motion depend on the relative sizes of a 2 and a 3 . Case 1: a 3 > a 2 . Then L 1 + a 3 > 0 and we must have a 2 1 + 4H(L 1 + a 3 ) ≥ 0 for trajectories. Unbounded trajectories correspond to H ≥ 0, and for H > 0 perigee occurs at For H = 0 perigee occurs at Bounded orbits correspond to H < 0, in which case perigee occurs at and apogee at We must have a 2 2 + 4L 1 (L 1 + a 3 ) ≥ 0 for trajectories. There is an equilibrium point at Period of bounded orbit. Since dr dt = 2 H − a 1 r − L 1 +a 3 r 2 , we have r dr 2 √ Hr 2 −a 1 r−(L 1 +a 3 ) = dt. Thus the time needed to go from perigee r 1 to apogee r 2 is Thus the period depends only on the energy.
The trajectories. Eliminating p θ , p r in the expressions for H, L 1 , L 2 we find the orbit equations From structure equation (3.4) and the fact that R 2 > 0 we see that the 2nd factor in S is a nonnegative constant. Further, from (3.4) we see that the denominator does not vanish for H < 0, that it has a single root sin θ = −2L 2 /a 1 for H = 0, and a double root for H > 0.
For the single root case with H = 0, assuming the plus sign in (3.5), we set z = sin θ + 2L 2 /a 1 in (3.5) and expand r in a Laurent series series to find the asymptotic expression If the minus sign is assumed in (3.5), then we have the finite limit lim ]. Parametrization of case 1 trajectories. Recall that a 1 , a 2 , a 3 are fixed constants with a 1 < 0, a 3 > a 2 > 0. To parametrize the trajectories we use perigee coordinates 1 , t p , r p where

Note that
For H < 0, apogee occurs for radial distance r a = a 1 r 2 p −rp|2 1 +a 1 rp| 2( 1 +a 1 rp) . For wedge boundaries S 1 , S 2 (radial lines containing a point on the trajectory where S vanishes) to exist at all, we must have a 2 2 + 4( 1 − a 3 ) 1 ≥ 0. Then the righthand boundary will always exist, but the left hand boundary will exist only if 2 1 ≥ a 2 . (However, from the R 2 > 0 inequality below, the left hand boundary must always exist.) The location of the wedge boundaries is defined by t = sin θ, where In order for the system to correspond to a trajectory, we must have R 2 ≥ 0. We will first require that R > 0, the generic case. We see that this is violated if either 1) 2 1 + a 1 r p = 0, or 2) t p lies outside the open interval (t 1 , t 2 ). Thus we have the conditions For H > 0 the zeros of the denominator occur at the values t ± of sin θ such that For H > 0, if the numerator of the trajectory equation (3.5) vanishes at sin(t) = t n then t n = t + or t − and t 1 < t n < t 2 . Since every zero of the denominator is a zero of one of the numerators, we must have t 1 < t − < t + < t 2 . Further, if 2 1 < a 2 then D has at most one zero. Thus, if this case occurs then there is a single asymptote and a single wedge boundary. Summary of necessary Case 1 conditions for trajectories with R 2 > 0.
For examples, see Figs. 28, 29, 30, 31. Case 1a: R 2 = 0. In this case equation (3.6) for the trajectories simplifies since S ≡ 0. We have N 2 = (a 2 2 + 4( 1 − a 3 ) 1 )D so Because N is linear in t = sin θ, it is easily seen that the trajectories are segments of ordinary conic sections: lines, ellipses, circles, hyperbolas and parabolas. Now R 2 = (    Case 1aI): Circular arc. 2 1 + a 1 r p = 0. We can express the constants of the motion in terms of r p , t p alone: The equation for the trajectories simplifies to r(θ) = a 2 2 +4( 1 −a 3 ) 1 N = r p , so the trajectory is a segment of a circle centered at the origin with radius r p . The wedge boundaries are t 1 = −a 2 − √ a 2 2 +a 2 1 r 2 p +2a 1 a 3 rp a 1 rp , t 2 = −a 2 + √ a 2 2 +a 2 1 r 2 p +2a 1 a 3 rp a 1 rp . Thus for this case to occur, r p must satisfy a 2 2 + a 2 1 r 2 p + 2a 1 a 3 r p ≥ 0. We conclude that this case always occurs for r p sufficiently small or sufficiently large, but that there is a forbidden intermediate band for which no circular trajectory exists.
Case 1aII): Conic section arc. 1 − a 3 − a 2 t p − 1 t 2 p = 0. Again we can express the constants of the motion in terms of r p , t p alone: The equation for the trajectories becomes , and the sign is chosen so that r > 0. Setting y = r sin θ, x = r cos θ we see that the trajectory is a segment of the curve αy + βr = ±1 or βr = −αy ± 1, so β 2 (x 2 + y 2 ) = α 2 y 2 ± 2αy + 1, α 2 + β 2 > 0. We have the following possibilities for the curve segments: The energy is always negative. Only strictly negative energies are possible for the elliptical arc.
In general, orbits are possible for H < 0 if and only if the constants of the motion satisfy the inequality a 3 + a 2 t p /(1 − t 2 p )(−a 1 ) < r p . Trajectories are possible in the positive energy case if and only if r p < a 3 + a 2 t p /(1 − t 2 p )(−a 1 ). Case 2: a 3 < a 2 . Now there are no restrictions on the sign of 1 = L 1 + a 3 . The force in the angular direction is always negative. In the first quadrant the force in the y direction is always negative. It follows from this that there are no periodic orbits contained entirely in the first quadrant. The positive y-axis is repelling, but the negative y-axis is attracting.
We assume, initially, that R 2 > 0. Equations (3.5) for the trajectories still hold. The condition a 2 2 + 4 1 ( 1 − a 3 ) > 0 is now satisfied automatically. There is at most a single wedge boundary The following general types of behavior occur: 1. If 1 < 0 the force in the radial direction is always negative. The trajectories all impact the lower y-axis and are perpendicular to the axis at impact. Sinceṗ r is always strictly negative along a trajectory, perigee must occur at y p = −r p on the negative y-axis. Note that R 2 is linear in H, always with positive linear coefficient. The critical value of H such that R 2 = 0 The qualitative behavior of the trajectories divides into two basic classes, depending on the sign of N (sin θ) for θ = − π 2 . Note that Case 2a: N (−1) > 0. Then if we increase H from E 0 to E 0 + for arbitrarily small > 0 there will be two points of intersection of the trajectory with the negative y-axis, one on the curve r = (N + 2 √ S)/D and one on the curve r = (N − 2 √ S)/D. The trajectory will remain bounded and there will be a wedge boundary. If we further increase H, leaving L 1 , L 2 unchanged, this behavior will persist until the critical value Here the intersection point with the negative y-axis and the curve r = (N −2 √ S)/D has moved to y = −∞ whereas the intersection point with the curve r = (N + 2 √ S)/D remains finite, As H is further increased the behavior of the trajectory is that it is unbounded and asymptotic to a radial line in the 2nd quadrant, that there is a wedge boundary and a single intersection point with the negative y-axis.
Case 2b: N (−1) < 0. Again, we increase H from E 0 to E 0 + for arbitrarily small > 0. Now there are no points of intersection with the negative y-axis and, in fact, the trajectory is non-physical until H is increased to beyond the critical value E 1 . The curve is entirely described by the function r = (N + 2 √ S)/D. The trajectory is unbounded, asymptotic to a radial line in the 2nd quadrant and intersects the negative y-axis a single time. There is no wedge boundary. Note that for this case, the region {N (−1) < 0, H ≤ E 1 } is nonphysical. (We do not have a complete proof of this but strong numerical evidence.) 2. If 1 > 0 the behavior is different. The trajectories do not necessarily intersect the negative y-axis. Some trajectories are disjoint. The qualitative behavior of the trajectories divides into two basic classes, depending on the sign of N (sin θ) for θ = − π 2 . Note again that Case 2c: N (−1) > 0. If we increase H from E 0 to E 0 + for arbitrarily small > 0 there will be two points of intersection of the trajectory with the negative y-axis, one on the curve r = (N + 2 √ S)/D and one on the curve r = (N − 2 √ S)/D. The trajectory will remain bounded if H < 0 (with a wedge boundary) and be unbounded if H > 0. The unbounded trajectory will divide into two disconnected parts. If we further increase H, leaving L 1 , L 2 unchanged, this behavior will persist until the critical value Here the intersection point with the negative y-axis and the curve r = (N −2 √ S)/D has moved to y = −∞ whereas the intersection point with the curve r = (N + 2 √ S)/D remains finite, As H is further increased the behavior of the trajectory is that it is unbounded and asymptotic to a radial line in the 2nd quadrant, and a single intersection point with the negative y-axis. There is a wedge boundary.
Case 2d: N (−1) < 0. Again, we increase H from E 0 to E 0 + for arbitrarily small > 0. There are no points of intersection with the negative y-axis and, the trajectory is non-physical for H < 0, but physical for the H > 0 interval until H is increased to beyond the critical value E 1 . Then the trajectory is unbounded and asymptotic to a radial line in the 2nd quadrant. Physical trajectories have a wedge boundary. See Fig. 32.
Case 2f: horizontal line segment, a 2 L 2 + a 1 1 = a 1 a 3 , y 0 = a 2 a 1 . The energy is H = If the energy is negative there is a wedge boundary at t 0 = sin < 0, and the trajectory impacts the negative y-axis. If the energy is positive, the trajectory is unbounded and there is no wedge boundary, but the trajectory still impacts the negative y-axis.
Case 2g: parabolic arc. If α = β we have (2 1 − a 2 ) 2 and the energy is positive. If α = −β we have (2 1 + a 2 ) 2 and the energy is negative. There is a wedge boundary at t 0 = sin , and the trajectory impacts the negative y-axis. Case 2h: elliptic arc.
, the general expression. There is a wedge boundary at t 0 = sin , and the trajectory impacts the negative y-axis.
Case 2i: hyperbolic arc. Here again, , the general expression. There is a wedge boundary at t 0 = sin , and the trajectory impacts the negative y-axis.
There are no radial line trajectories in this case.

The Post-Winternitz superintegrable systems
These systems are especially interesting as a pair of classical and quantum superintegrable systems that do not admit separation of variables, so they cannot be solved by traditional methods. We treat only the classical case here and show that its superintegrability allows us to determine the trajectories exactly. The quantum case can also be solved.

The classical system
The classical Post-Winternitz system is defined by the Hamiltonian [36] The generating constants of the motion are of 3rd and 4th degree: cyp y x 2/3 , Here {L 1 , L 2 } = −188c 3 , so the constants of the motion generate a Heisenberg algebra. Since there is no 2nd degree constant other than H, separation of the Hamilton-Jacobi equations is not possible in any coordinate system, so separation of variables methods cannot be used to compute the trajectories. However, we can find a parametric description of the trajectories. Choosing constants H = E, L 1 = 1 and L 2 = 2 , we can verify the two identities −p 3 y + 9cx 1/3 p x − 1 + 6Ep y = 0, From these results we can solve for p x and x, y as functions of the parameter p y : 5 Classical elliptic superintegrable systems

Elliptic superintegrability
Here we extend the ideas in [40,41] to a system that separates in Jacobi elliptic coordinates [35]. We start with a simple example, the Higgs oscillator S3 (2.22) [23], expressed in elliptic coordinates. .
The Hamilton-Jacobi equation can be taken as H = E for p α = ∂W ∂α , p β = ∂W ∂β . The separation equations in the coordinates α, β are Here, H and L 1 = −p 2 α + sn 2 (α, k)H − A sn 2 (α,k) are constants of the motion, and L 1 = Λ 1 on a coordinate hypersurface. To find another constant of the motion via the action-angle variable construction we introduce M and N which are computed according to where p α , p β are expressed in terms of coordinates α, β by (5.1). Applying action-angle theory, see, e.g., [20,21], we see that M − N is a constant of the motion. The integral M can be written as with a similar integral for N . This can be calculated as an elliptic integral by choosing the new variable ρ = sn 2 (α, k). We then have the relations The change of variables is then tantamount to the computation of where a = 1, b = 1 k 2 and c, d are roots of ρ 2 − Λ 1 E ρ − A E . A similar calculation applies to the substitution µ = sn 2 (β, k) with exactly the same a, b, c and d. A convenient choice of integral is where p and q are arbitrary nonzero numbers. The separation equations are Here, H and L 1 = p 2 α − sn 2 (pα, k)H + with a similar integral for N . This can be calculated as an elliptic integral by choosing the new variable ρ = sn 2 (pα, k). We then have the relations dρ dα = 2p sn(pα, k)cn(pα, k)dn(pα, k), The change of variables is then tantamount to the computation of whereα = pα, a = 1, b = 1 k 2 and c, d are roots of ρ 2 − Λ 1 E ρ − A E . A similar calculation applies to the substitution µ = sn 2 (qβ, k) with exactly the same a, b, c and d. A convenient choice of integral is This integral can be used for t = ρ or µ. If we now form sn(2k E(a − c)(b − d)pq(M − N )) and use the addition formula sn(qu + pv, κ) = sn(qu, κ)cn(pv, κ)dn(pv, κ) + sn(pv, κ)cn(qu, κ)dn(qu, κ) 1 − κ 2 sn 2 (qu, κ)sn 2 (pv, κ) we obtain a constant of the motion. Thus computation yields a rational constant of the motion whenever p/q is a rational number. Let us restrict ourselves to p = 1, q = 2. Then u = α/2 and v = β where We form , We also note that and ρ = sn 2 (α, k) and µ = sn 2 (2β, k). From this we deduce that where P is a constant of the motion rational in the momenta. We can calculate the components of this constant as follows: We could, of course, use the expressions above given in terms of c + d and cd.
For A = 0 the extra constant of the motion is G = N/D where Setting R = {L 1 , G} we obtain the polynomial structure relations Turning to the case when A = 0 we have G = N/D, where The structure relations are All of the examples that we have constructed using the addition theorem for elliptic functions obey polynomial structure equations, but we have no proof that this is true in general.

Example
, This is the special case of (5.3), (5.4) where we have replaced x, y by 1 − x 2 , 1 − y 2 , respectively, set k = 0 and rescaled. Note that H = 4p 2 We have {H, L 1 } = {H, L 2 } = 0, and with R = {L 1 , L 2 }, the polynomial structure relations are Using the three constants of the motion to eliminate the momenta terms we find and the trajectories must satisfy the implicit equation The left-hand side of this equation is quadratic in y, so we can solve for y as a function of x.
We obtain the solutions y = N ±2 Let (x 0 , y 0 ) be a point on the trajectory such that p x 0 = 0. A straightforward calculation gives Assumption: A > 0. Note that It follows that if we are in the region x > 1 and H > 0, then since L 2 1 + 4AH > |L 1 |, only the plus sign is possible in (5.5). We conclude that the trajectory must be unbounded to the right. Examples are Figs. 37 and 38. (The different colors in the figures correspond to different functions describing the trajectories. ) We investigate what happens if the trajectory intersects the lines y = ±x. If the intersection occurs for y = x then we find that x = −L 2 or p x = 2p y , or both. If the intersection occurs for y = −x then we find that x = L 2 , or p x = −2p y , or both. In both cases, we have p 2 Thus, in both cases if |L 2 | = 1, the trajectory intersects the lines y = ±x at finite momentum. However, the velocity is infinite. Indeedẋ = {H, x} = 2p x (1−x 2 )/(y 2 −x 2 ), which blows up unless p x = 0. This suggests that there are no bounded orbits in the square −1 < x, y < 1. However we find some bounded orbits in the region x > 1 for negative energies. Examples of these bounded orbits are Figs. 39 and 40.

Another version of the example
We set y = iz where z is real, multiply H by −1 and L 2 by −i to get the system , .
We will consider this system for A > 0 and resticted to the vertical strip of points (x, z) where −1 < x < 1. Then H must be > 0. The motion is unbounded, there are no bounded orbits, and the origin (0, 0) and the boundary lines x = ±1 are repulsive. The polynomial structure relations are R = {L 1 , L 2 } and Expressing the momenta in terms of the coordinates we have The equations for the trajectories turn out to be y = N ±2 The points (x 0 , z 0 , p x 0 , p z 0 ) on the trajectory where p x 0 = 0 satisfy Requiring first that R 2 > 0, we see that L 1 < 0 is a necessary condition for trajectories (otherwise S < 0 in the strip −1 < x < 1. If |L 1 |/2H > 1 then a further necessary condition for S > 0 in a strip is L 1 + H + A < 0. If |L 1 |/2H ≤ 1 then a necessary condition is L

A Review of some basic concepts in Hamiltonian mechanics
Hamiltonian mechanics is a reformulation of classical Newtonian mechanics alternate to Lagrangian mechanics. In the Hamiltonian formalism a physical system describing the motion of a particle at time t involves n generalized coordinates q j (t), and n generalized momenta p j (t). The phase space of the system is described by points (p j , q j ) ∈ R 2n . The generalized coordinates q j s and generalized momenta p j s are derived from the Lagrangian formulation. The dynamics of the system are given by Hamilton's equations [1,12] dq j dt = + ∂H ∂p j , dp j dt = − ∂H ∂q j , (A.1) where H = H(q 1 , . . . , q n , p 1 , . . . , p n , t) is the Hamiltonian, which in this paper corresponds to the total (time-independent) energy of the system: H = T + V , where T and V are kinetic and potential energy, respectively. Solutions of these equations give the trajectories of the system.
Here T is a function of q alone while V is a function of q alone. Explicitly, where g jk is a contravariant metric tensor on some real or complex Riemannian manifold. That is g −1 = det(g jk ) = 0, g jk = g kj and the metric on the manifold is given by ds 2 = n j,k=1 g jk dq j dq k , where (g jk ) is the covariant metric tensor, the matrix inverse to (g jk ). Under a local transformation q j = f j (q) the contravariant tensor and momenta transform according to so H is coordinate independent. Here m is a scaling parameter that can often be interpreted as the mass of the particle. For the Hamiltonian (A.2) the relation between momenta and the velocities is p j = m n =1 g j q , so that T = 1 2m j,k g jk (q)p j p k = m 2 n ,h=1 g h (q)q q h . Once the velocities are given, the momenta are scaled linearly in m. In mechanics the exact value of m may be important, but for our purposes can be scaled to any nonzero value using the above formulas.
To make direct contact with mechanics we will often set m = 1; for structure calculations we will usually set m = 1/2. The Poisson bracket of two arbitrary functions A(p, q), B(p, q) on the phase space is the function In terms of the Kronecker delta δ jk , coordinates (q, p) satisfy canonical relations {p j , p k } = {q j , q k } = 0, and {p j , q k } = δ jk . Using the Poisson bracket, we can rewrite Hamilton's equations (A.1) as For any function F(q, p), its dynamics along a trajectory q(t), p(t) is dF dt = {H, F}. Thus F(q, p) will be constant along a trajectory if and only {H, F} = 0. If {H, F} = 0, we say that F is a constant of the motion.

A.1 Classical integrability
A system with Hamiltonian H is integrable if it admits n constants of the motion L 1 = H, L 2 , . . . , L n that are in involution: and are functionally independent in the sense that det ∂L j ∂p k = 0. Suppose H is integrable with associated constants of the motion L j . Then by the inverse function theorem we can solve the n equations L j (q, p) = c j for the momenta to obtain p k = p k (q, c), k = 1, . . . , n, where c = (c 1 , . . . , c n ) is a vector of constants. For an integrable system, if a particle with position q lies on the common intersection of the hypersurfaces L j = c j for constants c j , then its momentum p is completely determined. Also, if a particle following a trajectory of an integrable system lies on the common intersection of the hypersurfaces L j = c j at time t 0 , where the L j are constants of the motion, then it lies on the same common intersection for all t near t 0 . Considering p j (q, c) and using the conditions (A.3) and the chain rule it is straightforward to verify ∂p j ∂q k = ∂p k ∂q j . Therefore, there exists a function u(q, c) such that p = ∂u ∂q , = 1, . . . , n. Note that P j (q, ∂u ∂q ) = c j , j = 1, . . . , n, and, in particular, u satisfies the Hamilton-Jacobi equation where E = c 1 . By construction det( ∂u ∂q j ∂c k ) = 0, and such a solution of the Hamilton-Jacobi equation depending nontrivially on n parameters c is called a complete integral. This argument is reversible: a complete integral of (A.4) determines n constants of the motion in involution, P 1 , . . . , P n .

Theorem A.1. A system is integrable if and only if (A.4) admits a complete integral.
A powerful method for demonstrating that a system is integrable is to exhibit a complete integral by using additive separation of variables.
It is a standard result in classical mechanics that for an integrable system one can integrate Hamilton's equations and obtain the trajectories [1,12]. The Hamiltonian formalism is is well suited to exploiting symmetries of the system and an important tool in laying the framework for quantum mechanics.

A.2 Classical superintegrability
Let F = (f 1 (q, p), . . . , f N (q, p)) be a set of N functions defined and locally analytic in some region of a 2n-dimensional phase space. We say F is functionally independent if the N × 2n matrix ∂f ∂q j , ∂f ∂p k has rank N throughout the region. Necessarily N ≤ 2n. The set is functionally dependent if the rank is strictly less than N on the region. In this case there is a nonzero analytic function F of N variables such that F (f 1 , . . . , f n ) = 0 identically on the region. Conversely, if F exists then the rank of the matrix is < N . We say that the Hamiltonian system H is (polynomially) integrable if there exists a set of N = n constants of the motion L 1 = H, . . . , L n , each polynomial in the momenta globally defined (except possibly for isolated singularities), that is functionally independent and in involution: {L j , L k }=0 for 1 ≤ j, k ≤ n. A classical Hamiltonian system in n dimensions is maximally (polynomially) superintegrable if it admits 2n − 1 functionally independent, globally defined constants, polynomial in the momenta (the maximum number possible [33]). At most n functionally independent constants of the motion can be in mutual involution [1]. However, several distinct n-subsets of the 2n − 1 polynomial constants of the motion for a superintegrable system could be in involution. In that case the system is multi-integrable. An important feature of superintegrable systems is that the orbits traced out by the trajectories can be determined algebraically, without the need for integration. Along any trajectory each of the symmetries is constant: L s = c s , s = 1, . . . , 2n − 1. Each equation L s (q, p) = c s determines a (2n − 1)-dimensional hypersurface in the 2n-dimensional phase space, and the trajectory must lie in that hypersurface. Thus the trajectory lies in the common intersection of 2n − 1 independent hypersurfaces; hence it must be a curve. Another important feature of the trajectories is that all bounded orbits are periodic [34].
The polynomial constants of the motion for a system with Hamiltonian H form the (polynomial) symmetry algebra S H of the system, closed under scalar multiplication and addition, multiplication and the Poisson bracket: Indeed if H is a Hamiltonian with constants of the motion L, K. Then αL + βK, LK and {L, K} are also constants of the motion. The n defining constants of the motion of a polynomially integrable system do not generate a very interesting symmetry algebra, because all Poisson brackets vanish. However, for the 2n − 1 generators of a polynomial superintegrable system the brackets cannot all vanish and the symmetry algebra has nontrivial structure. The degree (or order) O(L) of a polynomial constant of the motion L is its degree as a polynomial in the momenta. Here H has degree 2. The degree O(F k ) of a set of generators F k = {L 1 , . . . , L k }, is the maximum degree of the generators. Let S = S F k be a symmetry algebra of a Hamiltonian system generated by the set F k . Clearly, many different sets F k can generate the same symmetry algebra. Among all these there will be a set of generators F 0 k 0 for which = O(F 0 k 0 ) is a minimum. Here is unique, although F 0 k 0 is not. We define the degree (or order) of S to be . The theory of 2nd degree superintegrable systems has been developed in papers such as [8,16,17,18,19,22,28,37]; all such systems are known as are the structures of the symmetry algebras, all of which close at degree 6.