[current]500pt680pt

Polynomial algebras of superintegrable systems separating in Cartesian
coordinates from higher order ladder operators

Danilo Latini, Ian Marquette and Yao-Zhong Zhang

School of Mathematics and Physics, The University of Queensland, Brisbane, QLD
4072, Australia

<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>

abstract

###### Abstract

We introduce the general polynomial algebras characterizing a class of higher
order superintegrable systems that separate in Cartesian coordinates. The
construction relies on underlying polynomial Heisenberg algebras and their
defining higher order ladder operators. One feature of these algebras is that
they preserve by construction some aspects of the structure of the
$\mathfrak{gl}(n)$ Lie algebra. Among the classes of Hamiltonians arising in
this framework are various deformations of harmonic oscillator and singular
oscillator related to exceptional orthogonal polynomials and even Painlevé and
higher order Painlevé analogs. As an explicit example, we investigate a new
three-dimensional superintegrable system related to Hermite exceptional
orthogonal polynomials of type III. Among the main results is the
determination of the degeneracies of the model in terms of the finite-
dimensional irreducible representations of the polynomial algebra.

##  1 Introduction

Over the past years, it has been demonstrated how the study of superintegrable
systems is intimately connected with the study of various algebraic structures
beyond Lie algebras [1, 2, 3, 4], which take the form of quadratic [5, 6, 7,
8] and more generally polynomial algebras [9, 10, 11, 12, 13, 14]. This is due
to the fact that integrals are partial differential operators of second or
higher order.

Recently it has been shown that quadratic algebras associated with
$n$-dimensional systems are in general of higher rank [15, 16, 17, 18]. These
algebraic structures allow one to obtain useful information on quantum systems
and their degenerate spectrum [19]. In most cases, they do not display an
obvious basis and thus the construction of their representations is difficult.

Similarly to the two-dimensional anisotropic oscillator [20] and related
symmetry algebra, many papers were devoted to two-dimensional superintegrable
systems [21, 22, 23, 24] with deformations for which the wavefunctions were
written in terms of exceptional orthogonal polynomials [25, 26, 27, 28, 29,
30, 31, 32, 33, 34, 35, 36, 37, 38]. Since of their first appearance,
exceptional orthogonal polynomials have found application in many different
physical contexts [39, 40, 41, 42, 43, 44, 45]. One class of Hamiltonians with
richer structures such as different “regimes”for degeneracies were discovered,
leading to the path of different types of spectral design for the quantum
systems. It was found how polynomial algebras of arbitrary degree with three
generators allow complicated patterns for the irreducible representations
(irreps) connected with the degeneracies [22, 23]. These works also allowed
one to describe an important feature of the irreps of the polynomial algebras
for which, so far, only very limited results have been obtained.

This paper is devoted to introduce the general polynomial algebras for models
allowing separation of variables in Cartesian coordinates whose ladder
operators can be of arbitrary order. These polynomial algebras, whose
structures share some similarities with the well-known $\mathfrak{gl}(n)$ Lie
algebra, are very interesting from the point of view of representations, as
they allow one to rely on analog of using Cartan-like and ladder-like
generators.

We want to address that these quantum models with underlying polynomial
Heisenberg algebras [46, 47], even though they possess separation of variables
in Cartesian coordinates, have in general complicated expressions for their
wavefunctions beyond hypergeometric type of orthogonal polynomials. More
importantly, the patterns for their degenerate spectrum are non-trivial and
the separated equations do not provide information on the structures of
multiplets. In this paper, we obtain interesting patterns for finite-
dimensional unitary representations of the polynomial algebra and use them to
provide further insights into the degeneracies of the models. In particular,
for the system associated with type III Hermite exceptional orthogonal
polynomials, we find that the representations of the polynomial algebra
display similarities with the $\mathfrak{su}(3)$ representations, which we
will refer to as $\mathfrak{su}(3)$-like. We obtain direct sums of such
representations and the dimensions and multiplicities of the representations.

The paper is organized as follows. In Section 2 we construct the higher order
polynomial algebra characterizing $n$D superintegrable Hamiltonian systems
separating in Cartesian coordinates and whose constituents ladder operators,
together with their associated one-dimensional Hamiltonians, define polynomial
Heisenberg algebras. In Section 3 we restrict to the three-dimensional case.
Specifically, we study a new three-dimensional superintegrable Hamiltonian
system that saparates in Cartesian coordinates $(x_{1},x_{2},x_{3})$. This
model, that we introduce in Section 3.1, due to the presence of a rationally
extended oscillator in the $x_{1}$ axis, turns out to be characterized by a
separable solution arising as the product of two standard Hermite orthogonal
polynomials (in the $x_{2},x_{3}$ axes) and Hermite exceptional orthogonal
polynomials of type III (in the $x_{1}$ axis). As a direct consequence, the
model presents complicated patterns for the degeneracies of the energy levels,
that we aim to characterize in terms of finite-dimensional irreps of the
polynomial algebra introduced in the previous Section. Finally, Section 4 is
devoted to some concluding remarks.

##  2 From polynomial Heisenberg algebras to polynomial symmetry algebras

Let us consider a $n$-dimensional Hamiltonian system admitting separation of
variables in Cartesian coordinates:

$\hat{H}=\sum_{i=1}^{n}\hat{H}_{i}=\sum_{i=1}^{n}\bigl{(}-\partial_{x_{i}x_{i}}+V(x_{i})\bigl{)}$
(2.1)

and let us assume the existence of $n$ pairs of ladder operators
$(\hat{A}_{i},\hat{A}_{i}^{\dagger})$ of order $k_{i}$, $i=1,\dots,n$,
satisfying the following Polynomial Heisenberg Algebras (PHA) [46, 47]:

$\displaystyle[\hat{H}_{i},\hat{A}_{i}]$
$\displaystyle=-\alpha_{i}\hat{A}_{i}\qquad[\hat{H}_{i},\hat{A}^{\dagger}_{i}]=\alpha_{i}\hat{A}^{\dagger}_{i}$
(2.2) $\displaystyle\hat{A}_{i}^{\dagger}\hat{A}_{i}$
$\displaystyle=\phi_{i}(\hat{H}_{i})\qquad\quad\hat{A}_{i}\hat{A}_{i}^{\dagger}=\phi_{i}(\hat{H}_{i}+\alpha_{i})\,,$
(2.3)

where $\alpha_{i}$ are constants and $\phi_{i}(\hat{H}_{i})$ are $n$
polynomials of order $k_{i}$ in $\hat{H}_{i}$ . Clearly, by virtue of (2.3),
we can rewrite the commutators among the ladder operators as:

$[\hat{A}_{i},\hat{A}_{i}^{\dagger}]=\phi_{i}(\hat{H}_{i}+\alpha_{i})-\phi_{i}(\hat{H}_{i})\equiv
P_{i}(\hat{H}_{i})\,,$ (2.4)

where $P_{i}(\hat{H}_{i})$ are $n$ polynomials of order $k_{i}-1$. We also
notice that (2.3) implies:

$\phi_{i}(\hat{H}_{i}+\alpha_{i})\hat{A}_{i}=\hat{A}_{i}\phi_{i}(\hat{H}_{i})\,,\qquad\phi_{i}(\hat{H}_{i})\hat{A}^{\dagger}_{i}=\hat{A}^{\dagger}_{i}\phi_{i}(\hat{H}_{i}+\alpha_{i})\,.$
(2.5)

Moreover, the following additional relations hold:

$\displaystyle[\hat{H_{i}},(\hat{A}_{i}^{\dagger})^{s_{i}}]$
$\displaystyle=\alpha_{i}s_{i}(\hat{A}_{i}^{\dagger})^{s_{i}}\,,\qquad[\hat{H_{i}},(\hat{A}_{i})^{s_{i}}]=-\alpha_{i}{s_{i}}(\hat{A}_{i})^{s_{i}}$
(2.6) $\displaystyle(\hat{A}_{i}^{\dagger})^{s_{i}}(\hat{A}_{i})^{s_{i}}$
$\displaystyle=\Phi_{i}(s_{i},\hat{H}_{i})\,,\quad(\hat{A}_{i})^{s_{i}}(\hat{A}_{i}^{\dagger})^{s_{i}}=\Phi_{i}({s_{i}},\hat{H}_{i}+\alpha_{i}s_{i})\,,$
(2.7)

for $s_{i}\in\mathds{N}^{*}$, where we defined:

$\Phi_{i}(s_{i},x_{i}):=\prod_{j=1}^{s_{i}}\phi_{i}\bigl{(}x_{i}-\alpha_{i}(s_{i}-j)\bigl{)}\,.$
(2.8)

Proofs can be obtained by induction. For example, if we take the first
equation in (2.7), then for $s_{i}=1$ is trivially satisfied as the product
collapses to $\phi_{i}(x_{i})$. Now, suppose that for $s_{i}-1\geq 1$ is
satisfied, then:

$\displaystyle(\hat{A}_{i}^{\dagger})^{s_{i}}(\hat{A}_{i})^{s_{i}}=\hat{A}_{i}^{\dagger}(\hat{A}_{i}^{\dagger})^{s_{i}-1}(\hat{A}_{i})^{s_{i}-1}\hat{A}_{i}=A^{\dagger}_{i}\biggl{(}\prod_{j=1}^{s_{i}-1}\phi_{i}\bigl{(}x_{i}-\alpha_{i}(s_{i}-j-1)\bigl{)}\biggl{)}A_{i}$
$\displaystyle\overset{\eqref{eq:imp}}{=}\biggl{(}\prod_{j=1}^{s_{i}-1}\phi_{i}\bigl{(}x_{i}-\alpha_{i}(s_{i}-j)\bigl{)}\biggl{)}A^{\dagger}_{i}A_{i}$
$\displaystyle\overset{\eqref{eq:ladderop2}}{=}\biggl{(}\prod_{j=1}^{s_{i}-1}\phi_{i}\bigl{(}x_{i}-\alpha_{i}(s_{i}-j)\bigl{)}\biggl{)}\phi_{i}(x_{i})$
$\displaystyle\,\,=\prod_{j=1}^{s_{i}}\phi_{i}\bigl{(}x_{i}-\alpha_{i}(s_{i}-j)\bigl{)}\,,$

which is the desired result. The other relation can be obtained in a similar
way. We notice that an equivalent form for the product (2.8) also exists, as
the following equality holds:

$\prod_{j=1}^{s_{i}}\phi_{i}\bigl{(}x_{i}-\alpha_{i}(s_{i}-j)\bigl{)}=\prod_{j=1}^{s_{i}}\phi_{i}\bigl{(}x_{i}-\alpha_{i}(j-1)\bigl{)}\,.$
(2.9)

This can be shown by expanding the two products:

$\phi_{i}(x_{i}-\alpha_{i}(s_{i}-1))\phi_{i}(x_{i}-\alpha_{i}(s_{i}-2))\dots\phi_{i}(x_{i}-\alpha_{i})\phi_{i}(x_{i})=\phi_{i}(x_{i})\phi_{i}(x_{i}-\alpha_{i})\dots\phi_{i}(x_{i}-\alpha_{i}(s_{i}-1))\,.$
(2.10)

At this point, let us introduce the following $n(n-1)$ operators of order
$s_{i}k_{i}+s_{j}k_{j}$:

$\displaystyle\hat{E}_{ij}$
$\displaystyle:=(\hat{A}^{\dagger}_{i})^{s_{i}}(\hat{A}_{j})^{s_{j}}\qquad(i\neq
j)\,.$ (2.11)

If we require that they commute with the Hamiltonian (2.1), we get:

$[\hat{H},\hat{E}_{ij}]=(\alpha_{i}s_{i}-\alpha_{j}s_{j})\hat{E}_{ij}\,.$
(2.12)

This leads us to the condition: $\alpha_{i}s_{i}=\alpha_{j}s_{j}$, with
$\alpha_{1}s_{1}=\alpha_{2}s_{2}=\dots=\alpha_{n}s_{n}=:\alpha$. Together with
the above operators, we consider the one-dimensional Hamiltonians
$\hat{H}_{i}$ ($i=1,\dots,n$). The set composed by the $n^{2}$ operators
$\\{\hat{E}_{ij},\hat{H}_{i}\\}$ close in the following higher rank polynomial
algebra:

$\displaystyle[\hat{E}_{ij},\hat{E}_{k\ell}]$
$\displaystyle=\delta_{jk}\bigl{(}\Phi_{k}({s_{k}},\hat{H}_{k}+\alpha)-\Phi_{i}({s_{k}},\hat{H}_{k})\bigl{)}\hat{E}_{i\ell}-\delta_{i\ell}\bigl{(}\Phi_{\ell}({s_{\ell}},\hat{H}_{i}+\alpha)-\Phi_{i}({s_{\ell}},\hat{H}_{\ell})\bigl{)}\hat{E}_{kj}$
(2.13) $\displaystyle[\hat{E}_{ij},\hat{E}_{ji}]$
$\displaystyle=\Phi_{i}(s_{i},\hat{H}_{i})\Phi_{j}(s_{j},\hat{H}_{j}+\alpha)-\Phi_{i}(s_{i},\hat{H}_{i}+\alpha)\Phi_{j}(s_{j},\hat{H}_{j})$
(2.14) $\displaystyle[\hat{E}_{ij},\hat{H}_{k}]$
$\displaystyle=\alpha(\delta_{jk}\hat{E}_{ik}-\delta_{ik}\hat{E}_{kj})\,,$
(2.15)

where the functions $\Phi(s,x)$ are the ones given in (2.8).

######  Remark 1.

When we restrict to the $n$-dimensional harmonic oscillator:

$\hat{H}=\frac{1}{2}\sum_{i=1}^{n}\bigl{(}-\partial_{x_{i}x_{i}}+x^{2}_{i}\bigl{)}\,,$
(2.16)

the ladder operators are of the first-order $($$k_{1}=k_{2}=\dots=k_{n}=1$$)$.
Explicitly:

$\quad\hat{A}_{i}=\frac{\partial_{x_{i}}+x_{i}}{\sqrt{2}}\,,\quad\hat{A}^{\dagger}_{i}=\frac{-\partial_{x_{i}}+x_{i}}{\sqrt{2}}$
(2.17)

and the following commutation relations are satisfied:

$\displaystyle[\hat{H}_{i},\hat{A}_{i}]$ $\displaystyle=-\hat{A}_{i}\hskip
42.67912pt[\hat{H}_{i},\hat{A}^{\dagger}_{i}]=\hat{A}^{\dagger}_{i}$ (2.18)
$\displaystyle\hat{A}_{i}^{\dagger}\hat{A}_{i}$
$\displaystyle=\hat{H}_{i}-1/2\qquad\quad\,\,\hat{A}_{i}\hat{A}_{i}^{\dagger}=\hat{H}_{i}+1/2\,.$
(2.19)

Thus, in this limiting case we have $s_{1}=s_{2}=\dots=s_{n}$, and for
$s_{1}=s_{2}=\dots=s_{n}=1$ we get:

$\displaystyle[\hat{E}_{ij},\hat{E}_{k\ell}]$
$\displaystyle=\delta_{jk}{E}_{i\ell}-\delta_{i\ell}\hat{E}_{kj}$ (2.20)
$\displaystyle[\hat{E}_{ij},\hat{E}_{ji}]$
$\displaystyle=\hat{H}_{i}-\hat{H}_{j}$ (2.21)
$\displaystyle[\hat{E}_{ij},\hat{H}_{k}]$
$\displaystyle=\delta_{jk}\hat{E}_{ik}-\delta_{ik}\hat{E}_{kj}\,,$ (2.22)

as $\Phi_{i}(1,\hat{H}_{i})\equiv\phi_{i}(\hat{H}_{i})=\hat{H}_{i}-1/2$ and
$\alpha=1\,,\forall\,i=1,\dots,n$. Notice also that if we rewrite
$\hat{H}_{k}\equiv\hat{E}_{kk}$ $($i.e. we allow equal indices for the
generators (2.11)$)$ we can recast the algebra in the compact form:

$\displaystyle[\hat{E}_{ij},\hat{E}_{k\ell}]$
$\displaystyle=\delta_{jk}\hat{E}_{i\ell}-\delta_{i\ell}\hat{E}_{kj}\quad(i,j=1,\dots,n)\,,$
(2.23)

which gives the usual realisation of the $\mathfrak{u}(n)$ Lie algebra for the
$n$-dimensional harmonic oscillator.

The $n=2$ case has been already studied in the literature (see [22, 23] and
references therein). The Hamiltonian is:

$\hat{H}=\hat{H}_{1}+\hat{H}_{2}$ (2.24)

and the four operators
$\\{\hat{H}_{1},\hat{H}_{2},\hat{E}_{12},\hat{E}_{21}\\}$ are all conserved
when $\alpha_{1}s_{1}=\alpha_{2}s_{2}=:\alpha$, with commutation relations:

$\displaystyle[\hat{E}_{12},\hat{E}_{21}]$
$\displaystyle=\Phi_{1}(s_{1},\hat{H}_{1})\Phi_{2}(s_{2},\hat{H}_{2}+\alpha)-\Phi_{1}(s_{1},\hat{H}_{1}+\alpha)\Phi_{2}(s_{2},\hat{H}_{2})$
(2.25) $\displaystyle[\hat{E}_{12},\hat{H}_{1}]$
$\displaystyle=-\alpha\hat{E}_{12}\qquad[\hat{E}_{12},\hat{H}_{2}]=\alpha\hat{E}_{12}$
(2.26) $\displaystyle[\hat{E}_{21},\hat{H}_{1}]$
$\displaystyle=\alpha\hat{E}_{21}\hskip
28.45274pt[\hat{E}_{21},\hat{H}_{2}]=-\alpha\hat{E}_{21}\,.$ (2.27)

