Tools for Verifying Classical and Quantum Superintegrability

Recently many new classes of integrable systems in n dimensions occurring in classical and quantum mechanics have been shown to admit a functionally independent set of 2n-1 symmetries polynomial in the canonical momenta, so that they are in fact superintegrable. These newly discovered systems are all separable in some coordinate system and, typically, they depend on one or more parameters in such a way that the system is superintegrable exactly when some of the parameters are rational numbers. Most of the constructions to date are for n=2 but cases where n>2 are multiplying rapidly. In this article we organize a large class of such systems, many new, and emphasize the underlying mechanisms which enable this phenomena to occur and to prove superintegrability. In addition to proofs of classical superintegrability we show that the 2D caged anisotropic oscillator and a Stackel transformed version on the 2-sheet hyperboloid are quantum superintegrable for all rational relative frequencies, and that a deformed 2D Kepler-Coulomb system is quantum superintegrable for all rational values of a parameter k in the potential.


Introduction
We define an n-dimensional classical superintegrable system to be an integrable Hamiltonian system that not only possesses n mutually Poisson -commuting constants of the motion, but in addition, the Hamiltonian Poisson-commutes with 2n − 1 functions on the phase space that are globally defined and polynomial in the momenta. Similarly, we define a quantum superintegrable system to be a quantum Hamiltonian which is one of a set of n independent mutually commuting differential operators, and that commutes with a set of 2n − 1 independent differential operators of finite order. We restrict to classical systems of the form H = n i,j=1 g ij p i p j + V and corresponding quantum systems H = ∆ n +Ṽ . These systems, including the classical Kepler and anisotropic oscillator problems and the quantum anisotropic oscillator and hydrogen atom have great historical importance, due to their remarkable properties, [1,2]. The order of a classical superintegrable system is the maximum order of the generating constants of the motion (with the Hamiltonian excluded) as a polynomial in the momenta, and the maximum order of the quantum symmetries as differential operators. Systems of 2nd order have been well studied and there is now a structure and classification theory [3,4,5,6,7,8]. For 3rd and higher order superintegrable systems much less is known. In particular there have been relatively few examples and there is almost no structure theory. However, within the last three years there has been a dramatic increase in discovery of new families of possible higher order superintegrable classical and quantum systems [9,10,11,12,13,14,15,16]. The authors and collaborators have developed methods for verifying superintegrability of these proposed systems, [17,18,19,20]. In the cited papers the emphasis was on particular systems of special importance and recent interest. Here, however, the emphasis is on the methods themselves.
In Section 2 we review a method for constructing classical constants of the motion of all orders for n-dimensional Hamiltonians that admit a separation of variables in some orthogonal coordinate system. Then we apply it to the case n = 2 to show how to derive superintegrable systems for all values of a rational parameter k in the potential. Many of our examples are new. In Section 3 we review our method for establishing a canonical form for quantum symmetry operators of all orders for 2-dimensional Schrödinger operators such that the Schrödinger eigenvalue equation admits a separation of variables in some orthogonal coordinate system. Then we apply this method to establish the quantum superintegrability of the caged anisotropic oscillator for all frequencies that are rationally related, and for a Stäckel transformed version of this system on the 2-sheet hyperboloid. We give a second proof for the caged oscillator in all dimensions n that relies on recurrence relations for associated Laguerre polynomials. We also apply the canonical equations to establish the quantum superintegrability of a deformed Kepler-Coulomb system for all rational values of a parameter k in the potential.

