Structure Theory for Extended Kepler-Coulomb 3D Classical Superintegrable Systems

The classical Kepler-Coulomb system in 3 dimensions is well known to be 2nd order superintegrable, with a symmetry algebra that closes polynomially under Poisson brackets. This polynomial closure is typical for 2nd order superintegrable systems in 2D and for 2nd order systems in 3D with nondegenerate (4-parameter) potentials. However the degenerate 3-parameter potential for the 3D extended Kepler-Coulomb system (also 2nd order superintegrable) is an exception, as its quadratic symmetry algebra doesn't close polynomially. The 3D 4-parameter potential for the extended Kepler-Coulomb system is not even 2nd order superintegrable. However, Verrier and Evans (2008) showed it was 4th order superintegrable, and Tanoudis and Daskaloyannis (2011) showed that in the quantum case, if a second 4th order symmetry is added to the generators, the double commutators in the symmetry algebra close polynomially. Here, based on the Tremblay, Turbiner and Winternitz construction, we consider an infinite class of classical extended Kepler-Coulomb 3- and 4-parameter systems indexed by a pair of rational numbers $(k_1,k_2)$ and reducing to the usual systems when $k_1=k_2=1$. We show these systems to be superintegrable of arbitrarily high order and work out explicitly the structure of the symmetry algebras determined by the 5 basis generators we have constructed. We demonstrate that the symmetry algebras close rationally; only for systems admitting extra discrete symmetries is polynomial closure achieved. Underlying the structure theory is the existence of raising and lowering constants of the motion, not themselves polynomials in the momenta, that can be employed to construct the polynomial symmetries and their structure relations.

