Superintegrability and deformed oscillator realizations of quantum TTW Hamiltonians on constant-curvature manifolds and with reflections in a plane

We extend the method for constructing symmetry operators of higher order for two-dimensional quantum Hamiltonians by Kalnins et al (2010 J. Phys. A: Math. Theor. 43 265205). This expansion method expresses the integral in a finite power series in terms of lower degree integrals so as to exhibit it as a first-order differential operators. One advantage of this approach is that it does not require the a priori knowledge of the explicit eigenfunctions of the Hamiltonian nor the action of their raising and lowering operators as in their recurrence approach (Kalnins et al 2011 SIGMA 7 031). We obtain insight into the two-dimensional Hamiltonians of radial oscillator type with general second-order differential operators for the angular variable. We then re-examine the Hamiltonian of Tremblay et al (2009 J. Phys. A: Math. Theor. 42 242001) as well as a deformation discovered by Post et al (2011 J. Phys. A: Math. Theor. 44 505201) which possesses reflection operators. We will extend the analysis to spaces of constant curvature. We present explicit formulas for the integrals and the symmetry algebra, the Casimir invariant and oscillator realizations with finite-dimensional irreps which fill a gap in the literature.


Introduction
Let H be a quantum Hamiltonian on an N -dimensional Riemannian manifold: where ∆ is the Laplace-Beltrami operator.An integral (of H) is a differential operator with meromorphic coefficients which commutes with H.A family of integrals I 1 , I 2 , . . . is independent if there is a Jordan polynomial identity P (I 1 , I 2 , . ..) = 0 that resolves into a sum The Hamiltonian is a superintegrable system if there are more than N independent integrals (including the Hamiltonian) in which N of them pairwise commute.The system is of order d if the maximal order of the integrals, as differential operators, is d.In classical mechanics, the number of independent integrals is at most 2N − 1.This is typically the case for quantum systems also, and a model which reaches this limit in number of integrals is said to be maximally superintegrable.This family of independent integrals forms a generating set for a subspace of all the integrals of the Hamiltonian, called the symmetry algebra.The symmetry algebra contains elements of arbitrarily large order but is finitely generated.These algebras have steadily acquired their own status among the different algebraic structures of interest for application in mathematical physics or other areas, such as finite and infinite-dimensional Lie algebras and super Lie algebras.Repeated commutators of the generating integrals may produce additional integrals which are not independent.Thus, the symmetry algebra is a polynomial Lie algebra under the commutator operation [•, •] with additional relations.These polynomials are of Jordan-type and cross-multiplication can be concisely expressed with the anticommutator {•, •}.
Some systems of third-, fourth-and fifth-order have been studied [16,18,20,22].Searching for potentials such that resulting Hamiltonian is superintegrable involves non-linear compatibility equations, so that higher order systems are difficult to study.Indeed, solutions for these modest orders yield potentials of Painlevé, Chazy class and even beyond.
Several construction procedures have been developed for finding higher order integrals for a potential specified up to a few rational or integral parameters [19,20,21,27,28,29,30,31,32,33].Beyond a direct expansion in powers of the momenta [19,20,21], we have for separable systems that the symmetry algebra corresponds to the dynamic algebra of a simpler Hamiltonian, usually taking the form of generalized Heisenberg algebras [30,31] but obtaining explicit forms is a non-trivial problem.In the quantum case, these algebras are built up from ladder operators and the recurrence relations [32] of orthogonal polynomials and special functions from the Askey-Wilson scheme.
The expansion of integrals as a finite power series of the Hamiltonian and other known integrals so as to represent the whole as a differential operator of reduced order has proved effective [27,28,29].This technique, which we shall designate the expansion method, has been applied to the resolution of third-order integrals on Riemannian manifolds and to obtain various exotic potentials which depend upon Painlevé transcendents.The expansion method is powerful in that no knowledge of the eigenfunctions are required.Its application to the resolution of higher order integrals [27] was met with difficulties in producing explicit results, because the equations involved were of fourth degree.
It is our purpose here to apply the expansion and ladder operator methods to the construction of the symmetry algebra of certain separable Hamiltonians.We shall, throughout, work over a two-dimensional space of constant curvature κ.Positive curvature corresponds to elliptical space, negative curvature to hyperbolic space and no curvature to flat Euclidean space.The Hamiltonian families of interest here are of radial harmonic oscillator type where: ω > 0; M 2 is a second-order differential operator which depends on the angular variable θ and a rational parameter are the trigonometric, flat or hyperbolic sine and cosine.The expansion method, as applied to this Hamiltonian, consists in writing We may thus consolidate any differential operator in the variables r, θ into the form (1.2) where I 0 and I 1 are meromorphic in the first two arguments and polynomial in the last two, arranged so that the H-terms are ordered to the right of the variables and operators.Here, ρ denotes multiplication from the right, i.e. ρ(H)∂ r = ∂ r H.This can be extended to a large class of superintegrable systems.
The paper is organized in the following way.The main results is that we determine the full symmetry algebra which is not normally achieved through a straightforward application of the ladder operator method (see, e.g., Theorem 4.9).In §2, we introduce the two classes of models which we are of focus for this paper, the TTW and PVZ models, and discuss the different forms that their symmetries will take.The approach we shall take in finding their integrals I is to expand as a power series in s.In §3, we find that an infinite expansion corresponds to an infinite-order differential operator.The finite expansion enforces upon the operator the constraints that it is a linear expansion of formal operators L i satisfying where k = p/q in reduced form.The ladder operators L i decompose into a product of an angular operator J and radial operator K, which we determine by similar expansion techniques in §4.
In §5, the generators of the symmetry algebra are presented, their commutation and closure relations are given and finite-dimensional irreducible representations are derived.The latter will give us the spectrum of the Hamiltonian, obtained through purely algebraic means.

