Two-Variable Wilson Polynomials and the Generic Superintegrable System on the 3-Sphere

We show that the symmetry operators for the quantum superintegrable system on the 3-sphere with generic 4-parameter potential form a closed quadratic algebra with 6 linearly independent generators that closes at order 6 (as differential operators). Further there is an algebraic relation at order 8 expressing the fact that there are only 5 algebraically independent generators. We work out the details of modeling physically relevant irreducible representations of the quadratic algebra in terms of divided difference operators in two variables. We determine several ON bases for this model including spherical and cylindrical bases. These bases are expressed in terms of two variable Wilson and Racah polynomials with arbitrary parameters, as defined by Tratnik. The generators for the quadratic algebra are expressed in terms of recurrence operators for the one-variable Wilson polynomials. The quadratic algebra structure breaks the degeneracy of the space of these polynomials. In an earlier paper the authors found a similar characterization of one variable Wilson and Racah polynomials in terms of irreducible representations of the quadratic algebra for the quantum superintegrable system on the 2-sphere with generic 3-parameter potential. This indicates a general relationship between 2nd order superintegrable systems and discrete orthogonal polynomials.

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 quantum systems H = ∆ n +Ṽ . These systems, including the classical Kepler [1] and anisotropic oscillator systems and the quantum anisotropic oscillator and hydrogen atom have great historical importance, due to their remarkable properties [2,3,4,5,6]. One modern practical application among many is the Hohmann transfer, a fundamental tool for the positioning of earth satellites and for celestial navigation in general, which is based on the superintegrability of the Kepler system [7]. 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 order of a quantum superintegrable system is the maximum order of the quantum symmetries as differential operators.
The potential V corresponding to a 2nd order superintegrable system, classical or quantum, on an n-dimensional conformally flat manifold depends linearly on several parameters in general and can be shown to generate a vector space of dimension ≤ n + 2. (One dimension corresponds to the trivial addition of a constant to the potential and usually isn't included in a parameter count.) If the maximum is achieved, the potential is called nondegenerate. There is an invertible mapping between superintegrable systems on different manifolds, called the Stäckel transform, which preserves the structure of the algebra generated by the symmetries. In the cases n = 2, 3 it is known that all nondegenerate 2nd order superintegrable systems are Stäckel equivalent to a system on a constant curvature space [30,31]. An important fact for 2D systems is that all systems can be obtained from one generic superintegrable system on the complex 2-sphere by appropriately chosen limit processes, e.g. [32,33]. The use of these processes in separation of variables methods for wave and Helmholtz equations in n dimensions was pioneered by Bôcher [34]. For n = 3 it appears that all nondegenerate 3D systems can be obtained from one generic superintegrable system on the complex 3-sphere by similar limiting processes, but the proof is not yet complete [11,35].
For n = 2 we define the generic sphere system by embedding of the unit 2-sphere x 2 1 +x 2 2 +x 2 3 = 1 in three dimensional flat space. Then the Hamiltonian operator is The 3 operators that generate the symmetries are L 1 = L 12 , L 2 = L 13 , L 3 = L 23 where Here .
From the general structure theory for 2D 2nd order superintegrable systems with nondegenerate potential we know that the 3 defining symmetries will generate a symmetry algebra (a quadratic algebra) by taking operator commutators, which closes at order 6, [36]. That is, all possible symmetries can be written as symmetrized operator polynomials in the basis generators and in the 3rd order commutator R, where R occurs at most linearly. In particular, the dimension of the space of truly 2nd order symmetries for the Hamiltonian operator is 3, for the 3rd order symmetries it is 1, for the 4th order symmetries it is 6, and for the 6th order symmetries it is 10. For the generic 2-sphere quantum system the structure equations can be put in the symmetric form [12] ǫ ijk [L i , R] = 4{L i , L k } − 4{L i , L j } − (8 + 16a j )L j + (8 + 16a k )L k + 8(a j − a k ), (1.1) Here ǫ ijk is the pure skew-symmetric tensor, R = [L 1 , L 2 ] and {L i , L j } = L i L j + L j L i with an analogous definition of {L 1 , L 2 , L 3 } as a symmetrized sum of 6 terms. In practice we will substitute L 3 = H − L 1 − L 2 − a 1 − a 2 − a 3 into these equations.
In [12] we started from first principles and worked out some families of finite and infinite dimensional irreducible representations of the quadratic algebra with structure relations (1.1), (1.2), including those that corresponded to the bounded states of the associated quantum mechanical problem on the 2-sphere. Then we found 1-variable models of these representations in which the generators L i acted as divided difference operators in the variable t on a space of polynomials in t 2 . The eigenfunctions of one of the operators L i turned out to be the Wilson and Racah polynomials in their full generality. In essence, this described an isomorphism between the quadratic algebra of the generic quantum superintegrable system on the 2-sphere and the quadratic algebra generated by the Wilson polynomials.
The present paper is concerned with the extension of these results to the 3-sphere, where the situation is much more complicated. From the general structure theory for 3D 2nd order superintegrable systems with nondegenerate potential we know that although there are 2n − 1 = 5 algebraically independent 2nd order generators, there must exist a 6th 2nd order symmetry such that the 6 symmetries are linearly independent and generate a quadratic algebra that closes at order 6 [37]. (We call this the 5 =⇒ 6 Theorem.) Thus, all possible symmetries can be written as symmetrized operator polynomials in the basis generators and in the four 3rd order commutators R i , where the R i occur at most linearly. In particular, the dimension of the space of truly 2nd order symmetries is 3, for the 3rd order symmetries is 4, for the 4th order symmetries it is 21, and for the 6th order symmetries it is 56. In 3D there are 5 algebraically independent, but 6 linearly independent, generators. The algebra again closes at 6th order, but in addition there is an identity at 8th order that relates the 6 algebraically dependent generators. The representation theory of such quadratic algebras is much more complicated and we work out a very important instance of it here. In this case we will find an intimate relationship between these representations and Tratnik's 2-variable Wilson and Racah polynomials in their full generality [38,39,40].
For nD nondegenerate systems there are 2n − 1 functionally independent but n(n + 1)/2 linearly independent generators for the quadratic algebra. We expect that the relationships developed here will extend to n-spheres although the results will be of increasing complexity.
2 The quantum superintegrable system on the 3-sphere We define the Hamiltonian operator via the embedding of the unit 3-sphere A basis for the second order constants of the motion is In the following i, j, k, ℓ are pairwise distinct integers such that 1 ≤ i, j, k, ℓ ≤ 4, and ǫ ijk is the completely skew-symmetric tensor such that ǫ ijk = 1 if i < j < k. There are 4 linearly independent commutators of the second order symmetries (no sum on repeated indices): This implies, for example, that Here we define the commutator of linear operators F , G by [F, G] = F G − GF . The structure equations can be worked out via a relatively straightforward but tedious process. We get the following results.
The fourth order structure equations are Here {F, G} = F G + GF . The fifth order structure equations (obtainable directly from the fourth order equations and the Jacobi identity) are The sixth order structure equations are Here {A, B, C} = ABC + ACB + BAC + BCA + CAB + CBA.
The eighth order functional relation is Here {A, B, C, D} is the 24 term symmetrizer of 4 operators and the sum is taken over all pairwise distinct i, j, k, ℓ. For the purposes of the representation, it is useful to redefine the constants as a i = b 2 i − 1 4 . We note that the algebra described above contains several copies of the algebra generated by the corresponding potential on the two-sphere. Namely, let us define A to be the algebra generated by the set {L ij , I} for all i, j = 1, . . . , 4 where I is the identity operator. Then, we can see that there exist subalgebras A k generated by the set {L ij , I} for i, j = k and that these algebras are exactly those associated to the 2D analog of this system. Furthermore, if we define then H k will commute with all the elements of A k and will represent the Hamiltonian for the associated system. For example, take A 4 to be the algebra generated by the set {L 12 , L 13 , L 23 , I}.
In this algebra, we have the operator 3 )I which is in the center of A 4 and which is the Hamiltonian for the associated system on the two sphere immersed in Next we construct families of finite dimensional and infinite dimensional bounded below irreducible representations of this algebra that include those that arise from the bound states of the associated quantum mechanical eigenvalue problem. At the same time we will construct models of these representations via divided difference operators in two variables s and t. Important tools for this construction are the results of [12] giving the representations of the A k 's and known recurrence relations for one-variable Wilson and Racah polynomials.

