$\mathfrak{gl}(3)$ Polynomial Integrable System: Different Faces

$\boldsymbol{\mathfrak{gl}(3)}$ Polynomial Integrable System: Different
Faces of the 3-Body/$\boldsymbol{{\mathcal{A}}_{2}}$ Elliptic Calogero
Model††This paper is a contribution to the Special Issue on Symmetry,
Invariants, and their Applications in honor of Peter J. Olver. The full
collection is available at https://www.emis.de/journals/SIGMA/Olver.html

Alexander V. TURBINER, Juan Carlos LOPEZ VIEYRA
and Miguel A. GUADARRAMA-AYALA

A.V. Turbiner, J.C. Lopez Vieyra and M.A. Guadarrama-Ayala

Instituto de Ciencias Nucleares, Universidad Nacional Autónoma de México,
Apartado Postal 70-543, 04510 Ciudad de México, Mexico
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
Received July 26, 2023, in final form January 22, 2024; Published online
February 03, 2024

It is shown that the $\mathfrak{gl}(3)$ polynomial integrable system,
introduced by Sokolov–Turbiner in [J. Phys. A 48 (2015), 155201, 15 pages,
arXiv:1409.7439], is equivalent to the $\mathfrak{gl}(3)$ quantum Euler–Arnold
top in a constant magnetic field. Their Hamiltonian as well as their third-
order integral can be rewritten in terms of $\mathfrak{gl}(3)$ algebra
generators. In turn, all these $\mathfrak{gl}(3)$ generators can be
represented by the non-linear elements of the universal enveloping algebra of
the 5-dimensional Heisenberg algebra
$\mathfrak{h}_{5}(\hat{p}_{1,2},\hat{q}_{1,2},I)$, thus, the Hamiltonian and
integral are two elements of the universal enveloping algebra
$U_{\mathfrak{h}_{5}}$. In this paper, four different representations of the
$\mathfrak{h}_{5}$ Heisenberg algebra are used: (I) by differential operators
in two real (complex) variables, (II) by finite-difference operators on
uniform or exponential lattices. We discovered the existence of two
2-parametric bilinear and trilinear elements (denoted $H$ and $I$,
respectively) of the universal enveloping algebra $U(\mathfrak{gl}(3))$ such
that their Lie bracket (commutator) can be written as a linear superposition
of nine so-called artifacts – the special bilinear elements of
$U(\mathfrak{gl}(3))$, which vanish once the representation of the
$\mathfrak{gl}(3)$-algebra generators is written in terms of the
$\mathfrak{h}_{5}(\hat{p}_{1,2},\hat{q}_{1,2},I)$-algebra generators. In this
representation all nine artifacts vanish, two of the above-mentioned elements
of $U(\mathfrak{gl}(3))$ (called the Hamiltonian $H$ and the integral $I$)
commute(!); in particular, they become the Hamiltonian and the integral of the
3-body elliptic Calogero model, if $(\hat{p},\hat{q})$ are written in the
standard coordinate-momentum representation. If $(\hat{p},\hat{q})$ are
represented by finite-difference/discrete operators on uniform or exponential
lattice, the Hamiltonian and the integral of the 3-body elliptic Calogero
model become the isospectral, finite-difference operators on uniform-uniform
or exponential-exponential lattices (or mixed) with polynomial coefficients.
If $(\hat{p},\hat{q})$ are written in complex $(z,\bar{z})$ variables the
Hamiltonian corresponds to a complexification of the 3-body elliptic Calogero
model on ${\mathbb{C}^{2}}$.

elliptic Calogero model; integrable systems; 3-body systems

81R12; 81S05; 17J35; 81U15

##  1 Introduction

Let us take a finite-dimensional Lie algebra $\mathfrak{g}$ spanned by the
generators $J_{i}$, $i=1,2,\dots,\dim\mathfrak{g}$. The second degree
polynomial in the $J$-generators,

$H(J)=\sum_{i,j=1}^{\dim\mathfrak{g}}a_{ij}\\{J_{i},J_{j}\\}+\sum_{i}^{\dim\mathfrak{g}}b_{i}J_{i},$

where $\\{A,B\\}=AB+BA$ is the anti-commutator and $\\{a\\}$, $\\{b\\}$ are
parameters, defines the Hamiltonian of the quantum Euler–Arnold top in a
constant magnetic field with components $b_{i}$,
$i=1,2,\dots,\dim\mathfrak{g}$. It is well known that the generators $J_{i}$
of any semi-simple Lie algebra can be written in terms of the generators
$(\hat{p},\hat{q})$ of a Heisenberg algebra, hence,
$J_{i}=J_{i}(\hat{p},\hat{q})$. We call such a system a
$\mathfrak{g}$-polynomial system if its Hamiltonian is defined as

$H(\hat{p},\hat{q})\equiv H(J(\hat{p},\hat{q})).$

A particular example of an $\mathfrak{sl}(2)$-polynomial system was studied in
details in [10, equation (13)], which is associated with the harmonic
oscillator,

$H=-\hat{q}\hat{p}^{2}+\bigl{(}\hat{q}-p-\tfrac{1}{2}\bigr{)}\hat{p}=-J^{0}J^{-}+J^{0}-\bigl{(}p+\tfrac{1}{2}\bigr{)}J^{-},$

where $p=0,1$ and

$J^{0}=\hat{q}\hat{p},\qquad J^{-}=\hat{p},$

are two $\mathfrak{sl}(2)$ generators, $\big{[}J^{0},J^{-}\big{]}=-J^{-}$, see
below. The general $\mathfrak{sl}(2)$-polynomial system is associated with the
Heun operator, which is equivalent to the $BC_{1}$ elliptic Calogero model
[11]. The present paper is aimed at constructing an analogous but
$\mathfrak{gl}(3)$-polynomial system starting from the quantum
${\mathcal{A}}_{2}$ elliptic (3-body Calogero) model.

Celebrated 3-body elliptic Calogero model or, stated differently, the
${\mathcal{A}}_{2}$ elliptic model (in the Hamiltonian reduction nomenclature,
see, e.g., [3]), describes three point-like one-dimensional particles of unit
masses on the real line with pairwise interaction given by the Weierstrass
$\wp$-function with rectangular fundamental domain $(\omega,{\rm
i}\omega_{1})$. It is characterized by the Hamiltonian

$\displaystyle{\mathcal{H}}^{(e)}_{{\mathcal{A}}_{2}}=-\frac{1}{2}\sum_{i=1}^{3}\frac{\partial^{2}}{\partial
x_{i}^{2}}+\nu(\nu-1)(\wp(x_{1}-x_{2})+\wp(x_{2}-x_{3})+\wp(x_{3}-x_{1}))$
$\displaystyle\hphantom{{\mathcal{H}}^{(e)}_{{\mathcal{A}}_{2}}}{}\equiv-\frac{1}{2}\Delta^{(3)}+V_{{\mathcal{A}}_{2}},$
(1.1)

where $x_{1}$, $x_{2}$, $x_{3}$ are the coordinates of the bodies,
$\Delta^{(3)}$ is three-dimensional flat Laplace operator, which represents
the kinetic energy, $\kappa\equiv\nu(\nu-1)$ is the coupling constant. The
Weierstrass function $\wp(x)\equiv\wp(x\mid g_{2},g_{3})$ (see, e.g., [14]) is
defined as the solution of the equation

$\bigl{(}\wp^{\prime}(x)\bigr{)}^{2}=4\wp^{3}(x)-g_{2}\wp(x)-g_{3}=4(\wp(x)-e_{1})(\wp(x)-e_{2})(\wp(x)-e_{3}),$

where $g_{2}$, $g_{3}$ are the so-called elliptic invariants, which can be
conveniently parameterized as follows:

$g_{2}=12\bigl{(}\tau^{2}-\mu\bigr{)},\qquad
g_{3}=4\tau\bigl{(}2\tau^{2}-3\mu\bigr{)},$ (1.2)

where $\tau$, $\mu$ are parameters, and $e_{1}$, $e_{2}$, $e_{3}$ are its
roots which are chosen, conventionally, to obey $e\equiv e_{1}+e_{2}+e_{3}=0$.
Since the Hamiltonian (1.1) is translation-invariant, $x\rightarrow x+a$, one
can introduce the center-of-mass and relative coordinates,

$Y=\sum_{1}^{3}x_{i},\qquad y_{i}=x_{i}-\frac{1}{3}Y,$ (1.3)

with the condition $\sum_{1}^{3}y_{i}=0$. The Laplacian
$\Delta^{(3)}\equiv\sum_{i=1}^{3}\frac{\partial^{2}}{\partial x_{i}^{2}}$ in
these coordinates takes the form,

$\Delta^{(3)}=3\partial_{Y}^{2}+\frac{2}{3}\biggl{(}\frac{\partial^{2}}{\partial
y_{1}^{2}}+\frac{\partial^{2}}{\partial
y_{2}^{2}}-\frac{\partial^{2}}{\partial y_{1}\partial y_{2}}\biggr{)}.$

Separating out the center-of-mass coordinate $Y$, the two-dimensional
Hamiltonian arises

$\displaystyle{\mathcal{H}}_{{\mathcal{A}}_{2}}=-\frac{1}{3}\biggl{(}\frac{\partial^{2}}{\partial
y_{1}^{2}}+\frac{\partial^{2}}{\partial
y_{2}^{2}}-\frac{\partial^{2}}{\partial y_{1}\partial y_{2}}\biggr{)}$
$\displaystyle\hphantom{{\mathcal{H}}_{{\mathcal{A}}_{2}}=}{}+\nu(\nu-1)(\wp(y_{1}-y_{2})+\wp(2y_{1}+y_{2})+\wp(y_{1}+2y_{2})),$
(1.4)

which seemingly was already known to Charles Hermite as a two-dimensional
generalization of the celebrated one-dimensional Lamé operator (following
Serguei P. Novikov’s studies of unpublished notes by Charles Hermite
communicated to one of the authors (Alexander V. Turbiner)),

${\mathcal{H}}^{(e)}_{{\mathcal{A}}_{1}}=-\frac{1}{2}\frac{\partial^{2}}{\partial
y^{2}}+\kappa\wp(y),$

which is also the Hamiltonian of the ${\mathcal{A}}_{1}$ elliptic model [3],
see also [8]. We will call the operator (1.4) the two-dimensional Lamé
operator. In general, the above procedure allows us to connect the quantum
dynamics in the relative space of the three-body problem with two-dimensional
quantum dynamics [12].

For many years the question of the existence of polynomial eigenfunctions of
the operator (1.4) was a challenge to answer. It was eventually solved in 2015
by Sokolov–Turbiner in [6]: the discrete values of the coupling constant were
found

$\kappa\equiv\nu(\nu-1)=\frac{n}{9}(n+3),\qquad n=0,1,2,\dots,$ (1.5)

for which the $\frac{(n+2)(n+1)}{2}$ polynomial eigenfunctions exist in the
variables

$x=\frac{f^{\prime}(y_{1})-f^{\prime}(y_{2})}{f(y_{1})f^{\prime}(y_{2})-f(y_{2})f^{\prime}(y_{1})},\qquad
y=\frac{2(f(y_{1})-f(y_{2}))}{f(y_{1})f^{\prime}(y_{2})-f(y_{2})f^{\prime}(y_{1})},$
(1.6)

where

$f(x)=\wp(x|g_{2},g_{3})+\tau,$

is the shifted Weierstrass function.

In very tedious and highly non-trivial calculations, performed in [6], it was
found that the ${\mathcal{A}}_{2}$ elliptic Calogero model potential
$V_{{\mathcal{A}}_{2}}$ (see (1.1) and (1.4)) in variables (1.6) takes the
form of ratio of polynomials,

$V_{{\mathcal{A}}_{2}}(x,y)={3\nu(\nu-1)}\frac{\bigl{(}x+2\tau x^{2}+\mu
x^{3}-6\bigl{(}\mu-\tau^{2}\bigr{)}y^{2}+3\mu\tau xy^{2}\bigr{)}^{2}}{4D},$

where the denominator

$\displaystyle
4D(x,y)=3\mu^{2}x^{4}y^{2}+18\tau\mu^{2}x^{2}y^{4}+9\mu^{2}(3\tau^{2}-4\mu)y^{6}-4\mu
x^{5}-24\tau\mu x^{3}y^{2}$
$\displaystyle\hphantom{4D(x,y)=}{}-36\mu(\tau^{2}-2\mu)xy^{4}-4\tau
x^{4}-6(4\tau^{2}+5\mu)x^{2}y^{2}-18\tau(2\tau^{2}-3\mu)y^{4}$
$\displaystyle\hphantom{4D(x,y)=}{}-36\tau xy^{2}-\frac{4}{3}x^{3}-27y^{2}$
(1.7)

was called the determinant or the elliptic discriminant. In rational limit
$\tau=\mu=0$ this is the square of Vandermonde determinant or the discriminant
of the cubic equation, in trigonometric limit $\mu=0$ it becomes the so-called
trigonometric Vandermonde determinant or a trigonometric discriminant, see for
details [6]. Furthermore, the two-dimensional flat Laplacian in (1.4) becomes
the Laplace–Beltrami operator in $(x,y)$-coordinates

$\displaystyle\Delta_{g}(x,y;\tau,\mu)=3\biggl{(}\frac{x}{3}+\tau x^{2}+\mu
x^{3}+\bigl{(}\mu-\tau^{2}\bigr{)}y^{2}-\mu\tau
xy^{2}-\mu^{2}x^{2}y^{2}\biggr{)}\frac{\partial^{2}}{\partial x^{2}}$
$\displaystyle\hphantom{\Delta_{g}(x,y;\tau,\mu)=}{}+y\bigl{(}3+8\tau x+7\mu
x^{2}-3\mu\tau y^{2}-6\mu^{2}xy^{2}\bigr{)}\frac{\partial^{2}}{\partial
x\partial y}$
$\displaystyle\hphantom{\Delta_{g}(x,y;\tau,\mu)=}{}+\biggl{(}-\frac{x^{2}}{3}+3\tau
y^{2}+4\mu xy^{2}-3\mu^{2}y^{4}\biggr{)}\frac{\partial^{2}}{\partial y^{2}}$
$\displaystyle\hphantom{\Delta_{g}(x,y;\tau,\mu)=}{}+\bigl{(}1+4\tau x+5\mu
x^{2}-3\mu\tau y^{2}-6\mu^{2}xy^{2}\bigr{)}\frac{\partial}{\partial
x}+2y\bigl{(}2\tau+3\mu x-3\mu^{2}y^{2}\bigr{)}\frac{\partial}{\partial y},$

