Doubly Exotic $N$th-Order Superintegrable Classical Systems Separating

Doubly Exotic $\bm{N}$th-Order Superintegrable Classical
Systems Separating in Cartesian Coordinates

İsmet YURDUŞEN a, Adrián Mauricio ESCOBAR-RUIZ b
and Irlanda PALMA Y MEZA MONTOYA b

İ. Yurduşen, A.M. Escobar-Ruiz and I. Palma y Meza Montoya

a) Department of Mathematics, Hacettepe University, 06800 Beytepe, Ankara,
Turkey<EMAIL_ADDRESS>
b) Departamento de Física, Universidad Autónoma Metropolitana-Iztapalapa,
b) San Rafael Atlixco 186, México, CDMX, 09340 México<EMAIL_ADDRESS><EMAIL_ADDRESS>
Received December 18, 2021, in final form May 16, 2022; Published online May
27, 2022

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.

integrability in classical mechanics; higher-order superintegrability;
separation of variables; exotic potentials

70H06; 70H33; 70H50

##  1 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].

The systematic investigation of superintegrability has been initiated by Pavel
Winternitz and his collaborators in 1965 [20]. They first considered quadratic
superintegrability in Euclidean spaces and the subject has been subsequently
developed into many directions by several authors since then. For example, its
close relation with multiseparability was studied in detail in the references
[17, 18, 20, 31, 42, 49], the search for superintegrable systems in $2$\- and
$3$-dimensional spaces of constant and nonconstant curvature has been carried
out in the works [25, 26, 27, 31, 32, 33, 35, 36, 37, 50] and their
generalizations to $n$-dimensions have been analyzed in the papers [38, 39,
59].

Another important research direction in this field is to consider the
Hamiltonians with magnetic field and/or spin. Superintegrability with magnetic
field was first explored in the articles [5, 10] and much recently developed
in the articles [4, 43, 44]. The systematic investigation of integrability and
superintegrability for systems involving particles with spin was initiated in
the reference [68] and subsequently all the rotationally invariant
superintegrable systems in $E_{3}$ were classified in the articles [9, 69,
70]. On the other hand, spin dependent superintegrable systems were studied in
the works [53, 54] for matrix potentials simulating charged or neutral
fermions with non-trivial dipole moment in the presence of an electric field.

Still another interesting direction is to go beyond quadratic
superintegrability, i.e., the general theory of higher-order
superintegrability. Initial pioneering works were the articles of Drach [11,
12], where $10$ potentials allowing third-order integrals of motion were
announced. However, much later it was shown that $7$ of these potentials are
actually reducible, the third-order integral is the Poisson commutator of two
second-order integrals [58, 65]. Once again, the systematic investigation of
higher-order superintegrability, in particular the third-order one has been
initiated by Pavel Winternitz and his collaborators in the articles [22, 23,
45, 48, 55, 64]. Almost around the same time higher-order symmetry operators
were calculated for the Schrödinger operator and the determining equations for
the corresponding integrals of motion appeared in [52]. Nevertheless, it was
soon realized that the analysis became very complicated and some new ways of
approaching to the problem of higher-order superintegrability have to be
considered.

After the publication of the seminal paper “An infinite family of solvable and
integrable quantum systems on a plane” [62], the direction of the research has
been thoroughly shifted to higher-order integrability/superintegrability [56,
57, 63]. Moreover, in order to make them more easily tractable, new techniques
and methods have been implemented in the study of higher-order integrable and
superintegrable systems [6, 7].

From our point of view one of the main issues on higher-order
superintegrability is the classification of the superintegrable potentials. In
the case of 2D separable potentials in Cartesian coordinates, an $N$th-order
superintegrable system appeared for the first time in [61], where the
existence of nonlinear equations for the potential which makes the general
problem much more complicated was stressed as well. Recently in 2018, by means
of a systematic study an infinite 2-parametric family of superintegrable
potentials embracing those found in [28, 61] was presented in the paper [24]
by Grigoriev and Tsiganov. Their key element to construct polynomial integrals
of motion is the addition theorems for the action-angle variables, especially
the Chebyshev theorem applied to integrals on differential binomials (see also
[66, 67]). Such an elegant approach has the advantage that it uses the action-
angle variables which play a fundamental role in classical mechanics.

However, unlike the present direct approach some systems can be missed and the
generalization to the quantum case is not straightforward. The explicit list
of all $N$th-order 2D (polynomial) superintegrable potential separating in
Cartesian coordinates is far to be complete.

Throughout many recent research activity in the field of higher order
superintegrability, it has been clarified that three types of potentials can
occur, namely the standard, the doubly exotic and singly exotic potentials.
Standard ones are solutions of a linear compatibility condition for the
determining equations that govern the existence of a higher-order polynomial
integral 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., $\hbar\rightarrow 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

$\displaystyle{\mathcal{H}}={\mathcal{H}}_{1}(x)+{\mathcal{H}}_{2}(y)\equiv\frac{1}{2}\big{(}p_{1}^{2}+p_{2}^{2}\big{)}+V_{1}(x)+V_{2}(y),$
(1.1)

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

$\displaystyle{\mathcal{X}}={\mathcal{H}}_{1}(x)-{\mathcal{H}}_{2}(y)=\frac{1}{2}\big{(}p_{1}^{2}-p_{2}^{2}\big{)}+V_{1}(x)-V_{2}(y),$
(1.2)