Review of Wilson polynomials
Before we proceed to the model, we us present a basic overview of some of the characteristics of the Wilson polynomials [41] that we plan to employ in the creation of our model. The polynomials are given by the expression w n t 2 ≡ w n t 2 , α, β, γ, δ = (α + β) n (α + γ) n (α + δ) n where (a) n is the Pochhammer symbol and 4 F 3 (1) is a generalized hypergeometric function of unit argument. The polynomial w n (t 2 ) is symmetric in α, β, γ, δ. The Wilson polynomials are eigenfunctions of a divided difference operator given as where See [42] for a simple derivation. The Wilson polynomials Φ n (t 2 ) ≡ Φ (α,β,γ,δ) n (t 2 ), satisfy the three term recurrence formula
Finally, the weight function of the model will be based on a two dimensional generalization of the weight function of the Wilson polynomials.
For fixed α, β, γ, δ > 0 (or if they occur in complex conjugate pairs with positive real parts) [41], the Wilson polynomials are orthogonal with respect to the inner product When m is a nonnegative integer then α + β = −m < 0 so that the above continuous Wilson orthogonality does not apply. The representation becomes finite dimensional and the orthogonality is a finite sum Thus, the spectrum of the multiplication operator t 2 is the set {(α + k) 2 : k = 0, . . . , m}. Now, we are ready to determine the model.