with naturally-defined flat contravariant metric $g^{ij}$, $i,j=1,2$ with
polynomial entries. It can be easily checked that, remarkably, expression
(1.7) is equal to the determinant of this contravariant metric,

$D={\rm Det}\bigl{(}g^{ij}\bigr{)},$

which explains the name determinant, used in [6].

Surprisingly, the gauge rotation of the 2-dimensional Lamé operator (1.4) with
the determinant $D$ (1.7) to the power $\nu/2$ as a gauge factor transforms
operator (1.4) into the algebraic operator(!) with polynomial coefficients,

$\displaystyle
h_{{\mathcal{A}}_{2}}(x,y)=-3D^{-\frac{\nu}{2}}({\mathcal{H}}_{\rm{\mathcal{A}}_{2}}-3\nu(3\nu+1)\tau)D^{\frac{\nu}{2}}$
$\displaystyle\hphantom{h_{{\mathcal{A}}_{2}}(x,y)}{}=\bigl{(}x+3\tau
x^{2}+3\mu x^{3}+3\bigl{(}\mu-\tau^{2}\bigr{)}y^{2}-3\mu\tau
xy^{2}-3\mu^{2}x^{2}y^{2}\bigr{)}\frac{\partial^{2}}{\partial x^{2}}$
$\displaystyle\hphantom{h_{{\mathcal{A}}_{2}}(x,y)=}{}+y\bigl{(}3+8\tau x+7\mu
x^{2}-3\mu\tau y^{2}-6\mu^{2}xy^{2}\bigr{)}\frac{\partial^{2}}{\partial
x\partial y}$
$\displaystyle\hphantom{h_{{\mathcal{A}}_{2}}(x,y)=}{}+\frac{1}{3}\bigl{(}-x^{2}+9\tau
y^{2}+12\mu xy^{2}-9\mu^{2}y^{4}\bigr{)}\frac{\partial^{2}}{\partial y^{2}}$
$\displaystyle\hphantom{h_{{\mathcal{A}}_{2}}(x,y)=}{}+(1+3\nu)\bigl{(}1+4\tau
x+5\mu x^{2}-3\mu\tau y^{2}-6\mu^{2}xy^{2}\bigr{)}\frac{\partial}{\partial x}$
$\displaystyle\hphantom{h_{{\mathcal{A}}_{2}}(x,y)=}{}+2(1+3\nu)y\bigl{(}2\tau+3\mu
x-3\mu^{2}y^{2}\bigr{)}\frac{\partial}{\partial
y}+3\nu(1+3\nu)\mu\bigl{(}2x-3\mu y^{2}\bigr{)}.$ (1.8)

This was one of the principal results obtained in the article [6], which will
be essential in the present article. Let us emphasize that the operator
$h_{{\mathcal{A}}_{2}}(x,y)$ looks like the two-dimensional generalization of
the (algebraic) Heun operator, see, e.g., [11].

It was also found in [6] that the second-order algebraic differential operator
$h_{{\mathcal{A}}_{2}}(x,y)$ commutes with a non-trivial third-order algebraic
differential operator $k_{\rm{\mathcal{A}}_{2}}$ with polynomial coefficients,

$[h_{{\mathcal{A}}_{2}}(x,y),k_{\rm{\mathcal{A}}_{2}}(x,y)]=0,$

where

$\displaystyle k_{\rm{\mathcal{A}}_{2}}(x,y)=2\nu(1+3\nu)(2+3\nu)\mu
y\bigl{(}2\tau+3\mu x-3\mu^{2}y^{2}\bigr{)}$
$\displaystyle\hphantom{k_{\rm{\mathcal{A}}_{2}}(x,y)=}{}+\frac{1}{3}(1+3\nu)(2+3\nu)y\bigl{(}\mu+8\tau^{2}+28\mu\tau
x+21\mu^{2}x^{2}-9\mu^{2}\tau
y^{2}-18\mu^{3}xy^{2}\bigr{)}\frac{\partial}{\partial x}$
$\displaystyle\hphantom{k_{\rm{\mathcal{A}}_{2}}(x,y)=}{}-\frac{2}{9}(1+3\nu)(2+3\nu)\bigl{(}1+4\tau
x+6\mu x^{2}-24\mu\tau
y^{2}-36\mu^{2}x{y}^{2}+27\mu^{3}y^{4}\bigr{)}\frac{\partial}{\partial y}$
$\displaystyle\hphantom{k_{\rm{\mathcal{A}}_{2}}(x,y)=}{}+(2+3\nu)y\bigl{(}3\tau+4\bigl{(}2\tau^{2}+\mu\bigr{)}x+17\mu\tau
x^{2}+8\mu^{2}x^{3}$
$\displaystyle\hphantom{k_{\rm{\mathcal{A}}_{2}}(x,y)=}{}\quad{}+3\mu({\tau}^{2}-2\mu)y^{2}-6\mu^{2}\tau
xy^{2}-6{\mu}^{3}x^{2}y^{2}\bigr{)}\frac{\partial^{2}}{\partial x^{2}}$
$\displaystyle\hphantom{k_{\rm{\mathcal{A}}_{2}}(x,y)=}{}-\frac{2}{3}(2+3\nu)\bigl{(}x+4\tau
x^{2}+5\mu x^{3}+3\bigl{(}\mu-4\tau^{2}\bigr{)}y^{2}-27\mu^{2}x^{2}y^{2}$
$\displaystyle\hphantom{k_{\rm{\mathcal{A}}_{2}}(x,y)=}{}\quad{}-33\mu\tau
xy^{2}+9\mu^{2}\tau y^{4}+18\mu^{3}xy^{4}\bigr{)}\frac{\partial^{2}}{\partial
x\partial y}$
$\displaystyle\hphantom{k_{\rm{\mathcal{A}}_{2}}(x,y)=}{}-(2+3\nu)y\biggl{(}1+\frac{8}{3}\tau
x+3\mu x^{2}-7\mu\tau
y^{2}-10\mu^{2}xy^{2}+6\mu^{3}y^{4}\biggr{)}\frac{\partial^{2}}{\partial
y^{2}}$
$\displaystyle\hphantom{k_{\rm{\mathcal{A}}_{2}}(x,y)=}{}+y\bigl{(}1+5\tau
x+2\bigl{(}2\mu+3\tau^{2}\bigr{)}x^{2}+3\mu(\tau^{2}-2\mu)xy^{2}+9\mu\tau
x^{3}$
$\displaystyle\hphantom{k_{\rm{\mathcal{A}}_{2}}(x,y)=}{}\quad{}-\tau\bigl{(}3\mu-2\tau^{2}\bigr{)}y^{2}+3\mu^{2}x^{4}-3\mu^{2}\tau
x^{2}y^{2}-2\mu^{3}x^{3}y^{2}\bigr{)}\frac{\partial^{3}}{\partial x^{3}}$
$\displaystyle\hphantom{k_{\rm{\mathcal{A}}_{2}}(x,y)=}{}+\biggl{(}-\frac{2}{3}x^{2}+2\bigl{(}5{\tau}^{2}+\mu\bigr{)}xy^{2}-2\tau
x^{3}+3\tau y^{2}-2\mu x^{4}+3\mu\bigl{(}\tau^{2}-2\mu\bigr{)}y^{4}$
$\displaystyle\hphantom{k_{\rm{\mathcal{A}}_{2}}(x,y)=}{}\quad{}+19\mu\tau
x^{2}y^{2}-6{\mu}^{3}x^{2}y^{4}+10\mu^{2}x^{3}y^{2}-6\mu^{2}\tau
xy^{4}\biggr{)}\frac{\partial^{3}}{\partial x^{2}\partial y}$
$\displaystyle\hphantom{k_{\rm{\mathcal{A}}_{2}}(x,y)=}{}-y\biggl{(}x+\frac{10}{3}\tau{x}^{2}+\frac{11}{3}\mu{x}^{3}-13\mu\tau
xy^{2}+3\bigl{(}\mu-2\tau^{2}\bigr{)}y^{2}-11\mu^{2}x^{2}y^{2}$
$\displaystyle\hphantom{k_{\rm{\mathcal{A}}_{2}}(x,y)=}{}\quad{}+3\mu^{2}\tau
y^{4}+6\mu^{3}xy^{4}\biggr{)}\frac{\partial^{3}}{\partial x\partial y^{2}}$
$\displaystyle\hphantom{k_{\rm{\mathcal{A}}_{2}}(x,y)=}{}-\biggl{(}y^{2}+\frac{2}{27}x^{3}+2\tau
xy^{2}-3\mu\tau y^{4}+\frac{5}{3}\mu
x^{2}y^{2}-4{\mu}^{2}xy^{4}+2{\mu}^{3}y^{6}\biggr{)}\frac{\partial^{3}}{\partial
y^{3}}.$ (1.9)

Hence, $h_{{\mathcal{A}}_{2}}(x,y)$ and $k_{{\mathcal{A}}_{2}}(x,y)$ span the
two-dimensional commutative algebra of the differential operators in two
variables, which depend on three free parameters $\nu$, $\mu$, $\tau$. This is
the first non-trivial example of this. Naturally, the third-order differential
operator $k_{{\mathcal{A}}_{2}}(x,y)$ can be called the integral. By making
the inverse gauge rotation of the integral $k_{{\mathcal{A}}_{2}}(x,y)$,

$D^{\frac{\nu}{2}}k_{{\mathcal{A}}_{2}}(x,y)D^{-\frac{\nu}{2}},$

with the determinant $D$ (1.7) to the power $(-\nu/2)$ as the gauge factor and
changing variables $(x,y)\rightarrow(y_{1},y_{2})$ (1.6), we should arrive at
the third-order integral of the quantum 3-body elliptic Calogero model in the
form of the third-order differential operator with elliptic coefficients found
by Oshima [4]. This demonstrates explicitly the integrability of the original
3-body elliptic Calogero model written in $y_{1}$, $y_{2}$ coordinates.

It was concluded in [6] that the 3-body elliptic Calogero model defines a
polynomial integrable model with the algebraic Hamiltonian (1.8) and the
algebraic integral (1.9) with $\mu,\tau,\nu$-dependent parametric
coefficients. This model has $\mathfrak{sl}(3)$ hidden algebra in the
representation $(-3\nu,0)$. As a result the $\mathfrak{sl}(3)$ quantum
Euler–Arnold top in a constant magnetic field occurs. Note that for discrete
values of the coupling constant $\kappa$: $n=-3\nu$, $n=0,1,2,\dots$, the
$\mathfrak{sl}(3)$ hidden algebra emerges in the finite-dimensional
representation, thus, the top has a common finite-dimensional invariant
subspace for both $h$ and $k$.

The goal of this article is two-fold. First of all, the above-mentioned
polynomial integrable model, realized in terms of differential operators, will
be rewritten in terms of the generators of the Heisenberg algebra
$\mathfrak{h}_{5}$. Hence, its Hamiltonian will appear as an element of the
universal enveloping algebra $U_{\mathfrak{h}_{5}}$. Then we project it into
the translation-invariant or dilatation-invariant operators defining two
families of 3-parametric $\mu$, $\tau$, $\nu$ isospectral polynomial
integrable models on two-dimensional uniform or exponential lattices,
respectively, and two additional families on mixed two-dimensional
translation-invariant and dilatation-invariant lattices. All four families
admit 2-parametric $\mu$, $\tau$ polynomial eigenfunctions for certain
discrete values of the coupling constant $\nu$. An extra polynomial integrable
model occurs as a result of a special complexification of ${\mathbb{R}}^{2}$
to ${\mathbb{C}}^{2}$ via the Heisenberg algebra $\mathfrak{h}_{5}$ generators
acting on the two-dimensional Hilbert space with the Gaussian measure. The
spectrum of this model is characterized by infinite multiplicity and for
certain discrete values of the coupling constant $\kappa$ (1.5) the
eigenfunctions are poly-analytic functions in two complex variables of the
fixed degree. Second of all, it will be shown that $\mathfrak{gl}(3)$
polynomial integrable model, defined in the Fock space, is related with
special bilinear and trilinear, 2-parametric elements of the universal
enveloping algebra of the algebra $\mathfrak{gl}(3)$. It turns out that these
non-linear elements commute once they are written in terms of any concrete
realization of the algebra $\mathfrak{gl}(3)$ by elements of the universal
enveloping algebra $U_{\mathfrak{h}_{5}}$.

The article is organized with introduction, Sections 2–7, conclusions and
three Appendices A, B, C. In Section 2, the 3-body elliptic Calogero model in
algebraic form is reformulated in Fock space and its
$\mathfrak{gl}(3)$-polynomial integrable model is defined. Section 3 contains
four lattice versions of the 3-body elliptic Calogero model. Section 4 is
dedicated to complexification of the $\mathfrak{gl}(3)$-polynomial integrable
model into ${\mathbb{C}}^{2}$. In Section 5, all nine artifacts of the
$\mathfrak{gl}(3)$ algebra are presented as bilinear combinations of the
$\mathfrak{gl}(3)$ generators and Theorem is proved that all of them vanish if
the $\mathfrak{gl}(3)$ generators are written as non-linear elements of the
universal enveloping algebra $U_{\mathfrak{h}_{5}}$. Section 6 contains the
explicit expressions of the Hamiltonian, the cubic integral and their
commutator in terms of the $\mathfrak{gl}(3)$-algebra generators. In Section
7, the $G_{2}$/3-body (with pairwise and 3-body interactions) elliptic problem
is briefly discussed and the Fock space representation of the $G_{2}$ elliptic
3-body Hamiltonian is constructed.

Throughout the remaining text the hats in $p,q$’s will be dropped:
$({\hat{p}},{\hat{q}})\rightarrow(p,q)$.

##  2 3-body elliptic Calogero model in the Fock space

Let us take 5-dimensional Heisenberg algebra $\mathfrak{h}_{5}$ spanned by the
generators $p_{x}$, $p_{y}$, $q_{x}$, $q_{x}$, $I$, which obey the commutation
relations,