Models, Eigenfunctions and Spectra
In this section, we discuss the TTW and PVZ models which are of focus for this paper and give their spectral decomposition, which will be used for comparison with the representation theory of their symmetry algebras to be developed in §5.We shall discuss the differences between how their symmetry algebras will be constructed.
The celebrated TTW model [34,35] has the Posch-Teller operator as angular component: (2.1) where α, β > −1.The Hamiltonian (1.1) is separable in polar coordinates.Therefore the eigenfunctions can be determined analytically by simultaneously diagonalizing H and M 2 .The eigenfunctions of (2.1) are where P is a Jacobi polynomial.The corresponding eigenvalues are k 2 (2ℓ + α + β + 1) 2 .This model has attracted intensive interest [36,37,38,39,40,41,42,43,44] as a higherorder superintegrable system distinct from the anisotropic oscillator.The TTW model includes, for k = 1, 2, 3, the Smorodinsky-Winternitz model [45], the BC 2 system of Type V in the classification of Olshanetsky and Perelomov [46], and the three-body model introduced by Wolfes [47] respectively.The TTW systems in their full generality were introduced in 2009 by Tremblay, Turbiner and Winternitz where they verified the model was superintegrable for k = 4 and conjectured superintegrability for all rational k.Kalnins, Kress and Miller Jr. [32] and Quesne [37] proposed various approaches to demonstrating superintegrability.Studies of this model have generated new problem-solving techniques [27,32] which have been applied to other models that are structurally very similar to the TTW model [48].The ladder operator method of Kalnins et al. involved decomposing the integral into formal components L ±1 of the form 1  2 I I ± 1 2 I II M , I I , I II being differential operators.Then /M are integrals of H.Even for k = 1, this does not yield the lowest degree operators.Indeed, there is a differential operator I III and a polynomial P such that (2.2) The existence of the integral I III was inferred by Kalnins et al. from algebraic considerations and computed an expression for P (H).That M 2 , H, I III , H generate by multiplication the full symmetry algebra for all k is not evident from their analysis.That this is the case we demonstrate in §5.3.Moreover, we derive equation (2.2) from the explicit differential operator form.We also provide an algebraic derivation of the spectrum and a generalization of the results for curved spaces.The classical analogue [35] is well-understood and indeed has been shown to be superintegrable on curved space [42,49].
We shall also consider a model discovered by Post, Vinet and Zhedanov [50] (PVZ) which has angular momentum component where α, β > −1.Here, R is a reflection operator on functions of θ, i.e. (Rf The PVZ model is similar to the TTW model but exhibits its own peculiar properties.In particular, M 2 is a perfect square: (2.4) .
A symmetry algebra had been obtained for this system utilizing the ladder operator method.The ladder operators satisfy generalized Heisenberg algebra relations: for some polynomial P .The algebra stated by Post et al. was generated by M, H, L 1 , L −1 for q even and M, H, L 2 , L −2 for q odd.However, for q odd, we have the property that where ĨI , ĨII are differential-difference operators and ξ ± are constants.Then the full symmetry algebra is given by M, H, ĨI , ĨII , which we show in detail in §5.2.The model for even q does not possess this property, and forms only a simple polynomial algebra satisfying the generalized Heisenberg relations.
The reduction properties associated with the TTW model and the PVZ model for q even, do not appear to necessarily be specific to these models alone but a consequence of the commutation identities associated with the ladder operators, which are common to a variety of systems.
The eigenfunctions for both systems can be expressed as Ψ m,ℓ (r, θ) = S m,ℓ (r)Θ ℓ (θ) where the radial part is found to be The spectrum of H therefore splits into finite degeneracy levels {E m,ℓ } qm+pℓ=n .That is, the eigenvalues are invariant under the transformation (m, ℓ) → (m ± p, ℓ ∓ q) which is the implicit evidence for integral operators being built up out of ladder operators, that raise or lower the quantum numbers of the eigenfunctions.

General Form of the Symmetry
In this section, we demonstrate from the expansion method that every integral of the TTW and PVZ systems must be a linear combination of a product of ladder operators J u K ku where In the case of the TTW model, each equation is really two, with M being a formal operator.
For an operator I consolidated in the manner of (1.2), we write I ,r , I ,θ , . . .as its partial derivatives with respect to just the variables r, θ ignoring any of the operator terms.Then [I, H] = 0 generates the following system: For κ = 0, taking a power series expansion in r to solve the system (3.2a)-(3.2b)naturally suggests itself.For general κ, we find that the appropriate Ansatz is: We obtain the recurrence equations: Combining (3.3a) with (3.3b) in order to eliminate the (I 0,i )-sequence, we get It is not possible for this sequence to admit infinitely many non-zero values without I having infinite degree in H, M 2 .Therefore, we must truncate the series. If so m is 2 or 4. Furthermore, I 0,i = I 1,i = 0 for odd i.
In order to cut off the upper bound on the coefficients, we require that for some sufficiently large integer µ and, from (3.3b), we have with µ even.From here we wish to show I is a linear combination of operators L i satisfying (1.3) where i is any integer with magnitude less than or equal to for a non-zero integer i.Then satisfies (1.3).Conversely, any L i satisfying this commutation relation also satisfies where i = j, then L resolves into homogeneous solutions L ±i plus the particular solution This achieves the desired decomposition.It is evident that the {L i } i are linearly independent as operators with coefficients formally rational in M, H.It follows that each L i must commute with H also. Furthermore, there is a finite basis of operators L i satisfying (1.3) and we may choose a basis J i/p (θ, ∂ θ , R, M ) consisting of operators with no dependence upon r, ∂ r .Then ).These must satisfy (3.1).The basis J i/p can be further decomposed into operators satisfying either one of the equations However, for the TTW model, this is merely an ambiguity on the sign of M acting on the eigenfunctions of M .We take its action to always be of the same sign, i.e., we formally set M J u = J u (M + 2ku).The operators J, K modify respectively the quantum numbers ℓ, m of the eigenfunctions of H.