1 Introduction with S 1 = H and [H, S j ] ≡ HS j − S j H = 0, apparently the maximum number possible. The system is of order if the maximum order of the symmetry operators, other than the Schrödinger operator, is . Similarly, a classical superintegrable system with Hamiltonian on phase space with local coordinates x j , p j , where ds 2 = g jk dx j dx k is an integrable system such that there are 2n − 1 functionally independent functions polynomial in momenta, (easily provable to be the maximum number possible): S j (p, x), j = 1, . . . , 2n − 1, with S 1 = H, and globally defined such that {S j , H} = 0, where is the Poisson bracket. The system is of order if the maximum order of the generating constants of the motion is . As has been pointed out many times [2,22] such systems are of enormous historic and present day practical importance. In essence, superintegrable systems are those Hamiltonian systems that can be "solved" exactly, analytically and algebraically, without requiring numerical approximation. Superintegrable systems are used as the basis for exact and perturbation methods that underlie planetary motion determination, orbital maneuvering, the periodic table of the elements, the boson calculus and much of special function theory, [14,15].
The key property that makes a system "superintegrable" is that, in contrast to merely integrable systems, the symmetry algebra generated by the basis symmetries is nonabelian. This nonabelian structure can be analyzed and used to deduce properties of the system. Thus in the quantum case the irreducible representations of the symmetry algebra determine the multiplicities of the degenerate energy eigenspaces and permit algebraic computation of the eigenvalues.
The classical Kepler-Coulomb system in 3 dimensions is well known to be 2nd order superintegrable, with a symmetry algebra that closes polynomially in an so(4)-like structure, e.g. [1]. This polynomial closure (though not usually a Lie algebra) is typical for 2nd order superintegrable systems in 2D [8,11] and for 2nd order systems in 3D with nondegenerate (4-parameter) potentials. However the degenerate 3-parameter potential for the 3D extended Kepler-Coulomb system (also 2nd order superintegrable) is an exception, as its quadratic symmetry algebra doesn't close polynomially [7]. We write it in the form where x, y, z are the usual Cartesian coordinates with conjugate momenta p x , p y , p z in phase space and r = x 2 + y 2 + z 2 . The 3D 4-parameter extended Kepler-Coulomb system is not even 2nd order superintegrable. We write it in the form However, Verrier and Evans [28] showed it was 4th order superintegrable, and Tanoudis and Daskaloyannis [21] showed in the quantum case that, if a second 4th order symmetry is added to the generators, the symmetry algebra closes polynomially in the sense that all second commutators of the generators can be expresses as symmetrized polynomials in the generators. (Note that the 3-parameter potential is not just a restriction of the 4-parameter case, because it admits symmetries that are not inherited from the 4-parameter symmetries.) Here we introduce an analog of the TTW construction [23,24] and consider an infinite class of classical extended Kepler-Coulomb 3-and 4-parameter systems indexed by a pair of rational numbers (k 1 , k 2 ) and reducing to the usual systems when k 1 = k 2 = 1. We construct explicitly a set of generators, show these systems to be superintegrable of arbitrarily high order and determine the structure of the generated symmetry algebras. We demonstrate that the symmetry algebras close rationally; only for systems admitting extra discrete symmetries is polynomial closure achieved. Much of the paper is quite technical but, as this is the first work devoted to uncovering the structure of the symmetry algebras of 3D superintegrable systems of arbitrary order, we think it is important to expose the details of computations and concepts that later may prove to be routine.
For the 4-parameter system in the case k 1 = k 2 = 1, where discrete symmetry is present and a 6th generator is needed to obtain polynomial closure, we work out the 12th order functional relationship between the 6 generators.
In Section 2 we review the action angle construction that we employ to show that our systems are superintegrable and to enable the determination of the structure of the symmetry algebra. In Section 3 we use the fact that the 3-parameter Kepler-Coulomb system (1.1) separates in spherical coordinates r, θ 1 , θ 2 and, by replacing the angles by k 1 θ 1 , k 2 θ 2 for k 1 , k 2 rational, define an infinite family of extended Kepler-Coulomb systems, no longer restricted to flat space. We demonstrate that each of these systems is superintegrable, but of arbitrarily high order. We use our method of raising and lowering symmetries to determine the structure of the symmetry algebras generated by these systems. The general construction does not yield generators of minimum order and we show in Section 3.2 how a limit argument exposes the minimum order generators. Underlying the structure theory is the existence of raising and lowering constants of the motion, not themselves polynomials in the momenta, that can be employed to construct the polynomial symmetries and their structure relations. We show that the symmetry algebra closes rationally, not polynomially.
In Section 4 we apply our method to the 4-parameter Kepler-Coulomb system (1.2) and use its separation in spherical coordinates r, θ 1 , θ 2 . By replacing the angles by k 1 θ 1 , k 2 θ 2 for k 1 , k 2 rational, we define an infinite family of extended Kepler-Coulomb systems, again not restricted to flat space. We demonstrate that each of these systems is superintegrable, but of arbitrarily high order. We use raising and lowering symmetries to determine the structure of the symmetry algebras generated by these systems. Again the general construction does not yield generators of minimum order and we show in Section 5 how a limit argument exposes the minimum order generators. In general, the symmetry algebra closes rationally, but not polynomially. However in two cases k 1 = k 2 = 1 and k 1 = k 2 = 1/2 the system admits additional symmetry: the 6-element permutation group S 3 . The general construction shows that these systems admit 5 generators of orders 2, 2, 2, 2, 4, but in Section 5.2 we show that the permutation symmetry implies the existence of a 6th symmetry of order 4 such that the 6 symmetries are linearly independent. Then we demonstrate that the algebra generated by these 6 symmetries closes polynomially, in analogy with the computation for the quantum case in [21]. We go further and work out the 12th order functional relation between the 6 generators.
In Section 6 we present our conclusions and prospects for additional research.

