Symmetry, Integrability and Geometry: Methods and Applications A Recurrence Relation Approach to Higher Order Quantum Superintegrability ⋆

We develop our method to prove quantum superintegrability of an integrable 2D system, based on recurrence relations obeyed by the eigenfunctions of the system with respect to separable coordinates. We show that the method provides rigorous proofs of superintegrability and explicit constructions of higher order generators for the symmetry algebra. We apply the method to 5 families of systems, each depending on a parameter $k$, including most notably the caged anisotropic oscillator, the Tremblay, Turbiner and Winternitz system and a deformed Kepler-Coulomb system, and we give proofs of quantum superintegrability for all rational values of $k$, new for 4 of these systems. In addition, we show that the explicit information supplied by the special function recurrence relations allows us to prove, for the first time in 4 cases, that the symmetry algebra generated by our lowest order symmetries closes and to determine the associated structure equations of the algebras for each $k$. We have no proof that our generating symmetries are of lowest possible order, but we have no counterexamples, and we are confident we can can always find any missing generators from our raising and lowering operator recurrences. We also get for free, one variable models of the action of the symmetry algebra in terms of difference operators. We describe how the St\"ackel transform acts and show that it preserves the structure equations.

1 Introduction excluded). However, we always choose the generators such that the maximum order is as small as possible. Systems of 2nd order have been well studied and there is now a structure and classification theory [3,4,5,6,7,8]. Until very recently, there were comparatively few known superintegrable systems of order ≥ 3 and virtually no structure theory for the algebra generated by the symmetries. In the last three years, however, 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, and others, have developed methods for verifying superintegrability of these proposed systems, [17,18,19,20,21]. The first method developed by us, to verify superintegrability for 2-dimensional quantum systems, was based on a canonical form for symmetry operators of arbitrary order. This approach succeeded for several important systems, such as the caged anisotropic oscillator and the TTW potential for all rational k = p/q, but it led to multi-term recurrence formulas for which very careful analysis was needed to verify finite dimensional solution spaces, [20,22]. Furthermore, the approach yielded no information about the structure of the algebras generated by the symmetries. In [22] the authors introduced a new method, based on the recurrence relations obeyed by the separated eigenfunctions of the Hamiltonian, and sketched its application to the caged anisotropic oscillator. In this paper we further develop and make rigorous the special function recurrence relation method, and apply it to additional cases, including the TTW system. We show that this new method enables us to prove for the first time that the symmetry algebra for the TTW potential is closed for all rational k, and to compute the structure relations for the algebra. (In a series of recent papers [23,24,25,26] Marquette has used a ladder operator method to determine higher order symmetry operators and structure equations for 2D superintegrable systems that separate in Cartesian coordinates. In particular, he has found the structure equations for the caged anharmonic oscillator. Our method is related to his for separation in Cartesian coordinates, but more general in that it applies to systems that separate in any coordinate system, e.g., polar coordinates.) The basic idea that motivates the method is that, since the formal eigenspaces of the Hamiltonian are invariant under action of any symmetry operator, the operator must induce recurrence relations for the basis of eigenfunctions corresponding to the associated coordinate system in which the eigenvalue equation separates. The recurrence relation method uses the known recurrence relations for hypergeometric functions in order to reverse this process and determine a symmetry operator from a suitable combination of recurrence relations. We can compute the symmetry operators and structure equations for the symmetry algebra by restricting ourselves to a formal "basis" of separated eigenfunctions. Then we appeal to our theory of canonical forms for symmetry operators to show that results obtained on restriction to a formal eigenbasis actually hold as true identities for purely differential operators defined independent of "basis" functions. (We have no proof that this reverse process will always work but conjecture that it will succeed whenever all the separated eigenfunctions are of hypergeometric type.) We start with a simple example on the 2-sphere, to introduce the theory. Then we consider an example in Minkowski space followed by a revisiting of the caged anisotropic oscillator. Finally we treat the TTW potential again. In four of these cases we give the first proofs of the closure and the structure of the symmetry algebras.