$\displaystyle[p_{x},q_{x}]=1,\qquad$ $\displaystyle[p_{y},q_{y}]=1,\qquad$
$\displaystyle[p_{x},q_{y}]=0,\qquad$ $\displaystyle[p_{y},q_{x}]=0,$
$\displaystyle[p_{x},p_{y}]=0,\qquad$ $\displaystyle[q_{x},q_{y}]=0,\qquad$
$\displaystyle[p_{x,y},I]=0,\qquad$ $\displaystyle[q_{x,y},I]=0,$ (2.1)

see Appendix A.3. The universal enveloping algebra of the algebra
$\mathfrak{h}_{5}$: $U_{\mathfrak{h}_{5}}$, is spanned by all ordered
monomials in $p_{x}$, $p_{y}$, $q_{x}$, $q_{y}$.

Now let us form in $U_{\mathfrak{h}_{5}}$ a second degree polynomial in
$p$-generators,

$\displaystyle h_{{\mathcal{A}}_{2}}(p_{x},q_{x},p_{y},q_{y})$
$\displaystyle{}=\bigl{(}q_{x}+3\tau q_{x}^{2}+3\mu
q_{x}^{3}+3\bigl{(}\mu-\tau^{2}\bigr{)}q_{y}^{2}-3\mu\tau
q_{x}q_{y}^{2}-3\mu^{2}q_{x}^{2}q_{y}^{2}\bigr{)}p_{x}^{2}$
$\displaystyle\quad{}+q_{y}\bigl{(}3+8\tau q_{x}+7\mu q_{x}^{2}-3\mu\tau
q_{y}^{2}-6\mu^{2}q_{x}q_{y}^{2}\bigr{)}p_{x}p_{y}$
$\displaystyle\quad{}+\frac{1}{3}\bigl{(}-q_{x}^{2}+9\tau q_{y}^{2}+12\mu
q_{x}q_{y}^{2}-9\mu^{2}q_{y}^{4}\bigr{)}p_{y}^{2}$
$\displaystyle\quad{}+(1+3\nu)\bigl{(}1+4\tau q_{x}+5\mu q_{x}^{2}-3\mu\tau
q_{y}^{2}-6\mu^{2}q_{x}q_{y}^{2}\bigr{)}p_{x}$
$\displaystyle\quad{}+2(1+3\nu)q_{y}\bigl{(}2\tau+3\mu
q_{x}-3\mu^{2}q_{y}^{2}\bigr{)}p_{y}+3\nu(1+3\nu)\mu\bigl{(}2q_{x}-3\mu
q_{y}^{2}\bigr{)}$
$\displaystyle{}\equiv\sum_{i,j=x,y}c_{ij}(q)p_{i}p_{j}+\sum_{i=x,y}c_{i}(q)p_{i}+c_{0}(q),$
(2.2)

where $\tau$, $\mu$, $\nu$ are parameters. Here the coefficients $c_{ij}$ are
the 4th degree polynomials in $q$-generators, $c_{i}$ are the 3rd degree ones
and $c_{0}$ is the 2nd degree polynomial. We also form another non-linear
combination in $p,q$-generators in the $U_{\mathfrak{h}_{5}}$,

$\displaystyle
k_{{\mathcal{A}}_{2}}(p_{x},q_{x},p_{y},q_{y})=2\nu(1+3\nu)(2+3\nu)\mu
q_{y}\bigl{(}2\tau+3\mu q_{x}-3\mu^{2}q_{y}^{2}\bigr{)}$
$\displaystyle\qquad{}+\frac{1}{3}(1+3\nu)(2+3\nu)q_{y}\bigl{(}\mu+8\tau^{2}+28\mu\tau
q_{x}+21\mu^{2}q_{x}^{2}-9\mu^{2}\tau
q_{y}^{2}-18\mu^{3}q_{x}q_{y}^{2}\bigr{)}p_{x}$
$\displaystyle\qquad{}-\frac{2}{9}(1+3\nu)(2+3\nu)\bigl{(}1+4\tau q_{x}+6\mu
q_{x}^{2}-24\mu\tau
q_{y}^{2}-36\mu^{2}q_{x}{q_{y}}^{2}+27\mu^{3}q_{y}^{4}\bigr{)}p_{y}$
$\displaystyle\qquad{}+(2+3\nu)q_{y}\bigl{(}3\tau+4\bigl{(}2\tau^{2}+\mu\bigr{)}q_{x}+17\mu\tau
q_{x}^{2}+8\mu^{2}q_{x}^{3}$
$\displaystyle\qquad\quad{}+3\mu\bigl{(}{\tau}^{2}-2\mu\bigr{)}q_{y}^{2}-6\mu^{2}\tau
q_{x}q_{y}^{2}-6{\mu}^{3}q_{x}^{2}q_{y}^{2}\bigr{)}p_{x}^{2}$
$\displaystyle\qquad{}-\frac{2}{3}(2+3\nu)\bigl{(}q_{x}+4\tau q_{x}^{2}+5\mu
q_{x}^{3}+3\bigl{(}\mu-4\tau^{2}\bigr{)}q_{y}^{2}-27\mu^{2}q_{x}^{2}q_{y}^{2}$
$\displaystyle\qquad\quad{}-33\mu\tau q_{x}q_{y}^{2}+9\mu^{2}\tau
q_{y}^{4}+18\mu^{3}q_{x}q_{y}^{4}\bigr{)}p_{x}p_{y}$
$\displaystyle\qquad{}-(2+3\nu)q_{y}\biggl{(}1+\frac{8}{3}\tau q_{x}+3\mu
q_{x}^{2}-7\mu\tau
q_{y}^{2}-10\mu^{2}q_{x}q_{y}^{2}+6\mu^{3}q_{y}^{4}\biggr{)}p_{y}^{2}$
$\displaystyle\qquad{}+q_{y}\bigl{(}1+5\tau
q_{x}+2\bigl{(}2\mu+3\tau^{2}\bigr{)}q_{x}^{2}+3\mu\bigl{(}\tau^{2}-2\mu\bigr{)}q_{x}q_{y}^{2}+9\mu\tau
q_{x}^{3}$
$\displaystyle\qquad\quad{}-\tau\bigl{(}3\mu-2\tau^{2}\bigr{)}q_{y}^{2}+3\mu^{2}q_{x}^{4}-3\mu^{2}\tau
q_{x}^{2}q_{y}^{2}-2\mu^{3}q_{x}^{3}q_{y}^{2}\bigr{)}p_{x}^{3}$
$\displaystyle\qquad{}+\biggl{(}-\frac{2}{3}q_{x}^{2}+2\bigl{(}5{\tau}^{2}+\mu\bigr{)}q_{x}q_{y}^{2}-2\tau
q_{x}^{3}+3\tau q_{y}^{2}-2\mu q_{x}^{4}+3\mu(\tau^{2}-2\mu)q_{y}^{4}$
$\displaystyle\qquad\quad{}+19\mu\tau
q_{x}^{2}q_{y}^{2}-6{\mu}^{3}q_{x}^{2}q_{y}^{4}+10\mu^{2}q_{x}^{3}q_{y}^{2}-6\mu^{2}\tau
q_{x}q_{y}^{4}\biggr{)}p_{x}^{2}p_{y}$
$\displaystyle\qquad{}-q_{y}\biggl{(}q_{x}+\frac{10}{3}\tau{q_{x}}^{2}+\frac{11}{3}\mu{q_{x}}^{3}-13\mu\tau
q_{x}q_{y}^{2}+3(\mu-2\tau^{2})q_{y}^{2}-11\mu^{2}q_{x}^{2}q_{y}^{2}$
$\displaystyle\qquad\quad{}+3\mu^{2}\tau
q_{y}^{4}+6\mu^{3}q_{x}q_{y}^{4}\biggr{)}p_{x}p_{y}^{2}$
$\displaystyle\qquad{}-\biggl{(}q_{y}^{2}+\frac{2}{27}q_{x}^{3}+2\tau
q_{x}q_{y}^{2}-3\mu\tau q_{y}^{4}+\frac{5}{3}\mu
q_{x}^{2}q_{y}^{2}-4{\mu}^{2}q_{x}q_{y}^{4}+2{\mu}^{3}q_{y}^{6}\biggr{)}p_{y}^{3}$
$\displaystyle\qquad{}\equiv\sum_{i,j,k=x,y}d_{ijk}(q)p_{i}p_{j}p_{k}+\sum_{i,j=x,y}d_{ij}(q)p_{i}p_{j}+\sum_{i=x,y}d_{i}(q)p_{i}+d_{0},$
(2.3)

where the coefficients $d_{ijk}$, $d_{ij}$, $d_{i}$, $d_{0}$ are polynomials
in $q$ of degrees 6, 5, 4 and 3, respectively.

######  Theorem 2.1.

The expressions (2.2) and (2.3) form the commutative pair,

$[h_{{\mathcal{A}}_{2}},k_{{\mathcal{A}}_{2}}]=0,$

for any values of parameters $\tau$, $\mu$, $\nu$.

###### Proof.

By direct calculation. ∎

Note, that it can be checked by direct calculation that
$h_{{\mathcal{A}}_{2}}$, $k_{{\mathcal{A}}_{2}}$ written in the (classical)
phase space $(p,q)$-variables, $\\{p,q\\}=1$, do not form the commutative pair
with respect to the Poisson bracket,
$\\{h_{{\mathcal{A}}_{2}},k_{{\mathcal{A}}_{2}}\\}\neq 0$, for any values of
the parameters $\tau$, $\mu$, $\nu$. The reason for it is conceptual. Since we
do not discuss in this paper the classical-quantum correspondence being fully
focused on the quantum case, this will be discussed elsewhere.

From one side, it can be easily checked by direct calculation when $(p,q)$
generators of $\mathfrak{h}_{5}$ are written in the coordinate-momentum
representation (A.4) the expressions (2.2) and (2.3) become (1.8) and (1.9),
respectively. From another side, the expressions (2.2) and (2.3) can be
written in terms of $\mathfrak{gl}(3)$ generators in $(-3\nu,0)$
representation (A.3) as bilinear and trilinear combinations with
$\nu$-dependent coefficients, respectively, cf. [6, equations (20) and (25)].
Hence, the expressions (2.2) and (2.3) define the integrable
$\mathfrak{gl}(3)$ Euler–Arnold quantum top, or, equivalently, the integrable
$\mathfrak{sl}(3)$ Euler–Arnold quantum top of spin $(-3\nu)$.

By introducing the vacuum $|0\rangle$ as an object annihilated by
$p$-operators:

$p_{x}\,|0\rangle=0,\qquad p_{y}\,|0\rangle=0,$

in addition to the universal enveloping algebra $U_{\mathfrak{h}_{5}}$, this
leads to definition of the Fock space. The formal eigenvalue problem in the
Fock space for the Hamiltonian $h_{{\mathcal{A}}_{2}}$ is as follows:

$h_{{\mathcal{A}}_{2}}(p_{x},q_{x},p_{y},q_{y})\phi^{(h)}(q_{x},q_{y})|0\rangle=\lambda^{(h)}\phi^{(h)}(q_{x},q_{y})|0\rangle,$
(2.4)

where $\phi^{(h)}(q)$ is the eigen-operator and $\lambda^{(h)}$ is the
eigenvalue (spectral parameter). Analogously,

$k_{{\mathcal{A}}_{2}}(p_{x},q_{x},p_{y},q_{y})\phi^{(k)}(q_{x},q_{y})|0\rangle=\lambda^{(k)}\phi^{(k)}(q_{x},q_{y})|0\rangle.$
(2.5)

Owing to Theorem 2.1 the eigenvalue problems (2.4) and (2.5) have common
eigen-operators $\phi$. If spin $-\nu=n/3$, $n=0,1,2,\dots$, which corresponds
to the $\mathfrak{gl}(3)$ finite-dimensional representation $(n,0)$, the
eigenvalue problems (2.4) and (2.5) have $\frac{(n+2)(n+1)}{2}$ polynomial
eigen-operators $\phi^{(h,k)}(q_{x},q_{y})$.

######  Examples 2.2.

  * •

For $n=0$ (thus, at zero coupling, $\kappa=0$),

$\lambda^{(h)}_{0,1}=0,\qquad\phi^{(h)}_{0,1}=1.$

  * •

For $n=1$ at coupling

$\kappa=\frac{4}{9},$

the operator $h_{A_{2}}$ has a three-dimensional kernel (three zero modes) of
the type

$(a_{1}q_{x}+a_{2}q_{y}+b)$

with coefficients $a_{1}$, $a_{2}$ which do not vanish simultaneously and,

$\lambda^{(h)}_{1,i}=0,\qquad i=1,2,3.$

  * •

The first non-zero eigenvalues appear for $n=2$, thus, at

$\kappa=\frac{10}{9}.$

In total, there exist six polynomial eigenstates. Eigenvalues are given by the
roots of the factorized algebraic equation of degree 6,

$\displaystyle\bigl{(}\bigl{(}\lambda^{(h)}\bigr{)}^{2}+4\tau\lambda^{(h)}+4\mu\bigr{)}\bigl{(}\bigl{(}\lambda^{(h)}\bigr{)}^{2}+8\tau\lambda^{(h)}+4\mu+12\tau^{2}\bigr{)}$
$\displaystyle\qquad{}\times\bigl{(}\bigl{(}\lambda^{(h)}\bigr{)}^{2}+12\tau\lambda^{(h)}+4\mu+16\tau^{2}\bigr{)}=0.$

Explicitly,

$\displaystyle\bigl{(}{\lambda^{(h)}}\bigr{)}_{\pm}^{(1)}=-2\bigl{(}\tau\pm\sqrt{\tau^{2}-\mu}\bigr{)},\qquad\bigl{(}{\lambda^{(h)}}\bigr{)}_{\pm}^{(2)}=-2\bigl{(}2\tau\pm\sqrt{\tau^{2}-\mu}\bigr{)},$
$\displaystyle\bigl{(}{\lambda^{(h)}}\bigr{)}_{\pm}^{(3)}=-2\bigl{(}3\tau\pm\sqrt{5\tau^{2}-\mu}\bigr{)},$

and the corresponding six eigen-operators are of the form