Review of the action-angle construction
In [6,9,10,12,27] it was described how to determine a complete set of 2n − 1 functionally independent constants of the motion for a classical Hamiltonian on an n-dimensional Riemannian or pseudo-Riemannian manifold whose Hamilton-Jacobi equation separates in an orthogonal sub-group coordinate system. In the special case n = 3 the defining equations for the Hamiltonian H, expressed in the separable coordinates q 1 , q 2 , q 3 , take the form The additional constants of the motion can be constructed as Here, and L 4 ≡ 0. With this construction the functions L 1 , L 2 , L 3 , L 1 , L 2 are functionally independent constants of the motion. The functions L 2 , L 3 are second order polynomials in the momenta and determine the separation of variables in coordinates q 1 , q 2 , q 3 . In general the functions L 1 , L 2 are only locally defined and are not polynomials. The system will be superintegrable only if we can supplement L 1 , L 2 , L 3 with two more polynomial functions such that the full set is functionally independent. This will be possible only for very special systems In the following we will look at several candidate systems for superintegrability, show how to construct polynomial constants of the motion from L 1 , L 2 and work out the structure of the symmetry algebra generated by these constants.
3 The classical 3D extended Kepler-Coulomb system with 3-parameter potential The extended Kepler-Coulomb Hamiltonian is Here, L 2 , L 3 are constants of the motion that determine additive separation of the Hamilton-Jacobi equation. Further {L 2 , L 3 } = 0 so L 2 and L 3 are in involution. Applying our action angle construction to get two independent constants of the motion we note that q 1 = r, q 2 = θ 1 , q 3 = θ 2 and to obtain functions A j , B j , j = 1, 2 such that Here, k 1 = p 1 /q 1 , k 2 = p 2 /q 2 where p 1 , q 1 are relatively prime positive integers and p 2 , q 2 are relatively prime positive integers. From our general theory, are two constants of the motion such that the full set of five constants of the motion is functionally independent.
We work with the exponential functions, [4], see also [25,26]. We have, for j = 1, 2, where X 1 = sin(k 1 θ 1 )p θ 1 − i L 2 cos(k 1 θ 1 ), X 1 = sin(k 1 θ 1 )p θ 1 + i L 2 cos(k 1 θ 1 ), (Here, X, Y are, in general, not the complex conjugates of X, Y , respectively, unless all of the coordinates are real.) Now note that e q 1 A 1 −p 1 B 1 and e −q 1 A 1 +p 1 B 1 are constants of the motion, where Moreover, the identity e q 1 A 1 −p 1 B 1 e −q 1 A 1 +p 1 B 1 = 1 can be expressed as where P 1 is a polynomial in H, L 2 and L 3 .
Similarly, e p 1 q 2 A 2 −p 2 q 1 B 2 and e −p 1 q 2 A 2 +p 2 q 1 B 2 are constants of the motion, where The identity e p 1 q 2 A 2 −p 2 q 1 B 2 e −p 1 q 2 A 2 +p 2 q 1 B 2 = 1 can be written as where P 2 is a polynomial in H, L 2 and L 3 . Let a, b, c, d be nonzero complex numbers and consider the binomial expansion Here, either k = 2 or k = 3 and p, q are positive integers. Suppose p + q is odd. Then it is easy to see that the sum (3.1) takes the form √ L k T odd (L k ) where T odd is a polynomial in L k . On the other hand, if p + q is even then the sum (3.1) takes the form T even (L k ) where T even is a polynomial in L k .
Similarly, consider the binomial expansion On the other hand, if p + q is even then the sum (3.2) takes the form A third possibility is Then we must have p even to get a polynomial in L k . If p is odd the sum takes the form Then we must have p odd to get a polynomial in L k . If p is even the sum takes the form We define basic raising and lowering symmetries At this point we restrict to the case where each of p 1 , q 1 , p 2 , q 2 is an odd integer. (The other cases are very similar.) Let Then we see from the explicit expressions for the symmetries and the preceding parity argument that J 1 , J 2 , K 1 , K 2 are constants of the motion, polynomial in the momenta. Moreover, the identities The following relations are straightforward to derive from the definition of the Poisson bracket: From these results, we find Thus we obtain Similarly, Further, Then, using (3.3), we conclude that Similarly, we have the K-related identities Commutators relating the J and K symmetries are somewhat more complicated to compute. We have Now we can determine the nonpolynomial constant of the motion {J + , K + }: where the last term in braces vanishes identically. We conclude that Once we have {J + , K + } explicitly, we can obtain the remaining Poisson relations between the J and K symmetries with little additional work. We use the fact that J + J − = P 1 and We can write this relation in the more compact form Similarly, we have Thus, In summary: These relations prove closure of the symmetry algebra in the space of functions polynomial in J ± , K ± , rational in L 2 , L 3 , H and at most linear in

