Doubly Exotic $N$th-order Order Superintegrable Classical Systems Separating in Cartesian Coordinates

Superintegrable classical Hamiltonian systems in two-dimensional Euclidean space $E_2$ are explored. The study is restricted to Hamiltonians allowing separation of variables $V(x,y)=V_1(x)+V_2(y)$ in Cartesian coordinates. In particular, the Hamiltonian $\mathcal H$ admits a polynomial integral of order $N>2$. Only doubly exotic potentials are considered. These are potentials where none of their separated parts obey any linear ordinary differential equation. An improved procedure to calculate these higher-order superintegrable systems is described in detail. The two basic building blocks of the formalism are non-linear compatibility conditions and the algebra of the integrals of motion. The case $N=5$, where doubly exotic confining potentials appear for the first time, is completely solved to illustrate the present approach. The general case $N>2$ and a formulation of inverse problem in superintegrability are briefly discussed as well.


Introduction
For a classical Hamiltonian system with n degrees of freedom, the existence of n integrals of motion in involution is required to make it integrable in the Liouville sense. These integrals must be well-defined functions in the phase space and functionally independent. On the other hand, a superintegrable system possesses k additional integrals of motion being k = n − 1 the maximum possible number. The concept of superintegrability can be defined both in classical and quantum mechanics and it has been studied extensively for a very long period of time.
The outcome of such a long period of research activity has far reaching consequences both in mathematical and physical points of view. There exist several exhaustive review articles in literature which describe the history and current status of this topic [34,51].
Starting with a spherically symmetric standard Hamiltonian (i.e., the potential V = V (r) being velocity and spin independent), there exist only two superintegrable systems, namely the Kepler-Coulomb and the harmonic oscillator. Actually, these two potentials are exactly the ones which appear in the celebrated Bertrand's theorem [3,21]. Superintegrability of the Kepler-Coulomb problem is due to the existence of the conserved Laplace-Runge-Lenz vector [21,40,41,60] whilst in the case of the harmonic oscillator is a consequence of the existence of the quadrupole Jauch-Fradkin tensor [19,30]. of motion. For doubly exotic potentials this linear compatibility condition is satisfied trivially, it is identically zero, and the potentials satisfy non-linear equations. These classes of potentials, appearing in classical and quantum superintegrable systems have been studied both in Cartesian and polar coordinates [1,14,15,16,47].
The aim of this work is to establish in detail general properties of N th-order superintegrable classical systems that allow separation of variables in Cartesian coordinates. It can be considered as the classical counterpart of the general study on quantum superintegrable systems treated in [13]. However, unlike the latter, in this work we also study the algebra of the integrals of motion and provide exhaustive results for the case N = 5. In particular, it contains the classical analogues of all the quantum doubly exotic potentials obtained in [1] explicitly. We emphasize that for doubly exotic potentials, unlike the doubly standard ones, the limit from the quantum to the corresponding classical system (i.e., → 0) is singular for all the cases studied in the present work. Thus, the corresponding quantum and classical solutions are not connected at all. The Painlevé property characterizing the relevant determining equations in the quantum systems is completely lost in the classical case. In addition, a formulation of inverse problem in superintegrability is briefly discussed as well.
In the present article we focus on 2D classical Hamiltonian systems that are separable in Cartesian variables (x, y) and they also admit an extra polynomial integral of order N > 2. The generic Hamiltonian is given by where p i , i = 1, 2, are the canonical momenta conjugate to x and y, respectively. It describes a two-dimensional particle with unit mass m = 1 moving in the potential V (x, y) = V 1 (x)+V 2 (y). Thus, the phase space is four-dimensional. These systems are trivially second-order integrable because in addition to the Hamiltonian (1.1) they admit, for any V 1 (x) and V 2 (y), another 2nd-order symmetry of the form which Poisson commutes (i.e., {H, X } PB = 0) with the Hamiltonian (1.1). The existence of an N th-order third integral Y, makes the system N th-order superintegrable (more integrals of motion than degrees of freedom). In this case, the system is maximally superintegrable. Notice that H (1.1) is S 2 -invariant under the permutation x ⇔ y whilst the integral X is anti-invariant. From a physical point of view we are looking for 2D potentials V (x, y) = V 1 (x) + V 2 (y) for which all the bounded trajectories are closed and periodic. It is worth mentioning that (1.1) can also be interpreted as the Hamiltonian of the relative motion of a two-body problem on the plane with translational invariance. In this case, (x, y) are nothing but the Cartesian coordinates of the relative vector r = r 1 − r 2 ≡ (x, y) between the two bodies.
The outline of the paper is as follows. In Section 2, for an arbitrary potential V (x, y) not necessarily separable in a coordinate system we revisited the so called determining equations governing the existence of a general N th-order polynomial integral of motion Y N . In particular, the dominant N th-order terms in Y N lie in the enveloping algebra of the Euclidean Lie algebra e (2). From the next leading terms in Y N , a linear compatibility condition (LCC) can be obtained for the potential V only. The case of a separable potential in Cartesian coordinates is then analyzed in Sections 3-4, where we show and describe a well of determining equations and derive the first non-linear compatibility condition for the potential alone. The general form of the potentials is determined by solving these compatibility conditions. Afterwards, the surviving determining equations become linear and can be solved. In Section 5, based on the LCC, we introduce the doubly exotic potentials. A general formula for the corresponding integral Y N 4İ. Yurduşen, A.M. Escobar-Ruiz and I. Palma y Meza Montoya x Figure 1. The Hamiltonian (1.1) describes a particle with unit mass m = 1 moving in a two-dimensional is given. In Section 6 we discuss the role of the algebra of the integrals of motion in the search of superintegrable potentials, and a formulation of inverse problem in superintegrability is commented. Section 7 is devoted to the known examples with N = 3, 4. Finally, in Sections 8 and 9 we consider the case N = 5 and derive in detail all possible doubly exotic potentials. For conclusions see Section 10.
2 Superintegrability: existence of an N th-order polynomial integral In the present article we are considering Hamiltonian systems separable in Cartesian coordinates and hence they are second order integrable by construction. To further search for the superintegrability, we need to give the conditions for the existence of an additional integral of motion, which is a polynomial of order N > 2 in variables p 1 , p 2 . Although the general ideas for the existence of a N th-order integral is given in [29,57], here we would like to summarize those results for the sake of completeness.