$\bigl{(}a_{1}q_{x}^{2}+a_{2}q_{x}q_{y}+a_{3}q_{y}^{2}b_{1}q_{x}+b_{2}q_{y}+c\bigr{)}$

with parameters $a_{1}$, $a_{2}$, $a_{3}$, which do not vanish simultaneously,
and $b_{1}$, $b_{2}$, $c$. In the limit $\tau=\mu=0$ (it corresponds to the
rational ${\mathcal{A}}_{2}$ integrable model without the harmonic oscillator
terms) all six eigenvalues are degenerate to zero.

##  3 $\boldsymbol{\mathfrak{gl}(3)}$-polynomial integrable model on a lattice

###  3.1 Uniform translation-invariant lattice

Let us introduce the shift operator,

$T_{\delta}f(x)=f(x+\delta),\qquad T_{\delta}={\rm e}^{\delta\partial_{x}},$

where $\delta\in{{\mathbb{R}}}$ is parameter, which, sometimes, is called
spacing, and construct a canonical pair of shift operators (see, e.g., [5])

$D_{\delta}=\frac{T_{\delta}-1}{\delta},\qquad
X_{\delta}=xT_{-\delta}=x(1-\delta{D}_{-\delta}),$ (3.1)

where the operator $D_{\delta}$ is defined as

$D_{\delta}f(x)=\frac{f(x+\delta)-f(x)}{\delta},$

sometimes, it is called the Norlund derivative. It is easy to check that
$[D_{\delta},X_{\delta}]=1$, hence, $D_{\delta}$, $X_{\delta}$ form the
canonical pair, both operators are non-local. In the limit $\delta\rightarrow
0$ this pair becomes the well-known coordinate-momentum representation
$(\partial_{x},x)$ of the Heisenberg algebra $\mathfrak{h}_{3}(p,q,I)$,

$[p,q]=1,\qquad[p,I]=[q,I]=0.$

For non-vanishing $\delta$, the canonical pair (3.1) belongs to the extended
universal enveloping algebra ${\hat{U}}_{\mathfrak{h}_{3}}$. These operators
act on infinite uniform lattice space with spacing $\delta$

$\\{\dots,x-2\delta,x-\delta,x,x+\delta,x+2\delta,\dots\\}$

marked by $x\in{{\mathbb{R}}}$ – a position of a central (or reference) point
of the lattice.

By taking $D_{\delta}$, $X_{\delta}$ (3.1) as basic elements, it can be shown
that algebra $\mathfrak{h}_{5}$ of finite-difference (shift) operators can be
formed:

$\displaystyle[D_{\delta_{1},x},X_{\delta_{1},x}]=1,$
$\displaystyle[D_{\delta_{2},y},X_{\delta_{2},y}]=1,$
$\displaystyle[D_{\delta_{1},x},X_{\delta_{2},y}]=0,$
$\displaystyle[D_{\delta_{2},y},X_{\delta_{1},x}]=0,$
$\displaystyle[D_{\delta_{1},x},D_{\delta_{2},y}]=0,\quad$
$\displaystyle[X_{\delta_{1},x},X_{\delta_{2},y}]=0,\quad$
$\displaystyle[D_{\delta_{1}(\delta_{2}),x(y)},I]=0,\quad$
$\displaystyle[X_{\delta_{1}(\delta_{2}),x(y)},I]=0.$ (3.2)

Evidently, the vacuum vector,

$|0\rangle=1,$

for any $(\delta_{1},\delta_{2})\in{{\mathbb{R}}}^{2}$.

This algebra acts on the rectangular uniform lattice with spacings
$(\delta_{1},\delta_{2})$. By identifying in (2.2) and (2.3) the variables
$(p,q)$ with $(D_{\delta},X_{\delta})$, we arrive at the Hamiltonian and the
integral of the polynomial integrable model on the two-dimensional uniform
lattice with spacings $(\delta_{1},\delta_{2})$,

$h^{(\delta_{1},\delta_{2})}=h_{A_{2}}(D_{\delta_{1},x},X_{\delta_{1},x},D_{\delta_{2},y},X_{\delta_{2},y}),$
(3.3)

and

$k^{(\delta_{1},\delta_{2})}=k_{A_{2}}(D_{\delta_{1},x},X_{\delta_{1},x},D_{\delta_{2},y},X_{\delta_{2},y}),$
(3.4)

If parameter $-\nu=n/3$, $n=0,1,2,\dots$, the eigenvalue problems for the
operators (3.3) and (3.4) have $\frac{(n+2)(n+1)}{2}$ common polynomial
eigenfunctions $\phi^{(h,k)}(x,y)$ in the form of triangular polynomials,

$\langle x^{m_{x}}y^{m_{y}}\mid 0\leq m_{x}+m_{y}\leq n\rangle.$

The first polynomial eigenfunctions for $n=0,1,2$ can be easily found by using
the results collected in Examples 2.2.

###  3.2 Exponential dilatation-invariant lattice

Let us introduce the dilation operator,

$T_{q}f(x)=f(qx),\qquad T_{q}=q^{A},\qquad A\equiv x\partial_{x},$

where $q\in{{\mathbb{C}}}$, and construct a canonical pair of dilatation
operators

$D_{q}=x^{-1}\frac{T_{q}-1}{q-1},\qquad X_{q}=\frac{A(q-1)}{T_{q}-1}x,$ (3.5)

see [1], where $[D_{q},X_{q}]=1$ for any $q$. It can be checked that their
product is $q$-independent,

$X_{q}D_{q}=x\partial_{x}=A\qquad\mbox{and}\qquad
D_{q}X_{q}=\partial_{x}x=A+1.$

The operator $D_{q}$ is called the Jackson symbol (or the Jackson derivative).
Both operators $X_{q}$, $D_{q}$ are pseudodifferential operators with action
on monomials as follows:

$D_{q}x^{n}=\\{n\\}_{q}x^{n-1},\qquad
X_{q}x^{n}=\frac{n+1}{\\{n+1\\}_{q}}x^{n+1},$

where $\\{n\\}_{q}=\frac{1-q^{n}}{1-q}$ is the so called $q$-number $n$.

By taking $D_{q}$, $X_{q}$ (3.5) as basic elements it can be shown that
algebra $\mathfrak{h}_{5}$ of discrete operators can be formed:

$\displaystyle[D_{q_{1},x},X_{q_{1},x}]=1,\quad$
$\displaystyle[D_{q_{2},y},X_{q_{2},y}]=1,\quad$
$\displaystyle[D_{q_{1},x},X_{q_{2},y}]=0,\quad$
$\displaystyle[D_{q_{2},y},X_{q_{1},x}]=0,$
$\displaystyle[D_{q_{1},x},D_{q_{2},y}]=0,\qquad$
$\displaystyle[X_{q_{1},x},X_{q_{2},y}]=0,\qquad$
$\displaystyle[D_{q_{1}(q_{2}),x(y)},I]=0,\qquad$
$\displaystyle[X_{q_{1}(q_{2}),x(y)},I]=0,$

cf. (3.2). Evidently, the vacuum vector,

$|0\rangle=1$

for any $(q_{1},q_{2})\in{{\mathbb{R}}}^{2}$.

This algebra acts on the exponential lattice with spacings $(q_{1},q_{2})$. By
identifying in (2.2) and (2.3) the variables $(p,q)$ with $(D_{q},X_{q})$ we
arrive at the Hamiltonian and the integral of the polynomial integrable model
on the two-dimensional exponential lattice with spacings $(q_{1},q_{2})$,

$h^{(q_{1},q_{2})}=h_{A_{2}}(D_{q_{1},x},X_{q_{1},x},D_{q_{2},y},X_{q_{2},y})$
(3.6)

and

$k^{(q_{1},q_{2})}=k_{A_{2}}(D_{q_{1},x},X_{q_{1},x},D_{q_{2},y},X_{q_{2},y}).$
(3.7)

If parameter $-\nu=n/3$, $n=0,1,2,\dots$, the eigenvalue problems for (3.6)
and (3.7) have $\frac{(n+2)(n+1)}{2}$ common polynomial eigenfunctions
$\phi^{(h,k)}(x,y)$ in the form of triangular polynomials,

$\langle x^{m_{x}}y^{m_{y}}\mid 0\leq m_{x}+m_{y}\leq n\rangle.$

The first polynomial eigenfunctions for $n=0,1,2$ can be easily found by using
the results collected in Examples 2.2.

###  3.3 Mixed translation-invariant and dilatation-invariant lattice

It is evident that one can construct the operators $h$, $k$ acting in
$x$-direction on the uniform lattice and in $y$-direction on the exponential
lattice and visa versa. Therefore, there are two ways to realize it by taking

$p_{x}=D_{\delta_{1},x},\qquad q_{x}=X_{\delta_{1},x},\qquad
p_{y}=D_{q_{1},y},\qquad q_{y}=X_{q_{1},y},$ (3.8)

or

$p_{x}=D_{q_{2},x},\qquad q_{x}=X_{q_{2},x},\qquad
p_{y}=D_{\delta_{2},y},\qquad q_{y}=X_{\delta_{2},y}.$ (3.9)

In both cases the vacuum vector remains the same,

$|0\rangle=1.$

In a straightforward way one can build the Hamiltonian and the integral

$h^{(\delta_{1},q_{1})}=h_{A_{2}}(D_{\delta_{1},x},X_{\delta_{1},x},D_{q_{1},y},X_{q_{1},y}),$
(3.10)

and

$k^{(\delta_{1},q_{1})}=k_{A_{2}}(D_{\delta_{1},x},X_{\delta_{1},x},D_{q_{1},y},X_{q_{1},y}),$
(3.11)

for (3.8) and

$h^{(q_{2},\delta_{2})}=h_{A_{2}}(D_{q_{2},x},X_{q_{2},x},D_{\delta_{2},y},X_{\delta_{2},y}),$
(3.12)

and

$k^{(q_{2},\delta_{2})}=k_{A_{2}}(D_{q_{2},x},X_{q_{2},x},D_{\delta_{2},y},X_{\delta_{2},y}),$
(3.13)

for (3.9). In similar way as for (1.8)–(1.9), (3.3)–(3.4), (3.6)–(3.7), if
parameter $\nu=n/3$, ${n=0,1,2,\dots}$, the eigenvalue problems for
(3.10)–(3.11) and (3.12)–(3.13) have $\frac{(n+2)(n+1)}{2}$ common polynomial
eigenfunctions $\phi^{(h,k)}(x,y)$ in the form of triangular polynomials,

$\langle x^{m_{x}}y^{m_{y}}\mid 0\leq m_{x}+m_{y}\leq n\rangle.$

The first polynomial eigenfunctions for $n=0,1,2$ can be easily found by using
the results collected in Examples 2.2.

Remarkably, all these five integrable models (1.8)–(1.9), (3.3)–(3.4),
(3.6)–(3.7) and (3.10)–(3.11), (3.12)–(3.13) are isospectral.

##  4 $\boldsymbol{\mathfrak{gl}(3)}$-polynomial integrable model in
$\boldsymbol{{\bf{\mathbb{C}}^{2}}}$

Introduce the five-dimensional Heisenberg algebra
$\mathfrak{h}_{5}\equiv\mathbb{H}_{5}=\big{\\{}\mathfrak{a}_{1},\mathfrak{a}^{\dagger}_{1},\mathfrak{a}_{2},\mathfrak{a}^{\dagger}_{2},I\big{\\}}$
with commutator
$\big{[}\mathfrak{a}_{i},\mathfrak{a}^{\dagger}_{j}\big{]}=\delta_{ij}I$,
$i,j=1,2$,
$[\mathfrak{a}_{i},\mathfrak{a}_{j}]=\big{[}\mathfrak{a}^{\dagger}_{i},\mathfrak{a}^{\dagger}_{j}\big{]}=0$
and $[\mathfrak{a}_{i},I]=\big{[}\mathfrak{a}^{\dagger}_{i},I\big{]}=0$ by
using a new, mathematics-oriented notations [13]. Its representation on the
standard Hilbert space,

$L_{2}\bigl{(}{\mathbb{C}}^{2},{\rm d}\mu_{2}\bigr{)}=L_{2}({\mathbb{C}},{\rm
d}\mu)\otimes L_{2}({\mathbb{C}},{\rm d}\mu),$

with the Gaussian measure,

${\rm d}\mu(z)=\pi^{-1}{\rm e}^{-z\cdot{\bar{z}}}{\rm d}v(z),$

where ${\rm d}v(z)={\rm d}x{\rm d}y$ is the Euclidean volume measure on
${\mathbb{C}}={\mathbb{R}}^{2}$, is given by two canonical pairs of raising
and lowering operators related to $z=(z_{1},z_{2})\in{\mathbb{C}}^{2}$:

$\mathfrak{a}^{\dagger}_{j}={\bar{z}}_{j}-\frac{\partial}{\partial{z}_{j}},\qquad\mathfrak{a}_{j}=\frac{\partial}{\partial{\bar{z}}_{j}},$

where $\mathfrak{a}^{\dagger}_{j}$ is adjoint to $\mathfrak{a}_{j}$ with
$\big{[}\mathfrak{a}_{j},\mathfrak{a}^{\dagger}_{j}\big{]}=I$, $j=1,2$, see
[13] for details and $I$ is the identity operator. The vacuum vector
$|0\rangle$, defined by

$\mathfrak{a}_{1}|0\rangle=0,\qquad\mathfrak{a}_{2}|0\rangle=0,$

is any two-dimensional analytic function.

Formally, by taking (2.2) and (2.3) one can build the Hamiltonian

$h^{({\mathbb{C}}^{2})}=h_{{\mathcal{A}}_{2}}\bigl{(}\mathfrak{a}_{1},\mathfrak{a}^{\dagger}_{1},\mathfrak{a}_{2},\mathfrak{a}^{\dagger}_{2}\bigr{)},$

and the integral

$k^{({\mathbb{C}}^{2})}=k_{{\mathcal{A}}_{2}}\bigl{(}\mathfrak{a}_{1},\mathfrak{a}^{\dagger}_{1},\mathfrak{a}_{2},\mathfrak{a}^{\dagger}_{2}\bigr{)}.$