Structure relations for polynomial constants of the motion
Since, we have

Similar computations yield
Note: It can be verified that the numerators are divisible by are true polynomial constants of the motion, although not polynomial in the generators.

Minimal order generators
The generators for the polynomial symmetry algebra that we have produced so far are not of minimal order. Note that L 1 , L 2 , L 3 are of order 2 and the orders of J 1 , K 1 are one less than the orders of J 2 , K 2 , respectively. We will construct a symmetry K 0 of order one less than K 1 . Note that the symmetry K 2 is a polynomial in L 3 . The constant term in this polynomial expansion , itself a constant of the motion. In the case we are considering p 1 , q 1 , p 2 , q 2 are each odd, so the constant term is is a polynomial symmetry of order two less than K 2 . We have the identity The same construction fails for J 2 . It is a polynomial in L 2 , but the constant term in the expansion is not a constant of the motion. Indeed, in the special case k 1 = k 2 = 1, the symmetry J 1 is of minimal order 2, so J 1 cannot be realized as a commutator. Now we choose L 1 , L 2 , L 3 , J 1 , K 0 as the generators of our algebra. We define the basic nonzero commutators as Then we have a polynomial in the generators. Further, which again can be verified to be a polynomial in the generators. Note, however, that R 1 R 2 = 8p 2 1 p 2 J 2 K 1 , a product of Poisson brackets of the generators, is not a polynomial in the generators, although (R 1 R 2 ) 2 is such a polynomial. Using the identity (3.4) and the expression for {J 1 , K 2 }, it is easy to see that R 3 is rationally related to R 2 . It is clear that all additional commutators can be expressed as rational functions of the constants of the motion already computed.
We conclude that the polynomial symmetry algebra generated by the 5 basic generators and their 3 commutators closes rationally, but not polynomially. 4 The classical 3D extended Kepler-Coulomb system with 4-parameter potential Now we consider the Hamiltonian , .
Here L 2 , L 3 are constants of the motion, in involution. They determine additive separation in the variables r, θ, φ. Applying our usual construction to get two independent constants of the motion we note that q 1 = r, q 2 = θ 1 , q 3 = θ 2 and to obtain functions A j , B j , j = 1, 2 such that Here, k 1 = p 1 /q 1 , k 2 = p 2 /q 2 where p 1 , q 1 are relatively prime positive integers and p 2 , q 2 are relatively prime positive integers. From our general theory, are two constants of the motion such that the full set of five constants of the motion is functionally independent. We have Here, e q 1 A 1 −2p 1 B 1 and e −q 1 A 1 +p 1 B 1 are constants of the motion, where The identity e q 1 A 1 −p 1 B 1 e −q 1 A 1 +p 1 B 1 = 1 can be expressed as where P 1 is a polynomial in H, L 2 and L 3 .
Similarly, e p 1 q 2 A 2 −p 2 q 1 B 2 , e −p 1 q 2 A 2 +p 2 q 1 B 2 are constants of the motion, where The identity e p 1 q 2 A 2 −p 2 q 1 B 2 e −p 1 q 2 A 2 +p 2 q 1 B 2 = 1 becomes where P 2 is a polynomial in H, L 2 and L 3 . We define basic raising and lowering symmetries Further, we restrict to the case where each of p 1 , q 1 , p 2 , q 2 is an odd integer. The other cases are similar. Let Then we see from the explicit expressions for the symmetries that J 1 , J 2 , K 1 , K 2 are constants of the motion, polynomial in the momenta. Moreover, the identities hold. Note that the definitions of the K-symmetries differ from those for the 3-parameter potential.
The following relations are straightforward to derive from the definition of the Poisson bracket: From these results, we find Further, Similarly Commutators relating the J and K symmetries are somewhat more complicated to compute. We have Now we can compute the nonpolynomial constant of the motion {J + , K + }: where the last term in braces vanishes identically. We conclude that Once we have {J + , K + } explicitly, we can obtain the remaining Poisson relations between the J and K symmetries with little additional work. We use the fact that J + J − = P 1 and K + K − = P 2 . Then we have We can write this relation in the more compact form Similarly, we have Thus, In summary: These relations prove closure of the symmetry algebra in the space of functions polynomial in J ± , K ± , rational in L 2 , L 3 , H and at most linear in √ L 2 , √ L 3 .