The construction tool for classical systems
There are far more verified superintegrable Hamiltonian systems in classical mechanics than was the case 3 years ago. The principal method for constructing and verifying these new systems requires that the system is already integrable, in particular, that it admits a separation of variables in some coordinate system. For a Hamiltonian system in 2n-dimensional phase space the separation gives us n second order constants of the motion in involution. In this paper we first review a general procedure, essentially the construction of action angle variables, which yields an additional n − 1 constants, such that the set of 2n − 1 is functionally independent. This is of little interest unless it is possible to extract n new constants of the motion that are polynomial in the momenta, so that the system is superintegrable. We will show how this can be done in many cases. We start first with the construction of action angle variables for Hamiltonian systems in n dimensions, and later specialize to the case n = 2 to verify superintegrability. Consider a classical system in n variables on a complex Riemannian manifold that admits separation of variables in orthogonal separable coordinates x 1 , . . . , x n . Then there is an n × n Stäckel matrix such that Φ = det S = 0 and the Hamiltonian is and δ ik is the Kronecker delta. Here, we must require Π n i=1 T 1i = 0. We define the quadratic constants of the motion L k , k = 1, . . . , n by As is well known, Here, Furthermore, by differentiating identity (1) with respect to x h we obtain hj . Now we define nonzero functions M kj (x j , p j , L 1 , . . . , L n ) on the manifold by the requirement It is straightforward to check that these conditions are equivalent to the differential equations We can use the equalities (2) to consider M kj either as a function of x j alone, so that d dx j M kj = S jk /2p j or as a function of p j alone. In this paper we will take the former point of view. Now define functions Then we have This shows that the 2n − 1 functions H = L 1 , L 2 , . . . , L n ,L 2 , . . . ,L n , are constants of the motion and, due to relations (3), they are functionally independent. Now let's consider how this construction works in n = 2 dimensions. By replacing each separable coordinate by a suitable function of itself and the constants of the motion by suitable linear combinations of themselves, if necessary, e.g. [21], we can always assume that the Stäckel matrix and its inverse are of the form where f j is a function of the variable x j alone. The constants of the motion L 1 = H and L 2 are given to us via variable separation. We want to compute a new constant of the motionL 2 functionally independent of L 1 , L 2 . Setting M 21 = M, M 22 = −N , we see that from which we can determine M , N . ThenL 2 = M − N is the constant of the motion that we seek. The treatment of subgroup separable superintegrable systems in n dimensions, [18], is also a special case of the above construction. Suppose the Stäckel matrix takes the form which leads exactly to the construction found in [18].

Application of the construction for n = 2
As we have seen, for n = 2 and separable coordinates u 1 = x, u 2 = y we have We will present strategies for determining functions f 1 , f 2 , v 1 , v 1 such that there exists a 3rd constant of the motion, polynomial in the momenta. The constant of the motionL 2 = M − N constructed by solving equations (4) is usually not a polynomial in the momenta, hence not directly useful in verifying superintegrability. We describe a procedure for obtaining a polynomial constant from M − N , based on the observation that the integrals can often be expressed in terms of multiples of the inverse hyperbolic sine or cosine (or the ordinary inverse sine or cosine), and the hyperbolic sine and cosine satisfy addition formulas. Thus we will search for functions f j , v j such that M and N possess this property. There is a larger class of prototypes for this construction, namely the second order superintegrable systems. These have already been classified for 2-dimensional constant curvature spaces [4] and, due to the fact that every superintegrable system on a Riemannian or pseudo-Riemannian space in two dimensions is Stäckel equivalent to a constant curvature superintegrable system [22,23], this list includes all cases. We will typically start our construction with one of these second order systems and add parameters to get a family of higher order superintegrable systems. A basic observation is that to get inverse trig functions for the integrals M , N we can choose the potential functions f (z), v(z) from the list sometimes restricting the parameters to special cases. Many of these cases actually occur in the lists [4] so we are guaranteed that it will be possible to construct at least one second order polynomial constant of the motion from such a selection. However, some of the cases lead only to higher order constants.