By considering the following combinations of the generators:

$\hat{T}_{3}:=\frac{\hat{H}_{1}-\hat{H}_{2}}{2\alpha}\qquad\hat{T}_{+}:=\hat{E}_{12}\qquad\hat{T}_{-}:=\hat{E}_{21}\,,$
(2.28)

together with $[\hat{H},\cdot]=0$ we get:

$\displaystyle[\hat{T}_{+},\hat{T}_{-}]$
$\displaystyle=\Phi_{1}(s_{1},\hat{H}/2+\alpha\hat{T}_{3})\Phi_{2}(s_{2},\hat{H}/2-\alpha\hat{T}_{3}+\alpha)-\Phi_{1}(s_{1},\hat{H}/2+\alpha\hat{T}_{3}+\alpha)\Phi_{2}(s_{2},\hat{H}/2-\alpha\hat{T}_{3})$
(2.29) $\displaystyle[\hat{T}_{3},\hat{T}_{\pm}]$ $\displaystyle=\pm
T_{\pm}\,.$ (2.30)

In this way, we restrict to the set of operators
$\\{\hat{H},\hat{T}_{3},\hat{T}_{\pm}\\}$. At this point, by introducing the
new function:

$\Xi^{(s_{1},s_{2})}_{12}(\hat{H},\hat{T}_{3}):=\Phi_{1}(s_{1},\hat{H}/2+\alpha\hat{T}_{3})\Phi_{2}(s_{2},\hat{H}/2-\alpha\hat{T}_{3}+\alpha)$
(2.31)

we can recast the polynomial algebra in the following form:

$\displaystyle[\hat{H},\cdot]$ $\displaystyle=0$ (2.32)
$\displaystyle[\hat{T}_{+},\hat{T}_{-}]$
$\displaystyle=\Xi_{12}^{(s_{1},s_{2})}(\hat{H},\hat{T}_{3})-\Xi_{12}^{(s_{1},s_{2})}(\hat{H},\hat{T}_{3}+1)$
(2.33) $\displaystyle[\hat{T}_{3},\hat{T}_{\pm}]$
$\displaystyle=\pm\hat{T}_{\pm}\,.$ (2.34)

This polynomial algebra is of order $s_{1}+s_{2}-1$ in the generators
($k_{1}s_{1}+k_{2}s_{2}-1$ in terms of differential operators).

######  Remark 2.

When we restrict to the $2$D harmonic oscillator:

$s_{1}=s_{2}=\alpha=1$
$\Xi^{(s_{1},s_{2})}_{12}(\hat{H},\hat{T}_{3})-\Xi_{12}^{(s_{1},s_{2})}(\hat{H},\hat{T}_{3}+1)=2\hat{T}_{3}$

the polynomial algebra collapses to:

$\displaystyle[\hat{H},\cdot]=0\qquad[\hat{T}_{+},\hat{T}_{-}]=2\hat{T}_{3}\qquad[\hat{T}_{3},\hat{T}_{\pm}]=\pm\hat{T}_{\pm}\,,$

with $\hat{T}_{+}=\hat{A}_{1}^{\dagger}\hat{A}_{2}$ and
$\hat{T}_{-}=\hat{A}_{2}^{\dagger}\hat{A}_{1}$, where
$(\hat{A}_{i},\hat{A}_{i}^{\dagger})$ are the usual first-order ladder
operators defined in (2.17).

The first non-trivial example is represented by the $2$D anisotropic
oscillator, which has been investigated in [20] (see [2] for the general $n$D
case). The Hamiltonian and ladder operators read [20]:

$\hat{H}=\hat{H}_{1}+\hat{H}_{2}=\frac{1}{2}\bigl{(}-\partial_{x_{1}x_{1}}+x^{2}_{1}/m_{1}^{2}\bigl{)}+\frac{1}{2}\bigl{(}-\partial_{x_{2}x_{2}}+x^{2}_{2}/m_{2}^{2}\bigl{)}\,,\quad\hat{A}^{\dagger}_{i}=\frac{-\partial_{x_{i}}+x_{i}/m_{i}}{\sqrt{2}}\,,\quad\hat{A}_{i}=\frac{\partial_{x_{i}}+x_{i}/m_{i}}{\sqrt{2}}$
(2.35)

where $m_{1},m_{2}$ are mutually prime natural numbers. The above operators
close the commutation relations ($i=1,2$):

$\displaystyle[\hat{H}_{i},\hat{A}_{i}]$
$\displaystyle=-\hat{A}_{i}/m_{i}\hskip
46.94687pt[\hat{H}_{i},\hat{A}^{\dagger}_{i}]=\hat{A}^{\dagger}_{i}/m_{i}$
(2.36) $\displaystyle\hat{A}_{i}^{\dagger}\hat{A}_{i}$
$\displaystyle=\hat{H}_{i}-1/(2m_{i})\qquad\quad\,\,\,\hat{A}_{i}\hat{A}_{i}^{\dagger}=\hat{H}_{i}+1/(2m_{i})\,.$
(2.37)

The two operators $\hat{E}_{12}$, $\hat{E}_{21}$ commute with the Hamiltonian
$\hat{H}$ if $s_{1}/m_{1}=s_{2}/m_{2}=\alpha$, i.e. for $s_{i}=m_{i}$ and
$\alpha=1$. Thus, the generators of the polynomial algebra (of order
$m_{1}+m_{2}-1$) turn out to be:

$\hat{H}=\hat{H}_{1}+\hat{H}_{2}\,,\quad\hat{T}_{3}:=\frac{\hat{H}_{1}-\hat{H}_{2}}{2}\,,\quad\hat{T}_{+}:=\hat{E}_{12}=(A_{1}^{\dagger})^{m_{1}}(A_{2})^{m_{2}}\,,\quad\hat{T}_{-}:=\hat{E}_{21}=(A_{2}^{\dagger})^{m_{2}}(A_{1})^{m_{1}}\,,$
(2.38)

and the higher order polynomial relations (2.33)-(2.34) are obtained through
the following structure function:

$\small\Xi_{12}^{(m_{1},m_{2})}(\hat{H},\hat{T}_{3}):=\Phi_{1}(m_{1},\hat{H}/2+\hat{T}_{3})\Phi_{2}(m_{2},\hat{H}/2-\hat{T}_{3}+1)=\prod_{i=1}^{m_{1}}\biggl{(}\hat{H}/2+\hat{T}_{3}-\frac{2i-1}{2m_{1}}\biggl{)}\prod_{j=1}^{m_{2}}\biggl{(}\hat{H}/2-\hat{T}_{3}+\frac{2j-1}{2m_{2}}\biggl{)}\,.$
(2.39)

When $m_{1}=m_{2}=1$ the $\mathfrak{u}(2)$ algebra is recovered. Thus, the
polynomial algebra (2.32)-(2.34) characterized by the function (2.39)
represents a non-linear generalization of the $\mathfrak{u}(2)$ algebra. The
finite-dimensional irreps of this algebra have been constructed by using a
deformed oscillator algebra approach [3, 4, 5]. The same strategy has been in
used in [23] to obtain finite-dimensional irreps of the polynomial algebras
describing 2D superintegrable systems separable in Cartesian coordinates
related to a rational extension of the harmonic oscillator with type III
Hermite exceptional orthogonal polynomials. Many other 2D superintegrable
systems have been studied within this framework (see [10, 12, 48] and
references therein).

##  3 Three-dimensional superintegrable Hamiltonian systems

In this Section 3, we are interested in constructing polynomial algebras for
$3$D superintegrable systems separating in Cartesian coordinates starting from
polynomial Heisenberg algebras (2.2)-(2.3). The $3$D Hamiltonian reads:

$\hat{H}=\hat{H}_{1}+\hat{H}_{2}+\hat{H}_{3}$ (3.1)

and we can introduce the following set composed by the nine operators:

$\\{\hat{H}_{1},\hat{H}_{2},\hat{H}_{3},\hat{E}_{12},\hat{E}_{13},\hat{E}_{21},\hat{E}_{23},\hat{E}_{31},\hat{E}_{32}\\}\,,$
(3.2)

satisfying the condition
$\alpha_{1}s_{1}=\alpha_{2}s_{2}=\alpha_{3}s_{3}=\alpha$. The nonzero
commutation relations are obtained from (2.13)-(2.15) and, in complete analogy
with the $2$D case, we consider the following combinations:

$\displaystyle\hat{T}_{3}:=\frac{\hat{H}_{1}-\hat{H}_{2}}{2\alpha}\qquad\hat{Y}:=\frac{\hat{H}_{1}+\hat{H}_{2}-2\hat{H}_{3}}{3\alpha}\,,$
(3.3)
$\displaystyle\hat{T}_{+}:=\hat{E}_{12}\,,\quad\hat{T}_{-}:=\hat{E}_{21}\,,\quad\hat{U}_{+}:=\hat{E}_{23}\,,\quad\hat{U}_{-}:=\hat{E}_{32}\quad\hat{V}_{+}:=\hat{E}_{13}\,,\quad\hat{V}_{-}:=\hat{E}_{31}\,.$
(3.4)

In this way, we restrict to the set of operators
$\\{\hat{H},\hat{T}_{3},\hat{Y},\hat{T}_{\pm},\hat{U}_{\pm},\hat{V}_{\pm}\\}$.
Besides $[\hat{H},\cdot]=0$, these operators close under commutation to give
the polynomial algebra:

$[\hat{T}_{3},\hat{Y}]=0$ (3.5) $\displaystyle[\hat{T}_{+},\hat{T}_{-}]$
$\displaystyle=\Xi^{(s_{1},s_{2})}_{12}(\hat{H},\hat{T}_{3},\hat{Y})-\Xi_{12}^{(s_{1},s_{2})}(\hat{H},\hat{T}_{3}+1,\hat{Y})\hskip
56.9055pt[\hat{T}_{3},\hat{T}_{\pm}]=\pm\hat{T}_{\pm}\hskip
34.14322pt[\hat{Y},\hat{T}_{\pm}]=0$ (3.6)
$\displaystyle[\hat{V}_{+},\hat{V}_{-}]$
$\displaystyle=\Xi_{13}^{(s_{1},s_{3})}(\hat{H},\hat{T}_{3},\hat{Y})-\Xi_{13}^{(s_{1},s_{3})}(\hat{H},\hat{T}_{3}+1/2,\hat{Y}+1)\hskip
29.87547pt[\hat{T}_{3},\hat{V}_{\pm}]=\pm\frac{1}{2}\hat{V}_{\pm}\hskip
27.03003pt[\hat{Y},\hat{V}_{\pm}]=\pm\hat{V}_{\pm}$ (3.7)
$\displaystyle[\hat{U}_{+},\hat{U}_{-}]$
$\displaystyle=\Xi_{23}^{(s_{2},s_{3})}(\hat{H},\hat{T}_{3},\hat{Y})-\Xi_{23}^{(s_{2},s_{3})}(\hat{H},\hat{T}_{3}-1/2,\hat{Y}+1)\hskip
29.87547pt[\hat{T}_{3},\hat{U}_{\pm}]=\mp\frac{1}{2}\hat{U}_{\pm}\hskip
24.18501pt[\hat{Y},\hat{U}_{\pm}]=\pm\hat{U}_{\pm}\,$ (3.8)
$\displaystyle[\hat{T}_{\pm},\hat{V}_{\mp}]=\mp\bigl{(}\Phi_{1}(s_{1},\hat{H}/3+\alpha\hat{T}_{3}+\alpha\hat{Y}/2+\alpha)-\Phi_{1}(s_{1};\hat{H}/3+\alpha\hat{T}_{3}+\alpha\hat{Y}/2)\bigl{)}\hat{U}_{\mp}$
(3.9)
$\displaystyle[\hat{T}_{\pm},\hat{U}_{\pm}]=\pm\bigl{(}\Phi_{2}(s_{2},\hat{H}/3-\alpha\hat{T}_{3}+\alpha\hat{Y}/2+\alpha)-\Phi_{2}(s_{2},\hat{H}/3-\alpha\hat{T}_{3}+\alpha\hat{Y}/2)\bigl{)}\hat{V}_{\pm}$
(3.10)
$\displaystyle[\hat{U}_{\pm},\hat{V}_{\mp}]=\pm\bigl{(}\Phi_{3}(s_{3},\hat{H}/3-\alpha\hat{Y}+\alpha)-\Phi_{3}(s_{3},\hat{H}/3-\alpha\hat{Y})\bigl{)}\hat{T}_{\mp}\,,$
(3.11)

where we introduced the three functions:

$\displaystyle\Xi_{12}^{(s_{1},s_{2})}(\hat{H},\hat{T}_{3},\hat{Y})$
$\displaystyle:=\Phi_{1}(s_{1},\hat{H}/3+\alpha\hat{T}_{3}+\alpha\hat{Y}/2)\Phi_{2}(s_{2},\hat{H}/3-\alpha\hat{T}_{3}+\alpha\hat{Y}/2+\alpha)$
(3.12) $\displaystyle\Xi_{13}^{(s_{1},s_{3})}(\hat{H},\hat{T}_{3},\hat{Y})$
$\displaystyle:=\Phi_{1}(s_{1},\hat{H}/3+\alpha\hat{T}_{3}+\alpha\hat{Y}/2)\Phi_{3}(s_{3},\hat{H}/3-\alpha\hat{Y}+\alpha)$
(3.13) $\displaystyle\Xi_{23}^{(s_{2},s_{3})}(\hat{H},\hat{T}_{3},\hat{Y})$
$\displaystyle:=\Phi_{2}(s_{2},\hat{H}/3-\alpha\hat{T}_{3}+\alpha\hat{Y}/2)\Phi_{3}(s_{3},\hat{H}/3-\alpha\hat{Y}+\alpha)\,,$
(3.14)

and
$(\hat{T}_{-})^{\dagger}=\hat{T}_{+}\,,(\hat{U}_{-})^{\dagger}=\hat{U}_{+}\,,(\hat{V}_{-})^{\dagger}=\hat{V}_{+}$
whereas
$\hat{H}^{\dagger}=\hat{H}\,,\,\hat{T}_{3}^{\dagger}=\hat{T}_{3}\,,\hat{Y}^{\dagger}=\hat{Y}$.

The specific choice (3.3) is dictated by the following reasoning. Consider the
Hamiltonian of the isotropic harmonic oscillator (2.16) with $n=3$, i.e.
$s_{1}=s_{2}=s_{3}=\alpha=1$. In this limiting case, the polynomial algebra
collapses to:

$[\hat{H},\cdot]=[\hat{T}_{3},\hat{Y}]=0$ (3.15)
$\displaystyle[\hat{T}_{+},\hat{T}_{-}]$ $\displaystyle=2\hat{T}_{3}\hskip
56.9055pt[\hat{T}_{3},\hat{T}_{\pm}]=\pm\hat{T}_{\pm}\hskip
45.52458pt[\hat{Y},\hat{T}_{\pm}]=0$ (3.16)
$\displaystyle[\hat{V}_{+},\hat{V}_{-}]$
$\displaystyle=\hat{T}_{3}+\frac{3}{2}\hat{Y}\hskip
34.8548pt[\hat{T}_{3},\hat{V}_{\pm}]=\pm\frac{1}{2}\hat{V}_{\pm}\hskip
36.98866pt[\hat{Y},\hat{V}_{\pm}]=\pm\hat{V}_{\pm}$ (3.17)
$\displaystyle[\hat{U}_{+},\hat{U}_{-}]$
$\displaystyle=-\hat{T}_{3}+\frac{3}{2}\hat{Y}\hskip
27.03003pt[\hat{T}_{3},\hat{U}_{\pm}]=\mp\frac{1}{2}\hat{U}_{\pm}\hskip
34.14322pt[\hat{Y},\hat{U}_{\pm}]=\pm\hat{U}_{\pm}\,$ (3.18)
$\displaystyle[\hat{T}_{\pm},\hat{V}_{\mp}]$
$\displaystyle=\mp\hat{U}_{\mp}\hskip
49.79231pt[\hat{T}_{\pm},\hat{U}_{\pm}]=\pm\hat{V}_{\pm}\hskip
36.98866pt[\hat{U}_{\pm},\hat{V}_{\mp}]=\pm\hat{T}_{\mp}\,.$ (3.19)

The commutation relations involving the generators
$\\{\hat{T}_{3},\hat{Y},\hat{T}_{\pm},\hat{U}_{\pm},\hat{V}_{\pm}\\}$ are
those of the $\mathfrak{su}(3)$ Lie algebra [49], ($\hat{T}_{3},\hat{Y}$)
being the Cartan generators, which is known to be the symmetry algebra of the
three-dimensional isotropic harmonic oscillator [50, 51]. In particular, the
change of basis:

$\displaystyle\hat{X}_{1}$
$\displaystyle=\frac{1}{2}(\hat{T}_{+}+\hat{T}_{-})\qquad\hat{X}_{4}=\frac{1}{2}\bigl{(}\hat{V}_{+}+\hat{V}_{-})\qquad\hat{X}_{6}=\frac{1}{2}(\hat{U}_{+}+\hat{U}_{-})\qquad\hat{X}_{3}=\hat{T}_{3}$
(3.20) $\displaystyle\hat{X}_{2}$
$\displaystyle=\frac{\imath}{2}(\hat{T}_{-}-\hat{T}_{+})\qquad\hat{X}_{5}=\frac{\imath}{2}(\hat{V}_{-}-\hat{V}_{+})\qquad\hat{X}_{7}=\frac{\imath}{2}(\hat{U}_{-}-\hat{U}_{+})\qquad\hat{X}_{8}=\frac{\sqrt{3}}{2}\hat{Y}\,,$
(3.21)