Construction of the operators for the model
To begin, we review some basic facts about the representation. The original quantum spectral problem for (2.1) was studied in [43] from an entirely different point of view. It follows from this study that for the finite dimensional irreducible representations of the quadratic algebra the multiplicity of each energy eigenspace is (M + 2)(M + 1)/2 and we have where I is the identity operator. Of course, for an irreducible representation, the Hamiltonian will have to be represented by a constant times the identity and initially for the construction of the model, we assume We will obtain the quantized values of E from the model.
We recall that each operator L ij is a member of the subalgebras A k for k = i, j. Thus, we can use the known representations of these algebras, and symmetry in the indices, to see that the eigenvalues of each operator will be associated with eigenfunctions φ h,m indexed by integers 0 ≤ h ≤ m so that 4.1 A basis for L 13 , L 12 + L 13 + L 23 As described above, we seek to construct a representation of A by extending the representations obtained for the subalgebras A k . The most important difference for our new representation is that the operator H 4 = L 12 + L 13 is in the center of A 4 but not A. Hence, it can no longer be represented as a constant. We can still use the information about its eigenvalues to make an informed choice for its realization.
Restricting to bounded below irreducible representations of the quadratic algebra initially, we see from the representations of A 4 that the possible eigenvalues of H 4 are given as in (4.3) and the eigenvalues of L 13 are given as in (4.2).
We can begin our construction of a two-variable model for the realization of these representations by choosing variables t and s, such that i.e., the action of these operators is multiplication by the associated transform variables. From the eigenvalues of the operators, we can see that the spectrum of In this basis, the eigenfunctions d ℓ,m for a finite dimensional representation are given by delta functions