A simple system
Here we construct a proof of quantum integrability using recurrence formula techniques. As a trial system we consider a quantum Hamiltonian on the two-sphere: where k = p q and p, q are relatively prime positive integers. (Here, all our considerations are local. We do not require that the potential is globally defined on the two-sphere, or that there are any boundary conditions. Global issues can be examined on a case-by-case basis. Also, all parameters and variables can be complex, except for k which is rational.) This system is clearly integrable, since it admits a 2nd order symmetry, responsible for separation of the eigenvalue equation for H in spherical coordinates. Our aim is to show that this Hamiltonian admits an additional quantum symmetry, so is superintegrable. This is relatively straightforward to do using recurrence formula techniques. The typical separable solution in spherical coordinates has the form where T, U = P or Q are the solutions of Legendre's equation, α = k 2 ( 1 4 − a 2 ) and ψ = kϕ. Indeed if we write Ψ = Θ(θ)Φ(ϕ) then these solutions satisfy the "eigenfunction" equations HΨ = −n(n + 1)Ψ where H = L 1 and L 2 Ψ = −k 2 (N + 1 2 ) 2 Ψ, or (Note that L 2 is the symmetry operator associated with variable separation in spherical coordinates.) Here no boundary conditions are implied and all parameters can be complex. To proceed further we observe the following recurrence formulas for Legendre functions [27]: These relations enable the shifting of the indices ν and µ by ± 1. There are clearly similar relations for the functions U µ ν (y). If we let N → N + q then Ψ → T k(N + 1 2 )+p n (cos θ)(cos ψ) 1/2 U a N +q (sin ψ).
We note that N can be quite arbitrary and now consider the function where x = cos θ and y = sin ψ. This function can be obtained from via raising operators: Here (α) ℓ = α(α+1) · · · (α+ℓ−1) is the Pochhammer symbol, and we assume that we have chosen bases for the separable solutions such that the same recurrence formulas and normalizations hold for all elements of the basis. Thus relation (2.2) holds for all four elements in the basis. Similarly, we look at the possibility that N → N −q. This can be obtained from Υ 0 via lowering operators: This result is independent of which basis we choose. Now consider the differential operator ∆ = ∆ + + ∆ − . It follows from the form of the C and D operators that this is an even polynomial function of N + 1 2 . This can be seen from the ν → −ν − 1 symmetry of the operators D and the µ → −µ symmetry of the operators C: sin 2 kϕ Ψ for any eigenfunction Ψ of L 2 . It follows that under the transformation N → −N −1 we have ∆ + → ∆ − and ∆ − → ∆ + . Thus, everywhere the term (N + 1 2 ) 2ℓ occurs in the expansion of ∆ (where ℓ is a positive integer) we can replace it by (k −2 L 2 ) ℓ and obtain a pure differential operator, independent of the parameters n, N . As a consequence we see that we have constructed a pure differential operator, which we also call ∆, and which preserves each eigenspace of H when acting on functions of the form Similarly, the operator ∆ + − ∆ − goes to its negative under the mapping N → −N − 1, so it is an odd function of N + 1 2 . This implies that∆ = (1/(N + 1 2 ))(∆ + −∆ − ) is an even function of N + 1 2 , so it can also be represented as a pure differential operator, independent of the parameters N , n and it preserves each eigenspace of H. We have constructed two partial differential operators ∆ and∆, each of which commutes with the Hamiltonian H on formal eigenspaces. Thus they act like symmetry operators. However, to prove this we must show that they commute with H when acting on any analytic functions, not just eigenfunctions. To establish this fact we use the canonical form for symmetry operators studied in our papers [20,22].
Although we give the reasoning for this special system, our argument is quite general and immediately applicable to other systems. The operator L 2 determines separable coordinates x, y for the system (in this case (x, ψ). Now consider the commutator [H, ∆]. When acting on formal eigenfunctions Υ 0 the commutator gives 0. We want to show that it vanishes identically. To do this we write [H, ∆] in canonical form by recursively replacing all second derivatives ∂ 2 x , ∂ 2 y in terms of H, L 2 to obtain the expression This expression has to be interpreted as in [20], i.e., the parameters H, L 2 must be moved to the right before being identified as operators. Now applying this operator to any eigenfunction we obtain for all choices of the parameters H, L 2 . Noting that we have 4 linearly independent choices for Υ 0 = P i (x)Q j (y), we can write (2.4) as a set of 4 homogeneous equations for the 4 unknowns A, B, C, D: It is convenient to introduce the determinant function We then note that for our particular system W (P 1 (x), P 2 (x), Q 1 (y), Q 2 (y)) = 0, except at isolated points, since the Wronskian of two independent solutions of a separated eigenfunction equation is nonzero. Thus we conclude that A = B = C = D = 0. Consequently [H, ∆] = 0 identically. This proves that ∆ = L 3 is a symmetry operator for the system. The same proof shows that∆ = L 4 is also a symmetry operator. Since both L 3 and L 4 fail to commute with L 2 , each must be algebraically independent of H = L 1 , L 2 . Thus this system is superintegrable for all rational k.
Example 1. If k = p = q = 1 we have the familiar superintegrable system on the 2-sphere where the corresponding Schrödinger operator is i.e., the Laplacian on the two sphere s 2 1 + s 2 2 + s 2 3 = 1 (in Cartesian coordinates s j ) plus the potential V (s) = α s 2

1
. This corresponds to system [S3] in the list [4]. Then we have L 3 = ∆ + +∆ − and L 4 = 1 where x = cos θ and y = sin ϕ. Note that in this example ψ = ϕ since k = 1. We choose the standard spherical coordinates on the sphere viz. s 1 = sin θ cos ϕ, s 2 = sin θ sin ϕ, s 3 = cos θ and the corresponding expressions The Hamiltonian can then be written as From these calculations we deduce that where the separation equation in polar coordinates is due to L 2 = J 2 3 + α/cos 2 ψ with eigenvalue −(N + We can compute the general structure relations for the symmetries of system (2.1). We obtain To compute [L 3 , L 4 ] we need some preliminary results. We make note of the identities Recall that E = −n(n+1) is the eigenvalue of the Hamiltonian and −k 2 (N + 1 2 ) 2 is the eigenvalue of L 2 corresponding to a basis function. Using the property is an even polynomial function in N + 1 2 and a polynomial function of n(n + 1), hence when acting on separated basis functions F + = P + (H, L 2 ) is a polynomial in H and L 2 . In fact, F + = P (H, L 2 ) as pure differential operators, independent of basis. The proof of this fact is analogous to that given above. We write the operator F + − P + in canonical form A∂ 2 xy + b∂ x + C∂ y + D. Then on an arbitrary eigenbasis we have It follows via the usual Wronskian argument that N )) is an even polynomial function in N + 1 2 and a polynomial function of n(n + 1), so Now it is straightforward to obtain There is of course an extra constraint. In fact , and, symmetrizing, we find Then we can rewrite the structure equations as This shows that the symmetry algebra is generated by the symmetries H, L 2 , L 4 and is closed under commutation.
Example 2. We consider our system for the case k = 1 2 , i.e., p = 1, q = 2: where we now use the fact that ψ = ϕ/2. This is a special case of system [S7] in [4]. We form the functions ∆ + and ∆ − as before:

The Hamiltonian becomes
We obtain The operator describing the separation in polar coordinates is which corresponds to eigenvalue − 1 4 (N + 1 2 ) 2 . The basic commutation relation is

Another system
We can further investigate these ideas for the case of the special potential In polar coordinates the Schrödinger equation has the form If we write E = −β 2 then typical separable solutions in polar coordinates are We now construct operators which induce the transformation Ω → Ω ± 3 on the basis functions Ψ Ω . To do this transparently we use the variable w = δe 3iθ . We make the observation that Similarly we note that Clearly Φ = Φ + + Φ − is an even function of Ω, hence interpretable as a pure differential operator. Similarly the operator Φ + − Φ − is an odd function of Ω, soΦ = 1 Ω (Φ + − Φ − ) is a pure differential operator. We deduce as previously that [H, Φ] = [H,Φ] = 0, hence we have constructed a quantum superintegrable system.
We can extend these ideas to consider the potential where k = p/q and α = −k 2 δ 2 . The solutions have the form In order to map solutions of a fixed β eigenspace into solutions we can use the transformations Ω → Ω ± p. These transformations can be performed by the differential operators where w = δ exp(ipθ/q). If we make the transformation Ω → −Ω then we see that Φ − (−Ω) = (−1) p+q Φ + (Ω). Consequently there are two cases to consider: (a) p + q even: Φ (+) = Φ + + Φ − is an even function of Ω, hence is a pure differential operator, and is also a pure differential operator. Thus we have superintegrability in both cases.
Further, we can prove the finite closure of the symmetry algebra and construct the structure equations. We write and let L 2 be the differential operator whose eigenvalue corresponds to Ω 2 , i.e.
Then a direct computation verifies that the symmetry algebra structure relations are Here, {A, B} = AB + BA, {A, B, C} is the analogous 6-term symmetrizer of 3 operators, and R = 2pL 3 + p 2 L 4 . These relations hold no matter whether p + q is even or odd. This shows that the symmetries H, L 2 , L 4 generate the symmetry algebra, and that it closes. and