Cartesian type systems
For simplicity, we begin with Cartesian coordinates in flat space. The systems of this type are related to oscillators and are associated with functions Z j for j = 1, 5, 6. We give two examples.
[E1]. For our first construction we give yet another verification that the extended caged harmonic oscillator is classically superintegrable. We modify the second order superintegrable system [E1], [4], by looking at the potential .) It corresponds to the choice of a function of the form Z 5 for each of v 1 , v 2 . Evaluating the integrals, we obtain the solutions to within arbitrary additive constants. We choose these constants so that M , N are proportional to inverse hyperbolic sines. Then, due to the formula cosh 2 u − sinh 2 u = 1, we can use the identities (2) to compute cosh A and cosh B: Now suppose that ω 1 /ω 2 = k is rational, i.e. k = p q where p, q are relatively prime integers. Then ω 1 = pω, ω 2 = qω and are both constants of the motion. Each of these will lead to a polynomial constant of the motion. Indeed, we can use the relations (cosh x ± sinh x) n = cosh nx ± sinh nx, cosh(x + y) = cosh x cosh y + sinh x sinh y, sinh(x + y) = cosh x sinh y + sinh x cosh y, recursively to express each constant as a polynomial in cosh A, sinh A, cosh B, sinh B. Then, writing each constant as a single fraction with denominator of the form we see that the numerator is a polynomial constant of the motion. Note that, by construction, both sinh(−4ipqω[M −N ]) and cosh(−4ipqω[M −N ]) will have nonzero Poisson brackets with L 2 , hence they are each functionally independent of H, L 2 . Since each of our polynomial constants of the motion differs from these by a factor that is a function of H, L 2 alone, each polynomial constant of the motion is also functionally independent of H, L 2 . Similar remarks apply to all of our examples.
[E2]. There are several proofs of superintegrability for this system, but we add another. Here, we make the choice v 1 = Z 1 , v 2 = Z 5 , corresponding to the second order superintegrable system [E2] in [4]. The potential is where L 2 = p 2 x + ω 2 1 x 2 + αx and the system is second order superintegrable for the case ω 1 = ω 1 . Applying our method we obtain Thus if ω 1 /2ω 2 is rational we obtain a constant of the motion which is polynomial in the momenta.

Polar type systems
Next we look at flat space systems that separate in polar coordinates. The Hamiltonian is of the form Cases for which the whole process works can now be evaluated. Possible choices of f and g are In each case if p is rational there is an extra constant of the motion that is polynomial in the canonical momenta.
[E1]. Case (1) is system [E1] for k = 1. For general k this is the TTW system [12], which we have shown to be supperintegrable for k rational.
Case (2). This case is not quadratic superintegrable, but as shown in [19] it is superintegrable for all rational k.
[E8]. For our next example we take Case (3) where z = x + iy: for arbitrary k. (If k = 2 this is the nondegenerate superintegrable system [E8] listed in [4].) In polar coordinates, with variables r = e R and z = e R+iθ . Then we have This system is superintegrable for all rational k.
[E17]. Taking Case (4) we have, for z = x + iy: (For k = 1 this is the superintegrable system [E17] in [4].) Then we have Applying our procedure we find the functions This demonstrates superintegrability for all rational k.
[E16]. Case (5) corresponds to [E16] for k = 1, and our method shows that it is superintegrable for all rational k.
Case (6). This case is not quadratic superintegrable, but it is superintegrable for all rational k.

Spherical type systems
These are systems that separate in spherical type coordinates on the complex 2-sphere. The Hamiltonian is of the form Embedded in complex Euclidean 3-space with Cartesian coordinates, such 2-sphere systems can be written in the form where It is convenient to choose spherical coordinates In terms of these coordinates the Hamiltonian has the form [4]. We extend this Hamiltonian via the replacement ϕ → kϕ and proceed with our method. Thus with H expressed as above. The functions that determine the extra constant are Thus this system is superintegrable for all rational k.
[S7]. This system corresponds to v 1 (φ) = Z 2 and v 2 (ψ) a variant of Z 3 . The system is second order superintegrable: Choosing the coordinates ψ and ϕ we find We make the transformation ϕ → kϕ and obtain with H as before. The functions that determine the extra constants are Thus this system is superintegrable for all rational k. [S4].
Here v 1 (ϕ) = Z 4 and v 2 (ψ) is a variant of Z 2 . This is another system on the sphere that is second order superintegrable and separates in polar coordinates: In terms of angular coordinates ψ, ϕ the Hamiltonian is After the substitution ϕ → kϕ we have The functions that determine the extra constants are so this systems is also superintegrable for all rational k. [S2].
Here v 1 (ϕ) = Z 4 and v 2 (ψ) is a variant of Z 3 : which is is second order superintegrable. After the substitution ϕ → kϕ we obtain the system Applying our procedure we find Thus the system is superintegrable for all rational k.