Structure relations for polynomial symmetries of the 4-parameter potential
Note that Thus we have Similarly, Straightforward computations yield

Minimal order generators
The generators for the polynomial symmetry algebra that we have produced so far are not of minimal order. Here L 1 , L 2 , L 3 are of order 2 and the orders of J 1 , K 1 are one less than the orders of J 2 , K 2 , respectively. We will construct symmetries J 0 , K 0 of order one less than J 1 , K 1 , respectively. (In the standard case k 1 = k 2 = 1 it is easy to see that J 1 is of order 5 and J 2 is of order 6, whereas K 1 is of order 3 and K 2 is of order 4. Then J 0 , K 0 will be of orders 4 and 2, respectively, which we know corresponds to the minimal generators of the symmetry algebra in this case, [21,28].) The symmetry J 2 is a polynomial in L 2 with constant term itself a constant of the motion. Thus is a polynomial symmetry of order two less than J 2 . We have the identity The same construction works for K 2 . It is a polynomial in L 3 , with constant term is a polynomial symmetry of order two less than K 2 and we have the identity Further, Now we choose L 1 , L 2 , L 3 , J 0 , K 0 as the generators of our algebra. We define the basic nonzero commutators as Then we have where the last term on the right is a polynomial in the generators H, L 2 , L 3 . Further, which again can be verified to be a polynomial in the generators. Note that this symmetry algebra cannot close polynomially in the usual sense. If it did close then the product R 1 R 2 would be expressible as a polynomial in the generators. The preceding two equations show that R 2 1 R 2 2 is so expressible, but that the resulting polynomial is not a perfect square. Thus the only possibility to obtain closure is to add new generators to the algebra, necessarily functionally dependent on the original set. We see that where A and B are polynomial in the generators, so a polynomial in the generators, Continuing in this way, it is straightforward to show that L 1 , L 2 , L 3 , K 0 , J 0 generate a symmetry algebra that closes rationally. In particular, each of the commutators R 1 , R 2 , R 3 satisfies an explicit polynomial equation in the generators.

Stäckel equivalence of Kepler-Coulomb and caged isotropic oscillator systems
Consider the Hamiltonian for the caged isotropic oscillator where , .
Here L 2 , L 3 are constants of the motion, in involution. They determine additive separation in the spherical coordinates R, φ 1 , φ 2 . Also, j 1 , j 2 are nonzero rational numbers. If j 1 = j 2 = 1, then in terms of Cartesian coordinates we have H = p 2 x + p 2 y + p 2 z + α R 2 + β /x 2 + γ /y 2 + δ /z 2 . Note that system (5.9) can be considered as the 3-variable analog of the TTW system [23,24]. (Note, however, that this is a flat space system only if j 1 = 1.) Now consider the Hamilton-Jacobi equation H = E and take the Stäckel transform that corresponds to dividing by R 2 . Then, making the change of variables r = R 2 , 2φ 1 = θ 1 , 2φ 2 = θ 2 , we obtain the new Hamilton-Jacobi equation In other words, we obtain the extended Kepler-Coulomb system. Since the Stäckel transform preserves the structure of the symmetry algebra of a superintegrable system [3,13,19,20], all of our structure results apply to the caged isotropic oscillator. Note, however, that the standard case k 1 = k 2 = 1 for Kepler-Coulomb corresponds to j 1 = j 2 = 2 for the oscillator. Further, only for the cases k 1 = 1 and j 1 = 1 is the manifold flat. The similar analysis in two dimensions [19] is always restricted to flat space, but here the manifolds depend on k 1 and j 1 .