A basis for
Next, we construct L 12 in the model. Let f n,m be a basis for the model corresponding to simultaneous eigenvalues of L 12 , L 12 + L 13 + L 23 . From the representations of A 4 [12], we know that the action of L 13 on this basis is given by .
We already know that the bounded below representations of A 4 are intimately connected with the Wilson polynomials. The connection between these polynomials and the representation theory is the three term recurrence formula (4.4) for the action of L 13 on an L 12 basis, where the coefficients are given by (4.5) and (4.6).
We define the operator L on the representation space of the superintegrable system by the action of the three term recurrence relations for the Wilson polynomials given by expansion coefficients (3.2)-(3.4), i.e.
Note that with the choices we have a perfect match with the action of L 13 as Thus, the action of L 13 is given by and so we see that the action of L 13 on an L 12 basis is exactly the action of the variable t 2 on a basis of Wilson polynomials. Hence, we hypothesize that L 12 takes the form of an eigenvalue operator for Wilson polynomials in the variable t where τ , τ * t are given as (3.1) with the choice of parameters as given in (4.7). Here the subscript t expresses the fact that this is a difference operator in the variable t, although the parameters depend on the variable s.
The basis functions corresponding to diagonalizing H 4 and L 12 can be taken, essentially, as the Wilson polynomials where s m = m + 1 + (b 1 + b 2 + b 3 )/2 as above. Note that w n (t 2 ) actually depends on m (or s 2 ) through the parameters α, δ. Also α + δ is independent of m. Written in terms of the variable s, the parameters are given by Note that when s is restricted to s m , these parameters agree with (4.7).
Since the w n are symmetric with respect to arbitrary permutations of α, β, γ, δ, we can transpose α and β and verify that w n is a polynomial of order n in s 2 .

A basis for L 13 , L 24
For now, let us assume that we have a finite dimensional irreducible representation such that the simultaneous eigenspaces of L 12 , L 12 + L 13 + L 23 are indexed by integers n, m, respectively, such that 0 ≤ n ≤ m ≤ M . Each simultaneous eigenspace is one-dimensional and the total dimension of the representation space is (M + 1)(M + 2)/2. Now we need to determine the action of the operators L 14 , L 24 , L 34 in the model.
A reasonable guess of the form of the operator L 24 is as a difference operator in s, since it commutes with L 13 . We hypothesize that it takes the form of an eigenvalue equation for the Wilson polynomials in the variable s. We require that it have eigenvalues of the form (4.2).
Note that when acting on the delta basis d ℓ,m , it produces a three-term recursion relation. For our representation, we require that that the representation cut off at the appropriate bounds. That is if we write the expansion coefficients of L 24 acting on d ℓ,m , as For L 24 we take Hereτ s is the difference operator in s where the parameters areα,β,γ,δ.
With the operator L 24 thus defined, the unnormalized eigenfunctions of the commuting operators L 13 , L 24 in the model take the form g n,k where For this choice of parameters, the functions (4.10) constitute an alternative basis for the representation space, consisting of polynomials in s 2 , t 2 multiplied by a delta function in s.

Completion of the model
In this section, we finalize the construction of our model by realizing the operator L 34 . The operator L 34 must commute with L 12 , so we hypothesize that it is of the form where S u f (s, t) = f (s + u, t), A, B, C, D are rational functions of s to be determined, and the operators L αβ , R αβ , L, R, etc. are defined in Appendix B. The subscript t denotes difference The parameters are (4.8). Here On the other hand, we can consider the action of L 34 on the basis (4.10). Considering L 34 primarily as an operator on s we hypothesize that it must be of the form where the difference operators are defined in Appendix B with subscript s denoting difference operators in s and κ is a constant.
Finally, we express the operator L 14 as By a long and tedious computation we can verify that the 3rd order structure equations are satisfied if and only if E takes the values and the functional coefficients for L 34 in (4.11), (4.12) take the following form : and κ = −4. The expression forẼ(t) takes the formẼ(t) = µ 1 + µ 2 /(4t 2 − 1) where µ 1 , µ 2 are constants, but we will not list it here in detail. For finite dimensional representations, we have the requirement that M be a positive integer so we obtain the quantization of the energy obtained previously (4.1).