Horospherical systems
In terms of horospherical coordinates on the complex sphere we can construct the Hamiltonian If ω 1 /ω 2 is rational then this system is superintegrable. However, there is no need to go into much detail for the construction because a Stäckel transform, essentially multiplication by 1/y 2 , takes this system to the flat space system generalizing [E2] and with the same symmetry algebra. There is a second Hamiltonian which separates on the complex 2 sphere It is superintegrable for ω 1 /ω 2 rational, as follows from the Stäckel transform 1/y 2 from the sphere to flat space. If ω 1 = ω 2 then we obtain the system [S2]. We note that each of the potentials associated with horospherical coordinates can quite generally be written in terms of s 1 , s 2 , s 3 coordinates using the relations

Generic systems
These are systems of the form where j, k can independently take the values 2, 3, 4 and Z j ,Ẑ j depend on distinct parameters. Superintegrability is possible because of the integrals Consider an example of this last type of system: The equation H = E admits a separation constant As a consequence of this we can find functions M (x, p x ) and N (y, p y ) where We see that if q k is rational then we can generate an extra constant which is polynomial in the momenta.
The superintegrable systems that have involved the rational functions Z 1 , Z 5 and Z 6 work because of the integrals 2 Quantum superintegrability

The canonical form for a symmetry operator
We give a brief review of the construction of the canonical form for a symmetry operator [19]. Consider a Schrödinger equation on a 2D real or complex Riemannian manifold with Laplace-Beltrami operator ∆ 2 and potential V : that also admits an orthogonal separation of variables. If {u 1 , u 2 } is the orthogonal separable coordinate system the corresponding Schrödinger operator has the form and, due to the separability, there is the second-order symmetry operator i.e., [L 2 , H] = 0, and the operator identities We look for a partial differential operatorL(H, L 2 , u 1 , u 2 ) that satisfies [H,L] = 0.
We require that the symmetry operator take the standard form We have shown that we can writẽ and considerL as an at most second-order order differential operator in u 1 , u 2 that is analytic in the parameters H, L 2 . Then the conditions for a symmetry can be written in the compact form We can further simplify this system by noting that there are two functions F (u 1 , u 2 , H, L 2 ), G(u 1 , u 2 , H, L 2 ) such that (10) is satisfied by Then the integrability condition for (11), (12) is (with the shorthand ∂ j F = F j , ∂ jℓ F = F jℓ , etc., for F and G), and equation (13) becomes Here, any solution of (15), (16) with A, B, C not identically 0 corresponds to a symmetry operator that does not commute with L 2 , hence is algebraically independent of the symmetries H, L 2 . (Informally, this follows from the construction and uniqueness of the canonical form of a symmetry operator. The operatorsL = g(H, L 2 ) algebraically dependent on H and L 2 are exactly those such that A = B = C = 0, D = g(H, L 2 ). A formal proof is technical.) Note also that solutions of the canonical equations, with H, L 2 treated as parameters, must be interpreted in the form (8) with H and L 2 on the right, to get the explicit symmetry operators.