The Hamiltonian is
This corresponds to system [E14] in [4]. We find r and the constant describing the separation of variables in polar coordinates is which corresponds to the eigenvalue Ω 2 . We also have the structure relation [L 2 , Example 4. Take k = p = 2, q = 1 and w = δe 2iθ . Then and

The Hamiltonian is
This corresponds to system [E8] in [4]. We find and the constant describing the separation of variables in polar coordinates is which corresponds to the eigenvalue Ω 2 . We also have the structure relation [L 2 , L 4 ] = 4(L 3 +L 4 ).

The caged anisotropic oscillator revisited
In [22] we introduced the recurrence relation method by sketching a proof that the caged anisotropic oscillator was quantum superintegrable. Here we will provide more details and show that the symmetry algebra always closes. This result is not new [23] but we include it here to illustrate explicitly how it is obtainable from recurrences obeyed by Laguerre functions.The system is where µ 1 = pµ and µ 2 = qµ and p, q are positive integers that we assume are relatively prime. We look for eigenfunctions for the equation HΨ = λΨ of the form Ψ = XY and find the normalized solutions where the L α n (x) are associated Laguerre functions [27]. Separation in Cartesian coordinates is determined by either of the symmetry operators where H = L 1 + L 2 . For the separated solutions given above we have the eigenvalue equations L 1 Ψ = λ x Ψ, L 2 Ψ = λ y Ψ, where λ x = −2µ 1 (2n + a 1 + 1), λ y = −2µ 2 (2m + a 2 + 1).
Thus HΨ = EΨ where the energy eigenvalue is In the foregoing we impose no boundary conditions and the Laguerre functions are stand-ins for either of the two linearly independent solutions of the second order ordinary eigenvalue equations. In particular, n, m are allowed to be complex. In order that the eigenspace of H with eigenvalue E remain invariant under the action of a recurrence operator that changes m and n we must keep pn + qm constant. One possibility that suggests itself is that n → n + q, m → m − p. A second possibility is n → n − q, m → m + p. Now note the recurrence formulas for Laguerre functions (or confluent hypergeometric functions) viz We apply these for the cases z = µ 1 x 2 , z = µ 2 y 2 and make use of the eigenvalue equations for λ x , λ y , respectively. Again, we can choose two linearly independent solutions for each eigenvalue equation, each of which satisfies this recurrence. Then, considering the symmetry operators as acting on basis functions Ψ n = X n Y n , we have the recurrences x 2 X n = −4µ 1 (n + 1)X n+1 , (4.1) for either basis solution. Writing u = m + kn where k = p/q, we have m = u − kn and we can characterize a formal eigenfunction corresponding to energy E = −2µ(qu By direct calculation, using recurrences (4.1)-(4.4), we verify the relations (Note that this system is very simple to analyze compared to the other systems we study in this paper because these two operators are defined independent of n. Hence by the arguments that we have given for previous examples, each is a symmetry for the caged oscillator system that is algebraically independent of the pair L 1 , L 2 . No symmetrization or antisymmetrization is needed.) Thus the system is superintegrable. Now we construct the operators Though we have indicated a dependence of operators Φ 1 , Φ 2 on n in order to compute their action on a formal eigenbasis, in fact we see from relations (4.1)-(4.4) that they are pure differential operators, independent of the parameter n. Further they commute with both L 1 and L 2 . Hence by an argument that we have given for a previous example, they must be polynomials in the symmetries L 1 and H. Making the replacements we find The constant terms in the expansions of these relations should be interpreted as the constants times the identity operator. Using these results and the eigenvalue formulas L 1 Ψ n = −2µ 1 (2n + a 1 + 1)Ψ n , L 2 Ψ n = −2µ 2 (2u − 2kn + a 2 + 1)Ψ n we can derive the structure equations for the symmetry algebra by acting on a formal eigenbasis. Then we can use our previous argument to show that the structure equations must hold independent of basis. Let L 3 = Φ + + Φ − and Thus the symmetry algebra closes. We can take the symmetries H, L 1 , L 3 as the generators with R = [L 1 , L 3 ] = −4µpqL 4 and rewrite the structure equations as Example 5. We take the case of equal frequencies: µ 1 = µ 2 = µ, so p = q = 1. The corresponding recurrence operators are To proceed we need to have available the following operators: where L 1 , L 2 come from our general theory and we recall that H = L 1 + L 2 . We calculate the symmetry operators that our method implies.
Then we have This method implies the existence of the symmetry operator M from the expression for L 3 .

The TTW system
A similar but more complicated procedure works for the quantum TTW system [12,13]. Here the Hamiltonian is where we take k = p q as before. The general solution of the eigenvalue problem HΨ = EΨ is Ψ = e − ω 2 r 2 r k(2n+a+b+1) L k(2n+a+b+1) m (ωr 2 )(sin(kθ)) a+ 1 2 (cos(kθ)) b+ 1 2 P a,b n (cos(2kθ)), (5.1) where we have taken α = k 2 ( 1 4 −a 2 ) and β = k 2 ( 1 4 −b 2 ). The L-functions are associated Laguerre and the P -functions are Jacobi, not polynomials in general, [27]. We consider the functions where x = cos(2kθ), X = P or Q and Y = S or T . (Here Π is obtained from Ψ by a gauge transformation to remove the angular factors (sin(kθ)) a+ 1 2 (cos(kθ)) b+ 1 2 .) By this we mean that P a,b n (x) is a Jacobi polynomial if n is an integer. Otherwise it is given by its hypergeometric expression. If X = Q then this denotes the associated second solution of the Jacobi differential equation. Similar remarks apply to the choice of Y = S. We have defined this function to be S A m (r) = e − ω 2 r 2 r k(2n+a+b+1) L k(2n+a+b+1) m ωr 2 and T A m (r) to be a second independent solution. The energy eigenvalue is given by E = −2ω (2(m + nk) + 1 + (a + b + 1)k) (5.2) and A = k(2n + a + b + 1). The separation equation for Θ(θ) is andL 2 is a symmetry operator for the system. Under the gauge transformationL 2 goes to a symmetry that we shall call L 2 and which has the same eigenvalues. We see from the expression for E that in order that an energy eigenvalue be unchanged for different values of m, n we must fix u = m + nk. The two transformations n → n + q, m → m − p and n → n − q, m → m + p will each achieve this. Consider the functions X a,b n (x). If we want to raise or lower the index n we can do so with the operators [27]

Similarly, for the functions Y
(Note that E = −2ω[2(m + nk) + 1 + (a + b + 1)k], A = k(2n + a + b + 1) and, effectively, the operator K + is lowering m by 1 and raising n by q/p, whereas the operator K − is raising m by 1 and lowering n by q/p. Here, E is fixed.) We now construct the two operators and When applied to a basis function Ψ n = Y A m (R)X a,b n (x) for fixed u = m + kn, (so m = u − kn) these operators raise and lower indices according to Ξ + Ψ n = 2 q (−1) p ω p (n + 1) q (n + a + b + 1) q Ψ n+q , (5.5) where (α) q = (α)(α + 1) · · · (α + q − 1) for nonnegative integer q, and we note the relation (−α) q = (−1) q (α−q +1) q . From the explicit expressions (5.3), (5.4) for these operators it is easy to verify that under the transformation n → −n−a−b−1 we have Ξ + → Ξ − and Ξ − → Ξ + . Thus Ξ = Ξ + + Ξ − as a polynomial in n and u is unchanged under this transformation. Therefore it is a polynomial in (2n+a+b+1) 2 and u. As a consequence of the relation λ = −k 2 (2n+a+b+1) 2 , in the expansion of Ξ in terms of powers of (2n+a+b+1) 2 and E we can replace (2n+a+b+1) 2 by L 2 /k 2 and E by H everywhere they occur, and express Ξ as a pure differential operator, independent of parameters. (Note that in the expansion of Ξ in terms of the parameters, a term W E is replaced by W H with the W operator on the left.) Clearly this operator, which we will also call Ξ commutes with H on the eigenspaces of H. However, the same argument as used in Example 1 shows that in fact Ξ commutes with H in general, thus it is a symmetry operator for H. We can also easily see that under the transformation n → −n − a − b − 1 the operator Ξ + − Ξ − changes sign, hence the operatorΞ = (1/(2n + a + b + 1))(Ξ + − Ξ − ) is unchanged under this transformation. Again, making the replacements (2n + a + b + 1) 2 by L 2 /k 2 and u by −2 − (H + 2ωk(a + b + 1))/4ω we can expressΞ as a pure differential operator, independent of parameters, and it is a symmetry operator for H. Each of these symmetries has a nonzero commutator with L 2 , so each is algebraically independent of the set H, L 2 . This proves that the TTW system is quantum superintegrable for all rational k. We set L 3 = Ξ, L 4 =Ξ.
Using the explicit relations (5.5), (5.6) for the action of the raising and lowering operators Ξ ± on a basis we can obtain very detailed information about the structure of the symmetry algebra generated by L 2 , L 3 , L 4 . Applying the raising operator to a basis function, followed by the lowering operator, we obtain the result Reversing the order we find Thus the action of the operator Ξ − (n + q)Ξ + (n) + Ξ + (n − q)Ξ − (n) on any basis function Ψ n is to multiply it by ξ n + η n . However, it is easy to check from expressions (5.5), (5.6) and from (5.8), (5.7) that under the transformation n → −n − a − b − 1 we have Ξ − (n + q)Ξ + (n) ↔ Ξ + (n − q)Ξ − (n) and ξ n ↔ η n . Thus Ξ (+) = Ξ − (n + q)Ξ + (n) + Ξ + (n − q)Ξ − (n) is an even polynomial operator in (2n + a + b + 1), polynomial in u, and ξ n + η n is an even polynomial function in (2n + a + b + 1), polynomial in u. Furthermore, each of Ξ − (n + q)Ξ + (n) and Ξ + (n − q)Ξ − (n) is unchanged under the transformation u −→ −u − (a + b + 1) − 1, hence each is a polynomial of order p in (2u + (a + b + 1)k + 1) 2 = E 2 /4ω 2 . Due to the multiplicative factor ω 2p in each of these expressions we conclude that Ξ (+) is a symmetry operator whose action on a basis is given by a polynomial operator P (+) (H 2 , L 2 , ω 2 , a, b). In fact, The proof is analogous to that given in Example 1. We write the operator Ξ (+) − P (+) in canonical form A∂ 2 xR + b∂ x + C∂ R + D. Then on an arbitrary basis we have It follows via the usual Wronskian argument that A = B = C = D = 0. Similarly we note that the operator Ξ (−) = (Ξ − (n+q)Ξ + (n)−Ξ + (n−q)Ξ − (n))/(2n+a+b+1) is an even polynomial in (2n + a + b + 1), as is (ξ n − η n )/(2n + a + b + 1). Also it is polynomial in E 2 and ω 2 . Thus Ξ (−) = P (−) (H 2 , L 2 , ω 2 , a, b) is a symmetry operator which is a polynomial function of all of its variables. Now we can compute the structure relations explicitly by evaluating the operators on an eigenfunction basis. As we have demonstrated, these relations must then hold everywhere. The results are: Here, {A, B} = AB + BA and {A, B, C} is the analogous 6-term symmetrizer of 3 operators. A more transparent realization, for R = −4k 2 qL 3 − 4k 2 q 2 L 4 , is From this realization we see that the symmetries H, L 2 , L 4 generate a closed symmetry algebra.
Example 6. We consider the TTW system with k = p = q = 1. This is essentially the same as Example 5 but with different assumptions. We choose the operators X = ∂ 2 y −ω 2 +A 2 /y 2 and L 2 = M from the previous example where now A 1 = 1 4 −a 2 , A 2 = 1 4 −b 2 . Then where A = 2n + a + b + 1, and we form the usual combinations and From these expressions we deduce that [L 2 , L 4 ] = −4(L 3 + L 4 ), as expected. Note that it follows from these equations that the second order operator X − Y is a symmetry. Indeed, taking the commutator of H with L 3 we find [H, X − Y ]L 2 + L 2 [H, X − Y ] = 0. From this and formal adjoint properties one can conclude that [H, X − Y ] = 0, so X − Y is a 2nd order symmetry. However, Y − Y doesn't belong to the algebra we have already found. In other words our standard procedure didn't find the lowest order generator for the symmetry algebra in this case; it determined a proper subalgebra of the full symmetry structure. To find the full algebra we need to check that the symmetry algebra generated by X − Y , L 2 , H closes at finite order, though in this particular example we already know this to be the case. In the next subsection we will give another approach to this problem that describes how the missing symmetry can always be found and shows that it expressible in terms of the fundamental raising and lowering operators.

The symmetry L 5
We investigate the fact that, as shown in Example 6, our method doesn't always give generators of minimal order. For the case k = 1 we know from the structure theory of 2nd order superintegrable 2D systems with nondegenerate potential [28,3] that the space of 3rd order symmetries is 1-dimensional. Thus, we know that there is a 2nd order symmetry operator L 5 for this case, independent of L 2 , H, such that [L 2 , L 5 ] = L 4 . We will show how to obtain this symmetry from the raising and lowering operators Ξ ± , without making use of multiseparability. Thus, for general rational k we look for a symmetry operator L 5 such that [L 2 , L 5 ] = L 4 . Applying this condition to a formal eigenbasis of functions Ψ n we obtain the result The general solution is where β n is a rational scalar function. It is easy to check that the quantity in parentheses is a rational scalar function of (2n + a + b + 1) 2 . Thus we will have a true constant of the motion, polynomial in the momenta, provided we can choose β n such that the full quantity (5.9) is polynomial in (2n + a + b + 1) 2 . To determine the possibilities we need to investigate the singularities of this solution at n = − a+b+1 2 , n = − ±q+a+b+1 2 . It is easy to check that the operator has a removable singularity at n = − a+b+1 2 . In the special case k = 1 the singularities at n = − ±1+a+b+1 2 are removable, provided we set β n = − H(a 2 −b 2 ) 32(2n+a+b+2)(2n+a+b) . Further, the expression for L 5 in this case k = 1 implies Similar calculations show that the symmetries L 2 , L 5 , H generate the full closed symmetry algebra, which contains our original symmetry algebra properly. For general rational k we can set and determine a polynomial Q(H) such that the operator has a removable singularity at n = −(q + a + b + 1)/2, i.e., such that the residue is 0. (Since the solution is a polynomial in (2n + a + b + 1) 2 we don't have to worry about the singularity at n = −(−q + a + b + 1)/2.) Fixing E we can consider the operators K ± as depending on A alone. Note that Setting µ = −(a + b + 1)/2, we see that we have to evaluate the product Ξ + Ψ −q/2−µ , which factors as Thus to compute the residue we can pair up terms using the following consequences of the recurrence relations above: The caged anisotropic oscillator. From expressions (4.5), (4.6) and the eigenvalue equations for H and L 1 it is easy to write down a function space model for irreducible representations of the symmetry algebra. Note that since H commutes with all elements of the algebra, it corresponds to multiplication by a constant E in the model. We let the complex variable t correspond to n. Then the action of the symmetry algebra on the space of polynomials f (t) is given by difference operators 3) These operators satisfy relations (Note that, by using a similarity transformation via the Mellin transform and its inverse, we could also transform the model (6.1)-(6.3) into a realization by differential operators in one variable.) For a finite dimensional irreducible representation the unnormalized eigenfunctions of L 2 are delta functions Ψ t N (t) = δ(t + t N ). To derive an inner product on the representation space we can require that the adjoint of Φ + is Φ − and that L 2 is self-adjoint.
We can use the model to find several families of finite dimensional and infinite dimensional irreducible representations of the symmetry algebra, only some of which correspond to those that arise from the quantum mechanical eigenvalue associated with real physical systems. As an example of the use of the model to construct finite dimensional representations let us assume that a 1 , a 2 are real and look for a realization such that t N = t 0 + N q where N is an integer and 0 ≤ t 0 < q is fixed. Let N 0 be the smallest integer such that Ψ N 0 = δ(t+t N 0 ) is an eigenfunction, and let N 1 be the largest such integer. Then we must have Φ − Ψ N 0 = 0, Φ + Ψ N 1 = 0. We see from the model that one way to accomplish this is to choose where p 0 , q 0 are integers such that 0 ≤ p 0 < p and 0 ≤ q 0 < q. Now set M = N 1 − N 0 , so that the dimension of the representation is M + 1. Simple algebra gives the eigenvalue for H and eigenvalues The TTW potential. We use expressions (2.2), (2.3) and the eigenvalue equations for H and L 2 to write define a function space model for irreducible of the symmetry algebra. The results can be presented in a simpler form if we use a gauge transformation Φ n = g(n)Ψ n where g(n) g(n − q) = (n − q + 1) q (−u − k(n + a + b + 1)) p (−n − a) q and introduce a complex variable s, corresponding to n + (a + b + 1)/2. Then the action of the symmetry algebra on the space of polynomials f (s) is given by difference operators From these expressions it is straightforward to show that the space of polynomials in the variable s 2 is invariant under the action of the one-variable difference operators. For a finite dimensional irreducible representation the unnormalized eigenfunctions of L 2 are delta functions Φ s N (t) = δ(s + s N ). To derive an inner product on the representation space we can require that L 2 and L 4 are self-adjoint. We can use the difference model to find several families of finite dimensional and infinite dimensional irreducible representations of the symmetry algebra, only some of which correspond to those that arise from the quantum mechanical eigenvalue associated with real physical systems. As an example of the use of the model to construct finite dimensional representations let us assume that a, b, ω are real and look for a realization such that s N = s 0 + N q where N is an integer and 0 ≤ s 0 < q is fixed. Let N 0 be the smallest integer such that Φ N 0 = δ(s + s N 0 ) is an eigenfunction, and let N 1 be the largest such integer. We see from the model that one way to accomplish this is to choose where p 0 , q 0 are integers such that 0 ≤ p 0 < p and 0 ≤ q 0 < q. Now set M = N 1 − N 0 , so that the dimension of the representation is M + 1. We find for H and eigenvalues for L 2 . The system (3.1). For completeness we give the simple one variable model for this case:

Stäckel transforms and the recurrence method
The theory of Stäckel transforms does not guarantee that if two classical 2D superintegrable systems are related by a Stäckel transform and if one system is quantum superintegrable, then the other system is also quantum superintegrable and the quantum systems are related by a Stáckel transform. Some additional conditions must be fulfilled, [18]. However, if one of the systems is known to be quantum superintegrable via our recurrence relation method, then it is automatic that the second quantum system is also superintegrable and a quantum Stäckel transform of the first. We will give a single example which makes clear the general proof. In [16] there was introduced a new family of Hamiltonians with a deformed Kepler-Coulomb potential dependent on an indexing parameter k which was shown to be related to the TTW oscillator system system via coupling constant metamorphosis. The authors showed that this deformed Kepler system is classically superintegrable for all rational k = p/q, and in [22] we used the canonical equations for higher order symmetry operators to show that that it is quantum superintegrable. Here we use the recurrence relations obeyed by the eigenfunctions of the separating symmetry operators to give a new proof of quantum superintegrability and, also, to obtain the structure equations.
As stated above, expressed in polar coordinates r, θ, the quantum TTW system is HΨ = EΨ or The deformed Kepler-Coulomb system, expressed in polar coordinates R, φ, is H ′ Ψ = EΨ or Now note that if we divide both sides of expression (7.1) by r 2 , rearrange terms, and make the change of variables R = r 2 , φ = 2θ, then this expression is identical to (7.2) with the identifications Thus (7.2) is a Stäckel transform of (7.1) and the two systems are Stäckel equivalent, [18,16,22]. The principal observation that we need to make is that both systems have exactly the same formal eigenfunctions (5.1), modulo variable substitution and identifications (7.3), and exactly the same separation equations and recurrence formulas. Thus, substituting into (5.2) we see that the energy levels for the deformed Kepler-Coulomb system are E = Z 2 (2(m + nk) + 1 + (a + b + 1)k) 2 .
In the expressions for the raising and lowering operators of the TTW system we replace the operator E ∼ H by the constant 4Z to get the the raising and lowering operators for the deformed Kepler-Coulomb system. Similarly the structure equations for the Kepler-Coulomb system are obtained by simple permutations We find Similarly the operator L 5 can be added to the system.
We have developed a new method for verifying quantum superintegrability for 2D systems and given several applications to families of such systems, most notably the caged anisotropic oscillator, the Tremblay-Turbiner-Winternitz system and the deformed Kepler-Coulomb system and given new proofs of superintegrability for all rational k. The method relies on the assumption that the system has a second order symmetry operator, so that the Schrödinger eigenvalue equation HΨ = EΨ separates in a set of orthogonal coordinates u 1 , u 2 determined by the symmetry and so that the separated eigenfunctions Ψ = U 1 (u 1 )U 2 (u 2 ) satisfy computable recurrence relations. In practice this means that the separated solutions need to be of hypergeometric type. Then one employs the recurrence relations to construct operators that commute with H on a formal eigenbasis. We used our earlier developed canonical form for higher (≥ 2) order symmetry operators to show that operators commuting with H on a formal eigenbasis must be actual symmetry operators, further that operator identities verified on formal eigenbases must hold identically. Using this approach one can obtain explicit, though complicated, expressions for the higher order symmetry operators. We saw that in the case of the TTW potential our method didn't lead immediately to the lowest order generators, but that they could be found and expressed in terms of our raising and lowering operators. We have no proof as yet that we have found the lowest order generator for all rational k but this is the case for all examples that we know. Provable determination of the maximal symmetry algebra is a topic for future research. By acting on formal eigenbases we were able to compute symmetry algebras for each of the systems we studied and show that these algebras were closed under commutation. One striking result of these computations was that the structure equations for the symmetries were quite explicit and much simpler that the expressions for the symmetry operators themselves. In essence, one can determine the structure equations without knowing the higher order generating symmetries! In each case the action of the symmetries on a formal eigenbasis led us to a simple model of the associated symmetry algebra and its representations, in terms of difference operators in one variable. This greatly simplifies the analysis of the structure equations, classification of irreducible representations of the symmetry algebra and determination of the spectral properties of the generating symmetries. Of course adding the generator L 5 to the algebra complicates these models, though by construction it can always be realized as a difference operator. This is an issue for future research. Such models have independent interest [7,8,29].
Finally, we gave an example of the use of the Stäckel transform to map one superintegrable family of systems to another, while preserving the structure equations. In this case we mapped the TTW system to a deformed Kepler-Coulomb system and determined, for the first time, the structure equation for the Kepler-Coulomb system.
It is interesting that this approach to quantum superintegrability with structure results for symmetry algebras has preceded the classical approach; usually the reverse is true. To our knowledge, the closure of classical symmetry algebras has not been proven for the systems considered here. We are actively investigating the classical analog of the quantum construction. It is clear that these methods have much greater applicability than just the examples treated here. Indeed, all of the systems studied in [22] could be so analyzed. Further all our methods can clearly be extended to systems in dimensions > 2. There appears to be no obstacle other than growing complexity.
When it applies, the recurrence relation method is much simpler than our earlier introduced canonical form method for verification of quantum superintegrability and it provides much more information, including the structure equations. However, the recurrence method requires a detailed knowledge of recurrence relations obeyed by the separated solutions and there are examples of superintegrable systems where no such relations appear to exist [30, Section 4.2].
Thus the canonical form method for higher order symmetry operators seems to be more general. Furthermore, we used the canonical form at a crucial point in the recurrence approach to show that computations valid on formal eigenbases actually held identically. The approaches are mutually complementary.