Ladder Operators
In the previous section we established that the integral must be built up from linear combinations of ladder operators.We shall in this section derive the ladder operators for the values of the exponents in which they exist.This is the integers.We consider the angular ladder operators J for the PVZ and TTW model and finally, the radial ladder operator K which is the same for both.The result of this section is that the form of the integral must be where we permit χ i to be rational functions insofar as I is a differential or differential-difference operator.Having calculated the ladder operators, we determine what forms of coefficients are acceptable, that being when their fractional residue in M, H (treated as formal operators) on the whole vanishes.The TTW model has the additional constraint that I must be even in M .
We use the expansion method to simplify the angular partial derivatives.For the TTW model, we write in order to realize J u as a first-order operator in the form (4.1) For the PVZ model, we do the same by writing Hence, J u can be written as We can employ an analogous reduction for K u : In [32], the reducibility of the integrals was induced from simple cases (i.e.k = 1), but there was no way to verify that the results obtained were in fact of minimal degree.This is then an important problem: to develop systematically a means to obtain the lowest-order integrals.
Here we will develop the method for the TTW and PVZ models on constant-curvature space based on the explicit results of § §4.1-3.
For polynomials in M, H, we consider their divisibility where M, H are taken to be formal operators.The greatest common divisors of a family of polynomials P = {p i (M, H)} i are the polynomials Q which divide P such that all the divisors of P are divisors of Q.For a given representative q ∈ Q, we shall write gcd P ∝ q.The units are the non-zero real numbers.A formal differential operator of the form i 0≤a,b≤1 f abi (r, θ)∂ a r ∂ b θ g abi (M, H) where g abi are polynomials has a gcd given by gcd{g abi }.For the PVZ model, we substitute R for ∂ θ .