The caged anisotropic oscillator
In [19] we used the canonical form for symmetry operators to demonstrate the quantum superintegrability of the TTW system [12,13] for all rational k, as well as a system on the 2-hyperboloid of two sheets [19]. Here we give another illustration of this construction by applying it to the 2D caged oscillator [11]. For p = q This is the second order superintegrable system [E1] on complex Euclidean space, as listed in [4]. Here, where in Cartesian coordinates. We take p, q to be relatively prime positive integers. Thus The 2nd order symmetry operator is and we have the operator identities Based on the results of [18] for the classical case, we postulate expansions of F , G in finite series The sum is taken over terms of the form a = a 0 + m, b = b 0 + n, and c = 0, 1, where m, n are integers. The point (a 0 , b 0 ) could in principle be any point in R 2 , Taking coefficients with respect to the basis (18) in each of equation (15) and (16) gives recurrence relations for these coefficients. The shifts in the indices of A and B are integers and so we can view this as an equation on a two-dimensional lattice with integer spacings. While the shifts in the indices are of integer size, we haven't required that the indices themselves be integers, although they will turn out to be so in our solution. The 2 recurrence relations are of a similar complexity, but rather than write them out separately, we will combine them into a matrix recurrence relation by defining · · · · · · · · · · • · · · · · · · · • · · · · · · • · · · · · · • · • · • · · · · • · • · • · · · · • · · · · · · · · · · · · We write the 2 recurrence relations in matrix form as where each M i,j is a 2 × 2 matrix given below.
It is useful to visualize the the set of points in the lattice which enter into this recurrence for a given choice of (a, b). These are represented in Fig. 1. Although the recurrence relates 10 distinct points, some major simplifications are immediately apparent. We say that the lattice point (a, b) has even parity if a + b is an even integer and odd parity if a + b is odd. Each recurrence relates only lattice points of the same parity. Because of this we can assume that the nonzero terms C a,b will occur for points of one parity, while only the zero vector will occur for points with the opposite parity. A second simplification results from the fact that the recurrence matrices M a+m,b+n are of two distinct types. Either m, n are both even, in which case M a+m,b+n is diagonal, or m, n are both odd, in which case the diagonal elements of M a+m,b+n are zero. Another simplification follows from the observation that it is only the ratio p/q = r that matters in our construction. We want to demonstrate that the caged anisotropic oscillator is operator superintegrable for any rational r. The construction of any symmetry operator independent of H and L 2 will suffice. By writing ω = ω ′ /2, p ′ = 2p, q ′ = 2q in H, we see that, without loss of generality, we can always assume that p, q are both even positive integers with a single 2 as their only common factor.
If for a particular choice of p, q we can find a solution of the recurrence relations with only a finite number of nonzero vectors C a,b then there will be a minimal lattice rectangle with vertical sides and horizontal top and bottom that encloses the corresponding lattice points. The top row a 0 of the rectangle will be the highest row in which nonzero vectors C a 0 ,b occur. The bottom row a 1 will be the lowest row in which nonzero vectors C a 1 ,b occur. Similarly, the minimal rectangle will have right column b 0 and left column b 1 . Now slide the template horizontally across the top row such that only the lowest point on the template lies in the top row. The recurrence gives (a 0 + 1)bA a 0 ,b = 0 for all columns b. Based on examples for the classical system, we expect to find solutions for the quantum system such that a 0 ≥ 0 and b ≥ 0. Thus we require A a 0 ,b = 0 along the top of the minimal lattice rectangle, so that all vectors in the top row take the form Note that the recurrence relations preserve the following structure which we will require: 1. There is only the zero vector at any lattice point (a, b) with a + b odd.
2. If the row and column are both even then C a,b = 0 B a,b .
3. If the row and column are both odd then C a,b = A a,b 0 .
This does not mean that all solutions have this form, only that we are searching for at least one such solution.
To finish determining the size of the minimum lattice rectangle we slide the template vertically downward such that only the right hand point on the template lies in column b 1 . The recurrence 2)B a,b 1 = 0 for all rows a. There are several possible solutions that are in accordance with our assumptions, the most conservative of which is b 1 = 0. In the following we assume only that B a 0 ,b 0 = 0, the structure laid out above, and the necessary implications of these that follow from a step by step application of the recurrences. If we find a vector at a lattice point that is not determined by the recurrences we shall assume it to be zero. Our aim is to find one nonzero solution with support in the minimum rectangle, not to classify the multiplicity of all such solutions. Now we carry out an iterative procedure that calculates the values of C a,b at points in the lattice using only other points where the values of C i,j are already known. We position the template such that the recurrence M a,b , second row from the bottom and on the left (of the template), lies above the lattice point (a 0 , b), a 0 = q, and slide it from right to left along the top row. The case b = b 0 = p has already been considered. For the remaining cases we have Only even values of b need be considered. Now we lower the template one row and again slide it from right to left along the row. The contribution of the lowest point on the template is 0 and we find where we have replaced b by b + 1 so that again only even values of b need be considered.
For the first step we take b = p − 2. Then equation (19) becomes and (20) is vacuous in this case, leaving A q−1,p−1 undetermined. Note that several of the terms lie outside the minimal rectangle. From these two equations we can solve for B q,p−2 in terms of our given B q,p , A q−1,p . Now we march across both rows from right to left. At each stage our two equations now allow us to solve uniquely for A q−1,b+1 and B q,b in terms of A's and B's to the right (which have already been computed). We continue this until we reach b = 0 and then stop. At this point the top two rows of the minimal rectangle have been determined by our choice of B q,p and A q−1,p−1 . We repeat this construction for the third and fourth rows down, then the next two rows down, etc. The recurrence relations grow more complicated as the higher rows of the template give nonzero contributions. However, at each step we have two linear relations where the right hand side is expressed in terms of A's and B's either above or on the same line but to the right of A a−1,b+1 , B a,b , hence already determined. Since the determinant of the 2 × 2 matrix is nonzero over the minimal rectangle, except for the upper right corner and lower left corner, at all but those those two points we can compute A a−1,b+1 , B a,b uniquely in terms of quantities already determined. This process stops with rows a = 1, 2. Row a = 0, the bottom row, needs special attention. We position the template such that the recurrence M a,b lies above the lattice point (0, b), and slide it from right to left along the bottom row. The bottom point in the template now contributes 0 so at each step we obtain an expression for B 0,b in terms of quantities A, B from rows either above row 0 or to the right of column b in row 0, hence already determined. Thus we can determine the entire bottom row, with the exception of the value at (0, 0), the lower left hand corner. When the point M a,b on the template is above (0, 0) the coefficient of B 0,0 vanishes, so the value of B 0,0 is irrelevant and we have a true condition on the remaining points under the template. However, this is a linear homogeneous equation in the parameters B q,p , A q−1,p−1 : χB q,p + ηA q−1,p−1 = 0 for constants χ, η. Hence if, for example, χ = 0 we can require B q,p = −(η/χ)A q−1,p−1 and satisfy this condition while still keeping a nonzero solution. Thus after satisfying this linear condition we still have at least a one parameter family of solutions, along with the arbitrary B 0,0 (which is clearly irrelevant since it just adds a constant to the function G). However, we need to check those recurrences where the point M a,b on the template slides along rows a = −1, −2, −3, −4 to verify that these relations are satisfied. We have already utilized the case a = −4, and for cases a = −3, −2, −1 it is easy to check that the relations are vacuous.
The last issue is the left hand boundary. We have determined all vectors in the minimal rectangle, but we must verify that those recurrences are satisfied where the point M a,b on the template slides along columns b = −1, −2, −3, −4. However again the recurrences are vacuous.
We conclude that there is a one-parameter (at least) family of solutions to the caged anisotropic oscillator recurrence with support in the minimal rectangle. By choosing the arbitrary parameters to be polynomials in H, L 2 we get a finite order constant of the motion. Thus the quantum caged anisotropic oscillator is superintegrable. We note that the canonical operator construction permits easy generation of explicit expressions for the defining operators in a large number of examples. Once the basic rectangle of nonzero solutions is determined it is easy to compute dozens of explicit examples via Maple and simple Gaussian elimination. The generation of explicit examples is easy; the proof that the method works for all orders is more challenging.
Example. Taking p = 6 and q = 4, we find (via Gaussian elimination) a solution of the recurrence with nonzero coefficients We can take A 1,5 = a 1 and A 3,3 = a 2 and then all other coefficients will depend linearly on a 1 and a 2 . Calculating the A, B, C and D from equations (11), (12), (14) we find the coefficient of H −1 has a factor of 2a 2 + 9a 1 and so we set a 1 = ω 6 and a 2 = −9/2a 1 to obtain the symmetry operatorL, (9), where u 1 , u 2 are Cartesian coordinates and The total energy is −λ x − λ y = E. Taking µ 1 = pµ and µ 2 = qµ where p and q are integers we find that the total energy is E = −2µ(pn + qm + pa 1 + p + qa 2 + q).
Therefore, in order that E remain fixed we can admit values of integers m and n such that pn + qm is a constant. One possibility that suggests itself is that n → n + q, m → m − p. To see that this is achievable via differential operators we need only consider We now note the recurrence formulas for Laguerre polynomials viz.
Because of the separation equations we can associate p with a differential operator for both m and n in the expression for Ψ. We do not change the coefficients of L α p±1 (x). Therefore we can raise or lower the indices m and n in Ψ using differential operators. In particular we can perform the transformation n → n + q, m → m − p and preserve the energy eigenvalue. To see how this works we observe the formulas In particular, if we make the choice µ 1 = 2µ, µ 2 = µ then the operator D + (2µ, x)D − (µ, y) 2 transforms X n Y m to −128µ 3 (n+1)(m+α 2 )(m−1+α 2 )X n+1 Y m−2 . We see that this preserves the energy eigenspace. Thus can easily be extended to the case when µ 1 = pµ and µ 2 = qµ for p and q integers. A suitable operator is D + (pµ, x) q D − (qµ, y) p . Hence we have constructed a differential operator that commutes with H! The caged oscillator is quantum superintegrable. This works to prove superintegrability in all dimensions as we need only take coordinates pairwise.
Remark. This second proof of superintegrability, using differential recurrence relations for Laguerre polynomials is much more transparent than the canonical operator proof, and it generalizes immediately to all dimensions. Unfortunately, this oscillator system is very simple and the recurrence relation approach is more difficult to implement for more complicated potentials. Special function recurrence relations have to be worked out. Also, we only verified explicitly the commutivity on an eigenbasis for a bound state.