allows us to recover the Hermitian generators $\hat{X}_{a}$ ($a=1,\dots,8$)
as:

$\hat{X}_{a}=\frac{1}{2}\hat{A}^{\dagger}_{i}[\lambda_{a}]_{ij}\hat{A}_{j}\,,$
(3.22)

where $\lambda_{a}$ are the eight Gell-mann matrices [52]. In this basis, the
commutation relations read:

$[\hat{X}_{a},\hat{X}_{b}]=\imath f_{abc}\hat{X}_{c}\,,$ (3.23)

with $f_{123}=1$, $f_{147}=f_{165}=f_{246}=f_{257}=f_{345}=f_{376}=1/2$,
$f_{458}=f_{678}=\sqrt{3}/2$. In the above realisation, the quantum integrals
of motion arise as the following combinations:

$\displaystyle\hat{X}_{1}$
$\displaystyle=\hat{F}_{12}\qquad\quad\,\hat{X}_{4}=\hat{F}_{13}\qquad\quad\,\,\,\hat{X}_{6}=\hat{F}_{23}\qquad\quad\,\,\hat{X}_{3}=(\hat{F}_{11}-\hat{F}_{22})/2$
(3.24) $\displaystyle\hat{X}_{2}$
$\displaystyle=\hat{L}_{12}/2\qquad\,\,\hat{X}_{5}=\hat{L}_{13}/2\qquad\,\,\hat{X}_{7}=\hat{L}_{23}/2\qquad\,\,\hat{X}_{8}=(\hat{F}_{11}+\hat{F}_{22}-2\hat{F}_{33})/2\sqrt{3}\,,$
(3.25)

where
$\hat{F}_{ij}:=\frac{1}{2}(\hat{p}_{i}\hat{p}_{j}+\hat{x}_{i}\hat{x}_{j})$ is
the Demkov-Fradkin tensor [50, 53] and
$\boldsymbol{\hat{L}}=(\hat{L}_{23},-\hat{L}_{13},\hat{L}_{12})$ with
$\hat{L}_{ij}=\hat{x}_{i}\hat{p}_{j}-\hat{x}_{j}\hat{p}_{i}$, is the angular
momentum vector. In this basis, the second-order and third-order Casimir
invariants read:

$\displaystyle\hat{C}^{(2)}_{\mathfrak{su}(3)}$
$\displaystyle=\hat{X}_{1}^{2}+\hat{X}_{2}^{2}+\hat{X}_{3}^{2}+\hat{X}_{4}^{2}+\hat{X}_{5}^{2}+\hat{X}_{6}^{2}+\hat{X}_{7}^{2}+\hat{X}_{8}^{2}$
(3.26) $\displaystyle\hat{C}_{\mathfrak{su}(3)}^{(3)}$
$\displaystyle=\frac{3}{2}\hat{X}_{3}-\frac{\sqrt{3}}{2}\hat{X}_{8}+\sqrt{3}(\hat{X}_{1}^{2}\hat{X}_{8}+\hat{X}_{2}^{2}\hat{X}_{8}+\hat{X}_{3}^{2}\hat{X}_{8})+3(\hat{X}_{1}\hat{X}_{4}\hat{X}_{6}+\hat{X}_{1}\hat{X}_{5}\hat{X}_{7}-\hat{X}_{2}\hat{X}_{4}\hat{X}_{7}+\hat{X}_{2}\hat{X}_{5}\hat{X}_{6})$
$\displaystyle\,+\frac{3}{2}(\hat{X}_{3}\hat{X}_{4}^{2}+\hat{X}_{3}\hat{X}_{5}^{2}-\hat{X}_{3}\hat{X}_{6}^{2}-\hat{X}_{3}\hat{X}_{7}^{2})-\frac{\sqrt{3}}{2}(\hat{X}_{4}^{2}\hat{X}_{8}+\hat{X}_{5}^{2}\hat{X}_{8}+\hat{X}_{6}^{2}\hat{X}_{8}+\hat{X}_{7}^{2}\hat{X}_{8})-\frac{1}{\sqrt{3}}\hat{X}_{8}^{3}\,,$
(3.27)

and result in the following expressions involving the Hamiltonian operator:

$\displaystyle\hat{C}^{(2)}_{\mathfrak{su}(3)}=\frac{1}{3}(\hat{H}^{2}-9/4)\,,\qquad\hat{C}_{\mathfrak{su}(3)}^{(3)}=\frac{1}{9}(\hat{H}^{2}-9/4)\hat{H}\,.$
(3.28)

The existence of Casimir invariants for the general polynomial algebra defined
above is an open problem.

######  Remark 3.

If we consider Hamiltonian systems such as for $m_{1}=m_{2}=m_{3}=1$ they
describe the standard harmonic oscillator, then the polynomial algebra we have
constructed can be considered as a polynomial deformation of the Lie algebra
$\mathfrak{u}(1)\oplus\mathfrak{su}(3)$ involving higher order powers of both
the Cartan generators and the Hamiltonian, the latter being a central element.
Let us consider, as an example, the three-dimensional Hamiltonian system:

$\hat{H}=\hat{H}_{1}+\hat{H}_{2}+\hat{H}_{3}=\frac{1}{2}\bigl{(}-\partial_{x_{1}x_{1}}+x^{2}_{1}/m_{1}^{2}\bigl{)}+\frac{1}{2}\bigl{(}-\partial_{x_{2}x_{2}}+x^{2}_{2}/m_{2}^{2}\bigl{)}+\frac{1}{2}\bigl{(}-\partial_{x_{3}x_{3}}+x^{3}_{3}/m_{3}^{2}\bigl{)}\,,$
(3.29)

where $m_{1},m_{2},m_{3}$ are mutually prime natural numbers. The ladder
operators are first-order:

$\hat{A}_{i}=\frac{\partial_{x_{i}}+x_{i}/m_{i}}{\sqrt{2}}\,,\qquad\hat{A}^{\dagger}_{i}=\frac{-\partial_{x_{i}}+x_{i}/m_{i}}{\sqrt{2}}\,,$
(3.30)

and satisfy the commutation relations:

$\displaystyle[\hat{H}_{i},\hat{A}_{i}]$
$\displaystyle=-\frac{1}{m_{i}}\hat{A}_{i}\hskip
46.94687pt[\hat{H}_{i},\hat{A}^{\dagger}_{i}]=\frac{1}{m_{i}}\hat{A}^{\dagger}_{i}$
(3.31) $\displaystyle\hat{A}_{i}^{\dagger}\hat{A}_{i}$
$\displaystyle=\hat{H}_{i}-1/(2m_{i})\qquad\quad\,\,\hat{A}_{i}\hat{A}_{i}^{\dagger}\\!=\hat{H}_{i}+1/(2m_{i})\,.$
(3.32)

Since $\alpha_{k}=1/m_{k}$, the condition
$\alpha_{i}s_{i}=\alpha_{j}s_{j}=\alpha$ is satisfied for $s_{k}=m_{k}$,
$k=1,2,3$ and $\alpha=1$. This means that the generators of the polynomial
algebra (3.6)-(3.11) take the form:

$\hat{H}=\hat{H}_{1}+\hat{H}_{2}+\hat{H}_{3}\,,\quad\hat{T}_{3}=\frac{\hat{H}_{1}-\hat{H}_{2}}{2}\,,\quad\hat{Y}=\frac{\hat{H}_{1}+\hat{H}_{2}-2\hat{H}_{3}}{3}\,,\quad\hat{E}_{ij}=(\hat{A}^{\dagger}_{i})^{m_{i}}(\hat{A}_{j})^{m_{j}}\quad(i\neq
j)\,.$ (3.33)

In this case, the structure functions are given by:

$\displaystyle\Xi_{12}^{(m_{1},m_{2})}(\hat{H},\hat{T}_{3},\hat{Y})$
$\displaystyle:=\Phi_{1}(m_{1},\hat{H}/3+\hat{T}_{3}+\hat{Y}/2)\Phi_{2}(m_{2},\hat{H}/3-\hat{T}_{3}+\hat{Y}/2+1)$
(3.34) $\displaystyle\Xi_{13}^{(m_{1},m_{3})}(\hat{H},\hat{T}_{3},\hat{Y})$
$\displaystyle:=\Phi_{1}(m_{1},\hat{H}/3+\hat{T}_{3}+\hat{Y}/2)\Phi_{3}(m_{3},\hat{H}/3-\hat{Y}+1)$
(3.35) $\displaystyle\Xi_{23}^{(m_{2},m_{3})}(\hat{H},\hat{T}_{3},\hat{Y})$
$\displaystyle:=\Phi_{2}(m_{2},\hat{H}/3-\hat{T}_{3}+\hat{Y}/2)\Phi_{3}(m_{3},\hat{H}/3-\hat{Y}+1)\,,$
(3.36)

with:

$\displaystyle\Phi_{1}(m_{1},\hat{H}/3+\hat{T}_{3}+\hat{Y}/2)$
$\displaystyle=\prod_{i=1}^{m_{1}}\biggl{(}\hat{H}/3+\hat{T}_{3}+\hat{Y}/2-\frac{2i-1}{2m_{1}}\biggl{)}$
(3.37) $\displaystyle\Phi_{2}(m_{2},\hat{H}/3-\hat{T}_{3}+\hat{Y}/2)$
$\displaystyle=\prod_{i=1}^{m_{2}}\biggl{(}\hat{H}/3-\hat{T}_{3}+\hat{Y}/2-\frac{2i-1}{2m_{2}}\biggl{)}$
(3.38) $\displaystyle\Phi_{3}(m_{3},\hat{H}/3-\hat{Y})$
$\displaystyle=\prod_{i=1}^{m_{3}}\biggl{(}\hat{H}/3-\hat{Y}-\frac{2i-1}{2m_{3}}\biggl{)}\,.$
(3.39)

For $m_{1}=m_{2}=m_{3}=1$ the three-dimensional isotropic harmonic oscillator
is recovered.

###  3.1 Rational extension of the harmonic oscillator associated with type
III Hermite exceptional orthogonal polynomials

What we are interested to investigate in this Section 3.1 is a $3$D rational
extension of the harmonic oscillator associated with Hermite exceptional
orthogonal polynomials (EOP) of type III [22, 31, 34, 35, 36, 37, 38]. The
two-dimensional version of this problem has been investigated in [23], where a
new pair of ladder operators for the one-dimensional rationally extended
harmonic oscillator have been introduced. Here, we deal with the three-
dimensional case and, in particular, we consider the Hamiltonian system in
(3.1) with:

$\hat{H}_{1}=-\partial_{x_{1}x_{1}}+x_{1}^{2}-2\biggl{[}\frac{\mathcal{H}^{\prime\prime}_{m}}{\mathcal{H}_{m}}-\biggl{(}\frac{\mathcal{H}^{\prime}_{m}}{\mathcal{H}_{m}}\biggl{)}^{2}+\,1\biggl{]}\,,\quad\hat{H}_{2}=-\partial_{x_{2}x_{2}}+x_{2}^{2}\,,\quad\hat{H}_{3}=-\partial_{x_{3}x_{3}}+x_{3}^{2}$
(3.40)

where $\mathcal{H}_{m}:=(-\imath)^{m}H_{m}(\imath x)$, $H_{m}(z)$ being the
usual Hermite orthogonal polynomials with $m$ even. The one-dimensional
constituents of this three-dimensional Hamiltonian system are two standard
harmonic oscillators in the $x_{2}$ and $x_{3}$ axes and a rationally extended
oscillator in the $x_{1}$ axis. In one dimension, the latter arises as a
superpartner of the standard Hermite harmonic oscillator in the framework of
Supersymmetric Quantum Mechanics (SUSYQM) [54, 55]. Let us consider the
following pair of SUSY partners:

$\hat{H}_{+}:=-\frac{\text{d}^{2}}{\text{d}x^{2}}+V_{+}(x)-E_{m}\,,\quad\hat{H}_{-}:=-\frac{\text{d}^{2}}{\text{d}x^{2}}+V_{-}(x)-E_{m}\,,$
(3.41)

where the two potentials $V_{\pm}(x):=W^{2}(x)\mp W^{\prime}(x)+E_{m}$ are
given in terms of the superpotential $W(x)=-\phi^{\prime}_{m}(x)/\phi_{m}(x)$.
Here $E_{m}:=-2m-1$, $\phi_{m}(x)=\mathcal{H}_{m}(x)\exp(x^{2}/2)$ is solution
of the equation $\hat{H}_{+}\phi_{m}(x)=0$ and the partner potentials read:

$V_{+}(x)=x^{2}\,,\quad
V_{-}(x)=x^{2}-2\biggl{(}\frac{\mathcal{H}_{m}^{\prime\prime}}{\mathcal{H}_{m}}-\biggl{(}\frac{\mathcal{H}_{m}^{\prime}}{\mathcal{H}_{m}}\biggl{)}^{2}\,+\,1\biggl{)}\,.$
(3.42)

These two Hamiltonians are connected through the action of the first-order
ladder operators:

$a:=\frac{\text{d}}{\text{d}x}+W(x)=\frac{\text{d}}{\text{d}x}-\biggl{(}x+\frac{\mathcal{H}^{\prime}_{m}}{\mathcal{H}_{m}}\biggl{)}\,,\quad
a^{\dagger}=-\frac{\text{d}}{\text{d}x}+W(x)=-\frac{\text{d}}{\text{d}x}-\biggl{(}x+\frac{\mathcal{H}^{\prime}_{m}}{\mathcal{H}_{m}}\biggl{)}\,,$
(3.43)

as they interwine as:

$\hat{a}\hat{H}_{+}=\hat{H}_{-}\hat{a}\,,\quad\hat{a}^{\dagger}\hat{H}_{-}=\hat{H}_{+}\hat{a}^{\dagger}\,.$
(3.44)

In particular, for even values of $m$, the eigenfunctions and eigenvalues of
the two superpartner $\hat{H}_{\pm}$ turn out to be:

$\begin{cases}\psi_{n}^{+}(x)=\mathcal{N}_{n}^{+}H_{n}(x)\exp(-x^{2}/2)\,,\qquad\hskip
4.26773ptE_{n}^{+}=2(n+m+1)\qquad(n=0,1,2,\dots)\\\
\psi_{n}^{-}(x)=\mathcal{N}_{n}^{-}\frac{y_{n+m+1}^{(m)}}{\mathcal{H}_{m}(x)}\exp(-x^{2}/2)\,,\qquad
E_{n}^{-}=2(n+m+1)\qquad(n=-m-1,0,1,2,\dots)\end{cases}$ (3.45)

where $y_{0}^{(m)}(x)=1$ and
$y_{k+m+1}^{(m)}(x)=-\mathcal{H}_{m}(x)H_{k+1}(x)-2m\mathcal{H}_{m-1}(x)H_{k}(x)\quad(k=0,1,2,\dots)$
are given in terms of the type III Hermite EOP $y_{\mathsf{n}}^{(m)}(x)$, with
$\mathsf{n}:=n+m+1$. The normalisation constants read:

$\displaystyle\mathcal{N}_{n}^{+}:=(\sqrt{\pi}2^{n}n!)^{-1/2}\quad\text{and}\quad$
$\displaystyle\mathcal{N}_{n}^{-}:=\begin{cases}(\frac{2^{m}m!}{\sqrt{\pi}})^{1/2}&\quad\text{if
}n=-m-1\\\ (\sqrt{\pi}2^{n+1}(n+m+1)n!)^{-1/2}&\quad\text{if
}n=0,1,2,\dots\,.\end{cases}$

This model is characterized by an additional bound state placed below the
oscillator spectrum and associated to the energy eigenvalue $E_{-m-1}=-2m-1$.
The eigenfunctions defined above are orthonormal in $L^{2}(\mathbb{R})$ w.r.t.
the positive definite measure
$\text{d}\mu(x):=\exp(-x^{2})/(\mathcal{H}_{m})^{2}\text{d}x$. In the work
[23], the authors showed the existence of a pair of lowering and raising
operators ($\hat{A}_{m}$, $\hat{A}_{m}^{\dagger}$) for the Hamiltonian
$\hat{H}_{-}$ of the rationally extended oscillator. These ladder operators
are constructed through the composite action of the first-order ladder
operators (3.43) and the additional $m$ first-order auxiliary operators:

$a_{i}:=\frac{\text{d}}{\text{d}x}+\biggl{(}x+\frac{\mathcal{H}^{\prime}_{i-1}}{\mathcal{H}_{i-1}}-\frac{\mathcal{H}^{\prime}_{i}}{\mathcal{H}_{i}}\biggl{)}\,,\quad
a_{i}^{\dagger}=-\frac{\text{d}}{\text{d}x}+\biggl{(}x+\frac{\mathcal{H}^{\prime}_{i-1}}{\mathcal{H}_{i-1}}-\frac{\mathcal{H}^{\prime}_{i}}{\mathcal{H}_{i}}\biggl{)}\qquad(i=1,\dots,m)\,,$
(3.46)

such as:

$\hat{a}_{i}^{\dagger}\hat{a}_{i}=\hat{H}_{i}\,,\quad\hat{a}_{i}\hat{a}_{i}^{\dagger}=\hat{H}_{i+1}+2\,,\quad\text{with}\quad\hat{H}_{i}=-\frac{\text{d}^{2}}{\text{d}x^{2}}+x^{2}-2\biggl{(}\frac{\mathcal{H}_{i-1}^{\prime\prime}}{\mathcal{H}_{i-1}}-\biggl{(}\frac{\mathcal{H}^{\prime}_{i-1}}{\mathcal{H}_{i-1}}\biggl{)}^{2}\,\biggl{)}-\,3\,.$
(3.47)