PVZ Angular Operators
We determine the coefficients in (4.2) by finding solutions to (3.6) which are polynomials in M .The resultant system is ), where v = ±1 is to be determined.To solve this differential-difference system, we develop the coefficients into a sinusoidal power series with a cosine gauge on J u 1 : Substituting these series into (4.3a)-(4.3b),we get the recurrence relations For finite-order in M , the system (4.4a)-(4.4b)requires J u 0,i = J u 1,i = 0 for all but finite i.The appropriate limits are u .Thus, {M, J u } = −2kuJ u when u is odd and [M, J u ] = 2kuJ u when u is even.Up to proportionality the solution is unique, so we designate our principle solutions to be: and define the rest by J ±2j = (J ∓1 J ±1 ) j and J ±2j±1 = J ±1 J ±2j = J ∓2j J ±1 .In particular, Note that in (4.5) but not in (4.6), all the coefficients except the constant term are divisible by a linear factor in M .This, in fact, occurs for every odd u.
Proof.This is clearly true for |u| = 1 and we assume the lemma holds for |u| = 2j − 1, where j is a positive integer.Using J ±2j±1 = J ±1 J ∓2j±1 J ±1 , we have which is a constant.
The combinations (J u ) 2 and J −u J u for u odd and even respectively, commute with M .We designate them both by Φ u (M ), and is appropriately called the structure function.By induction, we compute: For u even, we have Φ −u (M ) = Φ u (M − 2ku) and Φ u (M ) = Φ j (−M − 2ku) for u odd.With knowledge of the structure function, we may now use equation (4.7) to compute J ±(2j+1) 0,0 (∓k(2j + 1)).The result of the computation, consolidated into a single formula, is: where u is any odd integer.By the lemma, it follows that is polynomial in M .We now need to prove that there are no redundant factors in J u generally.First, we compute some explicit terms in J u so as to illustrate the phenomenon which we shall exploit.For a concise expression, we use exp(θR) = cos(θ) + sin(θ)R.
Then for j > 1: where the ellipses indicates trigonometric functions with lower frequency.We observe in the above that the first few terms in (4.8) and (4.9) consist of products of M + (2i − 1)k where the number of factors is decreasing along with the frequency of the e kℓθR terms which they multiply.In particular, the common divisor of the second and third terms is the same as that of the fourth and fifth terms.This is why J u for u even cannot be reduced after removing the constant term, but J u for u odd can be.Theorem 4.2.Suppose α + β, α − β are not even integers.Then for all odd integer u, and for all non-zero even integer u, Proof.From (4.8) and (4.9), so we need to show that there are non-zero terms in J u when M is evaluated at each of the roots.By writing (− sgn(u)(2j + 1)) = 0 by the hypothesis.Similarly, by writing Theorem 4.2 will be combined with similar results for the K operators to determine that only the factor M + ku will divide the non-constant terms of J u K ku for u odd but this is indivisible for u even.