The model and basis functions
We shall now review what we have constructed, up to this point. We realize the algebra A by the following operators where the parameters for the τ t operators are given in (4.8), the parameters for the operatorsτ s are given in (4.9) and the functional coefficients of L 34 are given in (4.13 We also have a nonorthogonal basis given by Recall that the spectrum of the variables s, t is given by We finish the construction of the model by computing normalizations for the basis f n,m , and g ℓ,m and the weight function.

The weight function and normalizations
We begin this section by determining the weight function and normalization of the basis functions in the finite dimensional representations. Later, we shall extend the system to the infinite dimensional bounded below case.

The weight function and normalization
of the basisd ℓ,m (s, t) = δ(t − t ℓ )δ(s − s m ) We consider the normalization for the d ℓ,m = δ(t − t ℓ )δ(s − s m ) basis for finite dimensional representations where In order to derive these results we use the requirement that the generating operators L ij are formally self-adjoint. Consider a weight function ω(t, s) so that f (t, s), g(t, s) = f (t, s)g(t, s)ω(t, s)dsdt, then we assume that the basis functions are orthonormal with which implies that c 2 ℓ,m ω(t ℓ , s m ) = 1. The adjoint properties of L 13 and L 24 provide recurrence relations on the c ℓ,m . That is Similarly, the self-adjoint property of L 24 implies the recurrence relation .  When evaluated at s = s m , the parameters are given by and satisfy α + β = −m < 0. Thus, the Wilson orthogonality is realized as a finite sum over the weights of t 2 . However, the weight of the variable t is given by t ℓ = ℓ + β and we must adjust the equation for the Wilson orthogonality (3.6) by permuting α and β. This is allowed since the polynomial and the requirement α + β = −m are symmetric in the two parameters. In this form the Wilson orthogonality is given over the spectrum of the multiplication operator t 2 as the set {(β + ℓ) 2 : ℓ = 0, . . . , m} In light of this orthogonality, we hypothesize that the weight function is given by and so we look for normalization constants so that The orthogonality (5.5) in terms of the choices of parameters (5.4) is given by The weight function (5.3) can be rewritten as We can now solve the equation (5.6) for k n,m by comparing (5.8) and (5.7) to obtain With this normalization the basis functionsf n,m (s, t) are orthonormal.

Normalization of the w k (s 2 )δ(t − t ℓ ) basis
Next, we use the orthogonality of the Wilson polynomials to find the normalization of the g n,k basis in the finite dimensional representation. We take the normalized basis functions to be given bŷ Again, we want to show that there exist normalization constants h ℓ,k so that the following holds: When restricted to t = t ℓ the parametersα,β,γ,δ becomẽ where the index being summed over is m = ℓ, . . . , M instead of m − ℓ = 0, . . . , M − ℓ.
Comparing this orthogonality with the weight function (5.3) written as , the normalization constants are determined by the requirement for l ℓ,k . The proper choice of normalization is .

The norm of 1
Throughout the previous analysis, the weight and normalization have been proportional to an arbitrary overall factor c 0,0 . We can fix this constant by requiring that the function 1, belonging to the basis h n,k (s, t) of monomials in s 2 , t 2 , be normalized to length 1. We compute this by using the Wilson orthogonality for the 0th order Wilson polynomials (i.e. (3.6) with k, k ′ = 0). The norm of 1 is given by Evaluating the double sum gives Thus, setting will make 1, 1 = 1.

Infinite dimensional representations
For infinite dimensional but bounded below representations with −m and −M + ℓ nonnegative integers and all b j real we take the inner product of two functions f (t 2 , s 2 ), g(t 2 , s 2 ) in the form To compute the measure ω(t, s) such that our operators L jk are formally self-adjoint we use the fact that we already know the restriction of the measure to the s-constant basis (4.15) and the t-constant basis (4.16). For consistency, we see that the weight function should be

s).
Then we can compute the norm square of the constant function f (t 2 , s 2 ) = 1 by using (3.5) twice to evaluate the iterated integral: Comparing the measures that we have derived for the infinite dimensional and finite dimensional cases with the 2-variable Wilson and Racah polynomials introduced by Tratnik [38,39,40] we see that they agree. Thus we have found two-variable Wilson and Racah polynomials in complete generality.