which Poisson commutes (i.e., $\\{{\mathcal{H}},{\mathcal{X}}\\}_{\rm PB}=0$)
with the Hamiltonian (1.1). The existence of an $N$th-order third integral
$\mathcal{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 ${\mathcal{H}}$ (1.1) is
${\mathcal{S}}_{2}$-invariant under the permutation $x\Leftrightarrow y$
whilst the integral $\mathcal{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 ${\bf r}={\bf r}_{1}-{\bf
r}_{2}\equiv(x,y)$ between the two bodies.

$V(x,y)=V_{1}(x)+V_{2}(y)$$m$$y$$x$

Figure 1: The Hamiltonian (1.1) describes a particle with unit mass $m=1$
moving in a two-dimensional potential $V(x,y)=V_{1}(x)+V_{2}(x)$.

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 ${\mathcal{Y}}_{N}$. In
particular, the dominant $N$th-order terms in ${\mathcal{Y}}_{N}$ lie in the
enveloping algebra of the Euclidean Lie algebra $e(2)$. From the next leading
terms in ${\mathcal{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
${\mathcal{Y}}_{N}$ 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 $\bm{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.

###  2.1 General form

The most general form of an $N$th-order polynomial integral
${\mathcal{Y}}_{N}$ is given by

$\displaystyle{\mathcal{Y}}_{N}=\sum_{\ell=0}^{[\frac{N}{2}]}\sum_{j=0}^{N-2\ell}f_{j,2\ell}\,p_{1}^{j}p_{2}^{N-j-2\ell},$
(2.1)

see [29, 57], where ${f}_{j,2\ell}={f}_{j,2\ell}(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 ${\mathcal{Y}}_{N}$ (2.1) can be conveniently rewritten as
follows

$\displaystyle{\mathcal{Y}}_{N}=W_{N}+\text{lower order terms},$ (2.2)

where the leading term $W_{N}$ in (2.2)

$\displaystyle W_{N}=\sum_{0\leq m+n\leq
N}A_{N-m-n,m,n}L_{z}^{N-m-n}p_{1}^{m}p_{2}^{n},$ (2.3)

plays a fundamental role since it governs the existence or non-existence of
the integral ${\mathcal{Y}}_{N}$, here $A_{N-m-n,m,n}$ are
$\frac{(N+1)(N+2)}{2}$ real parameters and $L_{z}=xp_{2}-yp_{1}$ is the
$z$-component of angular momentum. If the quantity ${\mathcal{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.

###  2.2 The determining equations

The Poisson bracket of ${\mathcal{Y}}_{N}$ (2.1) with the Hamiltonian
$\mathcal{H}$ (1.1) gives a polynomial, in $p_{1}$ and $p_{2}$, of degree
$(N+1)$. Explicitly, we have

$\displaystyle\\{{\mathcal{H}},{\mathcal{Y}}_{N}\\}_{\rm
PB}=\sum_{n_{1}+n_{2}=0}^{N+1}M_{n_{1},n_{2}}p_{1}^{n_{1}}p_{2}^{n_{2}},$

where the coefficients
${M}_{n_{1},n_{2}}={M}_{n_{1},n_{2}}(x,y;f_{j,2\ell},V,N)$ depend on the
variables $x$, $y$, the functions $f_{j,2\ell}$ appearing in the integral
${\mathcal{Y}}_{N}$, the potential $V(x,y)$ we are looking for, and they also
carry an $N$-dependence. Superintegrability requires

$\displaystyle{M}_{n_{1},n_{2}}=0,\qquad n_{1}+n_{2}=0,1,2,\dots,(N+1),$ (2.4)

($\\{{\mathcal{H}},{\mathcal{Y}}_{N}\\}=0$). For an arbitrary potential
$V(x,y)$ not necessarily separable, the system (2.4) is equivalent to the
following set of determining equations (DE):

$\displaystyle(\partial_{x}{f}_{j-1,2\ell}+\partial_{y}{f}_{j,2\ell})-\big{[}(j+1){f}_{j+1,2\ell-2}\big{]}\partial_{x}V-\big{[}(N-2\ell+2-j){f}_{j,2\ell-2}\big{]}\partial_{y}V=0,$
(2.5)

$\ell=0,1,2,\dots,\big{[}\frac{N}{2}\big{]}$, $j=0,1,2,\dots,(N-2\ell)$. In
(2.5), the real functions ${f}_{j,\ell}\equiv 0$ identically for $\ell<0$ and
$j<0$ as well as for $j>N-2\ell$ (further details can be found in [29, 57]).
The DE correspond to the vanishing of all the coefficients, in the Poisson
bracket $\\{{\mathcal{H}},{\mathcal{Y}}_{N}\\}_{\rm PB}$, multiplying the
momentum terms of order $n_{1}+n_{2}=k=N+1,N-1,N-3,\dots,(N+1-2\ell)$. In
particular, for odd $N$ the coefficient multiplying the zero order term is
simply $f_{1,N-1}V_{1}^{\prime}+f_{0,N-1}V_{2}^{\prime}=0,$ obtained from
(2.5) by making the replacement $\ell\rightarrow\ell+1$. The DE govern the
existence of the integral ${\mathcal{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\ell}$ in ${\mathcal{Y}}_{N}$ depend on the preceding
$f_{j,2k}$, $0\leq k<\ell$.

  * •

The bottom level of equations (2.5) corresponds to $\ell=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,\dots,N$. The explicit
expression for $f_{j,0}$ is given by

$\displaystyle
f_{j,0}=\sum_{n=0}^{N-j}\sum_{m=0}^{j}\binom{N-n-m}{j-m}A_{N-n-m,m,n}x^{N-j-n}(-y)^{j-m},$

see [29, 57]. Accordingly, the leading part (2.3) of ${\mathcal{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 $\ell=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]

$\displaystyle\sum_{j=0}^{N-1}{(-1)}^{j}\partial_{x}^{N-1-j}\partial_{y}^{j}\big{[}(j+1)f_{j+1,0}\partial_{x}V+(N-j)f_{j,0}\partial_{y}V\big{]}=0.$
(2.6)

This above equation is a necessary but not sufficient condition for
$\\{{\mathcal{H}},{\mathcal{Y}}_{N}\\}=0$. Also, in the quantum case the LCC
remains identical to (2.6). However, the corresponding DE do acquire
$\hbar$-dependent terms.

  * •

Beginning from $\ell=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
$\hbar$. Hence, the classical and quantum cases can greatly differ, and it
requires to treat them separately.

##  3 Superintegrable potentials separable in Cartesian coordinates

###  3.1 The linear compatibility condition

In the case of a separable potential the LCC (2.6) leads to the ordinary
differential equations for $V_{1}(x)$

$\displaystyle\sum_{j=0}^{N-1}(j+1)!\sum_{n=0}^{N-j-1}\binom{N-1-n}{j}A_{N-1-n,1,n}\bigg{(}\frac{\rm
d}{{\rm d}x}\bigg{)}^{N-j+1}\big{[}x^{N-j-n-1}V_{1}^{\prime}(x)\big{]}=0,$
(3.1)

and

$\displaystyle\sum_{j=0}^{N-1}(j\\!+1)(j\\!+1)!(-1)^{2j+1}\\!\sum_{n=0}^{N-j-1}\\!\\!\binom{N\\!-n}{j\\!+1}A_{N-n,0,n}\bigg{(}\frac{\rm
d}{{\rm
d}x}\bigg{)}^{N-j+1}\\!\\!\big{[}x^{N-j-n-1}V_{1}^{\prime}(x)\big{]}\\!=0.$
(3.2)

For superintegrability, $\\{{\mathcal{H}},{\mathcal{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\leftrightarrow y$, respectively.

##  4 The first nonlinear compatibility condition

In the case of an arbitrary odd $N\geq 3$ polynomial integral of motion
${\mathcal{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 $\ell=2$ provides the form of the doubly
exotic potentials, see below.

As a first step, one solve the DE (2.5) with $\ell=1$. These equations define
all the coefficient-functions ${f}_{j,2}$ appearing in the integral (2.1).

Secondly, from the DE (2.5) with $\ell=2$ we compute the ($N-3$) functions
${f}_{j,4}$ except those with $j=\frac{N-5}{2}$ and $j=\frac{N-3}{2}$.
Eventually, we arrive at the equations

$\displaystyle\partial_{y}{f}_{\frac{N-5}{2},4}=\tilde{F}_{\frac{N-5}{2}},\qquad\partial_{x}{f}_{\frac{N-5}{2},4}+\partial_{y}{f}_{\frac{N-3}{2},4}=\tilde{F}_{\frac{N-3}{2}},\qquad\partial_{x}{f}_{\frac{N-3}{2},4}=\tilde{F}_{\frac{N-1}{2}},$
(4.1)

here the $\tilde{F}$’s, by construction, are real functions that solely depend
on the potential $V$ (and its derivatives).

Finally, from (4.1) it follows the equation

$\displaystyle\partial^{2}_{x}\tilde{F}_{\frac{N-5}{2}}+\partial^{2}_{y}\tilde{F}_{\frac{N-1}{2}}-\partial^{2}_{x,y}\tilde{F}_{\frac{N-3}{2}}\equiv
0,$

which gives the aforementioned NLCC for the potential $V$. In the case of
arbitrary even $N\geq 4$, the steps are quite similar, see details in [13].

From (2.5), it follows that more NLCC occur for each value of
$\ell=3,4,\dots,\big{[}\frac{N}{2}\big{]}$. Nevertheless, these additional
equations will simply restrict the general solution of the potential $V$ found
from the previous NLCC with $\ell=2$.

Therefore, for a separable potential $V=V_{1}(x)+V_{2}(y)$ the set of DE with
$\ell=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,\dots,N$) appearing in
(2.1). Then, the next level of DE with $\ell=1$ provide a LCC for the
potential alone and they also determine the functions $f_{j,2}$
($j=0,1,2,\dots,N-2$). At all further levels $\ell\geqslant 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$.

##  5 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
---|---|---|---
$V=V_{1}(x)+V_{2}(y)$ | potentials | potentials | potentials
Classical | Both functions |  Both $V_{1}(x)$, $V_{2}(y)$ |  The $x$-component
superintegrable |  $V_{1}(x)$, $V_{2}(y)$ satisfy | obey a NLCC, |  $V_{1}(x)$ satisfies
systems | non-trivially | a non-linear ODE | a linear/non-linear
| the LCC, | which do not pass |  ODE whilst $y$-component
| a linear ODE | the Painlevé test. |  $V_{2}(y)$ obeys
|  | The LCC is identically zero | a non-linear/linear OD

From this point of view, Cases 1–3 of Proposition 1 presented in [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.

###  5.1 Integral $\bm{{\mathcal{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 ${\mathcal{Y}}_{N}$ (2.1)
are given by

$\displaystyle W_{N}={}$ $\displaystyle
A_{0,N,0}p_{1}^{N}+A_{0,0,N}p_{2}^{N}+A_{N-4,2,2}L_{z}^{N-4}p_{1}^{2}p_{2}^{2}$
$\displaystyle+\sum_{4<m+n<N;\,|m-n|<N-4}A_{N-m-n,m,n}L_{z}^{N-m-n}p_{1}^{m}p_{2}^{n}$
$\displaystyle+\sum_{0\leq m+n=N;\,|m-n|\leq N-4}A_{0,m,n}p_{1}^{m}p_{2}^{n}.$
(5.1)

Therefore, $W_{N}$ (5.1) carries $\frac{1}{2}\big{(}10-5N+N^{2}\big{)}$ less
constants $A_{N-m-n,m,n}$ than the generic term (2.1). For large $N$, the
number of constants $A_{N-m-n,m,n}$ in $W_{N}$ grows quadratically with $N$.
Notice that $L_{z}$ occurs in $W_{N}$ starting from $N=5$. For $N\geq 5$, it
can also contain terms like $L_{z},L_{z}^{2},L_{z}^{3},\dots,L_{z}^{N-4}$
only.

Let us give the most general (leading) term $W_{N}$ of the integral
${\mathcal{Y}}_{N,{\rm doubly\ exotic}}$ for $N=3,4,5$ explicitly

$\displaystyle W_{3}=A_{030}p_{1}^{3}+A_{003}p_{2}^{3},$ $\displaystyle
W_{4}=A_{040}p_{1}^{4}+A_{004}p_{2}^{4}+A_{022}p_{1}^{2}p_{2}^{2},$
$\displaystyle
W_{5}=A_{050}p_{1}^{5}+A_{005}p_{2}^{5}+A_{032}p_{1}^{3}p_{2}^{2}+A_{023}p_{1}^{2}p_{2}^{3}+A_{122}L_{z}p_{1}^{2}p_{2}^{2}.$

##  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 $\mathcal{X}$ and ${\mathcal{Y}}_{N}$, we introduce the quantity

$\displaystyle C\equiv\\{{\mathcal{Y}}_{N},{\mathcal{X}}\\}_{\rm PB},$ (6.1)

which is a polynomial function in $p_{1}$ and $p_{2}$ of degree ($N+1$). If
${\mathcal{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 $({\mathcal{H}},{\mathcal{X}},{\mathcal{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 ${\mathcal{A}}$ and
${\mathcal{B}}$ in momentum variables $(p_{1},p_{2})$ and construct the new
object ${\mathcal{C}}=\\{{\mathcal{A}},{\mathcal{B}}\\}_{\rm PB}.$ Now, let us
assume that the algebra generated by
$({\mathcal{A}},{\mathcal{B}},{\mathcal{C}})$ is a closed polynomial algebra
with polynomial coefficients in ${\mathcal{H}}$. The question is under what
conditions these closure relations imply that ${\mathcal{A}}$ and
${\mathcal{B}}$ are integrals, i.e., they Poisson commute with $\mathcal{H}$?

##  7 Lowest order cases $\bm{N=3}$ and $\bm{N=4}$:
doubly exotic potentials

###  7.1 Case $\bm{N=3}$

The general integral (2.1) at $N=3$ is given by

$\displaystyle{\mathcal{Y}}_{3}=f_{3,0}p_{1}^{3}+f_{0,0}p_{2}^{3}+{f}_{1,0}p_{1}p_{2}^{2}+{f}_{2,0}p_{2}p_{1}^{2}+{f}_{1,2}p_{1}+{f}_{0,2}p_{2}.$
(7.1)

The first set of DE (2.5) with $\ell=0$ corresponds to the vanishing of all
the coefficients, in the Poisson bracket
$\\{{\mathcal{H}},{\mathcal{Y}}_{3}\\}_{\rm PB}$, multiplying the highest
momentum terms of order $4$. They can be solved directly to give the functions
$f_{j,0}$. For doubly exotic potentials, they read

$\displaystyle f_{30}=A_{030},\qquad f_{20}=0,\qquad f_{10}=0,\qquad
f_{00}=A_{003}.$

Thus, (7.1) reduces to

$\displaystyle{\mathcal{Y}}_{3}=A_{030}p_{1}^{3}+A_{003}p_{2}^{3}+{f}_{1,2}p_{1}+{f}_{0,2}p_{2}.$
(7.2)

The next set of DE is obtained by setting $\ell=1$ in (2.5). They correspond
to the vanishing of all the coefficients, in the Poisson bracket
$\\{{\mathcal{H}},{\mathcal{Y}}_{3}\\}_{\rm PB}$, multiplying the (next-to-
leading) momentum terms of order $2$. These DE take the form

$\displaystyle f_{1,2}^{(1,0)}=3A_{030}V_{1}^{\prime},\qquad
f_{1,2}^{(0,1)}+f_{0,2}^{(1,0)}=0,\qquad
f_{0,2}^{(0,1)}=3A_{003}V_{2}^{\prime}.$ (7.3)

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}$

$\displaystyle f_{1,2}=3A_{030}V_{1}+u_{2}(y),\qquad
f_{0,2}=3A_{003}V_{2}+u_{1}(x),$ (7.4)

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}^{\prime}+u_{2}^{\prime}=0$. Therefore,

$\displaystyle u_{1}=\alpha_{1}+\beta x,\qquad u_{2}=\alpha_{2}-\beta y,$

here $\alpha_{1}$, $\alpha_{2}$ are constants of integrations whilst $\beta$
is a separation constant. Finally, the last determining equation corresponds
to the vanishing of the coefficient, in the Poisson bracket
$\\{{\mathcal{H}},{\mathcal{Y}}_{3}\\}_{\rm PB}$, of order zero in momentum
variables. Explicitly, it takes the form

$\displaystyle
3A_{030}V_{1}V_{1}^{\prime}+3A_{003}V_{2}V_{2}^{\prime}+\alpha_{2}V_{1}^{\prime}-\beta
yV_{1}^{\prime}+\beta xV_{2}^{\prime}+\alpha_{1}V_{2}^{\prime}=0.$ (7.5)

Non trivial solutions of (7.5) correspond to separation of variables, namely
$\beta=0$. In this case, (7.5) leads to the following uncoupled equations

$\displaystyle
3A_{030}V_{1}V_{1}^{\prime}+\alpha_{2}V_{1}^{\prime}=\lambda,\qquad
3A_{003}V_{2}V_{2}^{\prime}+\alpha_{1}V_{2}^{\prime}=-\lambda,$

being $\lambda\neq 0$ (otherwise the functions $V_{1,2}$ are just constants)
the corresponding separation constant. Eventually, we arrive to the solutions

$\displaystyle V_{1}=\sqrt{\frac{2\lambda}{3A_{030}}}\sqrt{x},\qquad
V_{2}=\sqrt{\frac{-2\lambda}{3A_{003}}}\sqrt{y}.$ (7.6)

In general, the Poisson bracket $C=\\{{\mathcal{Y}}_{3},{\mathcal{X}}\\}_{\rm
PB}$ between ${\mathcal{Y}}_{N}$ (2.1) with $N=3$ and $\mathcal{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\propto\lambda$. Hence, in this
case the algebra of the three integrals of motion
($C,{\mathcal{Y}}_{3},{\mathcal{X}}$) takes the form

$\displaystyle\\{C,{\mathcal{X}}\\}_{\rm
PB}=0,\qquad\\{C,{\mathcal{Y}}_{3}\\}_{\rm PB}=0.$

$V=V_{1}(x)+V_{2}(y)$$0$$0.5$$1.0$$1.5$$2.0$$0$$0.2$$0.4$$x$$0.6$$0.8$$1.0$$0$$0.2$$0.4$$y$$0.6$$0.8$$1.0$

Figure 2: The doubly exotic potential (7.6) corresponding to $N=3$. It admits
the third-order integral ${\mathcal{Y}}_{3}$ (7.2). The values
$A_{030}=-A_{003}=\frac{2\lambda}{3}$ were used.

It is worth mentioning that the family of superintegrable potentials with

${\mathcal{Y}}_{N}=A_{0N0}p_{1}^{N}+A_{00N}p_{2}^{N}+{\rm lower\ order\
terms},$

has been analyzed in [28] by means of Heisenberg-type higher order symmetries.
For this family all three conserved quantities
$({\mathcal{H}},{\mathcal{X}},{\mathcal{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
${\mathcal{Y}}_{N}$.

###  7.2 Case $\bm{N=4}$

In this case $N=4$, for a doubly exotic potential the most general expression
of the fourth-order integral ${\mathcal{Y}}_{4}$ reads

$\displaystyle{\mathcal{Y}}_{4}=A_{040}p_{1}^{4}+A_{004}p_{2}^{4}+A_{022}p_{1}^{2}p_{2}^{2}+{\rm
lower\ order\ terms},$

where $A_{040}$, $A_{004}$ and $A_{022}$ are real constants. It can
immediately be rewritten as follows

$\displaystyle{\mathcal{Y}}_{4}=A_{040}{({\mathcal{H}}+{\mathcal{X}})}^{2}+A_{004}{({\mathcal{H}}-{\mathcal{X}})}^{2}+A_{022}{({\mathcal{H}}+{\mathcal{X}})}{({\mathcal{H}}-{\mathcal{X}})}+{\rm
lower\ order\ terms}.$ (7.7)

Now, without losing generality, one can always add to (7.7) any arbitrary
function of the second order trivial integrals ${\mathcal{H}}$ and
$\mathcal{X}$. This implies that no _bona fide_ doubly exotic potentials, with
a non-trivial fourth order integral, exist.

##  8 Case $\bm{N=5}$: doubly exotic potentials

We can write the most general $5$th-order polynomial integral
${\mathcal{Y}}_{5}$ in the form

$\displaystyle{\mathcal{Y}}_{5}=\sum_{\ell=0}^{2}\sum_{j=0}^{5-2\ell}f_{j,2\ell}\,p_{1}^{j}p_{2}^{5-j-2\ell}.$
(8.1)

###  8.1 Determining equations

Putting $\ell=0$ in (2.5) corresponds to the vanishing of all the
coefficients, in the Poisson bracket
$\\{{\mathcal{H}},{\mathcal{Y}}_{5}\\}_{\rm PB}$, multiplying the highest
momentum terms of order $6$. They can be solved directly to give the functions
$f_{j,0}$

$\displaystyle f_{50}=A_{050},\qquad f_{40}=0,\qquad f_{30}=A_{032}-yA_{122},$
$\displaystyle f_{20}=A_{023}+xA_{122},\qquad f_{10}=0,\qquad f_{00}=A_{005},$

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 $\ell=1$:

$\displaystyle f_{0,2}{}^{(0,1)}=5A_{005}V_{2}^{\prime},$ $\displaystyle
f_{1,2}{}^{(0,1)}+f_{0,2}{}^{(1,0)}=2(A_{023}+xA_{122})V_{1}^{\prime},$
$\displaystyle
f_{1,2}{}^{(1,0)}+f_{2,2}{}^{(0,1)}=3(A_{032}-yA_{122})V_{1}^{\prime}+3(A_{023}+xA_{122})V_{2}^{\prime},$
$\displaystyle
f_{3,2}{}^{(0,1)}+f_{2,2}{}^{(1,0)}=2(A_{032}-yA_{122})V_{2}^{\prime},$
$\displaystyle f_{3,2}{}^{(1,0)}=5A_{050}V_{1}^{\prime}.$ (8.2)

Now, the three DE (2.5) with $\ell=2$ are given by

$\displaystyle
f_{1,4}{}^{(1,0)}=3f_{3,2}V_{1}^{\prime}+f_{2,2}V_{2}^{\prime},$
$\displaystyle
f_{1,4}{}^{(0,1)}+f_{0,4}{}^{(1,0)}=2(f_{2,2}V_{1}^{\prime}+f_{1,2}V_{2}^{\prime}),$
$\displaystyle
f_{0,4}{}^{(0,1)}=3f_{0,2}V_{2}^{\prime}+f_{1,2}V_{1}^{\prime}.$ (8.3)

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

$\displaystyle\partial^{2}_{x}f_{0,4}^{(0,1)}+\partial^{2}_{y}f_{1,4}^{(1,0)}-\partial_{x}\partial_{y}\big{(}f_{0,4}^{(1,0)}-f_{1,4}^{(0,1)}\big{)}=0.$
(8.4)

Finally, the last determining equation with $\ell=2$ reads

$\displaystyle f_{1,4}V_{1}^{\prime}+f_{0,4}V_{2}^{\prime}=0.$

###  8.2 The (first) NLCC

The DE with $\ell=1$ (8.2) define the four functions $f_{0,2}$, $f_{1,2}$,
$f_{2,2}$ and $f_{3,2}$ appearing in the integral ${\mathcal{Y}}_{5}$ (8.1) in
front of the cubic terms ($p_{1}^{i}p_{2}^{j}$ with $i+j=3$). Explicitly

$\displaystyle
f_{0,2}=2xA_{122}T_{1}^{\prime}(x)+2A_{023}T_{1}^{\prime}(x)+A_{122}T_{1}(x)+5A_{005}T_{2}^{\prime}(y)+\alpha_{1}-\beta_{4}x^{3}+\sigma_{3}x^{2}+\alpha_{2}x,$
$\displaystyle
f_{1,2}=y(-3A_{122}T_{1}^{\prime}(x)-\alpha_{2}+3\beta_{4}x^{2}-2\sigma_{3}x)+3A_{032}T_{1}^{\prime}(x)+\nu_{1}+\nu_{3}x^{2}-\sigma_{2}x,$
$\displaystyle
f_{2,2}=x(3A_{122}T_{2}^{\prime}(y)-\beta_{2}-3\beta_{4}y^{2}-2\nu_{3}y)+3A_{023}T_{2}^{\prime}(y)+\sigma_{1}+\sigma_{3}y^{2}+\sigma_{2}y,$
$\displaystyle
f_{3,2}=-2yA_{122}T_{2}^{\prime}(y)+2A_{032}T_{2}^{\prime}(y)-A_{122}T_{2}(y)+5A_{050}T_{1}^{\prime}(x)$
$\displaystyle\hphantom{f_{3,2}=}{}+\beta_{1}+\beta_{4}y^{3}+\nu_{3}y^{2}+\beta_{2}y,$
(8.5)

where

$\displaystyle V_{1}(x)\equiv T_{1}^{\prime}(x),\qquad V_{2}(y)\equiv
T_{2}^{\prime}(y).$

Next, substituting (8.3) and (8.5) into (8.4) we obtain the following non-
linear compatibility condition (NLCC)

$\displaystyle\text{NLCC}=T_{1}{}^{(4)}(3T_{1}^{\prime}(A_{032}-yA_{122})+\nu_{1}+3\beta_{4}x^{2}y+x(-\sigma_{2}+\nu_{3}x-2\sigma_{3}y)-\alpha_{2}y)$
$\displaystyle\hphantom{\text{NLCC}=}{}+T_{2}{}^{(4)}(3T_{2}^{\prime}(xA_{122}+A_{023})+\sigma_{1}-3\beta_{4}xy^{2}+y(\sigma_{2}+\sigma_{3}y-2\nu_{3}x)-\beta_{2}x)$
$\displaystyle\hphantom{\text{NLCC}=}{}+T_{1}{}^{(3)}(-9yA_{122}T_{1}^{\prime\prime}+9A_{032}T_{1}^{\prime\prime}-4\sigma_{2}+8\nu_{3}x+24\beta_{4}xy-8\sigma_{3}y)$
$\displaystyle\hphantom{\text{NLCC}=}{}+T_{2}{}^{(3)}(9xA_{122}T_{2}^{\prime\prime}+9A_{023}T_{2}^{\prime\prime}+4\sigma_{2}-8\nu_{3}x-24\beta_{4}xy+8\sigma_{3}y)$
$\displaystyle\hphantom{\text{NLCC}=}{}+12\nu_{3}T_{1}^{\prime\prime}+12\sigma_{3}T_{2}^{\prime\prime}-36\beta_{4}xT_{2}^{\prime\prime}+36\beta_{4}yT_{1}^{\prime\prime}=0,$
(8.6)

where $\alpha_{1}$, $\beta^{\prime}$s, $\nu^{\prime}$s and $\sigma^{\prime}$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 ${\mathcal{Y}}_{5}$, we identify two cases for which
the above NLCC (8.6) admits separation of variables in Cartesian coordinates,
namely:

  * $(i)$

$A_{122}\neq 0$, $A_{023}=A_{032}=0$,

  * $(ii)$

$A_{023}^{2}+A_{032}^{2}\neq 0$, $A_{122}=0$,

with $A_{050}$ and $A_{005}$ arbitrary. These two cases are
${\mathcal{S}}_{2}$-invariant under the permutation $x\Leftrightarrow y$
(thus, $p_{1}\Leftrightarrow p_{2}$). Let us recall that the Hamiltonian
$\mathcal{H}$ and the integral $\mathcal{X}$ are ${\mathcal{S}}_{2}$-invariant
and ${\mathcal{S}}_{2}$-antiinvariant, respectively. Moreover, if

  * $(iii)$

$A_{122}=A_{023}=A_{032}=0,$

with $A^{2}_{050}+A^{2}_{005}\neq 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.

##  9 Results

###  9.1 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)$.

$\bullet$ The system $Q_{1}$: $A_{122}\neq 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

$\displaystyle(T_{1}^{\prime}){}^{2}-2\beta_{4}x^{2}T_{1}^{\prime}-4\beta_{4}T_{1}x+\beta_{4}^{2}x^{4}=0,$
$\displaystyle(T_{2}^{\prime}){}^{2}-2\beta_{4}y^{2}T_{2}^{\prime}-4\beta_{4}T_{2}y+\beta_{4}^{2}y^{4}=0,$
(9.1)

$\beta_{4}\neq 0$ is a real constant. The corresponding fifth-order integral
of motion is given by

$\displaystyle{\mathcal{Y}}_{5}^{(Q_{1})}=p_{1}^{2}p_{2}^{3}x-p_{1}^{3}p_{2}^{2}y+p_{1}^{3}\big{(}{-}2yT_{2}^{\prime}-T_{2}+\beta_{4}y^{3}\big{)}+p_{2}^{3}\big{(}2xT_{1}^{\prime}+T_{1}-\beta_{4}x^{3}\big{)}$
$\displaystyle\hphantom{{\mathcal{Y}}_{5}^{(Q_{1})}=}{}+p_{1}^{2}p_{2}x\big{(}3T_{2}^{\prime}-3\beta_{4}y^{2}\big{)}+p_{1}p_{2}^{2}y\big{(}3\beta_{4}x^{2}-3T_{1}^{\prime}\big{)}$
$\displaystyle\hphantom{{\mathcal{Y}}_{5}^{(Q_{1})}=}{}+p_{1}\bigg{(}{-}\frac{3}{2}\beta_{4}x^{2}y^{2}T_{2}^{\prime\prime}+3\beta_{4}y^{3}T_{1}^{\prime}+\frac{3}{2}x^{2}T_{2}^{\prime}T_{2}^{\prime\prime}-6yT_{1}^{\prime}T_{2}^{\prime}-3T_{2}T_{1}^{\prime}\bigg{)}$
$\displaystyle\hphantom{{\mathcal{Y}}_{5}^{(Q_{1})}=}{}+p_{2}\bigg{(}\frac{3}{2}\beta_{4}x^{2}y^{2}T_{1}^{\prime\prime}-3\beta_{4}x^{3}T_{2}^{\prime}-\frac{3}{2}y^{2}T_{1}^{\prime}T_{1}^{\prime\prime}+6xT_{1}^{\prime}T_{2}^{\prime}+3T_{1}T_{2}^{\prime}\bigg{)}.$
(9.2)

From $\mathcal{X}$ and ${\mathcal{Y}}_{5}^{(Q_{1})}$, we built the quantity

$\displaystyle
C\equiv\big{\\{}{\mathcal{Y}}_{5}^{(Q_{1})},{\mathcal{X}}\big{\\}}_{\rm PB},$

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 $\big{(}\mathcal{X}$,
${\mathcal{Y}}_{5}^{(Q_{1})},C\big{)}$ 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}^{\prime}$, $T_{2}^{\prime}$
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

$\displaystyle
3\beta_{4}T_{1}^{2}+8\beta_{4}^{3/2}x^{3/2}T_{1}^{3/2}+6\beta_{4}^{2}x^{3}T_{1}-\beta_{4}^{3}x^{6}-\delta=0,$
$\displaystyle
3\beta_{4}T_{2}^{2}+8\beta_{4}^{3/2}y^{3/2}T_{2}^{3/2}+6\beta_{4}^{2}y^{3}T_{2}-\beta_{4}^{3}y^{6}-\delta=0,$
(9.3)

where $\delta$ is an arbitrary constant. In the case $\delta=0$, we
immediately obtain the particular solutions

$T_{1}(x)=\beta_{4}x^{3},\ \frac{\beta_{4}}{9}x^{3},\qquad\text{and}\qquad
T_{2}(y)=\beta_{4}y^{3},\ \frac{\beta_{4}}{9}y^{3},$

which correspond to a well-known lower-order superintegrable system.

$V_{1}(x)$$3$$2$$1$$-1$$x$$-1.5$$-1.0$$-0.5$$0.5$$1.0$$1.5$

Figure 3: Case $N=5$: the $x$-component $V_{1}(x)$ of the doubly exotic
potential $V(x,y)=V_{1}(x)+V_{2}(y)$ of type $Q_{1}$. It corresponds to the
fifth-order integral ${\mathcal{Y}}_{5}$ (9.2). From the algebraic equations
(9.3) we obtain the four solutions $V_{1,i}(x)=T^{\prime}_{1,i}$, $i=1,2,3,4$,
displayed above. In the case $Q_{1}$ the $y$-component $V_{2}(y)$ is of the
same form with four similar solutions $V_{2,i}(y)$. The values
$\beta_{4}=\delta=1$ were used.

$V=V_{1}(x)+V_{2}(y)$$-1.25$$-0.1$$-0.75$$-0.5$$-0.25$$0$$0.25$$0$$0.2$$0.4$$x$$0.6$$0.8$$1.0$$0$$0.2$$0.4$$y$$0.6$$0.8$$1.0$

Figure 4: A doubly exotic potential $Q_{1}$ corresponding to $N=5$. It admits
the fifth-order integral ${\mathcal{Y}}_{5}$ (9.2). It also possesses bounded
trajectories which by construction are closed and periodic. The values
$\beta_{4}=\delta=1$ were used.

The algebra generated by the integrals takes the form

$\displaystyle\\{C,{\mathcal{X}}\\}_{\rm
PB}=-24\beta_{4}{\mathcal{Y}}_{5}^{(Q_{1})},$
$\displaystyle\big{\\{}C,{\mathcal{Y}}_{5}^{(Q_{1})}\big{\\}}_{\rm
PB}=12X\big{(}{\mathcal{H}}^{2}-{\mathcal{X}}^{2}\big{)}^{2}-48\delta{\mathcal{X}}{\mathcal{H}}.$

In the corresponding quantum system analyzed in [1], the case $(i)$ splits
into two subclasses of integrals ${\mathcal{Y}}_{5}$ that solely differ in
their lower order $\hbar$-dependable terms. Consequently, two systems called
$Q_{1}$ and $Q_{2}$ occur. However, in the classical limit $\hbar\rightarrow
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}\neq 0$, $A_{122}=A_{032}=0$) and $Q_{4}$
($A_{023}A_{005}\neq 0$, $A_{122}=A_{050}=A_{032}=0$) are _not_
superintegrable (like in the quantum case).

Case $(ii)$.

$\bullet$ The system $Q_{5}$: $A_{023}\neq 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}$

$\displaystyle
5A_{050}(T_{1}^{\prime})^{3}-12\tau^{2}xT_{1}^{\prime}+3\beta_{1}(T_{1}^{\prime})^{2}-12\tau^{2}T_{1}+\mu=0,$
(9.4)

$\tau\neq 0$, $\beta_{1}$ and $\mu$ are real constants, whereas

$\displaystyle V_{2}\equiv T_{2}^{\prime}=\pm 2\tau\sqrt{-y}.$

The corresponding highest-order integral of motion reads

$\displaystyle{\mathcal{Y}}_{5}^{(Q_{5})}=A_{050}p_{1}^{5}+p_{1}^{2}p_{2}^{3}+p_{1}^{3}(5A_{050}T_{1}^{\prime}+\beta_{1})+p_{1}\bigg{(}\frac{15}{2}A_{050}(T_{1}^{\prime}){}^{2}+3\beta_{1}T_{1}^{\prime}-6x\tau^{2}\bigg{)}$
$\displaystyle\hphantom{{\mathcal{Y}}_{5}^{(Q_{5})}=}{}+6\tau\sqrt{-y}p_{2}p_{1}^{2}+2T_{1}^{\prime}p_{2}^{3}+12\tau\sqrt{-y}T_{1}^{\prime}p_{2}.$
(9.5)

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}\rightarrow A_{005}$ and making the permutation
$x\Leftrightarrow y$ ($V_{1}\Leftrightarrow V_{2}$) in (9.4) and (9.5).

$V_{1}(x)$$4$$2$$-2$$-4$$-6$$x$$-10$$-5$$5$$10$

Figure 5: Case $N=5$: the $x$-component $V_{1}(x)$ of the doubly exotic
potential $V(x,y)=V_{1}(x)+V_{2}(y)$ of type $Q_{5}$. It admits the fifth-
order integral ${\mathcal{Y}}_{5}$ (9.5). From the equation (9.4) we obtain
the three numerical solutions $V_{1,i}(x)=T^{\prime}_{1,i}$, $i=1,2,3$,
displayed above. The values $A_{050}=\frac{1}{5}$, $\tau=\frac{1}{\sqrt{12}}$,
$\beta_{1}=1$ and $\mu=-3$ were used.

In this case, the algebra generated by the integrals
$C=\big{\\{}{\mathcal{Y}}_{5}^{(Q_{5})},{\mathcal{X}}\big{\\}}_{\rm PB}$,
${\mathcal{Y}}_{5}^{(Q_{5})}$ and $\mathcal{X}$ takes the form

$\displaystyle\\{C,{\mathcal{X}}\\}_{\rm PB}=0,$
$\displaystyle\big{\\{}C,{\mathcal{Y}}_{5}^{(Q_{5})}\big{\\}}_{\rm
PB}=-144\tau^{4}{({\mathcal{H}}+{\mathcal{X}})},$

and it does not provide further information about the solutions of (9.4).

$\bullet$ The system $Q_{6}$: $A_{023}A_{032}\neq 0$,
$A_{122}=A_{050}=A_{005}=0$. This system corresponds to $A_{023}A_{032}\neq
0$, all other $A_{ijk}=0$.

$\displaystyle
3A_{032}(T_{1}^{\prime}){}^{2}-\sigma_{2}xT_{1}^{\prime}-\sigma_{2}T_{1}=0,$
$\displaystyle
3A_{023}(T_{2}^{\prime}){}^{2}+\sigma_{2}yT_{2}^{\prime}+\sigma_{2}T_{2}=0,$
(9.6)

$\sigma_{2}\neq 0$ is a real constant.

$\displaystyle{\mathcal{Y}}_{5}^{(Q_{6})}=A_{032}p_{2}^{2}p_{1}^{3}+A_{023}p_{2}^{3}p_{1}^{2}+2p_{2}^{3}A_{023}T_{1}^{\prime}+2p_{1}^{3}A_{032}T_{2}^{\prime}+p_{1}p_{2}^{2}\bigg{(}3A_{032}T_{1}^{\prime}-\frac{\sigma_{2}x}{2}\bigg{)}$
$\displaystyle\hphantom{{\mathcal{Y}}_{5}^{(Q_{6})}=}{}+p_{1}^{2}p_{2}\bigg{(}3A_{023}T_{2}^{\prime}+\frac{\sigma_{2}y}{2}\bigg{)}+p_{1}\bigg{(}3xA_{023}T_{2}^{\prime}T_{2}^{\prime\prime}+6A_{032}T_{1}^{\prime}T_{2}^{\prime}+\frac{1}{2}\sigma_{2}xyT_{2}^{\prime\prime}\bigg{)}$
$\displaystyle\hphantom{{\mathcal{Y}}_{5}^{(Q_{6})}=}{}+p_{2}\bigg{(}3yA_{032}T_{1}^{\prime}T_{1}^{\prime\prime}+6A_{023}T_{1}^{\prime}T_{2}^{\prime}-\frac{1}{2}\sigma_{2}xyT_{1}^{\prime\prime}\bigg{)}.$
(9.7)

$V=V_{1}(x)+V_{2}(y)$$-0.4$$-0.3$$-0.2$$-0.1$$0$$0.1$$-5$$-4$$x$$-3$$-2$$-5$$-4$$-3$$-2$$y$

Figure 6: A doubly exotic potential $Q_{6}$ corresponding to $N=5$. It admits
the fifth-order integral ${\mathcal{Y}}_{5}$ (9.7). It also possesses bounded
trajectories which by construction are closed and periodic. The values
$A_{032}=-A_{023}=\frac{1}{12}$, $\sigma_{1}=1$ were used.

In this case, the algebra generated by the integrals
$C=\big{\\{}{\mathcal{Y}}_{5}^{(Q_{6})},{\mathcal{X}}\big{\\}}_{\rm PB}$,
${\mathcal{Y}}_{5}^{(Q_{6})}$ and $\mathcal{X}$ takes the form

$\displaystyle\\{C,{\mathcal{X}}\\}_{\rm PB}=0,$
$\displaystyle\big{\\{}C,{\mathcal{Y}}_{5}^{(Q_{6})}\big{\\}}_{\rm
PB}=2{\mathcal{X}}\sigma_{2}^{2}\big{(}{\mathcal{H}}^{2}-{\mathcal{X}}^{2}\big{)}^{2},$

and it does not provide further information about the solutions of (9.6).
However, it is easy to check that

$\displaystyle T_{1}(x)=\frac{W^{2}(x)-\sigma_{2}x^{2}}{12A_{032}},\qquad
T_{2}(y)=-\frac{W^{2}(y)-\sigma_{2}y^{2}}{12A_{023}},$

satisfy (9.6), where $W=W(z)$ is given by the following third order polynomial
equation

$\displaystyle\big{(}W-z\sqrt{\sigma_{2}}\big{)}\big{(}W+2z\sqrt{\sigma_{2}}\big{)}^{2}+\tau=0,$

here $\tau\neq 0$ is an integration constant.

$\bullet$ 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$.

$\displaystyle
3\beta_{1}(T_{1}^{\prime}){}^{2}-12\tau^{2}xT_{1}^{\prime}-12\tau^{2}T_{1}+\mu=0,\qquad
V_{2}\equiv T_{2}^{\prime}=\pm 2\tau\sqrt{-y},$

$\tau\neq 0$, $\beta_{1}$ and $\mu$ are real constants.

$\displaystyle{\mathcal{Y}}_{5}^{(Q_{7})}=p_{1}^{2}p_{2}^{3}+p_{1}^{3}\beta_{1}+p_{1}\big{(}3\beta_{1}T_{1}^{\prime}-6x\tau^{2}\big{)}+6\tau\sqrt{-y}p_{2}p_{1}^{2}+2T_{1}^{\prime}p_{2}^{3}+12\tau\sqrt{-y}T_{1}^{\prime}p_{2}.$

Case $(iii)$.

$\bullet$ The system $Q_{8}$: $A_{050}^{2}+A^{2}_{005}\neq 0$,
$A_{122}=A_{023}=A_{032}=0$. This system corresponds to the case
$A_{050}^{2}+A^{2}_{005}\neq 0$ and all other $A_{ijk}=0$.

$\displaystyle
A_{050}T_{1}^{\prime}{}^{3}+\beta_{1}T_{1}^{\prime}{}^{2}+\theta_{1}T_{1}^{\prime}-\Lambda
x+\kappa_{1}=0,$ $\displaystyle
A_{005}T_{2}^{\prime}{}^{3}+\alpha_{1}T_{2}^{\prime}{}^{2}+\phi_{1}T_{2}^{\prime}+\Lambda
y+\omega_{1}=0,$ (9.8)

$\Lambda\neq 0$, $\alpha_{1},\beta_{1},\theta_{1},\kappa_{1},\phi_{1}$ and
$\omega_{1}$ are real constants.

$\displaystyle{\mathcal{Y}}_{5}^{(Q_{8})}=A_{050}p_{1}^{5}+A_{005}p_{2}^{5}+p_{1}^{3}\bigg{(}5A_{050}T_{1}^{\prime}+\frac{5\beta_{1}}{3}\bigg{)}+p_{2}^{3}\bigg{(}5A_{005}T_{2}^{\prime}+\frac{5\alpha_{1}}{3}\bigg{)}$
$\displaystyle\hphantom{{\mathcal{Y}}_{5}^{(Q_{8})}=}{}+p_{1}\bigg{(}\frac{15}{2}A_{050}T_{1}^{\prime}{}^{2}+5\beta_{1}T_{1}^{\prime}+\frac{5\theta_{1}}{2}\bigg{)}+p_{2}\bigg{(}\frac{15}{2}A_{005}T_{2}^{\prime}{}^{2}+5\alpha_{1}T_{2}^{\prime}+\frac{5\phi_{1}}{2}\bigg{)}.$

In this case, for the solutions of (9.8) the Poisson bracket
$C=\big{\\{}{\mathcal{Y}}_{5}^{(Q_{8})},{\mathcal{X}}\big{\\}}_{\rm
PB}\propto\lambda$. Thus, the algebra of the integrals of motion takes the
form

$\\{C,{\mathcal{X}}\\}_{\rm
PB}=0,\qquad\big{\\{}C,{\mathcal{Y}}_{5}^{(Q_{8})}\big{\\}}_{\rm PB}=0.$

Again, the system $Q_{8}$ was found in [28] by means of Heisenberg-type higher
order symmetries. All three conserved quantities
$\big{(}{\mathcal{H}},{\mathcal{X}},{\mathcal{Y}}_{5}^{(Q_{8})}\big{)}$ admit
separation of variables in Cartesian coordinates. The particular case with
$\beta_{1}=\theta_{1}=\kappa_{1}=\alpha_{1}=\phi_{1}=\omega_{2}=0$, thus
$V_{1}\propto x^{\frac{1}{3}}$ and $V_{2}\propto y^{\frac{1}{3}}$, was studied
in [24] using action-angle variables.

$\bullet$ 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$.

$\displaystyle
A_{050}T_{1}^{\prime}{}^{3}+\beta_{1}T_{1}^{\prime}{}^{2}+\theta_{1}T_{1}^{\prime}-\Lambda
x+\kappa_{1}=0,$
$\displaystyle\alpha_{1}T_{2}^{\prime}{}^{2}+\phi_{1}T_{2}^{\prime}+\Lambda
y+\omega_{1}=0,$

$\Lambda\neq 0$, $\beta_{1}$, $\theta_{1}$, $\kappa_{1}$, $\phi_{1}$,
$\omega_{1}$ and $\alpha_{1}^{2}+\phi_{1}^{2}\neq 0$ are real constants.

$\displaystyle{\mathcal{Y}}_{5}^{(Q_{9})}=A_{050}p_{1}^{5}+p_{1}^{3}\bigg{(}5A_{050}T_{1}^{\prime}+\frac{5\beta_{1}}{3}\bigg{)}+\frac{5\alpha_{1}}{3}p_{2}^{3}+p_{1}\bigg{(}\frac{15}{2}A_{050}T_{1}^{\prime}{}^{2}+5\beta_{1}T_{1}^{\prime}+\frac{5\theta_{1}}{2}\bigg{)}$
$\displaystyle\hphantom{{\mathcal{Y}}_{5}^{(Q_{9})}=}{}+p_{2}\bigg{(}5\alpha_{1}T_{2}^{\prime}+\frac{5\phi_{1}}{2}\bigg{)}.$

##  10 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
(${\mathcal{H}},{\mathcal{X}},{\mathcal{Y}}_{N}$). Let us summarize the main
results reported in this paper:

  1. 1.

Higher-order ($N>2$) superintegrable classical systems

${\mathcal{H}}=\frac{1}{2}\big{(}p_{1}^{2}+p_{2}^{2}\big{)}+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 ${\mathcal{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. 2.

From the equation $\\{{\mathcal{Y}}_{N},{\mathcal{H}}\\}_{\rm 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}=\\{{\mathcal{X}},{\mathcal{Y}}_{N}\\}_{\rm
PB},{\mathcal{X}},{\mathcal{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. 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.

  4. 4.

The present study suggests to explore the inverse problem, namely we take two
polynomial functions $({\mathcal{A}}$ and ${\mathcal{B}})$ in momentum
variables ($p_{1},p_{2}$) and construct the new object
${\mathcal{C}}=\\{{\mathcal{A}},{\mathcal{B}}\\}_{\rm PB}$. If the algebra
generated by $({\mathcal{A}},{\mathcal{B}},{\mathcal{C}})$ is a closed
polynomial algebra with polynomial coefficients in ${\mathcal{H}}$, then under
what conditions this closure relations imply that ${\mathcal{A}}$ and
${\mathcal{B}}$ are integrals, i.e., they Poisson commute with $\mathcal{H}$?

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

$\displaystyle{\mathcal{Y}}_{N}^{(\rm Doubly\
exotic)}=L^{(N-4)}_{z}p_{1}^{2}p_{2}^{2}+{\rm(lower\ order\ terms)},\qquad
N\geq 5,$ $\displaystyle\big{\\{}{\mathcal{Y}}_{N}^{(\rm Doubly\
exotic)},{\mathcal{H}}\big{\\}}_{\rm PB}=0.$

The associated potential $V$ can be written as follows

$\displaystyle
V=V_{1}(x)+V_{2}(y)=\mathcal{G}^{\prime}(x;N)+\mathcal{G}^{\prime}(y;N),$

here $\mathcal{G}=\mathcal{G}(u;N)$ obeys a nonlinear first-order ODE of the
form

$\displaystyle{\mathcal{G}}^{\prime}\big{[}6u^{N-4}{\mathcal{G}}^{\prime}+4{(N-5)}u^{N-5}{\mathcal{G}}+F_{1}(u)+\sigma
u^{N-2}\big{]}$
$\displaystyle\qquad{}+{\mathcal{G}}\big{[}2{(N-5)}u^{N-6}{\mathcal{G}}+F_{2}(u)-2\sigma
u^{N-3}\big{]}+F_{3}(u)+bu^{N}=0.$ (10.1)

The three functions $F_{q}$’s in (10.1) are polynomials in the variable $u$ of
degree at most $(N-1)$, and $\sigma$, $b$ are real parameters as well. The
equation (10.1) is in complete agreement with the limit $\hbar\rightarrow 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.

### Acknowledgments

İY and AMER during a sabbatical leave and a postdoctoral academic stay at the
Centre de Recherches Mathématiques, Université de Montréal, respectively, were
introduced to the subject of higher-order superintegrability by Pavel
Winternitz. His enormous influence is present in this study as it does in the
whole subject. It is with admiration and great affection that we dedicate this
paper to his memory. We thank the anonymous referees and the editor for their
valuable comments and constructive suggestions on the manuscript.

## References

  * [1] Abouamal I., Winternitz P., Fifth-order superintegrable quantum systems separating in Cartesian coordinates: doubly exotic potentials, J. Math. Phys. 59 (2018), 022104, 27 pages, arXiv:1708.03379.
  * [2] Baldwin D., Hereman W., PainleveTestV4-2018.m: A Mathematica package for the Painlevé test of systems of nonlinear ordinary and partial differential equations, 2018, https://inside.mines.edu/~whereman/software/painleve/mathematica/.
  * [3] Bertrand J.L.F., Théorème relatif au mouvement d’un point attiré vers un centre fixe, C. R. Acad. Sci. Paris 77 (1873), 849–853.
  * [4] Bertrand S., Kubu O., Šnobl L., On superintegrability of 3D axially-symmetric non-subgroup-type systems with magnetic fields, J. Phys. A: Math. Theor. 54 (2021), 015201, 27 pages, arXiv:2008.01987.
  * [5] Bérubé J., Winternitz P., Integrable and superintegrable quantum systems in a magnetic field, J. Math. Phys. 45 (2004), 1959–1973, arXiv:math-ph/0311051.
  * [6] Chanu C.M., Degiovanni L., Rastelli G., The Tremblay–Turbiner–Winternitz system as extended Hamiltonian, J. Math. Phys. 55 (2014), 122701, 8 pages, arXiv:1404.4825.
  * [7] Chanu C.M., Rastelli G., Extended Hamiltonians and shift, ladder functions and operators, Ann. Physics 386 (2017), 254–274, arXiv:1705.09519.
  * [8] Daskaloyannis C., Generalized deformed oscillator and nonlinear algebras, J. Phys. A: Math. Gen. 24 (1991), L789–L794.
  * [9] Désilets J.F., Winternitz P., Yurduşen İ., Superintegrable systems with spin and second-order integrals of motion, J. Phys. A: Math. Theor. 45 (2012), 475201, 26 pages, arXiv:1208.2886.
  * [10] Dorizzi B., Grammaticos B., Ramani A., Winternitz P., Integrable Hamiltonian systems with velocity-dependent potentials, J. Math. Phys. 26 (1985), 3070–3079.
  * [11] Drach J., Sur l’integration logique des équations de la dynamique à deux variables: Forces constructives. Intégrales cubiques. Mouvements dans le plan, C. R. Acad. Sci. Paris 200 (1935), 22–26.
  * [12] Drach J., Sur l’integration logique et sur la tranformation des équations de la dynamique à deux variables: Forces conservatives. Intégrales cubiques. Mouvements dans le plan, C. R. Acad. Sci. Paris 200 (1935), 599–602.
  * [13] Escobar-Ruiz A.M., Linares R., Winternitz P., New infinite families of $N$th-order superintegrable systems separating in Cartesian coordinates, J. Phys. A: Math. Theor. 53 (2020), 445203, 26 pages, arXiv:2004.12173.
  * [14] Escobar-Ruiz A.M., López Vieyra J.C., Winternitz P., Fourth order superintegrable systems separating in polar coordinates. I. Exotic potentials, J. Phys. A: Math. Theor. 50 (2017), 495206, 35 pages, arXiv:1706.08655.
  * [15] Escobar-Ruiz A.M., López Vieyra J.C., Winternitz P., Yurduşen İ., Fourth-order superintegrable systems separating in polar coordinates. II. Standard potentials, J. Phys. A: Math. Theor. 51 (2018), 455202, 24 pages, arXiv:1804.05751.
  * [16] Escobar-Ruiz A.M., Winternitz P., Yurduşen İ., General $N$th-order superintegrable systems separating in polar coordinates, J. Phys. A: Math. Theor. 51 (2018), 40LT01, 12 pages, arXiv:1806.06849.
  * [17] Evans N.W., Super-integrability of the Winternitz system, Phys. Lett. A 147 (1990), 483–486.
  * [18] Evans N.W., Superintegrability in classical mechanics, Phys. Rev. A 41 (1990), 5666–5676.
  * [19] Fradkin D.M., Three-dimensional isotropic harmonic oscillator and ${\rm SU}_{3}$, Amer. J. Phys. 33 (1965), 207–211.
  * [20] Friš J., Mandrosov V., Smorodinsky Y.A., Uhlíř M., Winternitz P., On higher symmetries in quantum mechanics, Phys. Lett. 16 (1965), 354–356.
  * [21] Goldstein H., Classical mechanics, 2nd ed., Addison-Wesley Series in Physics, Addison-Wesley Publishing Co., Reading, Mass., 1980.
  * [22] Gravel S., Hamiltonians separable in Cartesian coordinates and third-order integrals of motion, J. Math. Phys. 45 (2004), 1003–1019, arXiv:math-ph/0302028.
  * [23] Gravel S., Winternitz P., Superintegrability with third-order integrals in quantum and classical mechanics, J. Math. Phys. 43 (2002), 5902–5912, arXiv:math-ph/0206046.
  * [24] Grigoriev Y.A., Tsiganov A.V., On superintegrable systems separable in Cartesian coordinates, Phys. Lett. A 382 (2018), 2092–2096, arXiv:1712.07321.
  * [25] Grosche C., Pogosyan G.S., Sissakian A.N., Path integral discussion for Smorodinsky–Winternitz potentials. II. The two- and three-dimensional sphere, Fortschr. Phys. 43 (1995), 523–563, arXiv:hep-th/9402121.
  * [26] Grosche C., Pogosyan G.S., Sissakian A.N., Path integral approach to superintegrable potentials. Two-dimensional hyperboloid, Phys. Part. Nuclei 27 (1996), 244–278.
  * [27] Grosche C., Pogosyan G.S., Sissakian A.N., Path integral discussion for superintegrable potentials. IV. Three-dimensional pseudosphere, Phys. Part. Nuclei 28 (1997), 486–519.
  * [28] Güngör F., Kuru S., Negro J., Nieto L.M., Heisenberg-type higher order symmetries of superintegrable systems separable in Cartesian coordinates, Nonlinearity 30 (2017), 1788–1808, arXiv:1411.6216.
  * [29] Hietarinta J., Direct methods for the search of the second invariant, Phys. Rep. 147 (1987), 87–154.
  * [30] Jauch J.M., Hill E.L., On the problem of degeneracy in quantum mechanics, Phys. Rev. 57 (1940), 641–645.
  * [31] Kalnins E.G., Separation of variables for Riemannian spaces of constant curvature, Pitman Monographs and Surveys in Pure and Applied Mathematics, Vol. 28, Longman Scientific & Technical, Harlow, John Wiley & Sons, Inc., New York, 1986.
  * [32] Kalnins E.G., Kress J.M., Miller Jr. W., Second-order superintegrable systems in conformally flat spaces. V. Two- and three-dimensional quantum systems, J. Math. Phys. 47 (2006), 093501, 25 pages.
  * [33] Kalnins E.G., Kress J.M., Miller Jr. W., Nondegenerate 2D complex Euclidean superintegrable systems and algebraic varieties, J. Phys. A: Math. Theor. 40 (2007), 3399–3411, arXiv:0708.3044.
  * [34] Kalnins E.G., Kress J.M., Miller Jr. W., Separation of variables and superintegrability. The symmetry of solvable systems, IOP Expanding Physics, IOP Publishing, Bristol, 2018.
  * [35] Kalnins E.G., Kress J.M., Miller Jr. W., Post S., Structure theory for second order 2D superintegrable systems with 1-parameter potentials, SIGMA 5 (2009), 008, 24 pages, arXiv:0901.3081.
  * [36] Kalnins E.G., Kress J.M., Miller Jr. W., Winternitz P., Superintegrable systems in Darboux spaces, J. Math. Phys. 44 (2003), 5811–5848, arXiv:math-ph/0307039.
  * [37] Kalnins E.G., Kress J.M., Winternitz P., Superintegrability in a two-dimensional space of nonconstant curvature, J. Math. Phys. 43 (2002), 970–983, arXiv:math-ph/0108015.
  * [38] Kalnins E.G., Miller Jr. W., Pogosyan G.S., Nondegenerate superintegrable systems in $n$-dimensional spaces of constant curvature, Phys. Atomic Nuclei 70 (2007), 545–553, arXiv:0912.2278.
  * [39] Kalnins E.G., Williams G.C., Miller Jr. W., Pogosyan G.S., On superintegrable symmetry-breaking potentials in $N$-dimensional Euclidean space, J. Phys. A: Math. Gen. 35 (2002), 4755–4773.
  * [40] Laplace P.S., Traité de mécanique céleste, Duprat, Paris, 1799.
  * [41] Lenz W., Über den Bewegungsverlauf und die Quantenzustände der gestörten Keplerbewegung, Z. Phys. 24 (1924), 197–207.
  * [42] Makarov A.A., Smorodinsky J.A., Valiev K., Winternitz P., A systematic search for nonrelativistic systems with dynamical symmetries, Nuovo Cimento A 52 (1967), 1061–1084.
  * [43] Marchesiello A., Šnobl L., Winternitz P., Three-dimensional superintegrable systems in a static electromagnetic field, J. Phys. A: Math. Theor. 48 (2015), 395206, 24 pages, arXiv:1507.04632.
  * [44] Marchesiello A., Šnobl L., Winternitz P., Spherical type integrable classical systems in a magnetic field, J. Phys. A: Math. Theor. 51 (2018), 135205, 24 pages.
  * [45] Marquette I., Superintegrability with third order integrals of motion, cubic algebras, and supersymmetric quantum mechanics. II. Painlevé transcendent potentials, J. Math. Phys. 50 (2009), 095202, 18, arXiv:0811.1568.
  * [46] Marquette I., Classical ladder operators, polynomial Poisson algebras, and classification of superintegrable systems, J. Math. Phys. 53 (2012), 012901, 12 pages, arXiv:1109.4471.
  * [47] Marquette I., Sajedi M., Winternitz P., Fourth order superintegrable systems separating in Cartesian coordinates I. Exotic quantum potentials, J. Phys. A: Math. Theor. 50 (2017), 315201, 29 pages, arXiv:1703.09751.
  * [48] Marquette I., Winternitz P., Polynomial Poisson algebras for classical superintegrable systems with a third-order integral of motion, J. Math. Phys. 48 (2007), 012902, 16 pages, arXiv:math-ph/0608021.
  * [49] Miller Jr. W., Symmetry and separation of variables, Encyclopedia of Mathematics and its Applications, Vol. 4, Addison-Wesley Publishing Co., Reading, Mass. – London – Amsterdam, 1977.
  * [50] Miller Jr. W., Second order superintegrable systems in three dimensions, SIGMA 1 (2005), 015, 17 pages, arXiv:nlin.SI/0511035.
  * [51] Miller Jr. W., Post S., Winternitz P., Classical and quantum superintegrability with applications, J. Phys. A: Math. Theor. 46 (2013), 423001, 97 pages, arXiv:1309.2694.
  * [52] Nikitin A.G., Higher-order symmetry operators for Schrödinger equation, in Superintegrability in Classical and Quantum Systems, CRM Proc. Lecture Notes, Vol. 37, Amer. Math. Soc., Providence, RI, 2004, 137–144, arXiv:1603.01715.
  * [53] Nikitin A.G., Matrix superpotentials and superintegrable systems for arbitrary spin, J. Phys. A: Math. Theor. 45 (2012), 225205, 13 pages, arXiv:1201.4929.
  * [54] Nikitin A.G., New exactly solvable systems with Fock symmetry, J. Phys. A: Math. Theor. 45 (2012), 485204, 9 pages, arXiv:1205.3094.
  * [55] Popper I., Post S., Winternitz P., Third-order superintegrable systems separable in parabolic coordinates, J. Math. Phys. 53 (2012), 062105, 20 pages, arXiv:1204.0700.
  * [56] Post S., Winternitz P., An infinite family of superintegrable deformations of the Coulomb potential, J. Phys. A: Math. Theor. 43 (2010), 222001, 11 pages, arXiv:1003.5230.
  * [57] Post S., Winternitz P., General $N$th order integrals of motion in the Euclidean plane, J. Phys. A: Math. Theor. 48 (2015), 405201, 24 pages, arXiv:1501.00471.
  * [58] Rañada M.F., Superintegrable $n=2$ systems, quadratic constants of motion, and potentials of Drach, J. Math. Phys. 38 (1997), 4165–4178.
  * [59] Rodriguez M.A., Winternitz P., Quantum superintegrability and exact solvability in $n$ dimensions, J. Math. Phys. 43 (2002), 1309–1322, arXiv:math-ph/0110018.
  * [60] Runge C., Vektoranalysis I, Hirzel, Leipzig, 1919.
  * [61] Thompson G., Polynomial constants of motion in flat space, J. Math. Phys. 25 (1984), 3474–3478.
  * [62] Tremblay F., Turbiner A.V., Winternitz P., An infinite family of solvable and integrable quantum systems on a plane, J. Phys. A: Math. Theor. 42 (2009), 242001, 10 pages, arXiv:0904.0738.
  * [63] Tremblay F., Turbiner A.V., Winternitz P., Periodic orbits for an infinite family of classical superintegrable systems, J. Phys. A: Math. Theor. 43 (2010), 015202, 14 pages, arXiv:0910.0299.
  * [64] Tremblay F., Winternitz P., Third-order superintegrable systems separating in polar coordinates, J. Phys. A: Math. Theor. 43 (2010), 175206, 17 pages, arXiv:1002.1989.
  * [65] Tsiganov A.V., The Drach superintegrable systems, J. Phys. A: Math. Gen. 33 (2000), 7407–7422, arXiv:nlin.SI/0001053.
  * [66] Tsiganov A.V., Addition theorems and the Drach superintegrable systems, J. Phys. A: Math. Theor. 41 (2008), 335204, 16 pages, arXiv:0805.3443.
  * [67] Tsiganov A.V., On maximally superintegrable systems, Regul. Chaotic Dyn. 13 (2008), 178–190, arXiv:0711.2225.
  * [68] Winternitz P., Yurduşen İ., Integrable and superintegrable systems with spin, J. Math. Phys. 47 (2006), 103509, 10 pages, arXiv:math-ph/0604050.
  * [69] Winternitz P., Yurduşen İ., Integrable and superintegrable systems with spin in three-dimensional Euclidean space, J. Phys. A: Math. Theor. 42 (2009), 385203, 20 pages, arXiv:0906.1021.
  * [70] Yurduşen İ., Tuncer O.O., Winternitz P., Superintegrable systems with spin and second-order tensor and pseudo-tensor integrals of motion, J. Phys. A: Math. Theor. 54 (2021), 305201, 32 pages, arXiv:2104.07436.


# AutoPINN: When AutoML Meets Physics-Informed Neural Networks

Xinle Wu1, Dalin Zhang1, Miao Zhang1, Chenjuan Guo1, Shuai Zhao2, Yi Zhang2,
Huai Wang2, Bin Yang1 1Department of Computer Science, Aalborg
University2Department of Energy Technology, Aalborg University xinlewu,
dalinz, miaoz, cguo<EMAIL_ADDRESS>szh, yiz<EMAIL_ADDRESS>
(2018)

###### Abstract.

Physics-Informed Neural Networks (PINNs) have recently been proposed to solve
scientific and engineering problems, where physical laws are introduced into
neural networks as prior knowledge. With the embedded physical laws, PINNs
enable the estimation of critical parameters, which are unobservable via
physical tools, through observable variables. For example, Power Electronic
Converters (PECs) are essential building blocks for the green energy
transition. PINNs have been applied to estimate the capacitance, which is
unobservable during PEC operations, using current and voltage, which can be
observed easily during operations. The estimated capacitance facilitates self-
diagnostics of PECs.

Existing PINNs are often manually designed, which is time-consuming and may
lead to suboptimal performance due to a large number of design choices for
neural network architectures and hyperparameters. In addition, PINNs are often
deployed on different physical devices, e.g., PECs, with limited and varying
resources. Therefore, it requires designing different PINN models under
different resource constraints, making it an even more challenging task for
manual design. To contend with the challenges, we propose Automated Physics-
Informed Neural Networks (AutoPINN), a framework that enables the automated
design of PINNs by combining AutoML and PINNs. Specifically, we first tailor a
search space that allows finding high-accuracy PINNs for PEC internal
parameter estimation. We then propose a resource-aware search strategy to
explore the search space to find the best PINN model under different resource
constraints. We experimentally demonstrate that AutoPINN is able to find more
accurate PINN models than human-designed, state-of-the-art PINN models using
fewer resources.

AutoML, Physics-informed Neural Network, Power Electronics

††copyright: acmcopyright††journalyear: 2018††doi:
XXXXXXX.XXXXXXX††conference: Make sure to enter the correct conference title
from your rights confirmation emai; June 03–05, 2018; Woodstock, NY††price:
15.00††isbn: 978-1-4503-XXXX-X/18/06††ccs: Computing methodologies Neural
networks††ccs: Applied computing Physical sciences and engineering

##  1\. Introduction

A variety of Physics-Informed Neural Networks (PINNs) has recently been
proposed to solve Scientific Machine Learning (SciML) problems, achieving
better performance than purely data-driven models in many engineering domains
(Karniadakis et al., 2021; Raissi et al., 2019; Mao et al., 2020; Cai et al.,
2021, 2022; Misyris et al., 2020). PINN models are typically built by
incorporating physical laws as hard or soft constraints into conventional
neural networks (NNs), which are then trained normally, e.g., using back-
propagation. For instance, Ji et al. (2021) transforms the evolution of the
species concentrations in chemical kinetics as the ordinary differential
equation, which is then encoded into the loss function to make the network
satisfy the governing equations. Since a PINN model embraces certain physical
laws, it can be also used to conduct parameter inference from observations,
e.g., estimating internal parameters of physical devices that are unobservable
via physical instrumentation. A more recent work (Zhao et al., 2022) builds a
PINN to estimate the internal parameters in Power Electronic Converters
(PECs), by incorporating the Buck converter physical knowledge into the
network training to estimate unknown model coefficients based on the
measurable data.

Although the PINN model has been applied to different industrial and physical
problems with achieving promising performance, it still has two main
limitations. First, existing physics-informed neural networks are manually
designed, while the inappropriate network architectures or hyperparameter
settings usually lead to the failure of learning relevant physical phenomena
for slightly more complex problems (Krishnapriyan et al., 2021). In other
words, different neural networks may favour different industrial or physical
problems. The range of hyperparameters and possible architectures of neural
networks constitute a huge design space so it is unrealistic for even experts
to find an optimal PINN model through extensive trials. Second, PINN models
are often deployed on different physical devices with limited computing
capabilities and storage, on e.g., different types of PECs. In addition,
different devices may also have quite diverse levels of resources. Thus, it is
very time-consuming to manually design PINN models that satisfy different
hardware constraints.

Recently, Automated Machine Learning (AutoML) has been extensively applied to
automatically design machine learning models, and the automatically designed
models usually achieve better performance than their manually designed
counterparts (Liu et al., 2018; Xu et al., 2019; Zoph and Le, 2017; Pham et
al., 2018; Wu et al., 2021, 2022). This naturally inspires to use AutoML to
solve the above two limitations. In this paper, we propose an Automated
Physics-Informed Machine Learning framework (AutoPINN) to design hardware-
aware PINN models, where predicting the internal parameters of PECs with our
AutoPINN is taken as an exemplary application. Specifically, we first tailor a
search space for this task, which contains a wide variety of design choices.
Next, we design a reinforcement learning-based search strategy to explore the
search space to find the best PINN model. Furthermore, we introduce a
hardware-aware reward into the search objective to search for high-accuracy
PINN models under given hardware constraints.

To the best of our knowledge, this is the first study integrating AutoML and
PINN. The contributions can be summarized as 1) a judiciously designed search
space for the task of estimating the internal parameters of PECs with the
physics-informed neural network; 2) a reinforcement learning-based search
strategy to explore the search space, and a refined search objective to search
for the best PINN model under given hardware constraints; 3) extensive
experiments on a public dataset, demonstrating that the proposed framework is
able to find more accurate and lightweight PINN models than manually designed
state-of-the-arts.

(a) PINN

(b) AutoPINN

Figure 1. PINN and AutoPINN. (a) shows the architecture of the existing PINN
model. (b) shows the overview of AutoPINN, where $n_{i}$ denotes the number of
units at the $i$-th layer, $a_{i}$ denotes the activation function at the
$i$-th layer.

##  2\. Methods

###  2.1. Physics-informed neural network

In this paper, we focus on the problem of estimating the internal parameters
of Power Electronic Converters (PECs) (Zhao et al., 2022), such as inductance
$L$, inductor resistance $RL$, and capacitance $C$, which are important
indicators for identifying the aging status of PECs. Accurately estimating
these values facilitates early replacement of the PECs that are about to fail.

For this task, a simple data-driven solution is to build a neural network to
predict the internal parameters of a PEC with its external observable values
(e.g., inductor current and capacitor voltage) as input. This solution demands
the measured internal parameters as labels for training the neural network.
However, during PEC operations, we cannot directly measure such interval
parameters, because measuring them requires to disassemble the PEC, thus
interrupting the normal operations of the PEC. This makes the simple data-
driven solution inapplicable.

PINNs are able to estimate the internal parameters without requiring observing
the internal parameters as labels. For example, recent work proposes an
autoencoder-paradigm by leveraging a physics-informed neural network to solve
this problem (Zhao et al., 2022), which is shown in Figure 1(a). It has two
components: a pure data-driven component that uses a Multilayer Perceptron
(MLP) neural network to encode input observations $X$ (e.g., inductor current
and capacitor voltage) into a latent representation $H$, and a physics
component that uses physical laws as the decoder $PI$ to decode the latent
representation into the evolved observations $\bar{X}$. Once the PINN model
converges, the internal parameters $\lambda$ of PECs can be estimated through
the learned parameters of the decoder part. More specifically, given the input
$X\in R^{F}$, where $F$ is the feature dimension, the PINN model first encode
it to a latent representation $H$ by using an MLP,

(1) $\displaystyle H=MLP(X;\phi),$

where $H\in\mathbb{R}^{q}$, $q$ is the dimension of the latent representation,
and $\phi$ is the parameters of the MLP network. Next, the PINN model can
decode the latent representation $H$ in the PEC system to an output $\bar{X}$
by using physical laws,

(2) $\displaystyle\bar{X}=PI(H;\lambda)$

where $PI(\cdot)$ is a physical function with the learned internal parameters
$\lambda$ based on the physical law in a PEC system. The optimization
objective is to minimize $|X-\bar{X}|$, formally:

(3)
$\displaystyle\lambda^{\ast}=argmin_{\lambda\in\mathbb{R}^{k},\phi}(|X-PI(MLP(X);\lambda)|)$

where $k$ is the number of parameters of the physical function $PI(\cdot)$,
and $\lambda^{\ast}$ is the parameters that lead to the minimal $|X-\bar{X}|$.
The PINN model is trained using backpropagation on the Mean Absolute Error
(MAE) loss, and the parameters $\phi$ and $\lambda$ are optimized jointly.

In the task of estimating the internal parameters of PECs, the input $X$ is
the vector composed of the inductor current and the capacitor voltage. And the
parameters of the physical function $PI(\cdot)$ are also the internal
parameters of PECs, that is, $\lambda$. Therefore, our goal is to obtain the
learned $\lambda$ after the PINN model is trained.

###  2.2. Automatically designed PINNs

A drawback of the PINN model in Section 2.1 is that we need to manually design
the architecture of the neural network. Due to a large number of choices, it
is often challenging in practice to find an appropriate physics-informed
neural network through manual network designing. In addition, real-world
industrial applications encounter different hardware constraints under
different scenarios, bringing more challenges to designing lightweight PINN
models to satisfy different memory and computational constraints.

To tackle this problem, we propose to combine AutoML methods and the PINN
model, and aim to automatically design the best PINN model under given
hardware constraints. In this section, we first introduce the search space
tailored for estimating the internal parameters of PECs, which contains
massive candidate PINN models. And then introduce how we use the reinforcement
learning based search strategy to explore the search space. Finally, we
introduce the hardware-aware search objective, which supports to search for
the best PINN model that satisfies given hardware constraints.

####  2.2.1. Search Space

We adopt MLP to implement the encoder in the PINN model and search for the
number of units $n_{i}$ and the activation function $a_{i}$ at the $i$-th
layer, where $i\leq T$ and $T$ is the maximum number of layers because
together they determine the architecture of the MLP and can lead to completely
different performance and memory consumption. The number of units can be
chosen from [0, 20, 30, 40, 50, 60], where $0$ means the current layer is
discarded. And the activation function can be chosen from [tanh, ReLu].
Therefore, for the $i$-th layer, each element of the possible choice
$(n_{i},a_{i})$ is chosen from the two lists above, respectively. For example,
(20, tanh) is a possible choice for a layer of the MLP, which means the layer
consists of 20 units and the activation function of this layer is tanh.

Overall, once the number of units and the activation function at each layer
are determined, a PINN model can be built.

####  2.2.2. Search Strategy

To explore the proposed search space, we apply a reinforcement learning based
search strategy. Figure 1(b) shows the overview of the search strategy.
Specifically, we first use a controller to generate the description of PINN
models. The controller is implemented by a Recurrent Neural Network (RNN),
which works in an autoregressive manner. At the $i$-th ($i\geq 1$) step, the
input of the RNN is the output from $(i-1)$-th step, except the first step,
whose input is a learned variable. The RNN encodes the input to a hidden state
and then inputs it to a softmax classifier to select the candidate PINN
models. At each step, the RNN outputs the number of units and the activation
function at each layer. After the RNN generates all the descriptions, a PINN
model is built based on it and trained to obtain the MAE.

Our goal is to learn the parameters of the RNN controller, denoted $\theta$,
to ensure that it can generate high-performance PINN models under given
hardware constraints. Therefore, we use the reciprocal of the MAE of the
generated PINN model as the reward to guide the optimization of the RNN.
Specifically, we maximize the following expected reward to learn the parameter
$\theta$ of the RNN controller,

(4) $\displaystyle J(\theta)=E_{P(c_{1:T};\theta)}[R]$

where $c_{1:T}$ is the output list of the RNN controller, and
$c_{1:T}=\\{(n_{1},a_{1}),(n_{2},a_{2}),(n_{3},a_{3}),(n_{4},a_{4}),(n_{5},a_{5})\\}$,
and $T=5$ in Figure 1(b). $P(c_{1:T};\theta)$ is the probability that the
controller generate the output $c_{1:T}$. $R=1/MAE$ is computed by training
the PINN model built based on $c_{1:T}$, and used as a reward to update the
parameter $\theta$ of the RNN controller. The objective encourages the RNN
controller to generate PINN models that yield small MAE. We follow Zoph et al.
(Zoph and Le, 2017) to use the REINFORCE rule (Williams, 1992) to compute the
gradient of $\theta$:

(5) $\displaystyle\bigtriangledown\theta
J(\theta)=\sum\limits_{t=1}^{T}\bigtriangledown_{\theta}logP(c_{t}|c_{(t-1):1};\theta)(1/MAE)$

####  2.2.3. Incorporating Hardware Constraints

Estimating the internal parameters of PECs is a multi-objective problem in
practice, as we usually need to design a PINN model with not only high
accuracy but also low memory usage. Given a storage constraint, our goal is to
the best PINN model under the constraint. To achieve this goal, we modify the
search objective in Equation 4 so that the search is hardware-aware, formally,

(6) $\displaystyle
J(\theta)=E_{P(c_{1:T};\theta)}(1/MAE)([\frac{Param(m)}{P_{0}}]^{w})$

where $Param(\cdot)$ is the function that counts the number of parameters of a
PINN model, which is used to approximate the memory usage. $P_{0}$ is the
given memory constraint, and $w$ is the weighting factor defined as:

$w=\begin{cases}\alpha,&if\ Param(m)\leq P_{0}\\\ \beta,&otherwise\end{cases}$

where $\alpha$ and $\beta$ are predefined constants. $m$ represents a PINN
model. In this work, we empirically set $\alpha=-0.02$ and $\beta=-0.1$. The
intuition behind this setting is that when the number of parameters of the
generated PINN model is less than $P0$, we increase $w$ because in this case,
we are more concerned with the accuracy of the searched PINN model.
Conversely, when the number of parameters of the generated PINN model is
larger than $P0$, we penalize the objective value with a smaller weight to
encourage the searched model to satisfy the given hardware constraint.

Table 1. Comparison Results w.r.t. MAE and # Param, which is the number of NN parameters. Comparison Models | $Average\>MAE$ | # Param | $L$ | $R_{L}$ | $C$ | $R_{C}$ | $R_{dson}$ | $R_{1}$ | $R_{2}$ | $R_{3}$ | $V_{in}$ | $V_{F}$
---|---|---|---|---|---|---|---|---|---|---|---|---
_M-PINN_ | 5.0 | 12,350 | 0.8 | 13.1 | 1.2 | 4.5 | 27.9 | 0.1 | 0.3 | 0.1 | 0.1 | 1.9
_R-PINN-1_ | 7.4 | 13,220 | 0.1 | 12.8 | 5.3 | 4.6 | 29.2 | 0.1 | 0.2 | 0.1 | 0.3 | 21.9
_R-PINN-2_ | 8.7 | 11,880 | 0.3 | 12.4 | 8.8 | 4.6 | 29.4 | 0.0 | 0.0 | 0.2 | 0.5 | 31.2
_R-PINN-3_ | 8.8 | 7,620 | 0.4 | 12.4 | 6.2 | 4.6 | 29.7 | 0.0 | 0.2 | 0.1 | 0.5 | 33.4
_AutoPINN-1_ | 4.7 | 15,480 | 0.9 | 13.0 | 0.8 | 4.4 | 27.5 | 0.1 | 0.2 | 0.2 | 0.1 | 0.2
_AutoPINN-2_ | 4.9 | 11,220 | 0.9 | 13.4 | 0.6 | 4.4 | 28.2 | 0.1 | 0.2 | 0.2 | 0.1 | 0.1
_AutoPINN-3_ | 5.0 | 9,220 | 0.8 | 12.9 | 1.5 | 4.4 | 27.4 | 0.1 | 0.2 | 0.1 | 0.0 | 2.9

##  3\. Experiments

###  3.1. Experimental Settings

Dataset.  To demonstrate the superiority of our proposed AutoPINN, we conduct
comparison experiments on a public PEC signal dataset (Zhao et al., 2022). The
data is collected by measuring the transient signals of the inductor current
$i$ and output voltage $u$ of a PEC. Only peaks of the inductor current and
output voltage are used. Therefore, the dataset $D$ contains a series of pairs
of $i$ and $u$: $D=\\{X_{k}\\}$, where $X_{k}=[i_{k},u_{k}]$. The input into
the achieved models by AutoPINN and comparison models is $X$, i.e., one pair
of $i$ and $u$. Overall, the dataset contains 360 pairs. For each pair
$X_{k}$, there is a corresponding set of internal parameters $\lambda$, which
is measured under laboratory conditions and used as the ground truth for
estimating the accuracy of all PINN models. Specifically,
$\lambda=\\{L,\,R_{L},\,C,\,R_{C},\,R_{dson},\,R_{1},\,R_{2},\,R_{3},\,V_{in},\,V_{F}\\}$,
where $L$ is inductance, $R_{L}$ is resistance of inductor, $C$ is
capacitance, $R_{C}$ is equivalent series resistance, $R_{dson}$ is onstate
resistance, $R_{1}$, $R_{2}$ and $R_{3}$ are three different loads, $V_{in}$
is input voltage, $V_{F}$ is diode forward voltage, respectively. Following
(Zhao et al., 2022), we use Mean Absolute Error (MAE) to evaluate the accuracy
of PINN models, where lower MAE values indicate higher accuracy.

Baselines.  We compare the AutoPINN with 1) the manually designed PINN model
proposed by Zhao et al. (Zhao et al., 2022), denoted as M-PINN, and 2) three
randomly selected PINN models from the search space, denoted as R-PINN-1,
R-PINN-2, and R-PINN-3, respectively. We reproduce M-PINN based on the source
code released by (Zhao et al., 2022). Note that M-PINN is also embraced in the
proposed search space. We run AutoPINN three times with different memory
constraints to search for three specific PINN models, denoted as AutoPINN-1,
AutoPINN-2, and AutoPINN-3, respectively.

Implementation details.  All experiments are conducted on an Nvidia Quadro RTX
8000 GPU. We adopt Adam (Kingma and Ba, 2014) followed by a full-batch L-BFGS
as the optimizer, with the learning rate $\eta$ = 0.001, and exponential decay
rates $\beta_{1}$ = 0.9 and $\beta_{2}$ = 0.999. The maximum number of MLP
layers $T$ is empirically set to 5. We train all PINN models for 2,000 epochs
to convergence and derive the unobservable parameters $\lambda$ from the well-
trained PINNs with the help of physical laws. For AutoPINN-1, AutoPINN-2, and
AutoPINN-3, we set their constraints on parameter numbers to 16,000, 12,000,
and 10,000, respectively.

####  3.1.1. Experimental Results

Table 1 summarizes the MAE performance and model size w.r.t number of
parameters of AutoPINN and all baselines. We use the bold font to highlight
the best-performed model in each column. Key observations are as follows.

First, the best-performed automatically designed PINN, AutoPINN-1, outperforms
the manually designed M-PINN on the average MAE (4.7 vs 5.0). When looking
into each parameter individually, we find that AutoPINN-1 performs better or
the same as M-PINN on eight out of ten parameters; for the remaining two
parameters (i.e., $L$ and $R_{3}$), AutoPINN-1 is only slightly worse than
M-PINN ($L$: 0.9 vs 0.8; $R_{3}$: 0.2 vs 0.1). In addition, the second-best
AutoPINN-2 not only has a lower average MAE than M-PINN (4.9 vs 5.0) but also
uses fewer model parameters (11,220 vs 12,350). This demonstrates that
AutoPINN can find better PINN models than manually designed ones with regard
to both accuracy and model size.

Second, when comparing the three automatically generated models, AutoPINN-1
outperforms AutoPINN-2 significantly but also contains extensively more
parameters. On the other hand, AutoPINN-3 performs slightly worse than
AutoPINN-2 but contains much fewer parameters. Moreover, AutoPINN-3 performs
the same as M-PINN regarding the average MAE. This demonstrates that AutoPINN
can effectively find well-performed but small PINN models; as the model size
increases, the PINN’s performance climbs as well.

Third, it is obvious that the randomly sampled PINN models perform extremely
poorly. For example, with much more parameters, R-PINN-1 performs far worse
than AutoPINN-3. This observation indicates that the proposed search strategy
is able to find high-accuracy PINN models from the search space.

##  4\. Related Work

###  4.1. Physics-informed Machine Learning

In recent years, there is a large body of work combining neural networks and
physical knowledge to solve scientific problems. A line of work introduces
physical knowledge into the traditional neural network by adding soft
constraints to the loss function (Krishnapriyan et al., 2021; Chen et al.,
2020; Geneva and Zabaras, 2020; Jin et al., 2021; Sahli Costabal et al.,
2020). However, this kind of methods cannot ensure the model strictly follow
the physical regulations, which may lead to unreasonable prediction results.
On the contrary, several existing works (Greydanus et al., 2019; Wang et al.,
2020; Beucler et al., 2019; Frerix et al., 2020; Cranmer et al., 2020) embed
specialized physical laws as hard constraints into neural networks, and
enforce the model to must satisfy physical constraints. Parameter estimation
is another application of PINN (Zhao et al., 2022; Regazzoni et al., 2021;
Lakshminarayana et al., 2022; Stiasny et al., 2021), which can conduct
parameter inference from the observations by leveraging linear or nonlinear
differential equations. The proposed automated framework is suitable for all
PINN models, and we introduce how we combine AutoML and PINN for parameter
estimation in this paper.

###  4.2. Automated Machine Learning

There is massive related work on automated machine learning, here we only
briefly discuss the most relevant work for this paper. Zoph and Le (2017)
propose to use an RNN-based controller to produce the description of a neural
network and train the controller with reinforcement learning to encourage it
to generate high-performance neural networks. Unlike this work, although we
also adopt the reinforcement learning based AutoML, we also consider the model
size of the produced neural network. Tan and Le (2019) search for small and
fast neural networks that can be deployed on mobile devices, and they
explicitly incorporate model latency into the search objective. In light of
the success of AutoML in various areas, we adopt this thrilling technique to
alleviate the costly human efforts in designing PINN models. To the best of
our knowledge, this paper is the first work to apply AutoML to Physics-
informed Machine Learning models.

##  5\. Conclusion

We present the AutoPINN to automatically design PINN models for estimating the
key but unobservable parameters of power electronic converters. In particular,
we first design a search space that contains massive candidate PINN models.
Then we apply a reinforcement learning based search strategy to explore the
search space to find the best PINN model. Furthermore, we introduce the model
size of the PINN model in the search objective to find the best PINN model
under given hardware constraints. The comparison between the proposed method
and baseline methods demonstrates the effectiveness of AutoPINN.

## References

  * (1)
  * Beucler et al. (2019) Tom Beucler, Stephan Rasp, Michael Pritchard, and Pierre Gentine. 2019.  Achieving conservation of energy in neural network emulators for climate modeling.  _arXiv preprint arXiv:1906.06622_ (2019).
  * Cai et al. (2022) Shengze Cai, Zhiping Mao, Zhicheng Wang, Minglang Yin, and George Em Karniadakis. 2022.  Physics-informed neural networks (PINNs) for fluid mechanics: A review.  _Acta Mechanica Sinica_ (2022), 1–12.
  * Cai et al. (2021) Shengze Cai, Zhicheng Wang, Sifan Wang, Paris Perdikaris, and George Em Karniadakis. 2021\.  Physics-informed neural networks for heat transfer problems.  _Journal of Heat Transfer_ 143, 6 (2021).
  * Chen et al. (2020) Yuyao Chen, Lu Lu, George Em Karniadakis, and Luca Dal Negro. 2020.  Physics-informed neural networks for inverse problems in nano-optics and metamaterials.  _Optics express_ 28, 8 (2020), 11618–11633.
  * Cranmer et al. (2020) Miles Cranmer, Sam Greydanus, Stephan Hoyer, Peter Battaglia, David Spergel, and Shirley Ho. 2020.  Lagrangian neural networks.  _arXiv preprint arXiv:2003.04630_ (2020).
  * Frerix et al. (2020) Thomas Frerix, Matthias Nießner, and Daniel Cremers. 2020\.  Homogeneous linear inequality constraints for neural network activations. In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition Workshops_. 748–749.
  * Geneva and Zabaras (2020) Nicholas Geneva and Nicholas Zabaras. 2020.  Modeling the dynamics of PDE systems with physics-constrained deep auto-regressive networks.  _J. Comput. Phys._ 403 (2020), 109056.
  * Greydanus et al. (2019) Samuel Greydanus, Misko Dzamba, and Jason Yosinski. 2019\.  Hamiltonian neural networks.  _Advances in Neural Information Processing Systems_ 32 (2019).
  * Ji et al. (2021) Weiqi Ji, Weilun Qiu, Zhiyu Shi, Shaowu Pan, and Sili Deng. 2021.  Stiff-pinn: Physics-informed neural network for stiff chemical kinetics.  _The Journal of Physical Chemistry A_ 125, 36 (2021), 8098–8106.
  * Jin et al. (2021) Xiaowei Jin, Shengze Cai, Hui Li, and George Em Karniadakis. 2021\.  NSFnets (Navier-Stokes flow nets): Physics-informed neural networks for the incompressible Navier-Stokes equations.  _J. Comput. Phys._ 426 (2021), 109951.
  * Karniadakis et al. (2021) George Em Karniadakis, Ioannis G Kevrekidis, Lu Lu, Paris Perdikaris, Sifan Wang, and Liu Yang. 2021.  Physics-informed machine learning.  _Nature Reviews Physics_ 3, 6 (2021), 422–440.
  * Kingma and Ba (2014) Diederik P Kingma and Jimmy Ba. 2014.  Adam: A method for stochastic optimization.  _arXiv preprint arXiv:1412.6980_ (2014).
  * Krishnapriyan et al. (2021) Aditi Krishnapriyan, Amir Gholami, Shandian Zhe, Robert Kirby, and Michael W Mahoney. 2021.  Characterizing possible failure modes in physics-informed neural networks.  _Advances in Neural Information Processing Systems_ 34 (2021).
  * Lakshminarayana et al. (2022) Subhash Lakshminarayana, Saurav Sthapit, and Carsten Maple. 2022\.  Application of Physics-Informed Machine Learning Techniques for Power Grid Parameter Estimation.  _Sustainability_ 14, 4 (2022), 2051.
  * Liu et al. (2018) Hanxiao Liu, Karen Simonyan, and Yiming Yang. 2018\.  DARTS: Differentiable Architecture Search. In _ICLR_.
  * Mao et al. (2020) Zhiping Mao, Ameya D Jagtap, and George Em Karniadakis. 2020\.  Physics-informed neural networks for high-speed flows.  _Computer Methods in Applied Mechanics and Engineering_ 360 (2020), 112789\.
  * Misyris et al. (2020) George S Misyris, Andreas Venzke, and Spyros Chatzivasileiadis. 2020\.  Physics-informed neural networks for power systems. In _2020 IEEE Power & Energy Society General Meeting (PESGM)_. IEEE, 1–5.
  * Pham et al. (2018) Hieu Pham, Melody Guan, Barret Zoph, Quoc Le, and Jeff Dean. 2018.  Efficient neural architecture search via parameters sharing. In _ICML_. 4095–4104.
  * Raissi et al. (2019) Maziar Raissi, Paris Perdikaris, and George E Karniadakis. 2019\.  Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations.  _Journal of Computational physics_ 378 (2019), 686–707.
  * Regazzoni et al. (2021) Francesco Regazzoni, Stefano Pagani, Cosenza Alessandro, Lombardi Alessandro, and Alfio Maria Quarteroni. 2021.  A physics-informed multi-fidelity approach for the estimation of differential equations parameters in low-data or large-noise regimes.  (2021).
  * Sahli Costabal et al. (2020) Francisco Sahli Costabal, Yibo Yang, Paris Perdikaris, Daniel E Hurtado, and Ellen Kuhl. 2020\.  Physics-informed neural networks for cardiac activation mapping.  _Frontiers in Physics_ 8 (2020), 42.
  * Stiasny et al. (2021) Jochen Stiasny, George S Misyris, and Spyros Chatzivasileiadis. 2021\.  Physics-informed neural networks for non-linear system identification for power system dynamics. In _2021 IEEE Madrid PowerTech_. IEEE, 1–6.
  * Tan and Le (2019) Mingxing Tan and Quoc V. Le. 2019.  EfficientNet: Rethinking Model Scaling for Convolutional Neural Networks. In _International Conference on Machine Learning 2019_. 6105–6114.
  * Wang et al. (2020) Bo Wang, Wenzhong Zhang, and Wei Cai. 2020.  Multi-scale deep neural network (mscalednn) methods for oscillatory stokes flows in complex domains.  _arXiv preprint arXiv:2009.12729_ (2020).
  * Williams (1992) Ronald J Williams. 1992\.  Simple statistical gradient-following algorithms for connectionist reinforcement learning.  _Machine learning_ 8, 3 (1992), 229–256.
  * Wu et al. (2021) Xinle Wu, Dalin Zhang, Chenjuan Guo, Chaoyang He, Bin Yang, and Christian S Jensen. 2021\.  AutoCTS: Automated correlated time series forecasting.  _Proceedings of the VLDB Endowment_ 15, 4 (2021), 971–983.
  * Wu et al. (2022) Xinle Wu, Dalin Zhang, Miao Zhang, Chenjuan Guo, Bin Yang, and Christian S Jensen. 2022\.  Joint Neural Architecture and Hyperparameter Search for Correlated Time Series Forecasting.  _arXiv preprint arXiv:2211.16126_ (2022).
  * Xu et al. (2019) Yuhui Xu, Lingxi Xie, Xiaopeng Zhang, Xin Chen, Guo-Jun Qi, Qi Tian, and Hongkai Xiong. 2019.  PC-DARTS: Partial Channel Connections for Memory-Efficient Architecture Search. In _ICLR_.
  * Zhao et al. (2022) Shuai Zhao, Yingzhou Peng, Yi Zhang, and Huai Wang. 2022\.  Parameter Estimation of Power Electronic Converters with Physics-informed Machine Learning.  _IEEE Transactions on Power Electronics_ (2022).
  * Zoph and Le (2017) Barret Zoph and Quoc V. Le. 2017.  Neural Architecture Search with Reinforcement Learning. In _International Conference on Learning Representation 2017_.