It is evident that they continue to commute. This procedure can be considered
as a complexification of the original polynomial model (1.8) and (1.9), which
emerged from the 3-body/${\mathcal{A}}_{2}$ elliptic Calogero model as its
algebraic version. Formally, the Hamiltonian is the sixth-order differential
operator in $\frac{\partial}{\partial{z}}$,
$\frac{\partial}{\partial{\bar{z}}}$. Note that the first polynomial
eigenfunctions of $h^{({\mathbb{C}}^{2})}$ for $n=0,1,2$ can be easily found
by using the results collected in Examples 2.2.

##  5 $\boldsymbol{\mathfrak{gl}(3)}$ algebra: artifacts

Long ago one of the authors (Alexander V. Turbiner) discovered in the algebra
$\mathfrak{gl}(3)$ with generators defined in (A.1) the existence of nine
bilinear combinations in generators with unusual property: all those bilinear
combinations vanish if the representation of $\mathfrak{gl}(3)$ generators by
the first-order differential operators (A.2) is taken! The explicit form of
the bilinear combinations is the following [9]:

$\displaystyle A_{1}=J_{8}J_{5}-J_{7}J_{6},\qquad$ $\displaystyle
A_{2}=J_{8}J_{3}-J_{7}J_{4},\qquad$ $\displaystyle
A_{3}=J_{7}J_{2}+J_{5}J_{0}+J_{5},$ $\displaystyle
A_{4}=J_{8}J_{1}+J_{4}J_{0}+J_{4},\qquad$ $\displaystyle
A_{5}=J_{7}J_{1}+J_{3}J_{0}+J_{3},\qquad$ $\displaystyle
A_{6}=J_{8}J_{2}+J_{6}J_{0}+J_{6},$ $\displaystyle
A_{7}=J_{6}J_{3}-J_{5}J_{4}+J_{3},\qquad$ $\displaystyle
A_{8}=J_{6}J_{1}-J_{4}J_{2},\qquad$ $\displaystyle
A_{9}=J_{5}J_{1}-J_{3}J_{2}.$ (5.1)

######  Theorem 5.1.

For the $\mathfrak{gl}(3)$ generators, written in terms of the Heisenberg
algebra $\mathfrak{h}_{5}$ generators (A.3), all nine artifacts (5.1) vanish

$A_{1,\dots,9}(p_{x,y},q_{x,y})=0.$

This theorem can be proved by direct calculation.

It can be checked that the commutators of the artifacts are of the form

$[A_{i},A_{J}]=\sum c_{ij}^{k}(J)A_{k},$

where $c_{ij}^{k}(J)$ for $i\neq j$ are linear combinations of the
$\mathfrak{gl}(3)$ generators. For example,

$[A_{1},A_{2}]=-(J_{8}A_{1}+J_{7}A_{2}),$

see Appendix C for the explicit calculation. Hence, $A_{i}$ do not span a Lie
algebra. Interesting question to ask is what would happen with the artifacts
$A$’s if instead of $\mathfrak{gl}(3)$ generators the $\mathfrak{gl}(3)$
Kac–Moody currents are taken. It will be addressed elsewhere.

##  6 The Hamiltonian and the integral in
$\boldsymbol{\mathfrak{gl}(3)}$-algebra generators

###  6.1 Hamiltonian

By taking the Hamiltonian (2.2), one can demonstrate that it can be rewritten
in the $\mathfrak{gl}(3)$ abstract generators, which obey formally the
commutation relations (A.1),

$\displaystyle h_{{\mathcal{A}}_{2}}(J)={}$ $\displaystyle
2J_{6}J_{1}-\frac{1}{3}J_{5}^{2}-J_{1}J_{0}+\mu(2J_{8}J_{5}+J_{7}J_{3}-2J_{7}J_{0}+3J_{4}^{2}+2J_{7})$
$\displaystyle{}+\tau(4J_{8}J_{2}+4J_{7}J_{1}-J_{6}^{2}-J_{3}^{2}+5J_{6}+5J_{3})-3\tau\mu
J_{8}J_{4}-3\mu^{2}J_{8}^{2}-3\tau^{2}J_{4}^{2},$ (6.1)

hence, in extremely compact form; here $\mu$, $\tau$ are parameters and the
dependence on $\nu$ can be included into the representation (into the
generators) and eventually is absent! Hence, (6.1) is two-parametric, bilinear
element of the universal enveloping algebra $U_{\mathfrak{gl}(3)}$. If
$\mu=\tau=0$ the element $h_{{\mathcal{A}}_{2}}$ (6.1) dramatically
simplifies,

$\displaystyle
h^{(\mu=\tau=0)}_{{\mathcal{A}}_{2}}(J)=2J_{6}J_{1}-\tfrac{1}{3}J_{5}^{2}-J_{1}J_{0}.$

By substituting the generators $J_{0,1,5,6}$ in the form of differential
operators (A.2), one can see that this element corresponds to the
3-body/${\mathcal{A}}_{2}$ rational Calogero model (without harmonic
oscillator term). Non-surprisingly, the raising generators $J_{7,8}$ are
absent in this case, as well as the generators $J_{4,3,2}$.

###  6.2 Cubic integral

In a similar way, as was done in order to construct (6.1), by taking the
integral (2.3) in the Fock space representation one can demonstrate that it
can be rewritten in the $\mathfrak{gl}(3)$ abstract generators which obey the
commutation relations (A.1),

$\displaystyle k_{{\mathcal{A}}_{2}}(J)={}$
$\displaystyle-\frac{2}{9}J_{6}^{2}J_{2}+\frac{2}{9}J_{6}J_{5}J_{1}+\frac{5}{9}J_{6}J_{2}J_{0}-\frac{2}{27}J_{5}^{3}+\frac{2}{9}J_{5}J_{1}J_{0}+J_{4}J_{1}^{2}$
$\displaystyle{}-\frac{2}{9}J_{3}^{2}J_{2}-\frac{2}{9}J_{2}J_{0}^{2}+\frac{2}{9}J_{6}J_{2}+\frac{2}{9}J_{5}J_{1}+\frac{2}{9}J_{2}J_{0}$
$\displaystyle{}-\tau\biggl{(}\frac{8}{9}J_{7}J_{6}J_{2}+\frac{8}{9}J_{7}J_{5}J_{1}-\frac{8}{9}J_{7}J_{2}J_{0}+\frac{2}{9}J_{6}J_{6}J_{5}-\frac{2}{9}J_{6}J_{5}J_{3}+\frac{2}{9}J_{5}J_{3}J_{3}$
$\displaystyle\quad{}-2J_{4}J_{3}J_{1}-3J_{8}J_{1}J_{1}+\frac{2}{3}J_{6}J_{5}+\frac{2}{3}J_{5}J_{3}-\frac{16}{9}J_{5}J_{0}-4J_{4}J_{1}\biggr{)}$
$\displaystyle{}+\tau^{2}\biggl{(}\frac{2}{3}J_{6}^{2}J_{4}-\frac{2}{3}J_{6}J_{4}J_{3}-\frac{8}{3}J_{6}J_{4}J_{0}+\frac{2}{3}J_{4}J_{3}^{2}-\frac{8}{3}J_{4}J_{3}J_{0}+\frac{8}{3}J_{4}J_{0}^{2}$
$\displaystyle\quad{}-\frac{4}{3}J_{6}J_{4}-\frac{4}{3}J_{4}J_{3}+\frac{8}{3}J_{4}J_{0}+2J_{4}\biggr{)}+2\tau^{3}J_{4}^{3}$
$\displaystyle{}-\mu\biggl{(}\frac{1}{3}J_{7}J_{6}J_{5}+\frac{2}{3}J_{7}J_{5}J_{3}-\frac{4}{3}J_{7}J_{5}J_{0}+\frac{2}{3}J_{6}^{2}J_{4}-\frac{2}{3}J_{6}J_{4}J_{3}-\frac{8}{3}J_{6}J_{4}J_{0}$
$\displaystyle\quad{}-\frac{1}{3}J_{4}J_{3}^{2}+\frac{10}{3}J_{4}J_{3}J_{0}-\frac{1}{3}J_{4}J_{0}^{2}+\frac{4}{3}J_{7}J_{5}-\frac{4}{3}J_{6}J_{4}+\frac{5}{3}J_{4}J_{3}-\frac{1}{3}J_{4}J_{0}\biggr{)}$
$\displaystyle{}-\mu\tau\biggl{(}4J_{8}J_{0}-\frac{1}{3}J_{8}J_{6}^{2}+\frac{28}{3}J_{8}J_{6}J_{3}+\frac{4}{3}J_{8}J_{6}J_{0}-\frac{7}{3}J_{8}J_{3}^{2}$
$\displaystyle\quad{}+\frac{16}{3}J_{8}J_{3}J_{0}-\frac{4}{3}J_{8}J_{0}^{2}-10J_{7}J_{6}J_{4}+3J_{4}^{3}-J_{8}J_{6}+7J_{8}J_{3}-\frac{8}{3}J_{8}\biggr{)}$
$\displaystyle{}+3\mu\tau^{2}J_{8}J_{4}^{2}-3\mu^{2}\tau J_{8}^{2}J_{4}$
$\displaystyle{}+\mu^{2}\bigl{(}2J_{8}J_{7}J_{6}+J_{8}J_{7}J_{3}-2J_{8}J_{7}J_{0}-6J_{8}J_{4}^{2}+4J_{8}J_{7}\bigr{)}-2\mu^{3}J_{8}^{3},$
(6.2)

where $\mu$, $\tau$ are parameters, see (1.2), and the explicit dependence on
$\nu$ is absent! Hence, it is two-parametric, trilinear element of the
universal enveloping algebra $U_{\mathfrak{gl}(3)}$. If $\mu=\tau=0$, the
element $k_{{\mathcal{A}}_{2}}$ (6.2) dramatically simplifies,

$\displaystyle
k^{(\mu=\tau=0)}_{{\mathcal{A}}_{2}}(J)=-\frac{2}{9}J_{6}^{2}J_{2}+\frac{2}{9}J_{6}J_{5}J_{1}+\frac{5}{9}J_{6}J_{2}J_{0}-\frac{2}{27}J_{5}^{3}+\frac{2}{9}J_{5}J_{1}J_{0}+J_{4}J_{1}^{2}-\frac{2}{9}J_{3}^{2}J_{2}$
$\displaystyle\hphantom{k^{(\mu=\tau=0)}_{{\mathcal{A}}_{2}}(J)=}{}-\frac{2}{9}J_{2}J_{0}^{2}+\frac{2}{9}J_{6}J_{2}+\frac{2}{9}J_{5}J_{1}+\frac{2}{9}J_{2}J_{0},$

it corresponds to the 3-body/${\mathcal{A}}_{2}$ rational integrable Calogero
model (without harmonic oscillator term). Since this is the exactly-solvable
problem, non-surprisingly, the raising generators $J_{7,8}$ are absent.

###  6.3 Commutator

By taking (6.1) and (6.2), one can make the extremely cumbersome (and very
lengthy) calculation of their Lie bracket (commutator) by using the specially
designed Maple 18 code. An example of the code is presented in Appendix C. It
was the main goal of the master thesis of one of the authors (Miguel A.
Guadarrama-Ayala). Eventually, it leads to the following statement:

######  Theorem 6.1.

The commutator of the expressions (6.1) and (6.2) is the linear superposition
of artifacts (5.1),

$[h_{{\mathcal{A}}_{2}}(J),k_{{\mathcal{A}}_{2}}(J)]=\sum_{i=1}^{9}c_{i}(J)A_{i},$
(6.3)

for any values of parameters $\tau$, $\mu$, where $c_{i}(J)$ are some
coefficient functions in $J$’s.

###### Proof.

By direct calculation by using the specially designed Maple 18 code. It was
carried out on a regular DELL desktop computer with 2.4 GhZ working frequency
and 6 GB RAM memory.

Intuitively, this result (6.3) is evident: in the Fock space representation,
where $h,k\in U_{\mathfrak{h}_{5}}$, the commutator should vanish, see
Theorems 2.1 and 5.1. Hence, the commutator should be a (non)-linear
combination of the artifacts. The fact that is a linear combination of the
artifacts is non-trivial.

Alternative way to represent the commutator (6.3) is as follows

$\displaystyle[h_{{\mathcal{A}}_{2}}(J),k_{{\mathcal{A}}_{2}}(J)]={}$
$\displaystyle
D_{1}+D_{2}\tau+D_{3}\mu+D_{4}\tau^{2}+D_{5}\tau\mu+D_{6}\mu^{2}+D_{7}\tau^{2}\mu+D_{8}\tau\mu^{2}$
$\displaystyle{}+D_{9}\mu^{3}+D_{10}\tau^{3}\mu+D_{11}\tau^{2}\mu^{2}+D_{12}\tau\mu^{3},$

where for the coefficients $D(J,A)$ are presented in Appendix B. ∎

##  7 $\boldsymbol{G_{2}}$ elliptic 3-body problem

By adding the 3-body interaction potential to the 3-body elliptic Calogero
Hamiltonian (1.4), we arrive at the 3-body Wolfes elliptic Hamiltonian in
$(y_{1},y_{2})$-coordinates (1.3),

$\displaystyle{\mathcal{H}}_{{G}_{2}}=-\frac{1}{3}\biggl{(}\frac{\partial^{2}}{\partial
y_{1}^{2}}+\frac{\partial^{2}}{\partial
y_{2}^{2}}-\frac{\partial^{2}}{\partial y_{1}\partial y_{2}}\biggr{)}$
$\displaystyle\hphantom{{\mathcal{H}}_{{G}_{2}}=}{}+(\nu-\lambda)(\nu-\lambda-1)(\wp(y_{1}-y_{2})+\wp(2y_{1}+y_{2})+\wp(y_{1}+2y_{2}))$
$\displaystyle\hphantom{{\mathcal{H}}_{{G}_{2}}=}{}+\lambda(3\lambda-1)(\wp(y_{1})+\wp(y_{2})+\wp(y_{1}+y_{2})),$
(7.1)