A proof of superintegrability for a deformed Kepler-Coulomb system
In [16] there is introduced a new family of Hamiltonians with a deformed Kepler-Coulomb potential dependent on an indexing parameter k which is shown to be related to the TTW oscillator system system via coupling constant metamorphosis. The authors showed that this system is classically superintegrable for all rational k. Here we demonstrate that this system is also quantum superintegrable. The proof follows easily from the canonical equations for the system. The quantum TTW system is HΨ = EΨ with H given by (7) where v 2 = β cos 2 (kθ) + γ sin 2 (kθ) = 2(γ + β) sin 2 (2kθ) + 2(γ − β) cos(2kθ) sin 2 (2kθ) .
(Setting r = e R we get the usual expression for this system in polar coordinates r, θ.) In our paper [19] we used the canonical form for symmetry operators to establish the superintegrability of this system. Our procedure was, based on the results of [20] for the classical case, to postulate expansions of F , G in finite series The sum is taken over terms of the form a = a 0 + m, b = b 0 + n, and c = 0, 1, where m, n are integers, a 0 is a positive integer and b 0 is a negative integer. Here F , G are the solutions of the canonical form equations (14), (15), (16) for the TTW system. We succeeded in finding finite series solutions for all rational k, and this proved superintegrability.
In [16] the authors point out that under the Stäckel transformation determined by the potential U = e 2R the TTW system transforms to the equivalent systemĤΨ = −αΨ, wherê Then, setting r = e 2R , φ = 2θ, we find the deformed Kepler-Coulomb system.
From the canonical equations, it is virtually immediate that system (22) is superintegrable. Indeed the only difference between the functions (21) defining the TTW system and the functions defining the Stäckel transformed system is that f 1 and v 1 are replaced by f 1 = exp(4R) and v 1 = −E exp(2R) (we can set E = H) and the former H is replaced by −α. From this we find that the only difference between the canonical equations for the TTW system and the canonical equations for the transformed system is that α becomes −Ĥ and H becomes −α. Since our proof of the superintegrability of the TTW system did not depend on the values of H and α, the same argument shows the superintegrability of the transformed system. Also, any solution (i.e., a second commuting operator) for the TTW system gives rise to a corresponding operator for the transformed system by replacing α and H with −Ĥ and −α, respectively.
It has been explicitly checked, using Maple, that the fifth order operator obtained from these expressions commutes withĤ.