General form
The most general form of an N th-order polynomial integral Y N is given by see [29,57], where f j,2 = f j,2 (x, y, V ) are assumed to be real functions which depend on the coordinates x and y and the potential V (x, y). The integral Y N (2.1) can be conveniently rewritten as follows Y N = W N + lower order terms, (2.2) where the leading term W N in (2.2) real parameters and L z = xp 2 − yp 1 is the z-component of angular momentum. If the quantity Y N Poisson commutes with the Hamiltonian (1.1) then the system becomes N th-order (N > 1) superintegrable. In fact, for 2D systems it would correspond to maximal superintegrability.

The determining equations
The Poisson bracket of Y N (2.1) with the Hamiltonian H (1.1) gives a polynomial, in p 1 and p 2 , of degree (N + 1). Explicitly, we have where the coefficients M n 1 ,n 2 = M n 1 ,n 2 (x, y; f j,2 , V, N ) depend on the variables x, y, the functions f j,2 appearing in the integral Y N , the potential V (x, y) we are looking for, and they also carry an N -dependence. Superintegrability requires M n 1 ,n 2 = 0, n 1 + n 2 = 0, 1, 2, . . . , (N + 1), . For an arbitrary potential V (x, y) not necessarily separable, the system (2.4) is equivalent to the following set of determining equations (DE): ( = 0, 1, 2, . . . , N 2 , j = 0, 1, 2, . . . , (N − 2 ). In (2.5), the real functions f j, ≡ 0 identically for < 0 and j < 0 as well as for j > N − 2 (further details can be found in [29,57]). The DE correspond to the vanishing of all the coefficients, in the Poisson bracket {H, Y N } PB , multiplying the momentum terms of order n 1 + n 2 = k = N + 1, N − 1, N − 3, . . . , (N + 1 − 2 ). In particular, for odd N the coefficient multiplying the zero order term is simply f 1,N −1 V 1 + f 0,N −1 V 2 = 0, obtained from (2.5) by making the replacement → + 1. The DE govern the existence of the integral Y N . In general, the system (2.5) is overdetermined. If the potential V (x, y) is not known a priori, then it must be calculated from the compatibility conditions of the DE.
The structure of the DE (2.5) can be summarized as follows: • The set of DE (2.5) can be seen as a well of recursive equations. The coefficients f j,2 in Y N depend on the preceding f j,2k , 0 ≤ k < .
• The bottom level of equations (2.5) corresponds to = 0. The associated DE do not depend on V , thus, allowing exact solvability. Indeed, they define the coefficient-functions f j,0 , j = 0, 1, 2, . . . , N . The explicit expression for f j,0 is given by see [29,57]. Accordingly, the leading part (2.3) of Y N is a polynomial of order N in the enveloping algebra of the Euclidean Lie algebra e(2) with basis {p 1 , p 2 , L z }.
• The 2nd level of DE (2.5) occurs at = 1. They provide a linear compatibility condition (LCC) for the potential V only. For arbitrary potential, this linear PDE can be written in the compact form [29,57] N −1 This above equation is a necessary but not sufficient condition for {H, Y N } = 0. Also, in the quantum case the LCC remains identical to (2.6). However, the corresponding DE do acquire -dependent terms.
• Beginning from = 2, the DE (2.5) will lead to nonlinear compatibility conditions (NLCC) for the potential V alone. We should remind here that in the quantum case these NLCC, unlike the LCC, do depend non-trivially on . Hence, the classical and quantum cases can greatly differ, and it requires to treat them separately.

The linear compatibility condition
In the case of a separable potential the LCC (2.6) leads to the ordinary differential equations For superintegrability, {H, Y N } = 0, these two linear equations (3.1) and (3.2) must be simultaneously satisfied. Similarly, for V 2 (y) there exist two ODEs which can be obtained from (3.1) and (3.2) using the symmetry x ↔ y, respectively.

The first nonlinear compatibility condition
In the case of an arbitrary odd N ≥ 3 polynomial integral of motion Y N , following the derivation presented in [13] we describe the procedure to construct the first NLCC in detail. In general, this equation obtained from the DE (2.5) with = 2 provides the form of the doubly exotic potentials, see below.
As a first step, one solve the DE (2.5) with = 1. These equations define all the coefficientfunctions f j,2 appearing in the integral (2.1).
Secondly, from the DE (2.5) with = 2 we compute the (N − 3) functions f j,4 except those with j = N −5 2 and j = N −3 2 . Eventually, we arrive at the equations here theF 's, by construction, are real functions that solely depend on the potential V (and its derivatives). Finally, from (4.1) it follows the equation which gives the aforementioned NLCC for the potential V . In the case of arbitrary even N ≥ 4, the steps are quite similar, see details in [13]. From (2.5), it follows that more NLCC occur for each value of = 3, 4, . . . , N 2 . Nevertheless, these additional equations will simply restrict the general solution of the potential V found from the previous NLCC with = 2.
Therefore, for a separable potential V = V 1 (x) + V 2 (y) the set of DE with = 0 are given by a system of ODEs which do not depend on V and they specify the coefficient-functions f j,0 (j = 0, 1, 2, . . . , N ) appearing in (2.1). Then, the next level of DE with = 1 provide a LCC for the potential alone and they also determine the functions f j,2 (j = 0, 1, 2, . . . , N − 2). At all further levels 2 the DE and their compatibility conditions are nonlinear ODEs for V . These compatibility conditions are instrumental to specify the general form of the potential V .

Doubly exotic potentials
Hereafter, we will restrict ourselves to the case of doubly exotic potentials. These potentials satisfy the LCC (2.6) trivially. In particular, the two linear ODEs (3.1) and (3.2) vanish identically for any V 1 (x). Hence, this LCC does not impose any constraint for V 1 (x) nor for V 2 (y). This situation occurs when the number of coefficients A N −m−n,m,n that figure in the LCC is less that those appearing in the integral Y N (2.1). In this case, we simply put equal to zero the coefficients A N −m−n,m,n in the LCC, thus, it vanishes identically, but still the integral Y N is of order N . In general, based on the LCC (2.6) one can classify the N th-order superintegrable systems into three major classes: doubly exotic potentials, singly exotic potentials and standard potentials (see [13]). This general classification is summarized in Table 1. Table 1. Classification of N th-order superintegrable classical systems (N > 2) separating in Cartesian coordinates. For a fixed value of N , there exist three generic types of potentials: doubly standard, doubly exotic and singly exotic potentials.

Potential
Doubly standart Doubly exotic Singly exotic  [24] are doubly exotic potentials for n 1 , n 2 > 1 whilst Cases 4-5 at n > 1 can not be doubly standard ones. Moreover, the aforementioned Cases 1-3 are nothing but particular solutions of the present direct approach. It is worth mentioning that we solely consider potentials where neither the x-part V 1 (x) nor the y-part V 2 (y) are constant functions.

Integral Y N for doubly exotic potentials
In the present work we will focus on doubly exotic potentials. In this case, the corresponding N th-order terms of the integral Y N (2.1) are given by Let us give the most general (leading) term W N of the integral Y N,doubly exotic for N = 3, 4, 5 explicitly 6 ODEs versus algebraic equations. Algebras of integrals of motion In the search of N th-order superintegrable potentials one faces the problem of solving an overdetermined system of ODEs where some of them are non-linear. Moreover, the number of involved equations increases with N . Therefore, the direct approach of solving all the DE (2.5) is far from being an efficient method. In order to simplify it, in the present consideration we propose to combine two basic elements, namely the non-linear compatibility conditions and the use of the algebra of the integrals (see below). As a result, in some cases the ODEs are reduced to pure algebraic equations. From X and Y N , we introduce the quantity which is a polynomial function in p 1 and p 2 of degree (N +1). If Y N is an integral of motion, then by construction C (6.1) is also conserved. The closure of the algebra generated by the integrals of motion (H, X , Y N , C) is guaranteed by the property of maximal superintegrability. The main question we aim to explore is the appearance and utility of a closed polynomial algebra.
It is important to mention that the study of the algebraic structure of the integrals of motion has been proven to be fruitful in the classification of higher-order superintegrable classical and quantum systems [8,46]. Also, the explicit results obtained in Section 9 suggests to explore the inverse problem, namely we take two polynomial functions A and B in momentum variables (p 1 , p 2 ) and construct the new object C = {A, B} PB . Now, let us assume that the algebra generated by (A, B, C) is a closed polynomial algebra with polynomial coefficients in H. The question is under what conditions these closure relations imply that A and B are integrals, i.e., they Poisson commute with H? 7 Lowest order cases N = 3 and N = 4: doubly exotic potentials The general integral (2.1) at N = 3 is given by Thus, (7.1) reduces to The next set of DE is obtained by setting = 1 in (2.5). They correspond to the vanishing of all the coefficients, in the Poisson bracket {H, Y 3 } PB , multiplying the (next-to-leading) momentum terms of order 2. These DE take the form The compatibility condition of the above system (7.3) does not provide further information on the potentials functions. However, the first and third equations can be solved immediately, they define the functions f 0,2 and f 1,2 where u 1 (x) and u 2 (y) are arbitrary functions of x and y, respectively. Substituting (7.4) into the second equation in (7.3) we obtain the equation u 1 + u 2 = 0. Therefore, here α 1 , α 2 are constants of integrations whilst β is a separation constant. Finally, the last determining equation corresponds to the vanishing of the coefficient, in the Poisson bracket {H, Y 3 } PB , of order zero in momentum variables. Explicitly, it takes the form Non trivial solutions of (7.5) correspond to separation of variables, namely β = 0. In this case, (7.5) leads to the following uncoupled equations being λ = 0 (otherwise the functions V 1,2 are just constants) the corresponding separation constant. Eventually, we arrive to the solutions In general, the Poisson bracket C = {Y 3 , X } PB between Y N (2.1) with N = 3 and X (1.2) is a polynomial in the variables p 1 and p 2 of degree four. However, for the potential functions (7.6), we obtain C ∝ λ. Hence, in this case the algebra of the three integrals of motion (C, Y 3 , X ) takes the form It is worth mentioning that the family of superintegrable potentials with has been analyzed in [28] by means of Heisenberg-type higher order symmetries. For this family all three conserved quantities (H, X , Y N ) admit separation of variables in Cartesian coordinates. However, such an approach does not allow us to obtain all the doubly exotic potentials with non-separable integrals Y N .

Case N = 4
In this case N = 4, for a doubly exotic potential the most general expression of the fourth-order integral Y 4 reads where A 040 , A 004 and A 022 are real constants. It can immediately be rewritten as follows Now, without losing generality, one can always add to (7.7) any arbitrary function of the second order trivial integrals H and X . This implies that no bona fide doubly exotic potentials, with a non-trivial fourth order integral, exist.

Case N = 5: doubly exotic potentials
We can write the most general 5th-order polynomial integral Y 5 in the form

Determining equations
Putting = 0 in (2.5) corresponds to the vanishing of all the coefficients, in the Poisson bracket {H, Y 5 } PB , multiplying the highest momentum terms of order 6. They can be solved directly to give the functions f j,0 where the condition that the LCC (2.6) is satisfied trivially was imposed, namely we consider doubly exotic potentials only. It implies that the existence or non-existence of fifth order doublyexotic potentials is governed by 5 parameters A 050 , A 005 , A 032 , A 023 , A 122 only. The next set of DE are obtained by setting = 1: Now, the three DE (2.5) with = 2 are given by Next, following the discussion of Section 4 we obtain from (8.2) the functions f 3,2 , f 2,2 , f 1,2 , f 0,2 in terms of V (see below). Afterwards, the r.h.s. in (8.3) would depend (non-linearly) on V and its derivatives alone. Consequently, (8.3) leads to the first NLCC in the form Finally, the last determining equation with = 2 reads

The (first) NLCC
The DE with = 1 (8.2) define the four functions f 0,2 , f 1,2 , f 2,2 and f 3,2 appearing in the integral Y 5 (8.1) in front of the cubic terms (p i 1 p j 2 with i + j = 3). Explicitly Next, substituting (8.3) and (8.5) into (8.4) we obtain the following non-linear compatibility condition (NLCC) where α 1 , β s, ν s and σ s are constants to be determined. We have the freedom to replace T 1(2) by T 1(2) + c for some real constant c to simplify the expressions. Also we can shift the variables x or y. Notice that the constants A 050 and A 005 do not appear in (8.6).
In terms of the parameters A 5−m−n,m,n that define the existence or non-existence of the integral Y 5 , we identify two cases for which the above NLCC (8.6) admits separation of variables in Cartesian coordinates, namely: with A 050 and A 005 arbitrary. These two cases are S 2 -invariant under the permutation x ⇔ y (thus, p 1 ⇔ p 2 ). Let us recall that the Hamiltonian H and the integral X are S 2 -invariant and S 2 -antiinvariant, respectively. Moreover, if (iii) A 122 = A 023 = A 032 = 0, with A 2 050 + A 2 005 = 0, the NLCC degenerates into a linear equation which must be identically zero for doubly exotic potentials. In such a case the NLCC does not provide any information on the potential. As a result of calculations, the cases (i), (ii) and (iii) are the only generic ones that satisfy all the DE.

Superintegrable potentials
Below, adopting the notation introduced in [1] we present the full list of doubly exotic fifth-order (N = 5) superintegrable potentials: Case (i).
• The system Q 1 : A 122 = 0, A 032 = A 023 = A 050 = A 005 = 0. This system corresponds to A 122 = 1, all other parameters A ijk = 0. In this case, by solving all the DE (8.2)-(8.4) we eventually arrive to the first-order nonlinear ODEs 4 y 4 = 0, (9.1) β 4 = 0 is a real constant. The corresponding fifth-order integral of motion is given by , we built the quantity which is a polynomial function in p 1 and p 2 of sixth degree. By construction, it is an integral when (9.1) are satisfied. Now, if we demand that the three elements X , Y (Q 1 ) 5 , C generate a closed polynomial algebra we eventually arrive to a nonlinear first-order differential equation for T 1 (x) and similarly for T 2 (y). Therefore, from these equations and (9.1) we can eliminate the first-derivative T 1 , T 2 terms and obtain an algebraic equation for both T 1 (x) and T 2 (y). The solutions of such algebraic equations turn out to be the general solutions of (9.1). Explicitly, these algebraic equations take the form where δ is an arbitrary constant. In the case δ = 0, we immediately obtain the particular solutions x 3 , and T 2 (y) = β 4 y 3 , β 4 9 y 3 , which correspond to a well-known lower-order superintegrable system. The algebra generated by the integrals takes the form In the corresponding quantum system analyzed in [1], the case (i) splits into two subclasses of integrals Y 5 that solely differ in their lower order -dependable terms. Consequently, two systems called Q 1 and Q 2 occur. However, in the classical limit → 0 the two systems Q 1 and Q 2 coincide. Next, within case (ii) the classical systems Q 3 (A 023 A 050 A 005 = 0, A 122 = A 032 = 0) and Q 4 (A 023 A 005 = 0, A 122 = A 050 = A 032 = 0) are not superintegrable (like in the quantum case).
• The system Q 5 : A 023 = 0, A 122 = A 032 = A 005 = 0, A 050 arbitrary. This system corresponds to A 023 = 1 and arbitrary A 050 , all other A ijk = 0. Again, by solving all the DE (8.2)-(8.4) we arrive to the first-order nonlinear ODE for T 1 (9.4) τ = 0, β 1 and µ are real constants, whereas The corresponding highest-order integral of motion reads Clearly, the case A 032 = 1 and arbitrary A 005 (all other A ijk = 0) also leads to a superintegrable potential. It can simply be obtained by replacing A 050 → A 005 and making the permutation x ⇔ y (V 1 ⇔ V 2 ) in (9.4) and (9.5). • The system Q 6 : A 023 A 032 = 0, A 122 = A 050 = A 005 = 0. This system corresponds to A 023 A 032 = 0, all other A ijk = 0. Figure 6. A doubly exotic potential Q 6 corresponding to N = 5. It admits the fifth-order integral Y 5 (9.7). It also possesses bounded trajectories which by construction are closed and periodic. The values A 032 = −A 023 = 1 12 , σ 1 = 1 were used.
In this case, the algebra generated by the integrals C = Y and X takes the form {C, X } PB = 0, C, Y (Q 6 ) 5 PB = 2X σ 2 2 H 2 − X 2 2 , and it does not provide further information about the solutions of (9.6). However, it is easy to check that satisfy (9.6), where W = W (z) is given by the following third order polynomial equation here τ = 0 is an integration constant.
• The system Q 7 . This system is a particular case of system Q 5 . It corresponds to the situation where A 023 = 1 and all other A ijk = 0.
In this case, for the solutions of (9.8) the Poisson bracket C = Y , X PB ∝ λ. Thus, the algebra of the integrals of motion takes the form PB = 0. Again, the system Q 8 was found in [28] by means of Heisenberg-type higher order symmetries. All three conserved quantities H, X , Y admit separation of variables in Cartesian coordinates. The particular case with β 1 = θ 1 = κ 1 = α 1 = φ 1 = ω 2 = 0, thus V 1 ∝ x 1 3 and V 2 ∝ y 1 3 , was studied in [24] using action-angle variables.
• The system Q 9 . This system is a particular case of system Q 8 . It corresponds to the situation where A 050 = 1 and all other A ijk = 0.

Conclusions
We considered N th-order superintegrable classical systems in a two-dimensional Euclidean space separating in Cartesian coordinates. They are characterized by three polynomial (in momentum variables) integrals of motion (H, X , Y N ). Let us summarize the main results reported in this paper: 1. Higher-order (N > 2) superintegrable classical systems H = 1 2 p 2 1 + p 2 2 + V 1 (x) + V 2 (y), can be classified into three classes: doubly standard, singly exotic and doubly exotic potentials. This classification is based on the nature of the equation that defines the most general form of the potential functions V 1 (x) and V 2 (x). For doubly standard potentials this equation is a linear compatibility condition necessary for the existence of the N th order integral of motion Y N (in general a PDE of order N in two variables) whilst in the case of doubly exotic potentials it is given by a nonlinear compatibility condition.
2. From the equation {Y N , H} PB = 0, we show in a systematic manner how to find and successively solve a "well" of NLCC separately for V 1 (x) and V 2 (y). It was also indicated that requiring the integrals of motion (C = {X , Y N } PB , X , Y N ) to span a closed polynomial algebra may help to simplify (reduce the order) of the DE and eventually to find the explicit solutions V 1 (x) and V 2 (y).
3. All fifth-order (N = 5) superintegrable doubly exotic potentials were derived explicitly by solving the set of DE. The DE lead to first order non-linear ODEs that define the functions V 1 and V 2 , respectively. Unlike the quantum case, these equations do not have the Painlevé property. This was verified either by finding their general solutions explicitly or by applying a standard test to them [2]. Interestingly, at N = 5 doubly exotic confining potentials appear for the first time. At N = 4 no doubly exotic potentials occur at all. Finally, a direct computation for the next two cases N = 6, 7 leads us to the following conjecture: Conjecture. There exists an infinite family of N th-order superintegrable systems with an integral The three functions F q 's in (10.1) are polynomials in the variable u of degree at most (N − 1), and σ, b are real parameters as well. The equation (10.1) is in complete agreement with the limit → 0 of its quantum analogue treated in [13]. In future work, we plan to establish in detail under what conditions the closed algebra of the integrals of motion is polynomial, and how to use it as a new systematic tool to solve the determining equations in a simpler and more efficient manner.