Specifically, by considering the recursive application of the $m+1$ first-
order operators (3.43) and (3.46), it is possible to define the higher order
lowering and raising operators:

$\hat{A}_{m}:=\biggl{(}\prod_{i=1}^{m}\hat{a}_{m-i+1}\biggl{)}\hat{a}^{\dagger}\qquad\hat{A}^{\dagger}_{m}:=\hat{a}\biggl{(}\prod_{i=1}^{m}\hat{a}^{\dagger}_{i}\biggl{)}\,,$
(3.48)

for which hold:

$[\hat{H}_{-},\hat{A}_{m}]=-(2m+2)\hat{A}_{m}\,\quad[\hat{H}_{-},\hat{A}^{\dagger}_{m}]=(2m+2)\hat{A}^{\dagger}_{m}\,\quad\,[\hat{A}_{m},\hat{A}^{\dagger}_{m}]=\phi(\hat{H}_{-}+2m+2)-\phi(\hat{H}_{-})\,,$
(3.49)

where the following polynomial in $\hat{H}_{-}$ has been introduced:

$\mathcal{\phi}(\hat{H}_{-}):=(\hat{H}_{-})\prod_{i=1}^{m}(\hat{H}_{-}-2m-2-2i)\,.$
(3.50)

In (3.48) we used the equality
$\bigl{(}\prod_{i=1}^{m}\hat{a}_{m-i+1}\bigl{)}^{\dagger}=\prod_{i=1}^{m}\hat{a}^{\dagger}_{i}$.
Then, this new pair of $(m+1)$th-order ladder operators, together with
$\hat{H}_{-}$, closes in a polynomial Heisenberg algebra of order $m$. These
are the ladder operators we will be interested in this paper in relation to
our polynomial algebra. We will specify later in the Section their explicit
action on the eigenstates of the rationally extended oscillator.

The 3D Hamiltonian that we have introduced arises in terms of the one-
dimensional superpartner Hamiltonians:

$\begin{cases}\hat{H}_{1}:=\hat{H}_{-}-2m-1\quad\text{in
the}\,\,x_{1}-\text{axis}\quad(\text{rationally extended oscillator})\\\
\hat{H}_{2}:=\hat{H}_{+}-2m-1\quad\text{in
the}\,\,x_{2}-\text{axis}\quad(\text{standard oscillator})\\\
\hat{H}_{3}:=\hat{H}_{+}-2m-1\quad\text{in
the}\,\,x_{3}-\text{axis}\quad(\text{standard oscillator})\,.\end{cases}$
(3.51)

Clearly, due to Cartesian separability, the eigenfunctions and eigenvalues
associated to the three-dimensional Hamiltonian operator
$\hat{H}=\hat{H}_{1}+\hat{H}_{2}+\hat{H}_{3}$ turns out to be:

$\displaystyle\Psi_{n_{1},n_{2},n_{3}}^{(m)}(x_{1},x_{2},x_{3})$
$\displaystyle=\mathcal{N}_{n_{1}}\mathcal{N}_{n_{2}}\mathcal{N}_{n_{3}}\frac{y^{(m)}_{n_{1}+m+1}(x_{1})}{\mathcal{H}_{m}(x_{1})}H_{n_{2}}(x_{2})H_{n_{3}}(x_{3})e^{-\frac{1}{2}(x_{1}^{2}+x_{2}^{2}+x_{3}^{2})}$
(3.52) $\displaystyle E_{n_{1},n_{2},n_{3}}$
$\displaystyle=2(n_{1}+n_{2}+n_{3})+3$ (3.53)

with $n_{1}=-m-1,0,1,2,\dots\,,n_{2}=0,1,2,\dots\,,n_{3}=0,1,2,\dots\,\,$.
Taking into account our previous discussion, let us introduce the following
pairs of ladder operators in the three directions $(x_{1},x_{2},x_{3})$:

$\displaystyle\hat{A}_{1;m}$
$\displaystyle=\hat{a}_{m}\dots\hat{a}_{1}\hat{a}^{\dagger}\,\qquad\quad\,\,\hat{A}_{1;m}^{\dagger}=\hat{a}\,\hat{a}_{1}^{\dagger}\dots\hat{a}_{m}^{\dagger}$
(3.54) $\displaystyle\hat{A}_{2}$
$\displaystyle=\partial_{x_{2}}+x_{2}\,\hskip
54.06006pt\hat{A}_{2}^{\dagger}=-\partial_{x_{2}}+x_{2}$ (3.55)
$\displaystyle\hat{A}_{3}$ $\displaystyle=\partial_{x_{3}}+x_{3}\,\hskip
54.06006pt\hat{A}_{3}^{\dagger}=-\partial_{x_{3}}+x_{3}\,,$ (3.56)

satisfying, together with their corresponding one-dimensional Hamiltonians,
the commutation relations:

$\displaystyle[\hat{H}_{1},\hat{A}_{1;m}]$
$\displaystyle=-(2m+2)\hat{A}_{1;m}\,\quad[\hat{H}_{1},\hat{A}^{\dagger}_{1;m}]=(2m+2)\hat{A}^{\dagger}_{1;m}\,\quad\,[\hat{A}_{1;m},\hat{A}^{\dagger}_{1;m}]=\phi_{1}(\hat{H}_{1}+2m+2)-\phi_{1}(\hat{H}_{1})$
(3.57) $\displaystyle[\hat{H}_{2},\hat{A}_{2}]$
$\displaystyle=-2\hat{A}_{2}\,\hskip
62.59596pt[\hat{H}_{2},\hat{A}^{\dagger}_{2}]=2\hat{A}^{\dagger}_{2}\,\hskip
73.26575pt[\hat{A}_{2},\hat{A}^{\dagger}_{2}]=2$ (3.58)
$\displaystyle[\hat{H}_{3},\hat{A}_{3}]$ $\displaystyle=-2\hat{A}_{3}\,\hskip
65.44142pt[\hat{H}_{3},\hat{a}^{\dagger}_{3}]=2\hat{A}^{\dagger}_{3}\,\hskip
73.26575pt[\hat{A}_{3},\hat{A}^{\dagger}_{3}]=2\,,$ (3.59)

where:

$\phi_{1}(\hat{H}_{1})=(\hat{H}_{1}+2m+1)\prod_{i=1}^{m}(\hat{H}_{1}-1-2i)\,.$
(3.60)

Given the above relations, we can see that the equalities
$2(m+1)s_{1}=2s_{2}=2s_{3}=\alpha$ are satisfied for the values
$s_{1}=1,s_{2}=s_{3}=m+1$ and $\alpha=2(m+1)$. Thus, from the general results,
we obtain that the following operators:

$\displaystyle\hat{T}_{3}:=\frac{\hat{H}_{1}-\hat{H}_{2}}{4(m+1)}\qquad\qquad\quad\hat{Y}:=\frac{\hat{H}_{1}+\hat{H}_{2}-2\hat{H}_{3}}{6(m+1)}\,,$
(3.61)

together with the higher order ones:

$\displaystyle\hat{T}_{+}:=\hat{A}^{\dagger}_{1;m}(\hat{A}_{2})^{m+1}\qquad\hat{U}_{+}:=(\hat{A}^{\dagger}_{2})^{m+1}(\hat{A}_{3})^{m+1}\qquad\hat{V}_{+}:=\hat{A}^{\dagger}_{1;m}(\hat{A}_{3})^{m+1}$
(3.62)
$\displaystyle\hat{T}_{-}:=(\hat{A}^{\dagger}_{2})^{m+1}\hat{A}_{1;m}\qquad\hat{U}_{-}:=(\hat{A}^{\dagger}_{3})^{m+1}(\hat{A}_{2})^{m+1}\qquad\hat{V}_{-}:=(\hat{A}^{\dagger}_{3})^{m+1}\hat{A}_{1;m}\,,$
(3.63)

commute with $\hat{H}$ and close to give the polynomial algebra (3.5)-(3.14)
with structure functions given by:

$\displaystyle\Xi_{12}^{(1,m+1)}(\hat{H},\hat{T}_{3},\hat{Y})$
$\displaystyle=\Phi_{1}(1,\hat{H}/3+2(m+1)\hat{T}_{3}+(m+1)\hat{Y})\Phi_{2}(m+1,\hat{H}/3-2(m+1)\hat{T}_{3}+(m+1)\hat{Y}+2(m+1))$
$\displaystyle\Xi_{13}^{(1,m+1)}(\hat{H},\hat{T}_{3},\hat{Y})$
$\displaystyle=\Phi_{1}(1,\hat{H}/3+2(m+1)\hat{T}_{3}+(m+1)\hat{Y})\Phi_{3}(m+1,\hat{H}/3-2(m+1)\hat{Y}+2(m+1))$
$\displaystyle\Xi_{23}^{(m+1,m+1)}(\hat{H},\hat{T}_{3},\hat{Y})$
$\displaystyle=\Phi_{2}(m+1,\hat{H}/3-2(m+1)\hat{T}_{3}+(m+1)\hat{Y})\Phi_{3}(m+1,\hat{H}/3-2(m+1)\hat{Y}+2(m+1))\,,$

where:

$\displaystyle\Phi_{1}(1,\hat{H}/3+2(m+1)\hat{T}_{3}+(m+1)\hat{Y})$
$\displaystyle=(\hat{H}/3+2(m+1)\hat{T}_{3}+(m+1)\hat{Y}+2m+1)$
$\displaystyle\hskip
22.76228pt\times\prod_{i=1}^{m}(\hat{H}/3+2(m+1)\hat{T}_{3}+(m+1)\hat{Y}-1-2i)$
(3.64) $\displaystyle\Phi_{2}(m+1,\hat{H}/3-2(m+1)\hat{T}_{3}+(m+1)\hat{Y})$
$\displaystyle=\prod_{i=1}^{m+1}\biggl{(}\hat{H}/3-2(m+1)\hat{T}_{3}+(m+1)\hat{Y}-2i+1\biggl{)}$
(3.65) $\displaystyle\Phi_{3}(m+1,\hat{H}/3-2(m+1)\hat{Y})$
$\displaystyle=\prod_{i=1}^{m+1}\biggl{(}\hat{H}/3-2(m+1)\hat{Y}-2i+1\biggl{)}\,.$
(3.66)

####  3.1.1 Degeneracies of energy levels and irreps of the polynomial algebra

The construction of finite-dimensional representations for the polynomial
algebras is in general a difficult task. It represents however an important
problem as the finite-dimensional representations, besides providing the
degeneracies, can also offer constraints on the spectrum of the models for
which other methods are not in general applicable, in particular when the
potential involves special functions defined only by non-linear ODEs. In such
cases, algebraic approaches offer powerful tools for providing insights into
the spectrum. Then, there is a need of developing further algebraic methods.

The construction of finite-dimensional representations of polynomial algebras
have been so far limited mostly to the context of two-dimensional quantum
superintegrable systems and their related rank-$1$ polynomial algebras with
three generators. For such cases the approach due to Daskaloyannis has been
applied widely. It consists of constructing the corresponding deformed
oscillator algebras [2, 3, 4, 5, 20, 56].

It was pointed out in [23] that the analysis of solutions by the Daskaloyannis
approach can be complicated because it requires solving the various algebraic
constraints and removing non physical solutions. In [24] It was shown that a
more direct approach of constructing finite-dimensional representations by
using zero modes of integrals allow one to treat systems with higher order
integrals and related algebras more efficiently. This is the approach we plan
to adapt and develop in this Section. The higher rank polynomial algebras we
have introduced display interesting properties that will make the construction
of representations feasible. In Section 3, we have been able to identify
integrals playing the role of Cartan, raising and lowering operators,
respectively. We will also develop ideas which have similarity with the
induced module construction. This will allow us to describe non-trivial
degeneracies in terms of states connected by the higher order quantum
integrals, the latter being referred to as multiplets. The actions of these
higher order operators on eigenstates can be found by combining the specific
actions of their constituent ladder operators, which reads:

$\displaystyle\hat{A}^{\dagger}_{1;m}\ket{n_{1}}_{x_{1}}$
$\displaystyle=\begin{cases}\sqrt{2^{m+1}m!}\ket{0}_{x_{1}}\hskip
199.16928pt\text{if}\quad n_{1}=-m-1\\\
-\sqrt{2^{m+1}(n_{1}+2m+2)\prod_{i=0}^{m-1}(n_{1}+m-i)}\ket{n_{1}+m+1}_{x_{1}}\hskip
11.38092pt\text{if}\quad n_{1}=0,1,2,\dots\end{cases}$
$\displaystyle\hat{A}_{1;m}\ket{n_{1}}_{x_{1}}$
$\displaystyle=\begin{cases}0\hskip 262.33444pt\text{if}\quad
n_{1}=-m-1,1,2\dots,m\\\ \sqrt{2^{m+1}m!}\ket{-m-1}_{x_{1}}\hskip
172.13925pt\text{if}\quad n_{1}=0\\\
-\sqrt{2^{m+1}(n_{1}+m+1)\prod_{i=0}^{m-1}(n_{1}-i-1)}\ket{n_{1}-m-1}_{x_{1}}\hskip
22.76228pt\text{if}\quad n_{1}=m+1,m+2,\dots\end{cases}$
$\displaystyle(\hat{A}_{2}^{\dagger})^{m+1}\ket{n_{2}}_{x_{2}}$
$\displaystyle=\sqrt{\prod_{j=1}^{m+1}(n_{2}+j)}\ket{n_{2}+m+1}_{x_{2}}\hskip
146.53151pt\text{if}\quad n_{2}=0,1,2,\dots$
$\displaystyle(\hat{A}_{2})^{m+1}\ket{n_{2}}_{x_{2}}$
$\displaystyle=\begin{cases}0\hskip 264.6107pt\text{if}\quad
n_{2}=0,1,2,\dots,m\\\
\sqrt{\prod_{j=0}^{m}(n_{2}-j)}\ket{n_{2}-m-1}_{x_{2}}\hskip
133.72786pt\text{if}\quad n_{2}=m+1,m+2,\dots\end{cases}$
$\displaystyle(\hat{A}_{3}^{\dagger})^{m+1}\ket{n_{3}}_{x_{3}}$
$\displaystyle=\sqrt{\prod_{j=1}^{m+1}(n_{3}+j)}\ket{n_{3}+m+1}_{x_{3}}\hskip
147.95424pt\text{if}\quad n_{3}=0,1,2,\dots$
$\displaystyle(\hat{A}_{3})^{m+1}\ket{n_{3}}_{x_{3}}$
$\displaystyle=\begin{cases}0\hskip 266.03343pt\text{if}\quad
n_{3}=0,1,2,\dots,m\\\
\sqrt{\prod_{j=0}^{m}(n_{3}-j)}\ket{n_{3}-m-1}_{x_{3}}\hskip
135.15059pt\text{if}\quad n_{3}=m+1,m+2,\dots\end{cases}$

Observing the action of the ladder operators in the $x_{1}$ axis, and
considering the standard action the ladder have in the other two axes, we see
that the (unnormalized) basis states are given by:

$\displaystyle\ket{n_{1}}_{x_{1}}\ket{n_{2}}_{x_{2}}\ket{n_{3}}_{x_{3}}=\begin{cases}[(\hat{A}_{1;m}^{\dagger})^{k}\ket{m-1}_{x_{1}}]\otimes[(\hat{A}_{2}^{\dagger})^{n_{2}}\ket{0}_{x_{2}}]\otimes[(\hat{A}_{3}^{\dagger})^{n_{3}}\ket{0}_{x_{3}}]\,\propto\,\\!\ket{(m+1)(k-1)}_{x_{1}}\ket{n_{2}}_{x_{2}}\ket{n_{3}}_{x_{3}}\\\
[(\hat{A}_{1;m}^{\dagger})^{k}\ket{1}_{x_{1}}]\otimes[(\hat{A}_{2}^{\dagger})^{n_{2}}\ket{0}_{x_{2}}]\otimes[(\hat{A}_{3}^{\dagger})^{n_{3}}\ket{0}_{x_{3}}]\,\,\qquad\propto\ket{1+k(m+1)}_{x_{1}}\ket{n_{2}}_{x_{2}}\ket{n_{3}}_{x_{3}}\\\
\vdots\hskip 221.93158pt\vdots\\\
[(\hat{A}_{1;m}^{\dagger})^{k}\ket{m}_{x_{1}}]\otimes[(\hat{A}_{2}^{\dagger})^{n_{2}}\ket{0}_{x_{2}}]\otimes[(\hat{A}_{3}^{\dagger})^{n_{3}}\ket{0}_{x_{3}}]\qquad\propto\ket{m+k(m+1)}_{x_{1}}\ket{n_{2}}_{x_{2}}\ket{n_{3}}_{x_{3}}\,,\end{cases}$

for $k=0,1\dots$, $n_{2}=0,1,\dots$ and $n_{3}=0,1,\dots\,\,.$

######  Remark 4.

This decomposition is due to the fact that we are using this specific type of
ladder operators in the $x_{1}$ axis, for which it is known the polynomial
Heisenberg algebra possesses $m+1$ infinite-dimensional unirreps spanned by
the states $\ket{j+(m+1)k}_{x_{1}}$, for $j=-m-1,1,\dots,m$ and
$k=0,1,\dots\,\,$. Other types of ladder operators exist in the literature for
the rationally extended Hermite oscillator associated with type III Hermite
EOP [57, 58, 59, 60], for example, the pair [22]:

$\hat{B}_{m}=\hat{a}\,\hat{b}\,\hat{a}^{\dagger}\,\quad\hat{B}^{\dagger}_{m}=\hat{a}\,\hat{b}^{\dagger}\,\hat{a}^{\dagger}\,,$
(3.67)