TTW Angular Operators
Here we will solve equation (3.1) for M 2 given by equation (2.1).Substituting (4.1) into (3.1)produces the system We develop the coefficients as a cosinusoidal power series to obtain the recurrence relations: Again, we require U, U , V, V to be finite.The necessary limits are U = |u|, U = |u| − 1, V = V = 0. We take as our principles and set In particular, we obtain a differential operator with no formal part if we consider symmetric combinations: where f is polynomial.We will need to determine, as with the PVZ model, to what extent f may be rational.Indeed, setting f (M ) = M −1 in (4.11) leaves no residue.
For j > 1, we have: where the ellipses denote trigonometric terms of lower frequency.The structure function Φ u (M ) is defined as J −u J u .For positive values, we determine: The divisibility property which held for odd integers in §4.1 this time holds for all integers.
Consequently, [J ±u − J ±u 0,0 (−ku)](M ± ku) −1 is of the form ±D 0 + D 1 M where D 0 , D 1 are differential operators.We wish to find an explicit formula for J u 0,0 (−ku).To this end, we express the structure function in terms of the coefficients: This gives J j 0,0 (−jk) up to a sign.By J j = J 1 J j−1 , we have: M ).This indicates that in powers of α: ; for all j ≥ 1.This enables us to choose the correct square-root in (4.13).At last, we calculate The result below shows that the coefficients, except in extraordinary cases, have no common factor.
Theorem 4.4.Suppose that α + β, α − β are not integers.Then gcd{J u 0,i , J u 1,i } i≥0 ∝ 1. Proof.We shall assume u is positive without loss of generality.From (4.12), we have . This shows that the coefficients cannot all be equal to zero simultaneously.

Radial Operators
For the determination of the radial operators K u , we resume the procedure of §3, namely by expanding as a power series in s: Equation (3.1) translates to the system We take as our principles The effect of curvature is manifested in the 1 2 κM (M ± 1) term, but otherwise the form of the operators is the same as that on flat space [27].Unlike in the computation of the eigenfunctions in §2, there is no need to divide into curved and non-curved space.This will serve only to increase the order of the symmetry algebra, but not its overall structure.Indeed, the calculations to compute the integrals remains the same since the K-operators only depend explicitly on integer powers of κ, c, s, M, H, whereas the Jacobi polynomials degenerate to the Laguerre polynomials only by a limit process.We set which gives the identity K −u = K u (r, ∂ r , −M ).Generally, we have for positive integer j.The operators K −u (r, ∂ r , M + 2u)K u (r, ∂ r , M ) commute with M, H and are thus polynomials in those operators.We designate this operator by the structure function Ψ u (M, H).We have the formula Observe that Ψ −j (M, H) = Ψ j (−M, H) = Ψ j (M − 2j, H).The base terms have the recurrence and a simple closed-form expression is readily obtainable: We see that K j 1,0 (−j, H) = 0 and the latter becomes Proposition 4.6.For all non-zero integers u, K u (−u, H) = K u 0,0 (−u, H).Proof.By equation (4.15), Now, suppose that for some 0 This equation is non-singular, so This proves the result.
Lemma 4.7.Let u be a non-zero integer, then must be polynomials in M .The denominators of z ± therefore must divide 2ju−1 i=0 (M ± i) respectively.If j = 1 then we are done, otherwise, we shall proceed inductively on j > 1.The common factor here between the two denominators is M , so we deduce the existence of polynomial operators That the leading order terms of Q ′ may be polynomial in M and ∂ r , the denominator of z ′ + must divide (M ± i).In particular, the denominators of y ′′′ − and z ′′′ + share no common factor.We may thus resolve Q ′′′ into a sum of two polynomial operators the result of substituting M − 2u for M in Q iv and denote by (•) vi all the coefficients therein.Then There are no common factors between this product and Φ −u (M, H).So the denominator of But the denominator of z ′′′ + = z iv + , by definition of Q ′′′ , cannot contain the factor M + (2j − 1)u.As , the result is proven.