which is also called the $G_{2}$ elliptic Hamiltonian in the Hamiltonian
reduction nomenclature [3]. It is characterized by two coupling constants
which can be parameterized conveniently as
$\kappa\equiv(\nu-\lambda)(\nu-\lambda-1)$ and
$\kappa_{2}\equiv\lambda(3\lambda-1)$. If $\kappa_{2}=0$ (or $\lambda=0,1/3$),
we return at the ${\mathcal{A}}_{2}$ elliptic model. It was shown in [6] that
by making the gauge rotation and changing variables to
$\bigl{(}u=x,v=y^{2}\bigr{)}$, see (1.6), the Hamiltonian (7.1) appears in the
form of the algebraic operator $h_{G_{2}}$ – the second-order differential
operator with polynomial coefficients,

$\displaystyle h_{G_{2}}(u,v)=\bigl{(}u+3\tau u^{2}+3\mu
u^{3}+3\bigl{(}\mu-\tau^{2}\bigr{)}v-3\mu\tau
uv-3\mu^{2}u^{2}v\bigr{)}\frac{\partial^{2}}{\partial u^{2}}$
$\displaystyle\hphantom{h_{G_{2}}(u,v)=}{}+2v\bigl{(}3+8\tau u+7\mu
u^{2}-3\mu\tau v-6\mu^{2}uv\bigr{)}\frac{\partial^{2}}{\partial u\partial v}$
$\displaystyle\hphantom{h_{G_{2}}(u,v)=}{}+4v\biggl{(}-\frac{u^{2}}{3}+3\tau
v+4\mu uv-3\mu^{2}v^{2}\biggr{)}\frac{\partial^{2}}{\partial v^{2}}$
$\displaystyle\hphantom{h_{G_{2}}(u,v)=}{}+(1+3\nu)\bigl{(}1+4\tau u+5\mu
u^{2}-3\mu\tau v-6\mu^{2}uv\bigr{)}\frac{\partial}{\partial u}$
$\displaystyle\hphantom{h_{G_{2}}(u,v)=}{}+2\biggl{(}-\frac{u^{2}}{3}+\tau(7+12\nu)v+2\mu(5+9\nu)uv-9\mu^{2}(1+2\nu)v^{2}\biggr{)}\frac{\partial}{\partial
v}$ $\displaystyle\hphantom{h_{G_{2}}(u,v)=}{}+3\nu(1+3\nu)\mu(2u-3\mu v)$
$\displaystyle\hphantom{h_{G_{2}}(u,v)=}{}+\lambda\biggl{(}6\bigl{(}1+2\tau
u+\mu u^{2}\bigr{)}\frac{\partial}{\partial u}+4\bigl{(}-u^{2}+3\tau v+3\mu
uv\bigr{)}\frac{\partial}{\partial v}+18\nu\mu u\biggr{)}.$ (7.2)

If $-3\nu=n$, where $n=0,1,2,3,\dots$, this operator has a finite-dimensional
triangular invariant subspace,

${\mathcal{P}}_{n}=\langle u^{n_{u}}v^{n_{v}}\mid 0\leq n_{u}+2n_{v}\leq
n\rangle.$

This space coincides with finite-dimensional representation space of the
algebra $\mathfrak{g}^{(2)}$: infinite-dimensional, eleven-generated algebra
of differential operators, see [6] and references therein, where it acts
irreducibly. It implies that $h_{G_{2}}$ can be rewritten in terms
$\mathfrak{g}^{(2)}$ generators (the Burnside theorem) [6].

After extremely tedious (and very lengthy) symbolic calculations by using the
Maple 18 code, see Appendix C for an example, it can be shown that the
existence of a differential operator $k_{m}(u,v)$ of degree five such that the
operator

$k_{\rm G_{2}}=k_{{\mathcal{A}}_{2}}^{2}(u,v)+\lambda k_{m}(u,v;\lambda),$

commutes with the $G_{2}$ elliptic Hamiltonian $h_{G_{2}}$ (7.2); $k_{m}$ has
the form of polynomial in $\lambda$ of finite degree. Note that in the
particular case of the $G_{2}$ rational Hamiltonian (see (7.2) at
$\mu=\tau=0$), this operator was calculated in [7] (where it corresponded to
the case $k=3$) in slightly different variables other than $u$, $v$: it is a
polynomial in $\lambda$ of degree four. In general, this operator will be
presented in its explicit form elsewhere. So far this operator is unknown in
the explicit form.

By taking the 5-dimensional Heisenberg algebra $\mathfrak{h}_{5}$ spanned by
the generators $p_{u}$, $p_{v}$, $q_{u}$, $q_{v}$, $I$, see (2.1), one can
form the following second degree polynomial in $p_{u}$, $p_{v}$:

$\displaystyle h_{G_{2}}(p_{u},p_{v},q_{u},q_{v})=\bigl{(}q_{u}+3\tau
q_{u}^{2}+3\mu q_{u}^{3}+3\bigl{(}\mu-\tau^{2}\bigr{)}q_{v}-3\mu\tau
q_{u}q_{v}-3\mu^{2}q_{u}^{2}q_{v}\bigr{)}p_{u}^{2}$
$\displaystyle\hphantom{h_{G_{2}}(p_{u},p_{v},q_{u},q_{v})=}{}+2q_{v}\bigl{(}3+8\tau
q_{u}+7\mu q_{u}^{2}-3\mu\tau q_{v}-6\mu^{2}q_{u}q_{v}\bigr{)}p_{u}p_{v}$
$\displaystyle\hphantom{h_{G_{2}}(p_{u},p_{v},q_{u},q_{v})=}{}+4q_{v}\biggl{(}-\frac{q_{u}^{2}}{3}+3\tau
q_{v}+4\mu q_{u}q_{v}-3\mu^{2}q_{v}^{2}\biggr{)}p_{v}^{2}$
$\displaystyle\hphantom{h_{G_{2}}(p_{u},p_{v},q_{u},q_{v})=}{}+(1+3\nu)\bigl{(}1+4\tau
q_{u}+5\mu q_{u}^{2}-3\mu\tau q_{v}-6\mu^{2}q_{u}q_{v}\bigr{)}p_{u}$
$\displaystyle\hphantom{h_{G_{2}}(p_{u},p_{v},q_{u},q_{v})=}{}+2\biggl{(}-\frac{q_{u}^{2}}{3}+\tau(7+12\nu)q_{v}+2\mu(5+9\nu)q_{u}q_{v}-9\mu^{2}(1+2\nu)q_{v}^{2}\biggr{)}p_{v}$
$\displaystyle\hphantom{h_{G_{2}}(p_{u},p_{v},q_{u},q_{v})=}{}+3\nu(1+3\nu)\mu\bigl{(}2q_{u}-3\mu
q_{v}\bigr{)}$ (7.3)
$\displaystyle\hphantom{h_{G_{2}}(p_{u},p_{v},q_{u},q_{v})=}{}+\lambda\bigl{(}6\bigl{(}1+2\tau
q_{u}+\mu q_{u}^{2}\bigr{)}p_{u}+4\bigl{(}-q_{u}^{2}+3\tau q_{v}+3\mu
q_{u}q_{v}\bigr{)}p_{v}+18\nu\mu q_{u}\bigr{)}.$

It is easy to check that if $(p,q)$-variables are taken in the coordinate-
momentum representation,

$p_{u}=\frac{\partial}{\partial u},\qquad p_{v}=\frac{\partial}{\partial
v},\qquad q_{u}=u,\qquad q_{v}=v,$

cf. (A.4), the expression (7.3) is reduced to the operator (7.2). The operator
$h_{G_{2}}(p_{u},p_{v},q_{u},q_{v})$ represents the $G_{2}$ elliptic model in
the Fock space.

By substituting into (7.3) the representations (A.5)–(A.6), we will arrive at
the $G_{2}$ elliptic lattice Hamiltonians defined on uniform-uniform, uniform-
exponential, exponential-uniform, exponential-exponential lattices in $(u,v)$
space as well as the complexified $G_{2}$ elliptic Hamiltonian in the
algebraic form.

##  8 Conclusions

In this paper, a polynomial integrable system, associated with the algebra
$U_{\mathfrak{h}_{5}}$ and inspired by the algebraic representation of the
${\mathcal{A}}_{2}$ elliptic model in Fock space is defined. It has the form
of a second degree polynomial in $p_{i}$, $i=1,2$,

$h_{{\mathcal{A}}_{2}}=c^{(2)}_{ij}p_{i}p_{j}+c^{(1)}_{i}p_{i}+c^{(0)},$ (8.1)

for the Hamiltonian and a 3rd degree polynomial in $p_{i}$, $i=1,2$,

$k_{{\mathcal{A}}_{2}}=d^{(3)}_{ijk}p_{i}p_{j}p_{k}+d^{(2)}_{ij}p_{i}p_{j}+d^{(1)}_{i}p_{i}+d^{(0)},$
(8.2)

for the integral, where the coefficients $\\{c\\}$ and $\\{d\\}$ are
polynomials in $q$ of a finite degrees, while $(p_{i},q_{i})$ form a canonical
pair. Overall, the operators $h_{{\mathcal{A}}_{2}}$ and
$k_{{\mathcal{A}}_{2}}$ depend on three free parameters $\mu$, $\tau$, $\nu$.
Remarkably, both operators $h_{{\mathcal{A}}_{2}}$ and $k_{{\mathcal{A}}_{2}}$
can be rewritten in terms of the $\mathfrak{sl}(3)$ generators
$J_{1,2,\dots,8}$ and they can be embedded into the $U_{\mathfrak{h}_{5}}$
algebra in the $(-3\nu,0)$ representation (A.3). Hence, $\nu$ corresponds to
the mark of the representation.

It can be conjectured that

######  Conjecture 8.1.

Up to canonical transformation

$p\rightarrow p+f(q),\qquad q\rightarrow q,$

there are no other non-trivial commuting operators in the
$U_{\mathfrak{h}_{5}}$ algebra of degree 2 and 3 in $p$ other than $h$ (8.1)
and $k$ (8.2).

The operators $h$ and $k$ can be rewritten in terms of the abstract
$\mathfrak{gl}(3)$ generators which obey the commutation relations (A.1) and
which give a non-vanishing commutator $[h,k]$. However, once the
$\mathfrak{gl}(3)$ generators are taken in the concrete representation (A.3)
the operators $h$ and $k$ becomes $h_{{\mathcal{A}}_{2}}$ (2.2) and
$k_{{\mathcal{A}}_{2}}$ (2.3), respectively, and their commutator
$[h_{{\mathcal{A}}_{2}},k_{{\mathcal{A}}_{2}}]=0$. The remarkable property of
the commutator $[h,k]$ is that it can be written as a linear superposition of
the artifacts $A_{1,2,\dots,9}$. We doubt there exist other elements of the
universal enveloping algebra $U_{\mathfrak{gl}(3)}$ (up to automorphisms) with
such a property.

Different realizations of $(p_{i},q_{i})$, $i=1,2$ by differential operators,
finite-difference operators, discrete operators, or the operators in $z$,
${\bar{z}}$ variables lead to a variety of concrete isospectral quantum
integrable polynomial systems in two continuous variables, on two-dimensional
uniform, exponential lattices or mixed ones, and on the ${{\mathbb{C}}}^{2}$
complex space. All these integrable models depend on the continuous parameter
$\nu$. If this parameter takes certain discrete values, all above-mentioned
integrable systems become quasi-exactly-solvable ones admitting a finite
number of polynomial eigenfunctions in the form of triangular polynomials.

##  Appendix A $\boldsymbol{\mathfrak{gl}(3)}$ algebra

The algebra $\mathfrak{gl}(3)$ is defined by nine generators $J_{i}$,
$i=0,1,2,\dots,8$, which obey the following commutation relations:

$\displaystyle[J_{0},J_{1}]=J_{1},\qquad$
$\displaystyle[J_{0},J_{2}]=J_{2},\qquad$
$\displaystyle[J_{0},J_{3}]=0,\qquad$ $\displaystyle[J_{0},J_{4}]=0,$
$\displaystyle[J_{0},J_{5}]=0,\qquad$ $\displaystyle[J_{0},J_{6}]=0,\qquad$
$\displaystyle[J_{0},J_{7}]=-J_{7},\qquad$
$\displaystyle[J_{0},J_{8}]=-J_{8},$ $\displaystyle[J_{1},J_{2}]=0,\qquad$
$\displaystyle[J_{1},J_{3}]=J_{1},\qquad$
$\displaystyle[J_{1},J_{4}]=0,\qquad$ $\displaystyle[J_{1},J_{5}]=J_{2},$
$\displaystyle[J_{1},J_{6}]=0,\qquad$
$\displaystyle[J_{1},J_{7}]=J_{3}-J_{0},\qquad$
$\displaystyle[J_{1},J_{8}]=J_{4},$ $\displaystyle[J_{2},J_{3}]=0,\qquad$
$\displaystyle[J_{2},J_{4}]=J_{1},\qquad$
$\displaystyle[J_{2},J_{5}]=0,\qquad$ $\displaystyle[J_{2},J_{6}]=J_{2},$
$\displaystyle[J_{2},J_{7}]=J_{5},\qquad$
$\displaystyle[J_{2},J_{8}]=J_{6}-J_{0},$
$\displaystyle[J_{3},J_{4}]=-J_{4},\qquad$
$\displaystyle[J_{3},J_{5}]=J_{5},\qquad$
$\displaystyle[J_{3},J_{6}]=0,\qquad$ $\displaystyle[J_{3},J_{7}]=J_{7},$
$\displaystyle[J_{3},J_{8}]=0,$
$\displaystyle[J_{4},J_{5}]=-J_{3}+J_{6},\qquad$
$\displaystyle[J_{4},J_{6}]=-J_{4},\qquad$
$\displaystyle[J_{4},J_{7}]=J_{8},\qquad$ $\displaystyle[J_{4},J_{8}]=0,$
$\displaystyle[J_{5},J_{6}]=J_{5},\qquad$
$\displaystyle[J_{5},J_{7}]=0,\qquad$ $\displaystyle[J_{5},J_{8}]=J_{7},$
$\displaystyle[J_{6},J_{7}]=0,\qquad$
$\displaystyle[J_{6},J_{8}]=J_{8},\qquad$ $\displaystyle[J_{7},J_{8}]=0.$
(A.1)

###  A.1 Structure constants

The commutation relations (A.1) of the $\mathfrak{gl}(3)$ algebra can be
represented as