Expansion coef f icients
We can easily determine the coefficients for the expansion of one of our bases in terms of another.
Here we write the expansion coefficients in terms of the unnormalized functions. The expansion of the d ℓ,m basis in terms of the f n,m basis is given by The expansion of the d ℓ,m basis in terms of the g ℓ,k basis is given by Finally, the expansion of the f n,m basis in terms of the g ℓ,k basis is given by In order to understand the significance of these expansions in quantum theory, it is useful to consider the results of [43]. There the Schrödinger eigenvalue problem for the generic potential on the n-sphere was considered, for general n. For n = 3 it was shown that all of the eigenfunctions of the pairs commuting operators treated in this paper separated in some version of either spherical or cylindrical coordinates and were expressible as continuous multivariable orthogonal polynomials orthogonal on a simplex. Thus the expansion coefficients derived here represent the expansion of one basis of solutions of the Schrödinger eigenvalue equation in terms of another.
6.1 A basis for L 12 , L 12 + L 14 + L 24 Now that we have computed the measures for our spaces of polynomials from first principles and established that they agree with those for the Tratnik generalization of Wilson and Racah polynomials to two variables [38,39], we can make use of known results for the Tratnik case to compute another ON basis for our spaces. In an appendix to [40] the authors show that the true 2-variable Racah polynomials defined by Tratnik are simultaneous eigenfunctions of two commuting difference operators L 1 , L 2 . We will identify these operators with our symmetry algebra and verify another eigenbasis for our representation space.
We construct the polynomials of Tratnik [38,39] and operators given in [40] via the definitions The original form of the polynomials were given in terms of the Racah polynomials which can be related to the Wilson polynomials via Then, the two-variable extension of the Wilson polynomials defined by Tratnik are given by equation (3.10) of [40] as with the requirements that 0 ≤ n 1 ≤ n 1 + n 2 ≤ M . We can express the 2-variable polynomial R 2 in terms of the Wilson polynomials using the original parameters and variables of the model as R 2 (n 1 , n 2 ; b i ; t, s; M ) = w n 1 α, β, γ, δ; t 2 (6.1) where as in (4.8) In particular, note that the parameters α, δ depend on s and so the polynomial w n 1 is a function of both s and t.
Note that it was already demonstrated in Section 4.2 that the polynomial, w n 1 , is an eigenfunction of L 12 . Furthermore, it is easy to see that w n 2 depends only on s and so will be left invariant by L 12 and so the 2-variable polynomial R 2 is an eigenfuction for L 12 . As was exhibited in [40], there is a set of two commuting difference operators whose simultaneous eigenfunctions are just these orthogonal polynomials. It is then natural to expect that these operators can be expressed in terms of the operators in our model which commute with L 12 , i.e., I, L 12 , L 34 , H 3 , and H 4 .
The commuting difference operators are given as follows, via [40]. Let I i be the operator which maps x i to −x i − β i and leaves fixed x j for j = i. Similarly, define E a x 1 as the operator which maps x i to x i + a and leaves fixed x j for j = i. We define functions B j,k i as and further extend these functions for k = −1, 0, 1 via Let ν be some multi-index ν = (ν 1 , ν 2 ) with ν i = −1, 0, 1 and µ be a single index µ = −1, 0, 1. Then the functions given by are enough to define the operators and describe the results of [40]. The operators commute and their eigenfunctions are given by R 2 (6.1). The eigenvalues of the operators L 1 , L 2 are given by so we hypothesize that L 1 is a linear combination of L 12 and the identity and L 2 is a linear combination of H 3 and the identity. In fact, it is straightforward to verify that The normalization of this basis can be found in [39] and [40]. Thus, we have shown that Tratnik's version of two-variable Racah polynomials corresponds to the L 12 , L 12 + L 14 + L 24 eigenbasis.