A proof of superintegrability for the caged oscillator on the hyperboloid
Using this same idea we can find a new superintegrable system on the 2-sheet hyperboloid by taking a Stäckel transformation of the quantum caged oscillator HΨ = EΨ, (17) in Cartesian coordinates u 1 , u 2 by multiplying with the potential U = 1/u 2 1 . The transformed system is We embed this system as the upper sheet s 0 > 0 of the 2-sheet hyperboloid s 2 0 − s 2 1 − s 2 2 = 1 in 3dimensional Minkowski space with Minkowski metric ds 2 = −ds 2 0 + ds 2 1 + ds 2 2 , via the coordinate transformation Then the potential for the transformed system is This is an extension of the complex sphere system [S2], distinct from the system (6) that we proved classically superintegrable. It is an easy consequence of the results of [17] that the classical version of system (23) is also superintegrable for all rationally related p, q. However, quantum superintegrability isn't obvious. However, from the results of Section 2.2 it is virtually immediate that this new system is quantum superintegrable for all relatively prime positive integers p, q. This follows from writing down the canonical equations (15), (16), first for the caged oscillator where and then for the Stäckel transformed system where The equations are identical except for the switches α 1 → −E, H → −α 1 . Thus our proof of quantum superintegrability for the caged oscillator carries over to show that the system on the hyperboloid is also quantum superintegrable.