$[J_{i},J_{j}]=c_{ij}^{k}J_{k},\qquad i,j,k=0,\dots,8,$

where $c_{ij}^{k}$ are the structure constants. The non-vanishing structure
constants are:

$\displaystyle c_{01}^{1}=1,\qquad$ $\displaystyle c_{02}^{2}=1,\qquad$
$\displaystyle c_{07}^{7}=-1,\qquad$ $\displaystyle c_{08}^{8}=-1,$
$\displaystyle c_{13}^{1}=1,\qquad$ $\displaystyle c_{15}^{2}=1,\qquad$
$\displaystyle c_{17}^{3}=1,\qquad$ $\displaystyle c_{17}^{0}=-1,\qquad$
$\displaystyle c_{18}^{4}=1,$ $\displaystyle c_{24}^{1}=1,\qquad$
$\displaystyle c_{26}^{2}=1,\qquad$ $\displaystyle c_{27}^{5}=1,\qquad$
$\displaystyle c_{28}^{6}=1,\qquad$ $\displaystyle c_{28}^{0}=-1,$
$\displaystyle c_{34}^{4}=-1,\qquad$ $\displaystyle c_{35}^{5}=1,\qquad$
$\displaystyle c_{37}^{7}=1,$ $\displaystyle c_{45}^{3}=-1,\qquad$
$\displaystyle c_{45}^{6}=1,\qquad$ $\displaystyle c_{46}^{4}=-1,\qquad$
$\displaystyle c_{47}^{8}=1,$ $\displaystyle c_{56}^{5}=1,\qquad$
$\displaystyle c_{58}^{7}=1,\qquad$ $\displaystyle c_{68}^{8}=1.$

###  A.2 Representation of $\boldsymbol{\mathfrak{gl}(3)}$ algebra in
differential operators

The algebra $\mathfrak{gl}(3)$ with commutation relations (A.1) can be
realized by the first-order differential operators in two variables,

$\displaystyle J_{1}=\frac{\partial}{\partial x},\qquad
J_{2}=\frac{\partial}{\partial y},\qquad J_{3}=x\frac{\partial}{\partial
x},\qquad J_{4}=y\frac{\partial}{\partial x},\qquad
J_{5}=x\frac{\partial}{\partial y},\qquad J_{6}=y\frac{\partial}{\partial y},$
$\displaystyle J_{7}=x\biggl{(}x\frac{\partial}{\partial
x}+y\frac{\partial}{\partial y}+3\nu\biggr{)},\qquad
J_{8}=y\biggl{(}x\frac{\partial}{\partial x}+y\frac{\partial}{\partial
y}+3\nu\biggr{)},$ (A.2)

and

$-J_{0}=x\frac{\partial}{\partial x}+y\frac{\partial}{\partial
y}+3\nu=J_{3}+J_{6}+3\nu,$

where $\nu$ is parameter. It corresponds to the irreducible representation of
the spin $(-3\nu,0)$. If $-3\nu=n$ is integer, the finite-dimensional
representation space which is spanned by triangular polynomials,

${\mathcal{P}}_{n}=\langle x^{p}y^{q}\mid 0\leq(p+q)\leq n\rangle,$

occurs.

###  A.3 Representation of $\boldsymbol{\mathfrak{gl}(3)}$ in
$\boldsymbol{(p,q)}$ space

Let us take 5-dimensional Heisenberg algebra $\mathfrak{h}_{5}$ spanned by the
generators $p_{x}$, $p_{y}$, $q_{x}$, $q_{x}$, $I$, which satisfy the
commutation relations,

$\displaystyle[p_{x},q_{x}]=1,\qquad$ $\displaystyle[p_{y},q_{y}]=1,\qquad$
$\displaystyle[p_{x},q_{y}]=0,\qquad$ $\displaystyle[p_{y},q_{x}]=0,$
$\displaystyle[p_{x},p_{y}]=0,\qquad$ $\displaystyle[q_{x},q_{y}]=0,\qquad$
$\displaystyle[p_{x,y},I]=0,\qquad$ $\displaystyle[q_{x,y},I]=0.$

Define its universal enveloping algebra $U_{\mathfrak{h}_{5}}$ as the algebra
of all ordered monomials
$\big{\\{}q_{x}^{i_{x}}q_{y}^{i_{y}}p_{x}^{j_{x}}p_{y}^{j_{y}}\big{\\}}$. It
is evident that the algebra $\mathfrak{gl}(3)$ realized as

$\displaystyle J_{1}=p_{x},\qquad J_{2}=p_{y},\qquad J_{3}=q_{x}p_{x},\qquad
J_{4}=q_{y}p_{x},\qquad J_{5}=q_{x}p_{y},\qquad J_{6}=q_{y}p_{y},$
$\displaystyle J_{7}=q_{x}(q_{x}p_{x}+q_{y}p_{y}+3\nu),\qquad
J_{8}=q_{y}(q_{x}p_{x}+q_{y}p_{y}+3\nu),$ (A.3)

and

$-J_{0}=q_{x}p_{x}+q_{y}p_{y}+3\nu=J_{3}+J_{6}+3\nu,$

is embedded into the universal enveloping algebra $U_{\mathfrak{h}_{5}}$.

Let us enlist four realizations of the commutation relation $[p_{x},q_{x}]=1$:

  * •

Continuous

$p_{x}=\frac{\partial}{\partial x}\equiv\partial_{x},\qquad q_{x}=x.$ (A.4)

It is well known, the so-called coordinate-momentum representation of the
$\mathfrak{h}_{3}$ Heisenberg algebra.

  * •

On uniform lattice

$p_{x}={\mathcal{D}}_{\delta},\qquad q_{x}=X_{\delta},$ (A.5)

with gap $\delta$, where ${\mathcal{D}}_{\delta}$ is the Norlund derivative
[5], it is the basis for the so-called umbral calculus.

  * •

On exponential lattice

$p_{x}={\mathcal{D}}_{q},\qquad q_{x}=X_{q},$

where ${\mathcal{D}}_{q}$ is the Jackson derivative, $q$ has the meaning of
the exponential spacing. It is described in details in [1].

  * •

Complex representation on ${\mathbb{C}}$

$\mathfrak{a}=\frac{\partial}{\partial{\bar{z}}},\qquad\mathfrak{a}^{\dagger}=-\frac{\partial}{\partial
z}+{\bar{z}},$ (A.6)

see [13] and references therein.

##  Appendix B Coefficients in the commutator (6.3)

The commutator between $h_{A_{2}}$ and $k_{A_{2}}$ can be written as the
polynomial in parameters $\tau$, $\mu$,

$\displaystyle[h_{A_{2}}(J),k_{A_{2}}(J)]={}$ $\displaystyle
D_{1}+D_{2}\tau+D_{3}\mu+D_{4}\tau^{2}+D_{5}\tau\mu+D_{6}\mu^{2}+D_{7}\tau^{2}\mu$
$\displaystyle{}+D_{8}\tau\mu^{2}+D_{9}\mu^{3}+D_{10}\tau^{3}\mu+D_{11}\tau^{2}\mu^{2}+D_{12}\tau\mu^{3},$

where the coefficients $D_{1,\dots,12}$ are presented by superposition of the
ordered polynomials in $\mathfrak{gl}(3)$-generators $J_{0,1,\dots,8}$
multiplied by the artifacts $A_{1,\dots,9}$ of the $\mathfrak{gl}(3)$ algebra,