with $(\hat{a},\hat{a}^{\dagger})$ given in (3.43) and
${b}=\frac{\text{d}}{\text{d}x}+x$,
${b}^{\dagger}=-\frac{\text{d}}{\text{d}x}+x$. However, as explained in [61],
the choice of higher order ladder operators (3.48) is the suitable one in the
study of higher order superintegrable models that are characterized by one-
dimensional potentials related to type III Hermite EOP.

######  Remark 5.

The set composed by the nine operators is not of minimal order for this model.
There are in fact two standard oscillators in the $x_{2}$ and $x_{3}$ axes and
therefore one generator of rotations is conserved. This would lead to a pair
of shift operators constructed with the product
$\hat{A}_{2}^{\dagger}\hat{A}_{3}$ and Hermitian conjugate.

If we set $N=n_{1}+n_{2}+n_{3}$, a direct computation shows that the
degeneracy of the energy level $E_{N}=2N+3$ is:

$\text{deg}(E_{N})=\begin{cases}1\hskip 113.81102pt\text{if}\quad N=-m-1\\\
2\hskip 113.81102pt\text{if}\quad N=-m\\\ \vdots\hskip 116.65646pt\vdots\\\
m+1\hskip 93.89418pt\text{if}\quad N=-1\\\ (N+2)(N+3)/2+m\hskip
21.33955pt\text{if}\quad N=0,1,2,\dots\,\,.\\\ \end{cases}$ (3.68)

At a fixed $N$, the basis states
$\ket{n_{1}}_{x_{1}}\ket{n_{2}}_{x_{2}}\ket{n_{3}}_{x_{3}}$ are connected by
the action of the operators $(\hat{T}_{\pm},\hat{U}_{\pm},\hat{V}_{\pm})$ and
find a useful graphical representation in the plane $(t_{3},y)$, the latter
being the eigenvalues associated to $\hat{T}_{3}$ and $\hat{Y}$ operators in
the above Cartesian basis, or equivalently, in the basis:
$\ket{N,n_{1},n_{2}}=\ket{n_{1}}_{x_{1}}\ket{n_{2}}_{x_{2}}\ket{N-(n_{1}+n_{2})}_{x_{3}}$.

In order to characterize degeneracies of the model in terms of finite-
dimensional irreps of the polynomial algebra, we denote the eigenstates of
$\hat{H}$ as $\ket{N,n_{1},n_{2}}$ and act on these physical states with shift
operators $\\{\hat{T}_{\pm},\hat{U}_{\pm},\hat{V}_{\pm}\\}$ to see how they
are connected for each $N$, i.e. for each energy level. Following up on the
approach developed in [24], we try to characterize some specific irreps of the
polynomial algebra by two integer numbers $(N,\mathbf{p})$ and their basis
states as $\ket{N,\tau,\mathbf{p};(\tilde{t}_{3},\tilde{y})}$ where
$\mathbf{p}:=(p+1)(p+2)/2$ specify the number of connected states, the pair
$(\tilde{t}_{3},\tilde{y})$ represents the eigenvalues of the operators
$\hat{T}_{3}+\alpha_{3}$ and $\hat{T}_{8}+\alpha_{8}$, where
$(\alpha_{3},\alpha_{8})$ are two representation-dependent constants and
$\tau$ distinguishes among repeated representations specified by the same
$\mathbf{p}$ value. Sets of connected states will be called “multiplets”and,
as we shall see in the following, they will be of triangular shape when
represented in the $(t_{3},y)$ plane. In particular, for each $N$ we require
the existence of a vector
$\ket{N,\tau,\mathbf{p};(\tilde{t}_{3}^{max},\tilde{y})}$ such as:

$\hat{T}_{+}\ket{N,\tau,\mathbf{p};(\tilde{t}_{3}^{max},\tilde{y})}=\hat{V}_{+}\ket{N,\tau,\mathbf{p};(\tilde{t}_{3}^{max},\tilde{y})}=0\,,$
(3.69)

and suppose that after $p+1$ applications of the operator $\hat{V}_{-}$ we
get:

$(\hat{V}_{-})^{p+1}\ket{N,\tau,\mathbf{p};(\tilde{t}_{3}^{max},\tilde{y})}=0\,.$
(3.70)

This operation completely specify the integer $p$ and, consequently, the
boundary of the multiplet $\mathbf{p}$. Then, to characterize all connected
states inside a given multiplet we act with:

$(\hat{T}_{-})^{a_{2}}(\hat{V}_{-})^{a_{1}}\ket{N,\tau,\mathbf{p};(\tilde{t}_{3}^{max},\tilde{y})}\,,$
(3.71)

for $0\leq a_{1}\leq p$ and $0\leq a_{2}\leq p-a_{1}$. To induce these irreps
we rely on the physical basis $\ket{N,n_{1},n_{2}}$ and on the explicit action
of the higher order integrals of motion on this basis. Specifically, for even
$m$, we adopt the following reparametrization:

$N=n_{1}+n_{2}+n_{3}=(m+1)\lambda+\mu\equiv N(m;\lambda,\mu)\,,$ (3.72)

for some appropriate values of $\lambda=-1,0,1,\dots$ and $\mu=0,1,\dots,m$.
This is helpful as in this way we are able to split the analysis of
degeneracies into four different regions of the energy spectrum specified by
given values of $\lambda$ and $\mu$. This turn out to be useful to simplify
the analysis of this particular problem. The energy spectrum reads:

$E_{N}|_{N=N(m;\lambda,\mu)}=2N(m;\lambda,\mu)+3=2((m+1)\lambda+\mu)+3\,,$
(3.73)

and it has to be intended, for fixed values of $\lambda=\lambda^{*}$ and
$\mu=\mu^{*}$, just as a linear function of the even integer $m$.

| $\lambda$ | $\mu$ | $N$ | $E_{N}$
---|---|---|---|---
I | $-1$ | $0,1,\dots,m$ | $-m-1,-m,-m+1,\dots,-1$ | $-2m+1,-2m+3,\dots,1$
II | $0$ | $0,1,\dots,m$ | $0,1,\dots,m$ | $3,5,7,\dots,2m+3$
III | $1$ | $0,1,\dots,m$ | $m+1,m+2,\dots,2m+1$ | $2m+5,2m+7,\dots,4m+5$
IV | $2,3,\dots$ | $0,1,\dots,m$ | $2m+2,2m+3,\dots,3m+2,3m+3,\dots,4m+3,\dots$ | $4m+7,4m+9,\dots,6m+7,6m+9,\dots$
Table 1: Values of $N$ ($\lambda$, $\mu$) specifying the four regions of the
energy spectrum.

Inside these four regions, that we report explicitly in Table 1, the eigenkets
will be then reparametrized in terms of $m$, $\lambda$ and $\mu$. The explicit
action of the higher order quantum integrals will then read:

$\displaystyle\hat{T}_{\pm}\ket{N,n_{1},n_{2}}$
$\displaystyle=t_{\pm}(m;N,n_{1},n_{2})\ket{N,n_{1}\pm(m+1),n_{2}\mp(m+1)}$
(3.74) $\displaystyle\hat{U}_{\pm}\ket{N,n_{1},n_{2}}$
$\displaystyle=u_{\pm}(m;N,n_{1},n_{2})\ket{N,n_{1},n_{2}\pm(m+1)}$ (3.75)
$\displaystyle\hat{V}_{\pm}\ket{N,n_{1},n_{2}}$
$\displaystyle=v_{\pm}(m;N,n_{1},n_{2})\ket{N,n_{1}\pm(m+1),n_{2}}\,,$ (3.76)

where
$\ket{N,n_{1},n_{2}}\equiv\ket{N(m;\lambda,\mu),n_{1}(m;\lambda,\mu),n_{2}(m;\lambda,\mu)}$
and the functions $(t_{\pm},u_{\pm},v_{\pm})$ will be directly computed from
the explicit action of the constituent ladder operators when specific values
of $\lambda$ and $\mu$ have been specified for the region we are considering.
The action of $\hat{T}_{3}$ and $\hat{Y}$ on this basis reads:

$\hat{T}_{3}\ket{N,n_{1},n_{2}}=\frac{n_{1}-n_{2}}{2(m+1)}\ket{N,n_{1},n_{2}}\qquad\hat{Y}\ket{N,n_{1},n_{2}}=\frac{3n_{1}+3n_{2}-2N}{3(m+1)}\ket{N,n_{1},n_{2}}\,.$
(3.77)

To have an idea on how higher order quantum integrals connect states, at fixed
energy level $N(m,\lambda^{*},\mu^{*})=(m+1)\lambda^{*}+\mu^{*}$, it is
convenient to rely on a graphical representation of their action in the plane
$(t_{3},y)$. In this plane physical states, as it can be understood from
(3.77), define triangular lattices. An illustrative example is reported in
Figure 1.

$-2$$-1$$1$$2$$-2$$-1$$1$$2$$t_{3}$$y$

$-2$$-1$$1$$2$$-2$$-1$$1$$2$$\text{\tiny{T}}_{\text{\raisebox{0.4pt}{\scalebox{0.4}{$\pm$}}}}$$\text{\tiny{V}}_{\text{\raisebox{0.4pt}{\scalebox{0.4}{$\pm$}}}}$$\text{\tiny{U}}_{\text{\raisebox{0.4pt}{\scalebox{0.4}{$\mp$}}}}$$t_{3}$$y$

Figure 1: On the left, a typical triangular lattice representing physical
states in the plane $(t_{3},y)$ is given, for a fixed value of
$N(m,\lambda^{*},\mu^{*})=N(2;0,2)=7$. On the right, we show an example of the
action of the higher order quantum integrals on a subset of connected states
(multiplet).

Throughout the paper, for the sake of clarity, we will be representing with
the same color physical states that are connected by the action of the higher
order quantum integrals. Then, the construction proceeds as follows:

  1. 1.

Inside a specific region of the energy spectrum, suppose a number $K=K(m;N)$
of basis eigenstates of $\hat{H}$ arise. We reparametrise these physical
states in the basis $\ket{N,n_{1},n_{2}}$ in terms of appropriate values of
$\lambda$ and $\mu$ valid inside the region we are considering. This specify
the following set composed by $K=K(m;N(m;\lambda,\mu))$ eigenkets:

$\small\mathcal{S}(m;\lambda,\mu):=\bigl{\\{}\ket{N(m;\lambda,\mu),n_{1}(m;\lambda,\mu),n_{2}(m;\lambda,\mu)}^{[1]},\dots,\ket{N(m;\lambda,\mu),n_{1}(m;\lambda,\mu),n_{2}(m;\lambda,\mu)}^{[\\#\mathcal{S}(m;\lambda,\mu)]}\bigl{\\}}\,,$
(3.78)

where we indicated $K\equiv\\#\mathcal{S}(m;\lambda,\mu)$, the cardinality of
the set (3.78).

  2. 2.

At fixed $\lambda=\lambda^{*}$, we decompose the above set into $m+1$ subsets
specified by values of $\mu=0,1,\dots,m$, i.e.:

$\mathcal{S}(m;\lambda^{*},\mu):=\mathcal{S}(m;\lambda^{*},0)\cup\mathcal{S}(m;\lambda^{*},1)\cup\dots\cup\mathcal{S}(m;\lambda^{*},m)\,.$
(3.79)

Each subset $\mathcal{S}(m;\lambda^{*},\mu^{*})$, for fixed $\mu=\mu^{*}$,
will then contain $\\#\mathcal{S}(m;\lambda^{*},\mu^{*})$ states and, of
course:

$\sum_{j=0}^{m}\\#\mathcal{S}(m;\lambda^{*},j)=\\#\mathcal{S}(m;\lambda^{*},\mu)\,.$

  3. 3.

We impose the conditions:

$\begin{cases}\hat{T}_{+}\ket{N(m;\lambda^{*},\mu),n_{1}(m;\lambda^{*},\mu),n_{2}(m;\lambda^{*},\mu)}^{[j]}=0\\\
\hat{V}_{+}\ket{N(m;\lambda^{*},\mu),n_{1}(m;\lambda^{*},\mu),n_{2}(m;\lambda^{*},\mu)}^{[j]}=0\end{cases}\qquad\text{for
each}\quad j=1,\dots,\\#\mathcal{S}(m;\lambda^{*},\mu)=0$ (3.80)

where
$\ket{N(m;\lambda^{*},\mu),n_{1}(m;\lambda^{*},\mu),n_{2}(m;\lambda^{*},\mu)}^{[j]}\in\mathcal{S}(m;\lambda^{*},\mu)$,
namely we look at those states that are simultaneously annihilated by the
action of two raising operators. These actions, that can be constructed from
the action of the constituent ladder operators in the three axes as specified
previously, when restricted to the specific regions ($\lambda^{*}$,
$\mu=0,1,\dots,m$), allow us to find the values $n_{1}=n_{1}(\lambda^{*},\mu)$
and $n_{2}=n_{2}(\lambda^{*},\mu)$ for which the system of equations is
satisfied inside each specific subset specified by $\mu=\mu^{*}$. This
operation leads us to identify, for each fixed
$N(m;\lambda^{*},\mu^{*})=(m+1)\lambda^{*}+\mu^{*}$, a number
$\\#\tilde{\mathcal{S}}(m;\lambda^{*},\mu^{*})$ of annihilated states:

$\footnotesize\tilde{\mathcal{S}}(m;\lambda^{*},\mu^{*})=\bigl{\\{}\ket{N(m;\lambda^{*},\mu^{*}),n_{1}(m;\lambda^{*},\mu^{*}),n_{2}(m;\lambda^{*},\mu^{*})}^{[1]},\dots,\ket{N(m;\lambda^{*},\mu^{*}),n_{1}(m;\lambda^{*},\mu^{*}),n_{2}(m;\lambda^{*},\mu^{*})}^{[\\#\tilde{\mathcal{S}}(m;\lambda^{*},\mu^{*})]}\bigl{\\}}\,.$
(3.81)

In the following, to simplify the notation, we indicate (3.80) as an operation
on sets, i.e.:

$(\hat{T}_{+},\hat{V}_{+})\rhd\mathcal{S}(m;\lambda^{*},\mu)=0\,.$

  4. 4.

We impose the condition:

$(\hat{V}_{-})^{p+1}\ket{N(m;\lambda^{*},\mu^{*}),n_{1}(m;\lambda^{*},\mu^{*}),n_{2}(m;\lambda^{*},\mu^{*})}^{[j]}=0\quad\text{for}\quad
j=1,\dots,\\#\tilde{\mathcal{S}}(m;\lambda^{*},\mu^{*})$ (3.82)

where
$\ket{N(m;\lambda^{*},\mu^{*}),n_{1}(m;\lambda^{*},\mu^{*}),n_{2}(m;\lambda^{*},\mu^{*})}^{[j]}\in\tilde{\mathcal{S}}(m;\lambda^{*},\mu^{*})$.
This allow us to determine the values $p_{k}$,
$k=1,\dots,\\#\tilde{\mathcal{S}}(m;\lambda^{*},\mu^{*})$, i.e. to determine
the boundary of the triangular multiplets composed by connected states. We
indicate the sets composed by these connected states as
$\mathcal{P}_{k}(m;\lambda^{*},\mu^{*})$. Each set
$\mathcal{P}_{k^{*}}(m;\lambda^{*},\mu^{*})$ turns out to have cardinality
$\\#\mathcal{P}_{k^{*}}(m;\lambda^{*},\mu^{*})=(p_{k^{*}}+1)(p_{k^{*}}+2)/2$.
In the following, we will specify these multiplets in term of their dimensions
$\boldsymbol{p_{k^{*}}}=(p_{k^{*}}+1)(p_{k^{*}}+2)/2$, so that depending on
the value assumed by $p_{k^{*}}$ we can have, besides singlets, $\mathbf{3}$,
$\mathbf{6}$, $\mathbf{10}$, $\dots$ multiplets of states connected by the
action of the higher order quantum integrals of motion. This means that, at
fixed energy, i.e. fixed $N=(m+1)\lambda^{*}+\mu^{*}$ the set of eigenstates:

$\footnotesize\mathcal{S}(m;\lambda^{*},\mu^{*})=\bigl{\\{}\ket{N(m;\lambda^{*},\mu^{*}),n_{1}(m;\lambda^{*},\mu^{*}),n_{2}(m;\lambda^{*},j^{*})}^{[1]},\dots,\ket{N(m;\lambda^{*},\mu^{*}),n_{1}(m;\lambda^{*},\mu^{*}),n_{2}(m;\lambda^{*},\mu^{*})}^{[\\#\mathcal{S}(m;\lambda^{*},\mu^{*})]}\bigl{\\}}$
(3.83)

will be organized into $\\#\tilde{\mathcal{S}}(m;\lambda^{*},\mu^{*})$
multiplets of dimension $(p_{k}+1)(p_{k}+2)/2$, for
$k=1,\dots,\\#\tilde{\mathcal{S}}(m;\lambda^{*},\mu^{*})$, namely:

$\mathcal{S}(m;\lambda^{*},\mu^{*})=\bigcup_{k=1}^{\\#\tilde{\mathcal{S}}(m;\lambda^{*},\mu^{*})}\mathcal{P}_{k}(m;\lambda^{*},\mu^{*})$

with:

$\\#\mathcal{S}(m;\lambda^{*},\mu^{*})=\sum_{k=1}^{\\#\tilde{\mathcal{S}}(m;\lambda^{*},\mu^{*})}(p_{k}+1)(p_{k}+2)/2\,.$

  5. 5.

Once multiplets have been identified, one can construct their “spectrum”by
applying:

$(\hat{T}_{-})^{a_{2}}(\hat{V}_{-})^{a_{1}}\ket{N(m;\lambda^{*},\mu^{*}),n_{1}(m;\lambda^{*},\mu^{*}),n_{2}(m;\lambda^{*},\mu^{*})}^{[j]}\quad\text{for}\quad
j=1,\dots,\\#\mathcal{P}_{k^{*}}(m;\lambda^{*},\mu^{*})$ (3.84)

with $0\leq a_{1}\leq p_{k^{*}}$ and $0\leq a_{2}\leq p_{k^{*}}-a_{1}$, where
$\ket{N(m;\lambda^{*},\mu^{*}),n_{1}(m;\lambda^{*},\mu^{*}),n_{2}(m;\lambda^{*},\mu^{*})}^{[j]}\in\mathcal{P}_{k^{*}}(m;\lambda^{*},\mu^{*})$.

At a fixed $N=(m+1)\lambda^{*}+\mu^{*}$ each state
$\ket{N(m;\lambda^{*},\mu^{*}),n_{1}(m;\lambda^{*},\mu^{*}),n_{2}(m;\lambda^{*},\mu^{*})}^{[j^{*}]}\in\mathcal{P}_{k^{*}}(m;\lambda^{*},\mu^{*})$
belonging to a given multiplet will be characterized by the values
($N,\boldsymbol{p_{k^{*}}}$) which provide information on the energy and on
the type of multiplet it belongs to, respectively. Then, the pair of values
$(t_{3},y)$ is used to specify their location inside a given multiplet, in
particular their rescaled values $(\tilde{t}_{3},\tilde{y})$ are considered to
set the center of gravity of the various triangles characterizing each irrep
in the origin of the $(t_{3},y)$ plane. Finally, when other multiplets are
considered for the same values of $N$, another label $\tau$ is necessary to
discriminate among repeated irreps arising with the same
$\boldsymbol{p_{k^{*}}}$. Clearly, normalizations of these states have to be
taken into account in the construction.

  6. 6.

The procedure terminates when for each specific region such that
$N=(m+1)\lambda+\mu$ all $\\#\mathcal{S}(m;\lambda,\mu)$ states in (3.78) are
organized into multiplets of connected states, so that irreps ($N,\mathbf{p}$)
and their basis states have been completely specified.

From now on, our analysis will be performed inside each of the four specific
regions we introduced. In particular, our main aim is restricted to the
characterization, for different values of $N$ (i.e. $\lambda$ and $\mu$), of
the explicit patterns of the $\mathbf{p}$ multiplets. This allows us to
provide the total number of irreps per given energy level and to characterize
the degeneracies of this model. We should again remark, however, that this is
not the algebra spanned by the minimal order generators for this model. This
means that states that are not connected by the action of the higher order
operators might be actually connected by the action of the existing first-
order ones. Therefore, even if this represents a good starting point to deal
with these kind of models in three dimensions, the polynomial algebra we
introduced arises more naturally when at least two rationally extended
oscillators (in two different axes) are considered, so that no residual polar
symmetry (leading to first-order quantum integrals) is preserved.

#### Region I

In region I, i.e. for $\lambda=-1$ and $\mu=0,1,\dots,m$, the basis states can
be parametrized as:

$\ket{N,n_{1},n_{2}}\qquad\text{with}\quad N=\mu-m-1,\quad
n_{1}=-m-1\quad\text{and}\quad 0\leq n_{2}\leq\mu\,.$

We thus have the following set of states:

$\mathcal{S}(m;-1,\mu)=\\{\ket{\mu-m-1,-m-1,n_{2}}\\}\,,\quad(0\leq
n_{2}\leq\mu)\,.$

If we decompose this set as
$\mathcal{S}(m;-1,\mu)=\mathcal{S}(m;-1,0)\cup\mathcal{S}(m;-1,1)\cup\dots\cup\mathcal{S}(m;-1,m)$,
namely:

$\small\mathcal{S}(m;-1,\mu)=\\{\ket{-m-1,-m-1,0}\\}\cup\\{\ket{-m,-m-1,0},\ket{-m,-m-1,1}\\}\cup\dots\cup\\{\ket{-1,-m-1,0},\dots,\ket{-1,-m-1,m}\\}$
(3.85)

its cardinality is easily computed to be
$\\#\mathcal{S}(m;-1,\mu)=1+2+\dots+m+1=\sum_{j=1}^{m+1}j=(m+1)(m+2)/2$. This
represents the total number of basis states arising in this region. All these
states are annihilated by the action of the higher order integrals. For
$N=-m-1$, i.e. $\mu=0$, there is just one state (a singlet) which we indicate
as $\mathbf{1}$. For $N=-m$, i.e. $\mu=1$, there are two singlets, so that we
can indicate them as $\mathbf{1}^{2}$. The process continues until we reach
$N=-1$, i.e. $\mu=m$, where we have $\mathbf{1}^{m+1}$. For each even $m$, the
eigenvalues of $(\hat{T}_{3},\hat{Y})$ operators associated to these states
turn out to be:

$t_{3}(m;n_{2})=-\frac{m+1+n_{2}}{2(m+1)}\,,\quad
y(m;\mu,n_{2})=\frac{3n_{2}-(2\mu+m+1)}{3(m+1)}\qquad(0\leq n_{2}\leq\mu)\,.$
(3.86)

These can be plotted in the plane $(t_{3},y)$ for different values of
$\mu=0,1,\dots,m$ at fixed $m$ (see Figure 2).

$-1$$-0.5$$0.5$$1$$-1$$-0.5$$0.5$$1$$t_{3}$$y$

$-1$$-0.5$$0.5$$1$$-1$$-0.5$$0.5$$1$$t_{3}$$y$

$-1$$-0.5$$0.5$$1$$-1$$-0.5$$0.5$$1$$t_{3}$$y$

Figure 2: Physical states arising for $N=-3$ (upper left), $N=-2$ (upper
right) and $N=-1$ ($m=2$).

For different values of $m$, the obtained results for region I are summarised
in Table 2.

$\lambda$ | $\mu$ |  $\mathbf{p}$ multiplets | Number of Irreps per level |  deg($E_{N}$)
---|---|---|---|---
$-1$ | $0,1,2,\dots,m$ | $\boldsymbol{1}^{\mu+1}$ | $\mu+1$ | $\mu+1$
Table 2: $\mathbf{p}$ multiplets with their number of occurrences, number of
irreps of the polynomial algebra and degeneracy per energy level
characterizing region I.

#### Region II

With the reparametrisation $N=(m+1)\lambda+\mu$, we see that states obtained
for $\lambda=0$ and $\mu=0,1\dots,m$ can be parametrized as:

$\ket{N,n_{1},n_{2}}\qquad\text{with}\quad N=\mu,\quad
n_{1}=-m-1,0,\dots,\mu\quad\text{and}\quad 0\leq n_{2}\leq\mu-
n_{1}\quad(\mu=0,1,\dots,m)\,.$ (3.87)

Explictly, we have the following set of states:

$\mathcal{S}(m;0,\mu)=\begin{Bmatrix}\ket{\mu,-m-1,0}&\ket{\mu,0,0}&\ket{\mu,1,0}&\dots&\ket{\mu,\mu-2,0}&\ket{\mu,\mu-1,0}&\ket{\mu,\mu,0}\\\
\ket{\mu,-m-1,1}&\ket{\mu,0,1}&\ket{\mu,1,1}&\dots&\ket{\mu,\mu-2,1}&\ket{\mu,\mu-1,1}&\\\
\ket{\mu,-m-1,2}&\ket{\mu,0,2}&\ket{\mu,1,2}&\dots&\ket{\mu,\mu-2,2}&&\\\
\vdots&\vdots&\vdots&&&&\\\
\ket{\mu,-m-1,\mu+m-1}&\ket{\mu,0,\mu+m-1}&\ket{\mu,1,\mu+m-1}&&&&\\\
\ket{\mu,-m-1,\mu+m}&\ket{\mu,0,\mu+m}&&&&&\\\
\ket{\mu,-m-1,\mu+m+1}&&&&&&\end{Bmatrix}$

At fixed $\mu=\mu^{*}$, the are a total number of $(\mu^{*}+2)(\mu^{*}+3)/2+m$
states. Thus, in this region, the total number of basis states is
$\sum_{j=0}^{m}\bigl{(}(j+2)(j+3)/2+m\bigl{)}=(m+1)(m^{2}+14m+18)/6$. Starting
from the set of states $\mathcal{S}(m;0,\mu)$, we first decompose it as:

$\mathcal{S}(m;0,\mu)=\mathcal{S}(m;0,0)\cup\mathcal{S}(m;0,1)\cup\dots\cup\mathcal{S}(m;0,m)\,,$

i.e.:

$\mathcal{S}(m;0,\mu)=\begin{Bmatrix}\ket{0,-m-1,0}&\ket{0,0,0}\\\
\ket{0,-m-1,1}&\\\ \ket{0,-m-1,2}&\\\ \vdots&\\\ \ket{0,-m-1,m-1}&\\\
\ket{0,-m-1,m}&\\\ \ket{0,-m-1,m+1}&\\\ \end{Bmatrix}$ $\cup$
$\begin{Bmatrix}\ket{1,-m-1,0}&\ket{1,0,0}&\ket{1,1,0}\\\
\ket{1,-m-1,1}&\ket{1,0,1}&\\\ \ket{1,-m-1,2}&&\\\ \vdots&&\\\
\ket{1,-m-1,m}&&\\\ \ket{1,-m-1,m+1}&&\\\ \ket{1,-m-1,m+2}&&\\\
\end{Bmatrix}\cup\dots\cup\begin{Bmatrix}\ket{m,-m-1,0}&\ket{m,0,0}&\ket{m,1,0}&\dots&\ket{m,m-2,0}&\ket{m,m-1,0}&\ket{m,m,0}\\\
\dots&\ket{m,m-2,1}&\ket{m,m-1,1}&\\\
\ket{m,-m-1,2}&\ket{m,0,2}&\ket{m,1,2}&\dots&\ket{m,m-2,2}&&\\\
\vdots&\vdots&\vdots&&&&\\\
\ket{m,-m-1,2m-1}&\ket{m,0,2m-1}&\ket{m,1,2m-1}&&&&\\\
\ket{m,-m-1,2m}&\ket{m,0,2m}&&&&&\\\ \ket{m,-m-1,2m+1}&&&&&&\\\
[1ex]\end{Bmatrix}$

Then, we impose:

$(\hat{T}_{+},\hat{V}_{+})\rhd\\{\mathcal{S}(m;0,0),\mathcal{S}(m;0,1),\dots,\mathcal{S}(m;0,m)\\}=0\,.$

In this way, we select the following subset composed by
$\\#\tilde{\mathcal{S}}(m;0,\mu)=(m+1)(m^{2}+8m+6)/6$ states:

$\tilde{\mathcal{S}}(m;0,\mu)=\begin{Bmatrix}&\ket{0,0,0}\\\
\ket{0,-m-1,1}&\\\ \ket{0,-m-1,2}&\\\ \vdots&\\\ \ket{0,-m-1,m-1}&\\\
\ket{0,-m-1,m}&\end{Bmatrix}$ $\cup$
$\begin{Bmatrix}&\ket{1,0,0}&\ket{1,1,0}\\\ &\ket{1,0,1}&\\\
\ket{1,-m-1,2}&&\\\ \vdots&&\\\ \ket{1,-m-1,m}&&\end{Bmatrix}$ $\cup\dots\cup$
$\begin{Bmatrix}&\ket{m,0,0}&\ket{m,1,0}&\dots&\ket{m,m-2,0}&\ket{m,m-1,0}&\ket{m,m,0}\\\
&\ket{m,0,1}&\ket{m,1,1}&\dots&\ket{m,m-2,1}&\ket{m,m-1,1}&\\\
&\ket{m,0,2}&\ket{m,1,2}&\dots&\ket{m,m-2,2}&&\\\
\vdots&\vdots&\vdots&\iddots&&&\\\ &\ket{m,0,m-1}&\ket{m,1,m-1}&&&&\\\
&\ket{m,0,m}&&&&&\end{Bmatrix}$ .

At this point, by imposing:

$(\hat{V}_{-})^{p+1}\rhd\\{\tilde{\mathcal{S}}(m;0,0),\tilde{\mathcal{S}}(m;0,1),\dots,\tilde{\mathcal{S}}(m;0,m)\\}=0\,,$

to each of these subsets, we find the values of $p$ associated to the various
states, so that we are finally able to organize them into $\mathbf{p}$
multiplets. For the sake of clarity, let us focus on some specific examples.
For $\mu=0$, the following $m+3$ states arise:

$\mathcal{S}(m;0,0)=\\{\ket{0,-m-1,0}\,\ket{0,-m-1,1}\,,\ket{0,-m-1,2},\dots,\ket{0,-m-1,m+1},\ket{0,0,0}\\}\,.$
(3.88)

Among these states, we select the ones annihilated by
$(\hat{T}_{+},\hat{V}_{+})$. The action of $\hat{T}_{+}$ specifies the
following values of $n_{2}$, $n_{2}=0,1,\dots,m$. Thus, the set of $m+3$
states reduces to the following set composed by $m+2$ states:

$\\{\ket{0,-m-1,0}\,\ket{0,-m-1,1}\,,\ket{0,-m-1,2},\dots,\ket{0,-m-1,m},\ket{0,0,0}\\}\,.$
(3.89)

Among these states, the following $m+1$ are also annihilated by the action of
$\hat{V}_{+}$:

$\tilde{\mathcal{S}}(m;0,0)=\\{\ket{0,-m-1,1}\,,\ket{0,-m-1,2},\dots,\ket{0,-m-1,m},\ket{0,0,0}\\}\,.$
(3.90)

This is because this operator annihilates states for which
$n_{1}+n_{2}=0,-1,-2,\dots,-m$ in this region for $\mu=0$. Therefore, the
states annihilated by the action of both operators turn out to be:

$\ket{0,n_{1},n_{2}}\quad\text{with}\quad
n_{1}=-m-1\,,n_{2}=1,\dots,m\quad\text{or}\quad n_{1}=n_{2}=0\,.$ (3.91)

######  Remark 6.

In general, for $N=\mu$, starting from the set of basis states (3.87), the
ones annihilated by $(\hat{T}_{+},\hat{V}_{+})$ are those for which the
following systems of equations is satisfied:

$\begin{cases}n_{1}+n_{2}=\mu,\mu-1,\dots,\mu-m\\\
n_{2}=0,1,\dots,m\,.\end{cases}$ (3.92)

This means that at fixed $\mu=\mu^{*}$, for different values of
$n_{1}=n_{1}^{*},n_{1}^{**},\dots$ we obtain the set of values of $n_{2}$ from
both equations. The intersection set of these values give us the solution we
are looking for. The are a total number of
$\\#\tilde{\mathcal{S}}(m;0,\mu)=(m+1)(5m^{2}+m+6)/6$ of these states.

At this point, by applying the operator $\hat{V}_{-}$ a repeated number of
times until we get zero, we can specify the multiplets composed by connected
states. In particular, for the above set we get:

$\begin{cases}\hat{V}_{-}\ket{0,-m-1,j}=0\hskip 28.45274pt\text{for}\quad
j=1,\dots,m\\\ (\hat{V}_{-})^{2}\ket{0,0,0}=0\,.\end{cases}$ (3.93)

This means that the $m+3$ states organize as one triplet $\mathbf{3}$ and $m$
singlets $\mathbf{1}^{m}$ under the action of the higher order integrals. If
$m=2$, for example, five physical states arise:

$\mathcal{S}(2;0,0)=\\{\ket{0,-3,0},\ket{0,-3,1},\ket{0,-3,2},\ket{0,-3,3},\ket{0,0,0}\\}\,.$

The states annihilated by the action of $(\hat{T}_{+},\hat{V}_{+})$ are
$\tilde{S}(2,0,0)=\\{\ket{0,-3,1},\ket{0,-3,2},\ket{0,0,0}\\}$. The operator
$\hat{V}_{-}$ annihilates the first two, so they can be considered as
singlets, whereas for the third one we have $(\hat{V}_{-})^{2}\ket{0,0,0}=0$.
At this point, by computing the explicit action of ladder operators on
$\ket{N,n_{1},n_{2}}=\ket{0,0,0}$ we get the path:

$\hat{T}_{+}\hat{U}_{+}\hat{V}_{-}\ket{0,0,0}\propto\ket{0,0,0}\,.$ (3.94)

We represent this situation in Figure 3, where we indicate in red states
belonging to the same $\mathbf{3}$, i.e. states connected by the action of the
higher order quantum integrals. The remaining states are the two singlets.

$-1$$-0.5$$0.5$$1$$-1$$-0.5$$0.5$$1$$t_{3}$$y$

$-1$$-0.5$$0.5$$1$$-1$$-0.5$$0.5$$1$$t_{3}$$y$

Figure 3: Physical states arising for $N=0$ represented in the plane
$(t_{3},y)$ for $m=2$. States connected by the action of higher order
integrals are indicated in red.

When $N=1$, i.e. $\lambda=0$ and $\mu=1$, for $m=2$ eight physical states
arise (see (3.68)), six of them organized as $\mathbf{3}^{2}$ and the
remaining two as $\mathbf{1}^{2}$. For $N=2$, i.e. $\lambda=0,\mu=2$ twelve
states arise, nine of them organized as $\mathbf{3}^{3}$ and the remaining
three as $\mathbf{1}^{3}$. We represent these situations in Figure 4.

$-1$$-0.5$$0.5$$1$$-1$$-0.5$$0.5$$1$$t_{3}$$y$

$-1$$-0.5$$0.5$$1$$-1$$-0.5$$0.5$$1$$t_{3}$$y$

Figure 4: Physical states for $N=1$ (left) and $N=2$ (right) organized in
$\mathbf{p}$ multiplets ($m=2$).

By repeating this procedure for different values of $m$ we have been able to
identify the pattern of the multiplets in region II. We report the obtained
results in Table 3.

$\lambda$ | $\mu$ |  $\mathbf{p}$ multiplets | Number of Irreps per level |  deg($E_{N}$)
---|---|---|---|---
$0$ | $0,1,2,\dots,m$ | $\begin{cases}\boldsymbol{3}^{\mu+1}\\\ \boldsymbol{1}^{m+\frac{\mu(\mu-1)}{2}}\end{cases}$ | $m+1+\frac{\mu(\mu+1)}{2}$ | $\frac{(\mu+2)(\mu+3)}{2}+m$
Table 3: $\mathbf{p}$ multiplets with their number of occurrences, number of
irreps of the polynomial algebra and degeneracy per energy level
characterizing region II.

#### Region III

For $N=m+1,m+2,\dots,2m+1$, i.e. $\lambda=1$, $\mu=0,1,2,\dots,m$ we enter
into a new regime, as a new larger $\mathbf{p}$ multiplet arise. This is
because we have sufficient room to accomodate connected physical states into a
$\mathbf{6}$ multiplet, i.e. we are allowed to ladder two times among physical
states with the higher order integrals (but no more than this). The set of
physical states in this region, for $\mu=0,1,\dots,m$, can be parametrized as:

$\ket{N,n_{1},n_{2}}\qquad\text{with}\quad N=\mu+m+1,\quad
n_{1}=-m-1,0,\dots,\mu+m+1\quad\text{and}\quad 0\leq
n_{2}\leq\mu+m+1-n_{1}\,.$ (3.95)

The are $\\#\mathcal{S}(m;1,\mu)=(m+1)(7m^{2}+38m+36)/6$ basis states in this
region. Among them, the ones for which:

$\begin{cases}n_{1}+n_{2}=\mu+m+1,\mu+m,\mu+m-1,\dots,\mu+1\\\
n_{2}=0,1,\dots,m\,,\end{cases}$ (3.96)

are all annihilated by the action of the $(\hat{T}_{+},\hat{V}_{+})$
operators. The are a total number
$\\#\tilde{\mathcal{S}}(m;1,\mu)=(m+1)(m+2)(5m+3)$ of these states. Let us
take the example $N=m+1$ (i.e. $\lambda=1$, $\mu=0$). We have a total number
of $\\#\mathcal{S}(m;1,0)=(m+3)(m+4)/2+m$ physical states:

$\mathcal{S}(m;1,0)=\begin{Bmatrix}\ket{m+1,-m-1,0}&\ket{m+1,0,0}&\ket{m+1,1,0}&\dots&\ket{m+1,m+1,0}\\\
\ket{m+1,-m-1,1}&\ket{m+1,0,1}&\ket{m+1,1,1}&&\\\ \vdots&\vdots&\vdots&&\\\
\ket{m+1,-m-1,m}&\ket{m+1,0,m}&\ket{m+1,1,m}&&\\\
\ket{m+1,-m-1,m+1}&\ket{m+1,0,m+1}&&&\\\ \vdots&&&&\\\
\ket{m+1,-m-1,2m+2}&&&&\end{Bmatrix}$

The ones annihilated by $(\hat{T}_{+},\hat{V}_{+})$ are the ones for which the
following system is satisfied:

$\begin{cases}n_{1}+n_{2}=m+1,m,m-1,\dots,1\\\
n_{2}=0,1,\dots,m\,.\end{cases}$ (3.97)

If $n_{1}=-m-1$, i.e. we consider the first column in $\mathcal{S}(m;1,0)$, we
get:

$\begin{cases}n_{2}=2m+2,2m+1,2m-1,\dots,m+2\\\
n_{2}=0,1,\dots,m\,.\end{cases}$ (3.98)

Thus, no common solutions exist. For $n_{1}=0,1,\dots,m+1$, we get:

$\begin{cases}n_{2}=m+1-n_{1},m-n_{1},m-1-n_{1},\dots,1-n_{1}\\\
n_{2}=0,1,\dots,m\,,\end{cases}$ (3.99)

and the solutions are given by:

$\begin{cases}n_{1}=0\qquad\quad\rightarrow\quad\quad n_{2}=1,\dots,m\\\
n_{1}=1\qquad\quad\rightarrow\quad\quad n_{2}=0,\dots,m\\\
n_{1}=2\qquad\quad\rightarrow\quad\quad n_{2}=0,\dots,m-1\\\
n_{1}=3\qquad\quad\rightarrow\quad\quad n_{2}=0,\dots,m-2\\\ \vdots\hskip
91.04872pt\vdots\\\ n_{1}=m+1\quad\rightarrow\quad\quad n_{2}=0\,.\end{cases}$
(3.100)

We thus have a total number of $m+(m+1)+(m+(m-1)+\dots+1)=(m+1)(m+2)/2+m$
states annihilated by $\hat{V}_{+}$ and $\hat{T}_{+}$. They can be all
parametrized as:

$\ket{m+1,0,n_{2}}\quad
n_{2}=1,\dots,m\quad\cup\quad\ket{m+1,n_{1},n_{2}}\quad\,n_{1}=1,\dots,m+1\,,\quad
0\leq n_{2}\leq m-n_{1}+1\,.$ (3.101)

If we now apply $(\hat{V}_{-})^{p+1}$ until we get zero on these states, we
can see how the total number of $(m+3)(m+4)/2+m$ states are organized. For
example, for $m=2$, the states with $N=3$ turn out to be connected as in
Figure 5.

$-1.5$$-1$$-0.5$$0.5$$1$$1.5$$-1$$1$$t_{3}$$y$ Figure 5: Physical states for
$N=3$ organized in $\mathbf{p}$ multiplets ($m=2$).

Observing Figure 5, we can see that the seventeen states organize as a
$\mathbf{6}$ multiplet (red), two $\mathbf{3}$ multiplets (blue, orange), and
five singlets under the action of the higher order integrals. By repeating the
construction for generic $N=m+1,m+2,\dots,2m+1$, i.e. $\lambda=1$,
$\mu=0,1,\dots,m$ we have been able to identify the pattern of multiplets in
region III. We report the obtained results in Table 4.

$\lambda$ | $\mu$ |  $\mathbf{p}$ multiplets | Number of Irreps per level |  deg($E_{N}$)
---|---|---|---|---
$1$ | $0,1,2,\dots,m$ | $\begin{cases}\boldsymbol{6}^{\mu+1}\\\ \boldsymbol{3}^{m+\frac{\mu(\mu-1)}{2}}\\\ \boldsymbol{1}^{\frac{m(m+3)}{2}+\mu(m-\mu-1)}\\\ \end{cases}$ | $\footnotesize\frac{(m+1)(m+2)}{2}+\frac{\mu(2m-\mu-1)}{2}+m$ | $\frac{(m+\mu+3)(m+\mu+4)}{2}+m$
Table 4: $\mathbf{p}$ multiplets with their number of occurrences, number of
irreps of the polynomial algebra and degeneracy per energy level
characterizing region III.

For $m=2$, states with $N=4,5$ also belong to this region. We report the
multiplets in Figure 6.

$-2$$-1$$1$$2$$-2$$-1$$1$$2$$t_{3}$$y$

$-2$$-1$$1$$2$$-2$$-1$$1$$2$$t_{3}$$y$

Figure 6: Physical states arising for $N=4$ (left) and $N=5$ (right) organized
in $\mathbf{p}$ multiplets ($m=2$).

#### Region IV

If we go ahead in the construction by considering higher values of $N$, we
realize that other regimes arise for $N=2m+2,2m+3,\dots,3m+2$,
$N=3m+3,3m+4,\dots,4m+2$ etc. In particular, for $\lambda=2,3\dots$ and
$\mu=0,1,\dots,m$, the set of basis states is the following:

$\ket{N,n_{1},n_{2}}\qquad\text{with}\quad N=(m+1)\lambda+\mu,\quad
n_{1}=-m-1,0,\dots,(m+1)\lambda+\mu\quad\text{and}\quad 0\leq
n_{2}\leq(m+1)\lambda+\mu-n_{1}\,.$

Let us focus on the first sub-region $\lambda=2$, $\mu=0,1,2,\dots,m$. We
have:

$\mathcal{S}(m;2,\mu)=$

$\begin{Bmatrix}\ket{2m+2,-m-1,0}&\ket{2m+2,0,0}&\ket{2m+2,1,0}&\dots&\ket{2m+2,2m+1,0}&\ket{2m+2,2m+2,0}\\\
\ket{2m+2,-m-1,1}&\ket{2m+2,0,1}&\ket{2m+2,1,1}&\dots&\ket{2m+2,2m+1,1}&\\\
\vdots&\vdots&\vdots&&&\\\
\ket{2m+2,-m-1,2m+1}&\ket{2m+2,0,2m+1}&\ket{2m+2,1,2m+1}&&&\\\
\ket{2m+2,-m-1,2m+2}&\ket{2m+2,0,2m+2}&&&&\\\ \vdots&&&&&\\\
\ket{2m+2,-m-1,3m+3}&&&&&\end{Bmatrix}$
$\cup\dots\cup$

$\begin{Bmatrix}\ket{3m+2,-m-1,0}&\ket{3m+2,0,0}&\ket{3m+2,1,0}&\dots&\ket{3m+2,3m+1,0}&\ket{3m+2,3m+2,0}\\\
\ket{3m+2,-m-1,1}&\ket{3m+2,0,1}&\ket{3m+2,1,1}&\dots&\ket{3m+2,3m+1,1}&\\\
\vdots&\vdots&\vdots&&&\\\
\ket{3m+2,-m-1,3m+1}&\ket{3m+2,0,3m+1}&\ket{3m+2,1,3m+1}&&&\\\
\ket{3m+2,-m-1,3m+2}&\ket{3m+2,0,3m+2}&&&&\\\ \vdots&&&&&\\\
\ket{3m+2,-m-1,4m+3}&&&&&\end{Bmatrix}$

Again, by imposing:

$(\hat{T}_{+},\hat{V}_{+})\rhd\\{\mathcal{S}(m;2,0),\mathcal{S}(m;2,1),\dots,\mathcal{S}(m;2,m)\\}=0\,,$

we select the subset $\tilde{\mathcal{S}}(m;2,\mu)$ composed by
$\\#\tilde{\mathcal{S}}(m;2,\mu)=(m+1)^{3}$ states annihilated by the above
operators. Once these states are found, we act with $(\hat{V}_{-})^{p+1}$ to
specify the boundary of the multiplets of connected states.

Let us perform the construction for $N=2m+2$, i.e. $\lambda=2,\mu=0$ to
provide some explicit example inside this region. For $\lambda=2,\mu=0$ the
following set of $\\#\mathcal{S}(m;2,0)=(2m+4)(2m+5)/2+m$ physical states
arise:

$\mathcal{S}(m;2,0)$ =
$\begin{Bmatrix}\ket{2m+2,-m-1,0}&\ket{2m+2,0,0}&\ket{2m+2,1,0}&\dots&\ket{2m+2,2m+1,0}&\ket{2m+2,2m+2,0}\\\
\ket{2m+2,-m-1,1}&\ket{2m+2,0,1}&\ket{2m+2,1,1}&\dots&\ket{2m+2,2m+1,1}&\\\
\vdots&\vdots&\vdots&&&\\\
\ket{2m+2,-m-1,2m+1}&\ket{2m+2,0,2m+1}&\ket{2m+2,1,2m+1}&&&\\\
\ket{2m+2,-m-1,2m+2}&\ket{2m+2,0,2m+2}&&&&\\\ \vdots&&&&&\\\
\ket{2m+2,-m-1,3m+3}&&&&&\end{Bmatrix}$

The ones annihilated by $(\hat{T}_{+},\hat{V}_{+})$, i.e. the subset
$\tilde{\mathcal{S}}(m;2,0)$, are those for which the following system is
satisfied:

$\begin{cases}n_{1}+n_{2}=2m+2,2m+1,\dots,m+3,m+2\\\
n_{2}=0,1,\dots,m\,,\end{cases}$ (3.102)

namely:

$\begin{cases}n_{1}=2\qquad\quad\,\,\rightarrow\quad\quad n_{2}=m\\\
n_{1}=3\qquad\quad\,\,\rightarrow\quad\quad n_{2}=m-1,m\\\ \vdots\hskip
93.89418pt\vdots\\\ n_{1}=m\qquad\,\,\,\,\,\,\rightarrow\quad\quad
n_{2}=2,\dots,m\\\ n_{1}=m+1\quad\,\,\rightarrow\quad\quad n_{2}=1,\dots,m\\\
n_{1}=m+2\quad\,\,\rightarrow\quad\quad n_{2}=0,\dots,m\\\ \vdots\hskip
93.89418pt\vdots\\\ n_{1}=2m+1\quad\rightarrow\quad\quad n_{2}=0,1\\\
n_{1}=2m+2\quad\rightarrow\quad\quad n_{2}=0\,.\end{cases}$ (3.103)

We thus have a total number of
$(1+2+\dots+m+1)+(m+(m-1)+\dots+1)=(m+1)(m+2)/2+m(m+1)/2=(m+1)^{2}$ states
annihilated by the action of both $\hat{V}_{+}$ and $\hat{T}_{+}$. Thus,
$\\#\tilde{\mathcal{S}}(m;2,0)=(m+1)^{2}$. They can be parametrized as:

$\ket{2m+2,n_{1},n_{2}}\quad n_{1}=2,\dots,m+2\,,\quad m+2-n_{1}\leq n_{2}\leq
m\quad\cup\quad n_{1}=m+3,\dots,2m+2\,,\quad 0\leq n_{2}\leq 2m+2-n_{1}\,.$

If we now apply $(\hat{V}_{-})^{p+1}$ until we get zero on these states, we
can see how the total number of $(2m+4)(2m+5)/2+m$ states are organized under
the action of the higher order integrals. The result we get is that they
decompose as $\mathbf{10}$, $\mathbf{6}^{m}$, $\mathbf{3}^{m(m+3)/2}$ and
$\mathbf{1}^{m(m-1)/2}$. Clearly, $10+6m+3m(m+3)/2+m(m-1)/2=(2m+4)(2m+5)+m$.
Now, by repeating this process for each $\lambda=3,4,\dots$, for
$\mu=0,1,\dots,m$ we have been able to identify the patterns for all the
multiplets.

We report the results that we have obtained in Table 5.

$\lambda$ | $\mu$ |  $\mathbf{p}$ multiplets | Number of Irreps per level |  deg($E_{N}$)
---|---|---|---|---
$2,3,\dots$ | $0,1,2,\dots,m$ | $\begin{cases}\boldsymbol{\frac{(\lambda+2)(\lambda+3)}{2}}^{\mu+1}\\\ \boldsymbol{\frac{(\lambda+1)(\lambda+2)}{2}}^{m+\frac{\mu(\mu-1)}{2}}\\\ \boldsymbol{\frac{(\lambda)(\lambda+1)}{2}}^{\frac{m(m+3)}{2}+\mu(m-\mu-1)}\\\ \boldsymbol{\frac{(\lambda-1)(\lambda)}{2}}^{\frac{(m-\mu)(m-\mu-1)}{2}}\end{cases}$ | $(m+1)^{2}$ | $\frac{((m+1)\lambda+\mu+2)((m+1)\lambda+\mu+3)}{2}+m$
Table 5: $\mathbf{p}$ multiplets with their number of occurrences, number of
irreps and degeneracy per energy level.

Let us fix, just as a clarifying example, $m=2$ and consider the case
$\lambda=2,\mu=0,1,2$ i.e. $N=6,7,8$. For $N=6$, $\\#\mathcal{S}(2;2,0)=38$
states $\mathcal{S}(2;2,0)$ arise. Among them, there are
$\\#\tilde{\mathcal{S}}(2;2,0)=9$ states belonging to the subset
$\tilde{\mathcal{S}}(2;2,0)$. For $N=7$, $\\#\mathcal{S}(2;2,1)=47$ states
$\mathcal{S}(2;2,0)$ arise. Among them, there are
$\\#\tilde{\mathcal{S}}(2;2,1)=9$ states belonging to the subset
$\tilde{\mathcal{S}}(2;2,1)$. Finally, for $N=8$, $\\#\mathcal{S}(2;2,2)=57$
states $\mathcal{S}(2;2,2)$ arise. Among them, there are again
$\\#\tilde{\mathcal{S}}(2;2,2)=9$ states belonging to the subset
$\tilde{\mathcal{S}}(2;2,2)$. These states, and the way they are connected
through the action of the higher order integrals is reported in Figure 7.

$-2$$-1$$1$$2$$-2$$-1$$1$$2$$t_{3}$$y$

$-2$$-1$$1$$2$$-2$$-1$$1$$2$$t_{3}$$y$

$-2$$-1$$1$$2$$-2$$-1$$1$$2$$t_{3}$$y$

Figure 7: Physical states for $N=6$ (upper left), $N=7$ (upper right) and
$N=8$ organized in $\mathbf{p}$ multiplets ($m=2$).

To summarize, for each energy level all multiplets composed by physical states
connected through the action of the higher order integrals have been
completely classified for all of the four regions into we splitted up our
initial problem. The general results we have obtained are reported in Table 6.

| $\lambda$ | $\mu$ |  $\mathbf{p}$ multiplets | Number of Irreps per level |  deg($E_{N}$)
---|---|---|---|---|---
I | $-1$ | $0,1,2,\dots,m$ | $\boldsymbol{1}^{\mu+1}$ | $\mu+1$ | $\mu+1$
II | $0$ | $0,1,2,\dots,m$ | $\begin{cases}\boldsymbol{3}^{\mu+1}\\\ \boldsymbol{1}^{m+\frac{\mu(\mu-1)}{2}}\end{cases}$ | $m+1+\frac{\mu(\mu+1)}{2}$ | $\frac{(\mu+2)(\mu+3)}{2}+m$
III | $1$ | $0,1,2,\dots,m$ | $\begin{cases}\boldsymbol{6}^{\mu+1}\\\ \boldsymbol{3}^{m+\frac{\mu(\mu-1)}{2}}\\\ \boldsymbol{1}^{\frac{m(m+3)}{2}+\mu(m-\mu-1)}\\\ \end{cases}$ | $\footnotesize\frac{(m+1)(m+2)}{2}+\frac{\mu(2m-\mu-1)}{2}+m$ | $\frac{(m+\mu+3)(m+\mu+4)}{2}+m$
IV | $2,3,\dots$ | $0,1,2,\dots,m$ | $\begin{cases}\boldsymbol{\frac{(\lambda+2)(\lambda+3)}{2}}^{\mu+1}\\\ \boldsymbol{\frac{(\lambda+1)(\lambda+2)}{2}}^{m+\frac{\mu(\mu-1)}{2}}\\\ \boldsymbol{\frac{(\lambda)(\lambda+1)}{2}}^{\frac{m(m+3)}{2}+\mu(m-\mu-1)}\\\ \boldsymbol{\frac{(\lambda-1)(\lambda)}{2}}^{\frac{(m-\mu)(m-\mu-1)}{2}}\end{cases}$ | $(m+1)^{2}$ | $\frac{((m+1)\lambda+\mu+2)((m+1)\lambda+\mu+3)}{2}+m$
Table 6: $\mathbf{p}$ multiplets with their number of occurrences, number of
irreps of the polynomial algebra and degeneracy per energy level for the 3D
rationally extended harmonic oscillator (3.40).

##  4 Concluding Remarks

In this work we have investigated the polynomial algebras associated with
$n$-dimensional superintegrable systems, with separation of variables in
Cartesian coordinates, that are endowed with higher order ladder operators
satisfying polynomial deformations of the Heisenberg algebra. These models are
part of a very large class of superintegrable systems (in fact the
classification of one-dimensional systems with operator algebra is still an
ongoing problem itself [48]). They have connection to Painlevé transcendents
[62, 63, 64] but in general they are associated with nonlinear equations of
Painlevé type. The study of these equations is itself a difficult problem as
shown by a series of papers by Cosgrove [65, 66]. However, there are families
of deformations that are of particular interest and that can be classified
completely. Among them is the case of $k$-step extension of the harmonic
oscillator [24], which is related to the multi-indexed Hermite exceptional
orthogonal polynomials.

We have introduced a new three-dimensional model involving Hermite exceptional
orthogonal polynomials of type III. This system exhibits very interesting
degeneracies for its energy spectrum. One of the key results is the
application of the underlying polynomial algebra structure to characterize
degeneracies in terms of its finite-dimensional representations. We have
performed calculations in order to obtain explicit expressions for a number of
distinct irreducible representations at a given energy level. Such
decomposition in terms of representations allows us to get further insights
into the structure of the states and further information on the degeneracies.

It is demonstrated that these representations share some similarities with
those of $\mathfrak{su}(3)$ and that states of given energy can be decomposed
into sum of irreducible representations which take the form of
$\mathfrak{su}(3)$-like triangular multiplets. In recent years, it was
discovered how quadratic, cubic and more generally polynomial algebras appear
in different contexts of mathematical physics [67, 68, 69, 70, 71, 72, 73,
74]. The class of polynomial algebras and related representations we have
obtained in this paper are of interest for further applications from both
mathematical and physical perspectives.

### Acknowledgements

This work was supported by the Future Fellowship FT180100099 and Discovery
Project DP190101529 from the Australian Research Council.

## References

  * [1] Ya. I. Granovskii, I. M. Lutzenko and A. S. Zhedanov, “Mutual integrability, quadratic algebras, and dynamical symmetry”. Ann. Phys. 217 (1992), 1.
  * [2] D. Bonatsos, C. Daskaloyannis, P. Kolokotronis and D. Lenis, “The symmetry algebra of the N-dimensional anisotropic quantum harmonic oscillator with rational ratios of frequencies and the Nilsson model”.  arXiv:hep-th/9411218 [hep-th].
  * [3] D. Bonatsos, C. Daskaloyannis and K. Kokkotas, “Deformed oscillator algebras for two-dimensional quantum superintegrable systems”. Phys. Rev. A 50 (1994), 3700.
  * [4] D. Bonatsos and C. Daskaloyannis, “Quantum groups and their applications in nuclear physics”. Progr. Part. Nucl. Phys. 43 (1999), 537.
  * [5] C. Daskaloyannis, “Quadratic Poisson algebras of two-dimensional classical superintegrable systems and quadratic associative algebras of quantum superintegrable systems”. J. Math. Phys. 42 (2001), 1100.
  * [6] S. Post, “Models of Quadratic Algebras Generated by Superintegrable Systems in 2D”. SIGMA 7 (2011), 036.
  * [7] W. Jr. Miller, S. Post and P. Winternitz, “Classical and Quantum Superintegrability with Applications”. J. Phys. A: Math. Theor. 46 (2013), 423001.
  * [8] E. G. Kalnins, “Contractions of 2D 2nd Order Quantum Superintegrable Systems and the Askey Scheme for Hypergeometric Orthogonal Polynomials”. SIGMA 9 (2013), 057.
  * [9] I. Marquette, “Superintegrability with third order integrals of motion, cubic algebras, and supersymmetric quantum mechanics. I. Rational function potentials”. J. Math. Phys. 50 (2009), 012101.
  * [10] I. Marquette, “Superintegrability and higher order polynomial algebras”. J. Phys. A: Math. Theor. 43 (2010), 135203.
  * [11] E. G. Kalnins, “A Recurrence Relation Approach to Higher Order Quantum Superintegrability”. SIGMA 7 (2011), 031.
  * [12] I. Marquette, “An infinite family of superintegrable systems from higher order ladder operators and supersymmetry”. J. Phys.: Conf. Ser. 284 (2011), 012047.
  * [13] A. M. Escobar-Ruiz, J. C. López Vieyra and P. Winternitz, “Fourth order superintegrable systems separating in polar coordinates. I. Exotic potentials”. J. Phys. A: Math. Theor. 50 (2017), 495206.
  * [14] I. Marquette, S. Post and L. Ritter, “A fourth-order superintegrable system with a rational potential related to Painlevé VI”. J. Phys. A: Math. Theor. 53 (2020), 50LT01.
  * [15] Md F. Hoque, I. Marquette and Y.-Z. Zhang, “Quadratic algebra structure and spectrum of a new superintegrable system in ${N}$-dimension”. J. Phys. A: Math. Theor. 48 (2015), 185201.
  * [16] Z. Chen, I. Marquette and Y.-Z. Zhang, “Superintegrable systems from block separation of variables and unified derivation of their quadratic algebras”. Ann. Phys. 411 (2019), 167970.
  * [17] H. De Bie, P. Iliev, W. van de Vijver and L. Vinet, “The Racah algebra: An overview and recent results”. Contemp. Math. 768 (2021), 3.
  * [18] D. Latini, I. Marquette and Y.-Z. Zhang, “Embedding of the Racah algebra $R(n)$ and superintegrability”. Ann. Phys. 426 (2021), 168397.
  * [19] F. Correa, Md F. Hoque, I. Marquette and Y.-Z. Zhang, “$N$-dimensional Smorodinsky–Winternitz model and related higher rank quadratic algebra $SW(N)$”. J. Phys. A: Math. Theor. 54 (2021), 395201.
  * [20] D. Bonatsos, C. Daskaloyannis, P. Kolokotronis and D. Lenis, “Symmetry algebra of the planar anisotropic quantum harmonic oscillator with rational ratio of frequencies”.  arXiv:hep-th/9402099 [hep-th].
  * [21] S. Post, S. Tsujimoto and L. Vinet, “Families of superintegrable Hamiltonians constructed from exceptional polynomials”. J. Phys. A: Math. Theor. 45 (2012), 405202.
  * [22] I. Marquette and C. Quesne, “New families of superintegrable systems from Hermite and Laguerre exceptional orthogonal polynomials”. J. Math. Phys. 54 (2013), 042102.
  * [23] I. Marquette and C. Quesne, “New ladder operators for a rational extension of the harmonic oscillator and superintegrability of some two-dimensional systems”. J. Math. Phys. 54 (2013), 102102.
  * [24] I. Marquette and C. Quesne, “Combined state-adding and state-deleting approaches to type III multi-step rationally extended potentials: Applications to ladder operators and superintegrability”. J. Math. Phys. 55 (2014), 112103.
  * [25] D. Gómez-Ullate, N. Kamran and R. Milson, “An extended class of orthogonal polynomials defined by a Sturm–Liouville problem”. J. Math. Anal. Appl. 359 (2009), 352.
  * [26] D. Gómez-Ullate, N. Kamran and R. Milson, “An extension of Bochner’s problem: Exceptional invariant subspaces”. J. Approx. Theory 162 (2010), 987.
  * [27] S. Odake and R. Sasaki, “Infinitely many shape invariant potentials and new orthogonal polynomials”. Phys. Lett. B 679 (2009), 414.
  * [28] S. Odake and R. Sasaki, “Another set of infinitely many exceptional ($X_{\ell}$) Laguerre polynomials”. Phys. Lett. B 684 (2010), 173.
  * [29] C. Quesne, “Exceptional orthogonal polynomials, exactly solvable potentials and supersymmetry”. J. Phys. A: Math. Theor. 41 (2008), 392001.
  * [30] C. Quesne, “Solvable Rational Potentials and Exceptional Orthogonal Polynomials in Supersymmetric Quantum Mechanics”. SIGMA 5 (2009), 084.
  * [31] C.-L. Ho, “Prepotential approach to solvable rational extensions of Harmonic Oscillator and Morse potentials”. J. Math. Phys. 52 (2011), 122107.
  * [32] Y. Grandati, “Solvable rational extensions of the isotonic oscillator”. Ann. Phys. 326 (2011), 2074.
  * [33] Y. Grandati and C. Quesne, “Disconjugacy, regularity of multi-indexed rationally extended potentials, and Laguerre exceptional polynomials”. J. Math. Phys. 54 (2013), 073512.
  * [34] D. Gómez-Ullate, N. Kamran and R. Milson, “The Darboux transformation and algebraic deformations of shape-invariant potentials”. J. Phys. A: Math. Gen. 37 (2004), 1789.
  * [35] J. F. Cariñena, A. M. Perelomov, M. F. Rañada and M. Santander, “A quantum exactly solvable nonlinear oscillator related to the isotonic oscillator”. J. Phys. A: Math. Theor. 41 (2008), 085301.
  * [36] J. M. Fellows and R. A. Smith, “Factorization solution of a family of quantum nonlinear oscillators”. J. Phys. A: Math. Theor. 42 (2009), 335303.
  * [37] Y. Grandati, “Solvable rational extensions of the Morse and Kepler-Coulomb potentials”. J. Math. Phys. 52 (2011), 103505.
  * [38] Y. Grandati, “Multistep DBT and regular rational extensions of the isotonic oscillator”. Ann. Phys. 327 (2012), 2411.
  * [39] S. Yu. Dubov, V. M. Eleonskii and N. E. Kulagin, “Equidistant spectra of anharmonic oscillators”. Chaos 4 (1994), 47.
  * [40] B. F. Samsonov, “New features in supersymmetry breakdown in quantum mechanics”. Mod. Phys. Lett. A 11 (1996), 1563.
  * [41] Md. F. Hoque, I. Marquette, S. Post and Y.-Z. Zhang, “Algebraic calculations for spectrum of superintegrable system from exceptional orthogonal polynomials”. Ann. Phys. 391 (2018), 203.
  * [42] S. E. Hoffmann., V. Hussin, I. Marquette and Y.-Z. Zhang, “Coherent states for ladder operators of general order related to exceptional orthogonal polynomials”. J. Phys. A: Math. Theor. 51 (2018), 315203.
  * [43] M. S. Plyushchay, Nonlinear Supersymmetry as a Hidden Symmetry, p. 163.  Springer International Publishing, Cham, 2019.
  * [44] L. Inzunza and M. S. Plyushchay, “Hidden symmetries of rationally deformed superconformal mechanics”. Phys. Rev. D 99 (2019), 025001.
  * [45] V. Chalifour and A. M. Grundland, “General solution of the exceptional Hermite differential equation and its minimal surface representation”. Ann. Henri Poincaré 21 (2020), 3341.
  * [46] D. J. Fernández C and V. Hussin, “Higher-order SUSY, linearized nonlinear Heisenberg algebras and coherent states”. J. Phys. A: Math. Gen. 32 (1999), 3603.
  * [47] J. M. Carballo, D. J. Fernández C, J. Negro and L. M. Nieto, “Polynomial Heisenberg algebras”. J. Phys. A: Math. Gen. 37 (2004), 10349.
  * [48] I. Marquette, M. Sajedi and P. Winternitz, “Two-dimensional superintegrable systems from operator algebras in one dimension”. J. Phys. A: Math. Theor. 52 (2019), 115202.
  * [49] W. Greiner, B. Müller, Quantum Mechanics: Symmetries.  Springer-Verlag, 1994.
  * [50] D. M. Fradkin, “Three-Dimensional Isotropic Harmonic Oscillator and ${SU}3$”. Amer. J. Phys. 33 (1965), 207.
  * [51] B. G. Wybourne, Classical Groups for Physicists.  Wiley, 1974.
  * [52] H. Georgi, Lie Algebras In Particle Physics: from Isospin To Unified Theories.  Avalon Publishing, 1999.
  * [53] Yu. N. Demkov, “Symmetry Group of the Isotropic Oscillator”. J. Exptl. Theoret. Phys. 36 (1959), 63.
  * [54] B. K. Bagchi, Supersymmetry in quantum and classical mechanics.  Chapman and Hall/CRC, 2001.
  * [55] F. Cooper, A. Khare and U. Sukhatme, “Supersymmetry and quantum mechanics”. Phys. Rep. 251 (1995), 267.
  * [56] Y. Tanoudis and C. Daskaloyannis, “Algebraic Calculation of the Energy Eigenvalues for the Nondegenerate Three-Dimensional Kepler-Coulomb Potential”. SIGMA 7 (2011), 054.
  * [57] D. J. Fernández C and N. Fernández‐García, “Higher‐order supersymmetric quantum mechanics”. AIP Conf. Proc. 744 (2004), 236.
  * [58] D. Bermúdez, “Supersymmetric Quantum Mechanics and Painlevé IV Equation”. SIGMA 7 (2011), 25.
  * [59] J. F. Cariñena and M. S. Plyushchay, “Ground-state isolation and discrete flows in a rationally extended quantum harmonic oscillator”. Phys. Rev. D 94 (2016), 105022.
  * [60] J. F. Cariñena and M. S. Plyushchay, “ABC of ladder operators for rationally extended quantum harmonic oscillator systems”. J. Phys. A: Math. Theor. 50 (2017), 275202.
  * [61] C. Quesne, “Ladder operators for solvable potentials connected with exceptional orthogonal polynomials”. J. Phys.: Conf. Ser. 597 (2015), 012064.
  * [62] I. Marquette, “Superintegrability with third order integrals of motion, cubic algebras, and supersymmetric quantum mechanics. II. Painlevé transcendent potentials”. J. Math. Phys. 50 (2009), 095202.
  * [63] I. Marquette, “Higher order superintegrability, Painlevé transcendents and representations of polynomial algebras”. J. Phys.: Conf. Ser. 1194 (2019), 012074.
  * [64] I. Marquette and P. Winternitz, Higher Order Quantum Superintegrability: A New “Painlevé Conjecture”, p. 103.  Springer International Publishing, Cham, 2019.
  * [65] C. M. Cosgrove, “Higher‐order Painlevé Equations in the Polynomial Class I. Bureau Symbol P2”. Stud. Appl. Math. 104 (2000), 1.
  * [66] C. M. Cosgrove, “Higher-Order Painlevé Equations in the Polynomial Class II: Bureau Symbol P1”. Stud. Appl. Math. 116 (2000), 321.
  * [67] A. Polishchuk and L. Positselski, Quadratic Algebras.  American Mathematical Society, 2005.
  * [68] H. Li, Gröbner Bases in Ring Theory.  World Scientific, 2011.
  * [69] P. D. Jarvis, G. Rudolph and L. A. Yates, “A class of quadratic deformations of Lie superalgebras”. J. Phys. A: Math. Theor. 44 (2011), 235205.
  * [70] H. Li, Noncommutative Polynomial Algebras of Solvable Type and Their Modules.  Chapman and Hall/CRC, 2021.
  * [71] L. A. Yates and P. D. Jarvis, “Hidden supersymmetry and quadratic deformations of the space-time conformal superalgebra”. J. Phys. A: Math. Theor. 51 (2018), 145203.
  * [72] F. A. Grünbaum, Luc Vinet and A. Zhedanov, “Algebraic Heun Operator and Band-Time Limiting”. Comm. Math. Phys. 364 (2018), 1041.
  * [73] N. Crampe, L. Poulain d’Andecy and L. Vinet, “A Calabi-Yau Algebra with E6 Symmetry and the Clebsch-Gordan Series of $\mathfrak{sl}$(3)”. J. Lie Theory 31 (2021), 1085.
  * [74] N. Crampé, L. Vinet and A. Zhedanov, “Heun algebras of Lie type”. Proc. Amer. Math. Soc. 148 (2019), 1079.