Conclusions and discussion
We have demonstrated explicitly the isomorphism between the quadratic algebra of the generic quantum superintegrable system on the 3-sphere and the quadratic algebra generated by the recurrence relations for two-variable Wilson [43].
Natural questions here are: what is the origin of these models of the symmetry algebra action and how can we determine when there is a differential operator model, a difference operator model or some other model? Clearly, the models are associated with the spectral resolutions of systems of commuting operators in the symmetry algebra. In [13] we showed how difference and differential operator models can be suggested by analysis of the corresponding classical systems, and these ideas are relevant here. In [44] we developed a recurrence relation approach for differential operators that allowed us to derive difference equation models for 2D quantum systems and, again, this approach should generalize to 3D quantum systems. Also, there is an obvious connection between the existence of models and bispectrality [40]. Another issue is that all models that we know of for quantum symmetry algebras of 2D and 3D superintegrable systems are associated with commuting operators whose simultaneous separated eigenfunctions are of hypergeometric type. Do there exist models with commuting operators whose simultaneous separated eigenfunctions are not hypergeometric?
It is suggested by our method that most of the quadratic algebras for all Stäckel equivalence classes of 3D second order quantum superintegrable systems on conformally flat spaces should be obtainable by appropriate limit processes from the quadratic algebra associated with the generic superintegrable system on the 3-sphere, namely that generated by the two-variable Wilson polynomials. However these limit processes are very intricate, see e.g. [33], and each equivalence class exhibits unique structure, so each class is important for study by itself. Moreover, within each class of Stäckel equivalent systems the structure of the quadratic algebra remains unchanged but the spectral analysis of the generators for the algebra can change. We conjecture that this limiting process for superintegrable quantum systems is analogous to the Askey scheme for obtaining various families of orthogonal polynomials as limits of Askey-Wilson polynomials.
As an example of this, in the paper [45] we studied the quadratic algebra associated with the quantum 3D caged isotropic oscillator. There, the Hamiltonian operator was and a basis for the second order constants of the motion was (with H = M 1 + M 2 + M 3 ) We found 3 two-variable models for physically relevant irreducible representations of the quadratic algebra. One was in terms of differential operators and led to monomial eigenfunctions for the generators that corresponded to separation of variables in Cartesian coordinates, one was in terms of mixed differential-difference operators and led to one-variable dual Hahn polynomial eigenfunctions for the generators that corresponded to separation of variables in cylindrical coordinates, and the third was in terms of pure difference operators and led to one-variable Wilson or Racah polynomial eigenfunctions for the generators that corresponded to separation of variables in spherical coordinates. It can be shown that the flat space caged isotropic oscillator system can be obtained as a limit of the generic system on the sphere, whereas at the quadratic algebra level one variable dual Hahn and Wilson polynomials can be obtained as limits of twovariable Wilson polynomials.
For nD nondegenerate superintegrable systems on conformally flat spaces there are 2n − 1 functionally independent but n(n + 1)/2 linearly independent generators for the quadratic algebra. It is reasonable to conjecture that the quadratic algebra of the generic potential on the n-sphere is uniquely associated with the (n − 1)-variable version of Tratnik's multivariable Wilson polynomials.
Finally, these results suggest the existence of a q version of superintegrability for quantum systems [46].

B Recurrence relations for construction of the spherical and cylindrical models
The spherical and cylindrical models are associated with the 8 basic raising and lowering operators for the Wilson polynomials, as well as the three term recurrence relation. We list these operators here and describe their actions on the basis polynomials Φ n ≡ Φ (α,β,γ,δ) n .
In addition we can use the three term recurrence (B.1) and multiplication by the operator y 2 to both raise and lower n by 1 while fixing the other parameters.