$\displaystyle
D_{1}=-\frac{2}{9}(8J_{4}J_{2}+3J_{3}J_{1})A_{9}-\frac{2}{9}(8J_{5}J_{1}-8J_{3}J_{2}-11J_{2}J_{0})A_{8}$
$\displaystyle\hphantom{D_{1}=}{}-\frac{4}{3}J_{2}J_{1}A_{7}-\frac{22}{9}J_{2}J_{1}A_{6}+\frac{4}{9}J_{2}J_{1}A_{5}+\frac{22}{9}J_{2}^{2}A_{4}-\frac{4}{9}J_{1}^{2}A_{3},$
$\displaystyle
D_{2}=\frac{2}{9}\bigl{(}-6J_{6}^{2}-6J_{5}J_{4}+3J_{3}J_{0}+4J_{0}^{2}-8J_{6}+3J_{3}+10J_{0}-14\bigr{)}A_{9}$
$\displaystyle\hphantom{D_{2}=}{}+\frac{8}{9}(3J_{6}J_{5}+9J_{4}J_{1}+4J_{5})A_{8}-\frac{2}{9}(12J_{5}J_{1}-13J_{2}J_{0})A_{7}-\frac{28}{9}J_{6}J_{2}A_{5}+\frac{28}{9}J_{6}J_{1}A_{3},$
$\displaystyle
D_{3}=\frac{2}{9}\bigl{(}2J_{8}J_{5}-4J_{7}J_{3}+3J_{7}J_{0}-36J_{4}^{2}+4J_{7}\bigr{)}A_{9}$
$\displaystyle\hphantom{D_{3}=}{}+\frac{1}{3}(2J_{8}J_{1}-7J_{7}J_{5}+24J_{4}J_{3}+30J_{4})A_{8}$
$\displaystyle\hphantom{D_{3}=}{}+\frac{1}{9}(5J_{7}J_{2}+12J_{6}J_{5}-12J_{5}J_{3}+36J_{5}J_{0}-10J_{5})A_{7}$
$\displaystyle\hphantom{D_{3}=}{}-\frac{4}{9}(3J_{5}J_{0}-4J_{5})A_{6}-\frac{1}{9}(36J_{6}J_{5}-16J_{5}J_{3}+12J_{5}J_{0}+63J_{4}J_{1})A_{5}$
$\displaystyle\hphantom{D_{3}=}{}+\frac{1}{3}(-8J_{6}J_{1}-10J_{4}J_{2}+3J_{3}J_{1}+6J_{1}J_{0}+17J_{1})A_{4}$
$\displaystyle\hphantom{D_{3}=}{}+\frac{1}{9}\bigl{(}4J_{6}^{2}-J_{6}J_{0}-4J_{5}J_{4}-19J_{6}+8J_{0}-12\bigr{)}A_{3}+\frac{4}{3}J_{5}J_{2}A_{2}+\frac{2}{3}J_{6}J_{2}A_{1},$
$\displaystyle
D_{4}=\frac{8}{3}(3J_{4}J_{3}-2J_{4}J_{0})A_{8}-4J_{4}J_{1}A_{7}-10J_{4}J_{1}A_{6}+10J_{4}J_{2}A_{4},$
$\displaystyle
D_{5}=\frac{1}{3}(9J_{8}J_{6}+48J_{8}J_{3}+14J_{7}J_{4}+71J_{8})A_{8}-\frac{2}{3}(2J_{7}J_{5}-3J_{4}J_{3}+20J_{4}J_{0})A_{7}$
$\displaystyle\hphantom{D_{5}=}{}-\frac{2}{3}(16J_{8}J_{1}-23J_{4}J_{3})A_{6}+\frac{1}{6}(83J_{8}J_{1}-78J_{4}J_{3}+219J_{4}J_{0}+242J_{4})A_{5}$
$\displaystyle\hphantom{D_{5}=}{}+\frac{1}{6}\bigl{(}64J_{8}J_{2}-83J_{7}J_{1}-124J_{6}^{2}+34J_{6}J_{0}-40J_{5}J_{4}+50J_{3}^{2}$
$\displaystyle\hphantom{D_{5}=+\frac{1}{6}\bigl{(}}{}-229J_{3}J_{0}+54J_{6}+32J_{0}^{2}-297J_{0}-66\bigr{)}A_{4}$
$\displaystyle\hphantom{D_{5}=}{}-\frac{2}{3}\bigl{(}41J_{6}J_{1}-13J_{5}^{2}+7J_{4}J_{2}\bigr{)}A_{2}-\frac{2}{3}(9J_{6}J_{5}+4J_{5}J_{0}-J_{5})A_{1},$
$\displaystyle
D_{6}=\frac{26}{3}J_{8}^{2}A_{9}+\frac{2}{3}(3J_{8}J_{6}+3J_{8}J_{3}-26J_{8}J_{0}-J_{8})A_{7}$
$\displaystyle\hphantom{D_{6}=}{}-6(J_{8}J_{6}-J_{8}J_{3}-J_{8}J_{0})A_{6}-\frac{1}{3}(7J_{8}J_{3}+10J_{8}J_{0}+20J_{8}-19J_{7}J_{4})A_{5}$
$\displaystyle\hphantom{D_{6}=}{}+\frac{1}{3}\bigl{(}36J_{7}J_{6}-19J_{7}J_{0}-90J_{4}^{2}+21J_{7}\bigr{)}A_{4}$
$\displaystyle\hphantom{D_{6}=}{}+\frac{1}{3}\bigl{(}19J_{7}J_{1}-8J_{6}^{2}-4J_{6}J_{3}+50J_{6}J_{0}-6J_{5}J_{4}+J_{3}^{2}+54J_{6}+20J_{3}+50\bigr{)}A_{2}$
$\displaystyle\hphantom{D_{6}=}{}-(3J_{8}J_{1}-J_{7}J_{5})A_{1},$
$\displaystyle
D_{7}=2(-9J_{8}J_{6}+4J_{8}J_{3}-3J_{8})A_{7}-8J_{7}J_{4}A_{6}+4(7J_{8}J_{6}-2J_{8}J_{3}+4J_{8}J_{0}+2J_{7}J_{4})A_{5}$
$\displaystyle\hphantom{D_{7}=}{}+4\bigl{(}2J_{8}J_{5}-9J_{7}J_{6}-4J_{7}J_{0}+6J_{4}^{2}+5J_{7}\bigr{)}A_{4}+8J_{8}J_{4}A_{3}$
$\displaystyle\hphantom{D_{7}=}{}+2(4J_{7}J_{1}-2J_{6}J_{3}-23J_{6}+4J_{0}+6J_{3}+23)A_{2}-2(4J_{8}J_{1}-9J_{6}J_{4})A_{1},$
$\displaystyle
D_{8}=-6J_{8}J_{4}A_{4}+\bigl{(}75J_{4}^{2}-27J_{7}J_{6}-2J_{7}J_{3}+4J_{7}J_{0}\bigr{)}A_{2}$
$\displaystyle\hphantom{D_{8}=}{}+(15J_{8}J_{6}-16J_{8}J_{3}+20J_{8}J_{0}+25J_{8})A_{1},$
$\displaystyle D_{9}=-18J_{8}^{2}A_{4}+18J_{8}J_{4}A_{2}-12J_{8}J_{7}A_{1},$
$\displaystyle D_{10}=-66J_{4}^{2}A_{2},\qquad
D_{11}=-48J_{8}J_{4}A_{2},\qquad D_{12}=-30J_{8}^{2}A_{2}.$

##  Appendix C Maple code examples

In this appendix, we will provide an example of Maple code which carry out two
test calculations: the commutator $[J_{4}J_{1},J_{5}J_{1}]$ and the commutator
of two artifacts $[A_{1},A_{2}]$. Ultimate goal is to get a combination of
lexicographically ordered monomials.

First upload Maple packages: Physics, Library and StringTools:

    > with(Physics);
     with(Library);
     with(StringTools);

Here we define the $\mathfrak{gl}(3)$ commutation relations among the
generators $j0,j1,\dots,j8$

    > Setup(noncommutativeprefix = j,
    %Commutator(j1, j0) = -j1,
    %Commutator(j2, j0) = -j2,
    %Commutator(j2, j1) = 0,
    %Commutator(j3, j0) = 0,
    %Commutator(j3, j1) = -j1,
    %Commutator(j3, j2) = 0,
    %Commutator(j2, j3) = 0,
    %Commutator(j4, j0) = 0,
    %Commutator(j4, j1) = 0,
    %Commutator(j4, j2) = -j1,
    %Commutator(j4, j3) = j4,
    %Commutator(j5, j0) = 0,
    %Commutator(j5, j1) = -j2,
    %Commutator(j5, j2) = 0,
    %Commutator(j5, j3) = -j5,
    %Commutator(j5, j4) = j3-j6,
    %Commutator(j5, j6) = j5,
    %Commutator(j5, j8) = j7,
    %Commutator(j6, j0) = 0,
    %Commutator(j6, j1) = 0,
    %Commutator(j6, j2) = -j2,
    %Commutator(j6, j3) = 0,
    %Commutator(j6, j4) = j4,
    %Commutator(j6, j5) = -j5,
    %Commutator(j6, j7) = 0,
    %Commutator(j7, j0) = j7,
    %Commutator(j7, j1) = j0-j3,
    %Commutator(j7, j2) = -j5,
    %Commutator(j7, j3) = -j7,
    %Commutator(j7, j4) = -j8,
    %Commutator(j7, j5) = 0,
    %Commutator(j8, j0) = j8,
    %Commutator(j8, j1) = -j4,
    %Commutator(j8, j2) = j0-j6,
    %Commutator(j8, j3) = 0,
    %Commutator(j8, j4) = 0,
    %Commutator(j8, j5) = -j7,
    %Commutator(j8, j6) = -j8,
    %Commutator(j8, j7) = 0);

$\displaystyle{\color[rgb]{0,0,1}[{\it
algebrarules}=\\{[j1,j0]_{-}=-j1,[j2,j0]_{-}=-j2,[j2,j1]_{-}=0,[j3,j0]_{-}=0,}$
$\displaystyle{\color[rgb]{0,0,1}[j3,j1]_{-}=-j1,[j3,j2]_{-}=0,[j4,j0]_{-}=0,[j4,j1]_{-}=0,[j4,j2]_{-}=-j1,[j4,j3]_{-}=j4,}$
$\displaystyle{\color[rgb]{0,0,1}[j5,j0]_{-}=0,[j5,j1]_{-}=-j2,[j5,j2]_{-}=0,[j5,j3]_{-}=-j5,[j5,j4]_{-}=j3-j6,}$
$\displaystyle{\color[rgb]{0,0,1}[j6,j0]_{-}=0,[j6,j1]_{-}=0,[j6,j2]_{-}=-j2,[j6,j3]_{-}=0,[j6,j4]_{-}=j4,[j6,j5]_{-}=-j5,}$
$\displaystyle{\color[rgb]{0,0,1}[j6,j7]_{-}=0,[j7,j0]_{-}=j7,[j7,j1]_{-}=j0-j3,[j7,j2]_{-}=-j5,[j7,j3]_{-}=-j7,}$
$\displaystyle{\color[rgb]{0,0,1}[j7,j4]_{-}=-j8,[j7,j5]_{-}=0,[j8,j0]_{-}=j8,[j8,j1]_{-}=-j4,[j8,j2]_{-}=j0-j6,}$
$\displaystyle{\color[rgb]{0,0,1}[j8,j3]_{-}=0,[j8,j4]_{-}=0,[j8,j5]_{-}=-j7,[j8,j6]_{-}=-j8,[j8,j7]_{-}=0\\},}$
$\displaystyle{\color[rgb]{0,0,1}noncommutativeprefix=\\{j\\}]}$

Here we define the lexicographic ordering of the $j$-generators via the list
“jorder”

    > jorder := [j8, j7, j6, j5, j4, j3, j2, j1, j0]:

Here we define the function OR which orders the $J$-generators in every term
of an expression in the lexicographic order defined above.

The function OR uses the command SortProducts which is a command of the
Physics Library which acts on a Maple expression containing noncommutative
products. It uses the list “jorder” which defines the ordering of the
$J$-generators. The option “usecommutator” indicates that commutators defined
previously.

    > OR := proc (x)
    Simplify(SortProducts(Expand(x), jorder, usecommutator))
    end proc;

Examples: product and commutator of two monomials

    > monomial1 := j4 * j1:
    > monomial2 := j5 * j1:

Products of monomials

    > productmonomials12 := monomial1 * monomial2;

${\color[rgb]{0,0,1}{\it productmonomials12:=}j4j1j5j1}$

    > productmonomials21 := monomial2 * monomial1;

${\color[rgb]{0,0,1}{\it productmonomials21:=}j5j1j4j1}$

We apply the function OR to reorder the product of monomials in lexicographic
order

    > OR(productmonomials12);

${\color[rgb]{0,0,1}j5j4j1^{2}+j4j1j2-j3j1^{2}+j1^{2}j6}$

When two generators commute, like $j1$ and $j2$, we need to use SortProducts
with the desired ordering list: $[j2,j1]$

    > orderedproductmonomials12:=SortProducts(SortProducts(OR(productmonomials12),
    [j2, j1], usecommutator), [j6, j1], usecommutator);

${\color[rgb]{0,0,1}orderedproductmonomials12:=j5j4j1^{2}+j4j2j1-j3j1^{2}+j6j1^{2}}$

    > orderedproductmonomials21:= OR(productmonomials21);

${\color[rgb]{0,0,1}{\it orderedproductmonomials21:=}j5j4j1^{2}}$

Commutator of monomials using the Physics Library command “Commutator”:

    > commutatormonomials12 := Commutator(monomial1, monomial2);

${\color[rgb]{0,0,1}{\it commutatormonomials12:=}j4j2j1+(-j3+j6)j1^{2}}$

Reordering of terms in commutator using function OR and SortProducts

    > commutatormonomials12 := SortProducts(SortProducts(OR(commutatormonomials12),
    [j2, j1], usecommutator), [j6, j1], usecommutator);

${\color[rgb]{0,0,1}{\it commutatormonomials12:=}j4j2j1-j3j1^{2}+j6j1^{2}}$

Defining artifacts

    > A[1] := j8*j5-j7*j6:
     A[2] := j8*j3-j7*j4:
     A[3] := j7*j2+j5*j0+j5:
     A[4] := j8*j1+j4*j0+j4:
     A[5] := j7*j1+j3*j0+j3:
     A[6] := j8*j2+j6*j0+j6:
     A[7] := j6*j3-j5*j4+j3:
     A[8] := j6*j1-j4*j2:
     A[9] := j5*j1-j3*j2:

Commutator of artifacts $A[1]$, $A[2]$:

    > CommutatorA1A2 := SortProducts(SortProducts(OR(Commutator(A[1], A[2])),
     [j8, j7], usecommutator), [j7, j4], usecommutator)

${\color[rgb]{0,0,1}{\it
CommutatorA1A2:=}-j8^{2}j5+j8j7j6-j8j7j3+j7j4j7-j8j7}$

     > AnsatzcommutatorA1A2:= ’-j8*A[1]-j7*A[2]’;

${\color[rgb]{0,0,1}{\it AnsatzcommutatorA1A2:=}-j8A_{1}-j7A_{2}}$

Testing commutator of artifacts $A[1]$, $A[2]$:

     > OR(Simplify(CommutatorA1A2 - AnsatzcommutatorA1A2));

${\color[rgb]{0,0,1}0}$

Note added in proof. We know that $A_{2}$ rational system, see (1.8) for the
Hamiltonian $h_{3}\equiv h_{{\cal A}_{2}}(x,y)$ and (1.9) for the cubic
integral $k_{3}\equiv k_{{\cal A}_{2}}(x,y)$ in the algebraic form at
$\mu=\tau=0$, is integrable, $[h_{3},k_{3}]=0$; however, this rational system
admits separation of variables in polar coordinates in the space of relative
motion and, hence, the existence of the additional quadratic integral of
motion $x_{3}$, see [7]: $[h_{3},x_{3}]=0$ with property
$I\equiv[x_{3},k_{3}]\neq 0$; hence, $[h_{3},I]=0$. Thus, the $A_{2}$ rational
integrable system is superintegrable. Recently, it was shown in [2] that
double commutators $[x_{3},I]=P_{3}(h_{3},k_{3},x_{3},I)$ and
$[k_{3},I]=Q_{3}(h_{3},k_{3},x_{3},I)$ are cubic polynomials in
$(h_{3},k_{3},x_{3},I)$. Hence, we arrive at 4-generated, infinite-dimensional
cubic polynomial algebra of integrals.

### Acknowledgments

A.V.T. is thankful to Willard Miller Jr.(1937–2023) and Peter Olver
(University of Minnesota, USA) for helpful discussions in different stages of
the project and the general encouragement to proceed and to complete this
work. Due to enormous computational complexity, this research was running for
many years, it was supported in part by the PAPIIT grants IN109512 and
IN108815 (Mexico) at the initial stage of the study and by the PAPIIT grant
IN113022 (Mexico) at its final stage. M.A.G.A. thanks the CONACyT grant for
master degree studies (Mexico) in 2016–2018, when the key calculations of the
commutator (6.3) were partially carried out.

All symbolic calculations were carried out by using the MAPLE-18 on a regular
DELL desktop computer with CPU processor Intel(R) Core (TM) i7-3770 $@$
3.40GHz with 6Gb RAM, although some pieces of calculations were done on a
regular PC laptop.

A.V.T. thanks PASPA-UNAM grant (Mexico) for its support during his sabbatical
stay in 2021–2022 at the University of Miami, where this work was mostly
completed. We thank all anonymous referees for careful reading of the text,
many inspiring comments, remarks and proposals, which improved significantly
the presentation.

This work is dedicated to the 70th birthday of Peter Olver to whom we always
had admiration as an exemplary mathematician and scientist.

## References

  * [1] Chryssomalakos C., Turbiner A.V., Canonical commutation relation preserving maps, J. Phys. A 34 (2001), 10475–10485, arXiv:math-ph/0104004.
  * [2] Lopez Vieyra J.C., Turbiner A.V., Wolfes model aka $G_{2}/I_{6}$-rational integrable model: $\mathfrak{g}^{(2)}$, $\mathfrak{g}^{(3)}$ hidden algebras and quartic polynomial algebra of integrals, arXiv:2310.20481.
  * [3] Olshanetsky M.A., Perelomov A.M., Quantum integrable systems related to Lie algebras, Phys. Rep. 94 (1983), 313–404.
  * [4] Oshima T., Completely integrable systems associated with classical root systems, SIGMA 3 (2007), 061, 50 pages, arXiv:math-ph/0502028.
  * [5] Smirnov Yu.F., Turbiner A.V., Lie algebraic discretization of differential equations, Modern Phys. Lett. A 10 (1995), 1795–1802, arXiv:funct-an/9501001.
  * [6] Sokolov V.V., Turbiner A.V., Quasi-exact-solvability of the $A_{2}/G_{2}$ elliptic model: algebraic forms, $\mathfrak{sl}(3)/\mathfrak{g}^{(2)}$ hidden algebra, polynomial eigenfunctions, J. Phys. A 48 (2015), 155201, 15 pages, Corrigendum, J. Phys. A 48 (2015), 359501, 2 pages, arXiv:1409.7439.
  * [7] Tremblay F., Turbiner A.V., Winternitz P., An infinite family of solvable and integrable quantum systems on a plane, J. Phys. A 42 (2009), 242001, 10 pages, arXiv:0904.0738.
  * [8] Turbiner A.V., Lamé equation, $\mathfrak{sl}(2)$ algebra and isospectral deformations, J. Phys. A 22 (1989), L1–L3.
  * [9] Turbiner A.V., Lie-algebras and linear operators with invariant subspaces, in Lie Algebras, Cohomology, and New Applications to Quantum Mechanics (Springfield, MO, 1992),Contemp. Math., Vol. 160, American Mathematical Society, Providence, RI, 1994, 263–310, arXiv:funct-an/9301001.
  * [10] Turbiner A.V., Different faces of harmonic oscillator, in SIDE III—Symmetries and Integrability of Difference Equations (Sabaudia, 1998), CRM Proc. Lecture Notes, Vol. 25, American Mathematical Society, Providence, RI, 2000, 407–414, arXiv:math-ph/9905006.
  * [11] Turbiner A.V., The Heun operator as a Hamiltonian, J. Phys. A 49 (2016), 26LT01, 8 pages, arXiv:1603.02053.
  * [12] Turbiner A.V., Miller Jr. W., Escobar-Ruiz M.A., From two-dimensional (super-integrable) quantum dynamics to (super-integrable) three-body dynamics, J. Phys. A 54 (2021), 015204, 10 pages, arXiv:1912.05726.
  * [13] Turbiner A.V., Vasilevski N., Poly-analytic functions and representation theory, Complex Anal. Oper. Theory 15 (2021), 110, 24 pages, arXiv:2103.12771.
  * [14] Whittaker E.T., Watson G.N., A course of modern analysis, Cambridge Math. Lib., Cambridge University Press, Cambridge, 1996.