Symmetry Algebra
The construction of the symmetry algebra for superintegrable systems allows us to describe the degeneracies of the Hamiltonian.Each degenerate family consolidates into the whole eigenspace of the Hamiltonian, and so gives us the whole spectrum.
For J u K ku to be well-defined, then u and ku must both be integers, so u is an integral multiple of q.Therefore the general form of an integral of H is where χ i are formal rational functions, with only finitely many that are non-zero, of M, H such that I is a polynomial in those operators, and, for the TTW model, even in M .Let χ 2 N + χ 2 −N = 0 and χ i = 0 for |i| > N .Then by Theorem 4.9, we require (5.1) 0,j (M ), J ±N q 1,j (M )} j≥0 to be a polynomial in M .By Theorems 4.2 and 4.4, the denominator in (5.1) is a unit.So, by linear independence of {J iq K ip } i , the integral must be of order ≥ 2N (p + q − 1) for the TTW model and ≥ N (2p + q − 2) for the PVZ model.Whether or not we can eliminate any of the factors in the product must depend, by Corollary 4.8, on So far the representation theory for the TTW and the PVZ model has not been developed.In this section we will also provide a description of the representations from the symmetry algebra in the context of deformed oscillator realizations.This is also interesting in its own right as the construction of representations has been mainly limited to quadratic algebras such as Racah or some cubic algebras as in the Heun-Lie algebras.
The deformed oscillator realization [52] was introduced within the context of quadratically superintegrable systems.Here we will rely on the analogous construction for polynomial algebras of three generators of arbitrarily large order.This realization consists in the system generated by the three operators N, b, b † satisfying the relations where Ξ is, for our purposes, an analytic function.The finite-dimensional irreducible representations of this algebra are, up to a translation of N, given by where {|i } n i=0 is the basis of the module.The constraints on finite dimension imply that Ξ(0) = Ξ(n) = 0 with no intermediate roots of Ξ.The central elements, namely R[H] act as constants on the representation, and we denote the eigenvalue of H as E. with no other roots of Ξ. From one of the two roots of Ξ, we may determine the eigenvalue E of H.The energy spectrum will include as a subset the physical spectrum (2.5) found by analysis.However, even the non-physical states can still be interesting in its links with finding new special functions.

PVZ q even
For q even, 2p + q ≥ N (2p + q − 2) for all N > 1.Thus we need only consider the form to find the smallest integral with φ 0 a polynomial and φ ±1 at most quadratic.Let us show that an integral of order 2p + q − 2 or 2p + q − 1 is impossible.We note that

Res
where the ellipses indicate lower frequency terms, cannot be zero unless M divides φ ±1 .Similarly,

Res
Therefore, H is a superintegrable Hamiltonian of order 2p + q with generators M, J q K p , J −q K p .We set and Λ ± = Φ ±q (2pA)Ψ ±p (2pA, H).The symmetry algebra relations are The realisation as a deformed oscillator algebra is straightforward: where u is an arbitrary constant.

PVZ q odd
For q odd, 2p + q − 1 ≥ N (2p + q − 2) for all N > 1.So as before, we consider By a similar argument to §5.1, φ 0 , φ ±1 must divide M .Lemma 4.1 guarantees that may take φ ±1 to be any real number proportion of M .Therefore, H is a superintegrable Hamiltonian of order 2p + q − 1.If we write then X, Y, M, H together generate all the integrals of H.The Hamiltonian is therefore superintegrable of order 2p + q − 1.We take and 2  .
The algebra relations are This is a symmetry algebra of order 4p + 2q − 3 if κ = 0 and of order 2p + 2q − 3 if κ = 0. We have the Casimir operator of the enveloping algebra The element D evaluates to zero in the differential-difference operator realisation.
When k = 1, we have For κ = 0, Λ + − Λ − is of first-order, and the system generated by A, B, C is a Bannai-Ito algebra over the ring R[H]: In order to produce a deformed oscillator realization, we require a 'fermionic' number operator in contrast to the 'bosonic' operator N that alternates the sign of b, b † for left-and right-multiplication.This is achieved if we take (−1) N = cos(πN), giving:

Conclusion
We have applied the expansion method to the TTW and PVZ model.This has proved successful on account of their separability, where we are able to obtain ODE in the radial variable with operator coefficients.The solution lead us to ladder operators connected with orthogonal polynomials which we computed and analyzed their properties in order to obtain the respective symmetry algebras.This shows how the expansion method can lead to the resolution of a range of superintegrable systems.The approach here has also proved more effective in analyzing the ladder operators as we have not needed to consider the action on a basis of eigenfunctions.Moreover, we have shown that the ladder operators do not give the generators directly but have given a method to reduce them down to basic operators.
Other models in two dimensions, such as the Post-Winternitz systems, could be considered as well as various Darboux deformations which provide isospectral or almost isospectral problems, which lead to exceptional Jacobi polynomial of Type I, II and III [53,54].In addition, we may also consider the case where the angular part consists of exotic non-elementary functions, since explicit knowledge of the eigenfunctions is not readily obtainable [18,25,26].
We have also found the degenerate energy spectrum of the two Hamiltonians via the algebraic method.That is, the three generator polynomial algebra formed by the integrals and the deformed oscillator realization.It is not known in general if the spectrum can always be recovered in this way for superintegrable systems, even on two-dimensional spaces.This paper has given examples where this is the case for a symmetry algebra which may be of arbitrarily large order.