Discussion
A basic issue in discovering and verifying higher order superintegrabily of classical and quantum systems is the difficulty of manipulating high order constants of the motion and, particularly, higher order partial differential operators. We have described several approaches to simplify such calculations. Although our primary emphasis in this paper was to develop tools for verifying classical and quantum superintegrabity at all orders, we have presented many new results. The classical Eucidean systems [E8], [E17] and most of the examples of superintegrability for spaces with non-zero scalar curvature are new. We have explored the limits of the construction of classical superintegrable systems via the methods of Section 1.1 for the case n = 2. Further use of this method will require looking at n > 2, where new types of behavior occur, such as appearance of superintegrable systems that are not conformally flat. We also developed the use of the canonical form for a symmetry operator to prove quantum superintegrability. We applied the method to the n = 2 caged anisotropic oscillator to give the first proof of quantum superintegrability for all rational k. We used the Stäckel transform together with the canonical equations to give the first proofs of quantum superintegrability for all rational k of a 2D deformed Kepler-Coulomb system, and of the caged anisotropic oscillator on the 2-hyperboloid. Then we introduced a new approach to proving quantum superintegrability via recurrence relations obeyed by the energy eigenfunctions of a quantum system, and gave an alternate proof of the superintegrabilty for all rational k of the caged anisotropic oscillator. This proof clearly extends to all n. The recurrences for the caged oscillator are particularly simple, but the method we presented shows great promise for broader application. Clearly, we are just at the beginning of the process of discovery and classification of higher order superintegrable systems in all dimensions n > 2.