5.2
The special case k 1 = k 2 = 1 In the case k 1 = k 2 = 1 we are in Euclidean space and our system has additional symmetry. In terms of Cartesian coordinates x = r sin θ 1 cos θ 2 , y = r sin θ 1 sin θ 2 , z = r cos θ 1 , the Hamiltonian is Note that any permutation of the ordered pairs (x, β), (y, γ), (z, δ) leaves the Hamiltonian unchanged. This leads to additional structure in the symmetry algebra. The basic symmetries are Note that the permutation symmetry of the Hamiltonian shows that I xz , I yz are also constants of the motion, and that where If δ = 0 then M 3 is the analog of the 3rd component of the Laplace vector and is itself a constant of the motion. The symmetries H, L 2 , L 3 , J 0 , K 0 form a generating (rational) basis for the constants of the motion. Under the transposition (x, β) ↔ (z, δ) this basis is mapped to an alternate basis H, All of the identities in Section 4.1 hold for the primed symmetries. It is easy to see that the K symmetries are simple polynomials in the L, K symmetries already constructed, e.g., However, the J symmetries are new. In particular, Note that the transposition (y, γ) ↔ (z, δ) does not lead to anything new. Indeed, under the symmetry we would obtain a constant of the motion but it is straightforward to check that so that the new constant depends linearly on the previous constants. For further use, note that under the transposition symmetry (x, β) ↔ (y, γ) the constants of the motion L 2 , L 3 , J 0 , J 1 are invariant, whereas K 0 , K 1 change sign. The action on J 0 is more complicated. Under the transposition J 0 and J 0 switch places. Thus from expression (5.11) we see that In the paper [21], Tanoudis and Daskaloyannis show that the quantum symmetry algebra generated by the 6 functionally dependent symmetries H, L 2 , L 3 , J 0 , K 0 and J 0 closes polynomially, in the sense that all double commutators of the generators are again expressible as polynomials in the generators, very strong evidence that the classical analog also closes polynomially. (However, they did not address the issue of determining the functional indenpendence explicitly.) Keys to understanding the polynomial closure are the rational structure equations (5.3)-(5.8) and the terms J 1 K 1 and If J 1 K 1 can be expressed as a polynomial in the generators then, with this result substituted in the rational structure equations, we will obtain polynomial structure equations. The requirement that the rational structure equations become polynomial is a strong restriction on the polynomial Let's focus on equation (5.8). Note that the left hand side of this equation is of order 6 in the momenta and J 1 K 1 is of order 8. What can we say about the polynomial P in order that its substitution into (5.8) turns the right hand side into a polynomial structure equation?
First note that the polynomial is not determined uniquely by this requirement. Indeed, if P 1 , P 2 are two solutions their difference is of the form SQ where S is a polynomial in the generators of order ≤ 4 in the momenta. Similarly, we can add such a SQ to any solution and get another solution. A straightforward computation using the polynomial and degree conditions alone yields the result where S = c 1 J 0 + c 2 J 0 + S 1 , c 1 and c 2 = 0 are parameters, and S 1 (H, L 2 , L 3 , K 0 ) is a polynomial of order ≤ 2 in its arguments. Under the transposition (x, β) ↔ (y, γ) the left hand sides of equations (5.12) and (5.13) change sign. Since Q is invariant, S must change sign. With all of these hints it is fairly easy to compute the exact result. It is in agreement with all our conditions. To determine the functional relationship between the 6 generators we can use the identity F ≡ J 2 1 K 2 1 − (J 1 K 1 ) 2 = 0, where the first term is obtained from (5.1) and (5.2) and the second term from (5.12). The result is a polynomial in the momenta of order 16, and of the simple form where A j = A j (H, L 2 , L 3 , K 0 ). Since Q factors out of the equation, the basic identity is of order 12: The leading coefficient is A 1 = −4Q, and By substituting (5.13) into expressions (4.1)-(4.4) and (5.4)-(5.7) for k 1 = k 2 = 1, and simplifying, we can recast each of these relations into polynomial structure equations in the generators. Further by squaring (5.3), substituting the structure equations for R 2 1 , R 2 2 and R 1 R 2 into the result and simplifying, we obtain the polynomial structure equations for R 2 3 . Multiplying (5.3) by R 1 substituting the structure equations for R 2 1 and R 1 R 2 into the result and simplifying, we obtain the polynomial structure equations for R 1 R 3 . Similarly, multiplying (5.3) by R 2 we can obtain the polynomial structure equations for R 2 R 3 . This entire construction carries over easily for the primed symmetries and their basic commutators (5.10). The only gap remaining to prove polynomial closure is consideration of the basic commutator R 0 = {J 0 , J 0 } and double commutators involving J 0 , J 0 simultaneously. An important observation here is that under any variable-parameter transposition R 0 changes sign and R 2 0 is invariant. Using this observation, we have verified that R 2 0 = 65536H 4 I xy I xz I yz − βI yz (I xy + I xz ) − γI xz (I xy + I yz ) − δI xy (I xz + I yz ) − βI 2 yz − γI 2 xz − δI 2 xy + (β(β + 3γ + 3δ) + γδ)I yz + (γ(γ + 3δ + 3β) + δβ)I xz + (δ(δ + 3β + 3γ) + βγ)I xy − 2(βγ 2 + β 2 γ + βδ 2 + β 2 δ + γδ 2 + γ 2 δ) .
A simple consequence is {L 2 , R 0 } = 0. The expression K 2 1 R 2 0 is a perfect square in the 6 generators and gives where the sign is determined by comparing the highest order terms. The expression J 2 1 R 2 0 is not a perfect square, but we can make it so by adding an appropriate, uniquely determined, multiple of the functional relation: (J 1 R 0 ) 2 = J 2 1 R 2 0 − 64 2 H 4 F . We then obtain where the sign is determined by comparing highest order terms. Note that The problem of computing the double commutator {J 0 , R 0 } is greatly simplified by noting that it changes sign under the transposition symmetry (x, β) ↔ (y, γ). We find {J 0 , R 0 } = 512H 2 (J 0 I yz − J 0 I xz ) + δ(J 0 − J 0 ) − γJ 0 + βJ 0 + 2α 2 (I xz − I yz ) − 2α 2 (β − γ) .
All other commutators and products can be obtained from the preceding results by use of the discrete transposition symmetries. Thus the symmetry algebra closes polynomially, and the 6 generators obey a functional identity of order 12.

Final comments
We have studied families of Hamiltonians that generalize the 3-and 4-parameter extended Kepler-Coulomb systems, found explicit generators for the symmetries that show these systems to be superintegrable, and worked out the explicit structure equations for the symmetry algebras determined by taking repeated Poisson brackets of the generators. We found it amazing that the structures could be computed exactly! This analysis strongly suggests that for higher order superintegrable systems in n ≥ 3 dimensions, polynomial closure of the symmetry algebra is relatively rare and dependent on additional discrete symmetry. Rational closure seems to be common. The structure analysis shows the fundamental importance of raising and lowering symmetries for these systems, [5,18]. These are nonpolynomial constants of the motion. However, all of the polynomial constants can be formed from them. We have no proof that the generators found by us are of minimal order in all cases. However, it is clear that all such generators must be expressible in terms of the basic rational raising and lowering symmetries.
An obvious issue is that of quantum analogs of the classical constructions. How is the problem to be quantized? What are the operator analogs of raising and lowering symmetries and of rational closure? We will give results on this in a second paper.
In recent papers [16,17] a new and very interesting approach to classical superintegrability has been developed, based on the Galois theory for differential equations. It remains to be understood how this approach relates to the techniques in the present paper.