Integrable and superintegrable quantum mechanical systems with position
dependent masses invariant with respect to one parametric Lie groups. 1.
Systems with cylindric symmetry.

A. G. Nikitin 111E-mail<EMAIL_ADDRESS>Institute of Mathematics,
National Academy of Sciences of Ukraine,
3 Tereshchenkivs’ka Street, Kyiv-4, Ukraine, 01024, and

Universit´a del Piemonte Orientale,

Dipartimento di Scienze e Innovazione Tecnologica,

viale T. Michel 11, 15121 Alessandria, Italy

Cylindrically symmetric quantum mechanical systems with position dependent
masses (PDM) admitting at least one second order integral of motion are
classified. It is proved that there exist 68 such systems which are
inequivalent. Among them there are twenty seven superintegrable and twelve
maximally superintegrable. The arbitrary elements of the correspondinding
Hamiltonians (i.e.,masses and potentials) are presented explicitly.

##  1 Introduction

Integrals of motion belong to very important part of classical and quantum
mechanics. Just the existence of the sufficient number of integrals of motion
for a Hamiltonian system makes it integrable or exactly solvable, and we are
not supposed to calculate its approximate solutions to describe its behavior.
Moreover, sometimes it is possible to describe the main physical properties of
the mentioned system using its integrals of motion an d ignoring the motion
equations. A classical example of such situation was presented by Pauli who
had calculated the energy spectrum of the Hydrogen atom using its integrals of
motion forming the Laplace-Runge-Lentz vector [1]. Ad it was done before the
discovery of Schrödinger equation!

The systematic search for integrals of motion admitted by the Schrödinger
equation equations started with papers [2], [3] were all inequivalent second
order integrals of motion for the 2d one particle quantum systems had been
classified. And it needed 24 years to extend this result to the 3d case [4],
[5].

Papers [2] and [3] were indeed seminal. They were followed by a great number
of research works, see, e.g. survey [6]. In particular, integrable and
superintegrable system with matrix potentials have been classified for both
spin-orbit [7, 8] and Pauli type interactions [9, 10]. One of the modern
trends is to search for the third and even arbitrary order integrals of motion
[11], see also [12] where the determining equations for such integrals were
presented.

The higher (in particular, second) order integrals of motion are requested for
description of systems admitting solutions in separated variables [13] just
such integrals of motion characterize integrable and superintegrable systems
[6]. Let us mention also the nice conjecture of Marquette and Winternitz [14]
which predicts a surprising connection of higher order superintegrability with
the soliton theory.

In any case the integrable and superintegrable systems of the standard quantum
mechanics belong to the well developed research field which, however, still
have some white spots. Among them is the classification of arbitrary order
symmetry operators of generic form, which still are described only at the
level of determining equations [12].

One more important research field which is closely related to the mentioned
one is the classification of second order integrals of motion for quantum
mechanical systems with position dependent mass. Such systems are used in many
branches of modern theoretical physics, whose list can be found, e.g., in [15,
16]. Their symmetries and integrals of motion were studied much less than
those ones for the standard SE. However, the classification of Lie symmetries
of the PDM Schrödinger equations with scalar potentials have been obtained
already [17, 18, 19].

Second order integrals of motion for 2d PDM SEs are perfectly classified [20,
21, 22, 23]. The majority of such systems admits also at least one continuous
Lie symmetry. The two dimensional second-order (maximally) superintegrable
systems for Euclidean 2-space had been classified even algebraic geometrically
[24].

The situation with the 3d systems is still indefinite. At the best of my
knowledge the completed classification results were presented only for the
maximally superintegrable (i.e., admitting the maximal possible number of
integrals of motion) systems [25, 26], and (or) for the system whose integrals
of motion are supposed to satisfy some special conditions like the
functionally linearly dependence [27]. The nondegenerate systems, i.e, those
ones which have have 5 linearly independent, contained in 6 linearly
independent (but functionally dependent) 2nd order integrals of motion are
known [28], see also [29] for the contemporary trends in this field. In
addition, a certain progress can be recognized in the classification of the so
called semidegenerate systems which admit five linearly independent integrals
of motion and whose potentials are linear combinations of three functionally
independent terms [30].

Surely, the maximally superintegrable systems are very important and
interesting. In particular, they admit solutions in multi separated
coordinates [31, 32, 33, 34]. On the other hand, there are no reasons to
ignore the PDM systems which admit second order integrals of motion but are
not necessary maximally superintegrable. And just such systems are the subject
of our study.

In view of the complexity of the total classification of integrals of motion
for 3d PDM systems it is reasonable to separate this generic problem to well
defined subproblems which can have their own values. The set of such
subproblems can be treated as optimal one if solving them step by step we can
obtain the complete classification of PDM systems admitting integrals of
motion. We choose the optimal set of subproblems in the following way.

The first subproblem consists in the classification of the PDM systems
admitting the first order integrals of motion. This problem is already solved,
refer to [17].

The first order integrals of motion are nothing but generators of Lie
symmetries. The important aspect of the results presented in [17] is the
complete description of possible Lie symmetry groups which can be admitted by
the stationary PDM Schrödinger equation. And this property, i.e., the presence
of Lie symmetry, can be effectively used to separate the the problem of the
classification of the PDM systems admitting second order integrals of motion
for PDM systems to a well defined subproblems corresponding to the fixed
symmetries.

It was shown in [17] the PDM Schrödinger equation can admit six, four, three,
two or one parametric Lie symmetry groups. In addition, there are also such
equations which have no Lie symmetry. In other words, there are six well
defined classes of such equations which admit n-parametric Lie groups with
$n=6,4,3,2,1$ or do not have any Lie symmetry. And it is a natural idea to
search for second order integrals of motion consequently for all these
classes.

The systems admitting six- or four-parametric Lie groups are not too
interesting since the related Hamiltonians cannot include non-trivial
potentials. That is why we started our research with the case of three-
parametric groups. The classification of the corresponding PDM systems
admitting second order integrals of motion was obtained in [35]. There were
specified 38 inequivalent PDM systems together with their integrals of motion.
The majority of them are new systems which are not not maximally
superintagrable.

Notice that the superintegrable 3d PDM systems invariant with respect to the
3d rotations have been classified a bit earlier in paper [36] where their
supersymmetric aspects were discussed also. For relativistic aspects of
superintegrability see [37], [38].

The systems admitting two-parametric Lie groups and second order integrals of
motion have been classified in [39]. We again find a number of new systems in
addition to the known maximally superintegrable ones.

The natural next step is to classify the systems which admit one parametric
Lie groups and second order integrals of motion. As it is shown in [17], up to
equivalence these groups are reduced to dilatations, shifts along the fixed
coordinate axis, rotations around this axis and some specific combinations of
the mentioned transformations. We will conventionally call them the natural
and combined symmetries respectively.

In the present paper we start the systematic search for integrable and
superintegrable PDM systems admitting one parametric Lie symmetry groups.
Namely, we present the classification of the mentioned systems which are
invariant w.r.t. the rotations around the fixed coordinate axis. The number of
such systems appears to be rather extended, and to keep the reasonable size of
paper we restrict ourselves to this particular symmetry, while the systems
with the other Lie symmetries will be presented in the following papers.
Notice that the number of inequivalent systems which admit the other one
parametric symmetry groups is much less extended than in the case of the
cylindric symmetry.

In spite of the fact that the usual strategy in studying of superintegrable
systems with PDM is to start with the classical Hamiltonian systems an then
quantize them if necessary, we deal directly with quantum mechanical systems.
This way is more complicated but it guaranties obtaining of all integrals of
motion including those ones which can disappear in the classical limit [40].

The main result of the present paper consist in the complete classification of
integrable, superintegrable and maximally superintegrable PDM systems with
cylindric symmetry. In addition, we optimise the algorithm of solution of the
related determining equations which can be used for a classification of other
PDM systems.

##  2 Formulation of the problem

We are studying the stationary Schrödinger equations with position dependent
mass of the following generic form:

$\displaystyle H\psi=E\psi,$ (1)

where

$\displaystyle H=p_{a}f({\bf x})p_{a}+V({\bf x}),$ (2)

${\bf x}=(x^{1},x^{2},x^{3}),$ $p_{a}=-i\partial_{a}$, $V({\bf x})$ and
$f({\bf x})=\frac{1}{2m({\bf x})}$ are functions associated with the potential
and inverse PDM, and summation from 1 to 3 is imposed over the repeating index
$a$.

The particular form (2) of the hamiltonian is convenient for study of its
symmetries and integrals of motion. Moreover, more generic formulations
including the arbitrary ambiguity parameters (refer, e.g. to [16]) are
mathematically equivalent to (2).

In paper [17] all equations (1) admitting at least one first order integral of
motion was found. Such integrals of motion generate Lie groups which leave the
equations invariant. In accordance with [17] there are six inequivalent Lie
symmetry groups which can be accepted by the PDM Schrödinger equations. They
include three ”natural” groups, rotation around the third coordinate axis,
shift along this axis and dilatation groups. In addition, we can fix three
combined symmetries which are superpositions of rotations and shifts,
rotations and dilatations, and shifts, rotations and conformal
transformations.

The generic form of the corresponding inverse masses $f$ and potentials $V$
can be represented by the following formulae [17]:

$\displaystyle f=F(\tilde{r},x_{3}),\quad V=G(\tilde{r},x_{3}),$ (3)
$\displaystyle f=F(x_{1},x_{2}),\quad V=V(x_{1},x_{2}),$ (4) $\displaystyle
f=r^{2}F(\varphi,\theta),\quad V=V(\varphi,\theta)$ (5)

where $F(.)$ and $V(.)$ are arbitrary functions whose arguments are fixed in
the brackets,

$\displaystyle
r=(x^{2}_{1}+x_{2}^{2}+x_{3}^{2})^{\frac{1}{2}},\quad\tilde{r}=(x^{2}_{1}+x_{2}^{2})^{\frac{1}{2}},\quad\varphi=\arctan\left(\frac{x_{2}}{x_{1}}\right),\
\ \theta=\arctan\left(\frac{\tilde{r}}{x_{3}}\right).$

In the present paper we classify the PDM systems which admit one of the
mentioned natural symmetries, namely the rotations around the third coordinate
axis, and, in addition, have at least one second order integral of motion. The
generic form of the corresponding inverse masses $f$ and potentials $V$ are
represented in (3).

Equations (1), (2) with arbitrary parameters presented in (3) admit the
following first order integrals of motion

$\displaystyle L_{3}=x_{1}p_{2}-x_{2}p_{1}$ (6)

which is nothing but the third component of the orbital momentum.

Our goal is to fix such systems (1), (3) which, in addition to their Lie
symmetries, admit second order integrals of motion whose generic form is:

$Q=\partial_{a}\mu^{ab}\partial_{b}+\eta$ (7)

where $\mu^{ab}=\mu^{ba}$ and $\eta$ are unknown functions of $\bf x$ and
summation from 1 to 3 is imposed over all repeating indices.

Operators (7) are formally hermitian. In addition, just representation (7)
leads to the most compact and simple systems of determining equations for
unknown parameters $\mu^{ab}$ and $\eta$.

By definition, operators $Q$ should commute with $H$:

$[H,Q]\equiv HQ-QH=0.$ (8)

Evaluating the commutator in (8) and equating to zero the coefficients for the
linearly independent differential operators
$\partial_{a}\partial_{b}\partial_{c}$ and $\partial_{a}$ we come to the
following determining equations

$\displaystyle
5\left(\mu^{ab}_{c}+\mu^{ac}_{b}+\mu^{bc}_{a}\right)=\delta^{ab}\left(\mu^{nn}_{c}+2\mu^{cn}_{n}\right)+\delta^{bc}\left(\mu^{nn}_{a}+2\mu^{an}_{n}\right)+\delta^{ac}\left(\mu^{nn}_{b}+2\mu^{bn}_{n}\right),$
(9) $\displaystyle\left(\mu^{nn}_{a}+2\mu^{na}_{n}\right)f-5\mu^{an}f_{n}=0,$
(10) $\displaystyle\mu^{ab}V_{b}-f\eta_{a}=0$ (11)

where $\delta^{bc}$ is the Kronecker delta, $f_{n}=\frac{\partial f}{\partial
x_{n}},\ \mu^{an}_{n}=\frac{\partial\mu^{an}}{\partial x_{n}}$, etc., and
summation is imposed over the repeating indices $n$ over the values $n=1,2,3$.

Equations (9) , (10) and (11) give the necessary and sufficient conditions for
commutativity of operators $H$ (2) and $Q$ (7) [39].

##  3 Evolution of the determining equations

A particular solution of equations (9) is $\mu^{ab}=\mu_{0}^{ab}$ where

$\displaystyle\mu_{0}^{ab}=\delta^{ab}g({\bf r})$ (12)

with arbitrary function $g({\bf r})$.

Whenever tensor $\mu_{0}^{ab}$ is nontrivial, the determining equations (10)
represent the coupled system of three nonlinear partial differential equation
equations for two unknowns $g(x)$ and $f({\bf x})$. Fortunately, this system
can be linearizing by introduction of the new dependent variables

$\displaystyle M=\frac{1}{f},\ N=\frac{g}{f}$ (13)

which reduces (10) to the following form:

$\displaystyle\left(\mu^{nn}_{a}+2\mu^{na}_{n}\right)M+5\mu^{an}M_{n}+N_{a}=0.$
(14)

We see that for variables (13) the determining equation (10) is linear.
Equation (11) in its turn can be effectively linearised by introducing the
following new dependent variables $\tilde{M}$ and $R$:

$\displaystyle\tilde{M}=MV,\ R=\frac{N\tilde{M}}{M}-\eta$ (15)

which reduce (11) to the following equation:

$\displaystyle\left(\mu^{nn}_{a}+2\mu^{na}_{n}\right)\tilde{M}+5\mu^{an}\tilde{M}_{n}+R_{a}=0$
(16)

which simple coincides with (14). Surely, it does not mean that $M$ and $N$
coincide with $\tilde{M}$ and $\tilde{N}$ respectively, since these functions
can include different arbitrary elements, say, different integration
constants. In accordance with (13) and (15) the related inverse mass potential
have the form:

$\displaystyle f=\frac{2}{M},\ V=\frac{\tilde{M}}{M}$ (17)

where $\tilde{M}$ and $M$ are different solutions of the same equation, while
the corresponding functions $g$ and $\eta$ are expressed via $M$ and
$\tilde{M}$ in the following manner:

$\displaystyle g=\frac{N}{M},\ \eta=-\frac{N\tilde{M}}{M}-R.$ (18)

We see that to find the admissible inverse mass and potential it is sufficient
to solve the only linear equation (14) and then find the desired functions
$f,V,g$ and $\eta$ using definitions (17) and (18).

Just linearised determining equation (13) together with the mentioned
definitions will be used in the following to solve our classification problem.

Let us represent generic integral of motion (7) in terms of new dependent
variables $M,N,\tilde{M}$ and $\tilde{N}$ (refer to (13) and (15))

$\displaystyle Q=P_{a}\mu^{ab}P_{b}+(N\cdot H)-R$ (19)

where we denote

$\displaystyle(N\cdot H)=P_{a}(Nf)P_{a}+NV.$ (20)

The latter definition includes a Hermitized product of function $N$ with
Hamiltonian (2).

It is necessary to note than whenever function $g({\bf r})$ is equal to zero,
i.e., tensor $\mu_{0}^{ab}$ is trivial, the above presented speculations are
forbidden. We still can deal with the determining equations (14), but we are
supposed to deal with the initial determining equations (11) instead of (16).

##  4 Equivalence relations

An important step of our classification problem is the definition of
equivalence relations which will be presented in this section.

Non degenerated changes of dependent and independent variables of equations
(1), (2) are called equivalence transformations provided they keep their
generic form up to the explicit expressions for the arbitrary elements $f$ and
$V$. They have the structure of a continuous group which however can be
extended by some discrete elements. Let us remind that a particular subset of
the equivalence transformations are invariance transformation which by
definition keep the mentioned arbitrary elements uncharged.

It was shown in [17] that the maximal continuous equivalence group of equation
(1) is the group of conformal transformations of the 3d Euclidean space which
we denote as C(3). The corresponding Lie algebra is a linear span of the
following first order differential operators [17]:

$\displaystyle\begin{split}&P^{a}=p^{a}=-i\frac{\partial}{\partial
x_{a}},\quad L^{a}=\varepsilon^{abc}x^{b}p^{c},\\\
&D=x_{n}p^{n}-\frac{3{\mathrm{i}}}{2},\quad
K^{a}=r^{2}p^{a}-2x^{a}D,\end{split}$ (21)

where $r^{2}=x_{1}^{2}+x_{2}^{2}+x_{3}^{2}$ and
$p_{a}=-i\frac{\partial}{\partial x_{a}}.$ Operators $P^{a},$ $L^{a},$ $D$ and
$K^{a}$ generate shifts, rotations, dilatations and pure conformal
transformations respectively.

In addition, equation (1) is form invariant with respect to the following
discrete transformations:

$\displaystyle x_{a}\to\tilde{x}_{a}=\frac{x_{a}}{r^{2}},\quad\psi({\bf
x})\to\tilde{x}^{3}\psi(\tilde{\bf x})$ (22)

where
$\tilde{x}=\sqrt{\tilde{x}_{1}^{2}+\tilde{x}_{2}^{2}+\tilde{x}_{3}^{2}}.$

Notice that the related Lie algebra c(3) is isomorphic to the algebra so(1,4)
whose basic elements $S_{\mu\nu}$ can be expressed via generators (21) as:

$\displaystyle S_{ab}=\varepsilon_{abc}L_{c},\quad
S_{4a}=\frac{1}{2}(K_{a}-P_{a}),\quad S_{0a}=\frac{1}{2}(K_{a}+P_{a}),\quad
S_{04}=D$ (23)

where $a,b=1,2,3.$ The related Lie group SO(1,4) is the Lorentz group in
(1+4)-dimensional space. The discrete transformation (22) anticommutes with
$S_{4a}$ and $S_{40}$ but commutes with the remaining generators (23). Thus
its action on operators (23) can be represented as follows:

$\displaystyle S_{4a}\to-S_{4a},\ S_{04}\to-S_{04},\ S_{0a}\to S_{0a},\
S_{ab}\to S_{ab}.$ (24)

The presented speculations are valid for an abstract system (1) which is free
of any additional constrains. However, for the systems whose arbitrary
elements satisfy condition (3) the equivalence group is reduced since it is
supposed that it does not change the invariance groups of these equations. It
means that the set of generators (23) should be reduced to such ones which
commute with $L_{3}$. There are four the generators satisfying this condition,
namely:

$\displaystyle\ P_{3},\ L_{3},\ K_{3},\ D.$ (25)

They generate the reduced equivalence algebra so(2,1)$\oplus$e(1) where e(1)
includes the only basis element $L_{3}$.

Thus, in comparison with the generic case, for special arbitrary elements
presented in (3) the admissible continuous equivalence transformations are
reduced to group SO(2,1) extended by rotations w.r.t. the third coordinate
axis. However, the admissible discrete transformations are extended by the
reflection of one out of two the independent variables $x_{1}$ and $x_{2}$,
i.e.,

$\displaystyle x_{1}\to-x_{1},\ x_{2}\to x_{2},\ x_{3}\to x_{3}$ (26)

or

$\displaystyle\ x_{1}\to x_{1},\ x_{2}\to-x_{2},\ x_{3}\to x_{3}$ (27)

which keep the related equation (1) invariant. These discrete transformations
can be added to the universal discrete transformation (22).

The presented equivalence relations will be used in the following to simplify
calculations and and to optimize the representation of the classification
results.

##  5 Identities in the extended enveloping algebra of c(3)

It was noted in [39] that integrals of motion (7) where $\mu^{ab}$ Killing
tensors, i.e., solutions of equations can be represented as bilinear
combinations of the basic elements of algebra c(3) (21) added by the special
term with $\mu^{ab}=\delta_{ab}g({\bf x})$ and potential term $\eta$. In other
words they admit the following representation:

$\displaystyle
Q=c^{\mu\nu,\lambda\sigma}\\{S_{\mu\nu},S_{\lambda\sigma}\\}+P_{a}g({\bf
x})P_{a}+\eta$ (28)

where $S_{\mu\nu}$ are generators (23) and $c^{\mu\nu,\lambda\sigma}$ are
numeric parameters.

In accordance with (28) the second order integrals of motion of the PDM
Schrödinger equation belong to the enveloping algebra of their equivalence
algebra, i.e., c(3), extended by special terms $p_{a}g({\bf x})p_{a}$. Notice
that the same is true for the any order integrals of motion with the
appropriate generalization of the extending term.

Representation (28) is very important. Being combined with the equivalence
relations discussed in the previous section it enables essentially simplify
both the calculations and the representation of their results. In addition to
the equivalence relations we will use numerous identities in the enveloping
algebra of algebra c(3) which take place for its particular realization (21).
These identities are presented in the following formulae:

$\displaystyle\\{P_{a},D\\}+\varepsilon_{abc}\\{P_{b},L_{c}\\}=2P_{c}x_{a}P_{c},$
(29) $\displaystyle L_{1}^{2}+L_{2}^{2}+L_{3}^{2}+D^{2}=P_{a}r^{2}P_{a},$ (30)
$\displaystyle\\{L_{a},L_{b}\\}+\\{P_{a},K_{b}\\}-\delta^{ab}(L_{1}^{2}+L_{2}^{2}+L_{3}^{2})=2Q^{ab}$
(31) $\displaystyle\begin{split}&P_{1}^{2}+P_{2}^{2}+P_{3}^{2}=P_{a}P_{a},\\\
&\\{P_{a},K_{b}\\}-\\{P_{b},K_{a}\\}=2\varepsilon_{abc}L_{c}D,\\\
&P_{1}L_{1}+P_{2}L_{2}+P_{3}L_{3}=0\end{split}$ (32)

where $Q^{ab}=-P_{c}x_{a}x_{b}P_{c}.$

The message given by relations (29)-(32) is that the terms in the l.h.s. can
be treated as linearly dependent whenever they are included into second order
integral of motion (28) since the latter one includes a yet indefinite term of
the kind presented in the r.h.s. of equations (29)-(32).

We will use relations (29)-(32) to produce maximally compact presentations for
the integrals of motion.

##  6 Solution of determining equations

The autonomous subsystem (9) defines the conformal Killing tensor which is the
fourth order polynomial in variables $x_{a}$ and includes an arbitrary
function which multiplies the Kronecker delta. The explicit expression for
this polynomial are presented in (12) and in the following formulae (see,
e.g., [41]) :

$\displaystyle\mu^{ab}_{1}=\lambda_{1}^{ab},$ (33)
$\displaystyle\mu^{ab}_{2}=\lambda_{2}^{a}x^{b}+\lambda_{2}^{b}x^{a}-2\delta^{ab}\lambda_{2}^{c}x^{c},$
(34)
$\displaystyle\mu^{ab}_{3}=(\varepsilon^{acd}\lambda_{3}^{cb}+\varepsilon^{bcd}\lambda_{3}^{ca})x^{d},$
(35)
$\displaystyle\mu^{ab}_{4}=(x^{a}\varepsilon^{bcd}+x^{b}\varepsilon^{acd})x^{c}\lambda^{d}_{4},$
(36)
$\displaystyle\mu^{ab}_{5}=\delta^{ab}r^{2}+k(x^{a}x^{b}-\delta^{ab}r^{2}),$
(37)
$\displaystyle\mu^{ab}_{6}=\lambda_{6}^{ab}r^{2}-(x^{a}\lambda_{6}^{bc}+x^{b}\lambda_{5}^{ac})x^{c}-\delta^{ab}\lambda_{6}^{cd}x^{c}x^{d},$
(38)
$\displaystyle\mu^{ab}_{7}=(x^{a}\lambda_{7}^{b}+x^{b}\lambda_{7}^{a})r^{2}-4x^{a}x^{b}\lambda_{7}^{c}x^{c}+\delta^{ab}\lambda_{7}^{c}x^{c}r^{2},$
(39)
$\displaystyle\mu^{ab}_{8}=2(x^{a}\varepsilon^{bcd}+x^{b}\varepsilon^{acd})\lambda_{8}^{dn}x^{c}x^{n}-(\varepsilon^{ack}\lambda_{8}^{bk}+\varepsilon^{bck}\lambda_{8}^{ak})x^{c}r^{2},$
(40)
$\displaystyle\mu^{ab}_{9}=\lambda_{9}^{ab}r^{4}-2(x^{a}\lambda_{9}^{bc}+x^{b}\lambda_{9}^{ac})x^{c}r^{2}+(4x^{a}x^{b}+\delta^{ab}r^{2})\lambda_{9}^{cd}x^{c}x^{d}+\delta^{ab}\lambda_{9}^{cd}x^{c}x^{d}r^{2}$
(41)

where $r=\sqrt{x_{1}^{2}+x_{2}^{2}+x_{3}^{2}}$, $\lambda_{n}^{a}$ and
$\lambda_{n}^{ab}$ are arbitrary parameters, satisfying the conditions
$\lambda_{n}^{ab}=\lambda_{n}^{ba},\lambda_{n}^{bb}=0$.

Thus we have to search for solutions of the determining equations (10) where
$\mu^{ab}$ are linear combinations of tensors (12), (33)-(41).

Formulae (12) and (33)-(41) include an arbitrary function $g({\bf x})$ and 35
arbitrary parameters $\lambda_{n}^{a}$ and $\lambda_{n}^{ab}$,
$a,b=1,2,...,9$. In addition, we have eight coefficients which appear in an
arbitrary linear combination of tensors (12), and three more unknown function,
i.e., $f$, $V$ and $\eta$. Thus the problem of the complete classification of
the 3d PDM systems admitting second order integrals of motion looks to be
huge. However, our strategy is to solve it step by step for the systems
admitting three, two and one parametric Lie groups, and, finally, for the
systems which do not admit any Lie symmetry. The first two steps have been
already done in papers [35] and [39]. The third step is the subject of the
current paper.

We are studying the PDM systems which are invariant with respect to rotations
around the fixed axis (say, the third one). In accordance with (3) the
corresponding Hamiltonian (2) is reduced to the following form:

$\displaystyle H=p_{a}f(\tilde{r},x_{3})p_{a}+V(\tilde{r},x_{3})$ (42)

where $\tilde{r}=\sqrt{x_{1}^{2}+x_{2}^{2}}$.

As it was noted in Section 3 the equivalence group of equation (1) is reduced
to SO(2,1)$\otimes$E(1) provided the related Hamiltonian (2) has the reduced
form (42). The corresponding infinitesimal operators are presented in (25). In
addition, there is the discrete equivalence transformations (22), (26) and
(27).

However, the symmetry fixed above makes it possible to decouple the second
order integrals of motion to the following three subclasses: scalars, vectors
and second rank tensors with respect to rotations around the third coordinate
axis. We will consider them consequently.

###  6.1 Scalar integrals of motion

Since the searched integrals of motion include bilinear combinations of the
generators $S_{\mu\nu}$ of the extended Lorentz group SO(1,4), it is
relatively easy to select such of them which commute with $S_{12}=L_{3}$ and
so are scalars with respect to the rotations. First, all $S_{\mu\nu}$ with
$\mu,\nu\neq$ 1 and 2 commute with $S_{12}$ and so their bilinear combinations
are scalars. There exactly three such $S_{\mu\nu}$, namely, $S_{04},S_{03}$
and $S_{43}$. Surely, this list can be added by $S_{12}$. Thus we have
indicated the following scalars:

$\displaystyle\begin{array}[]{c}S_{04}^{2},\ S_{03}^{2},\ S_{34}^{2},\
S_{12}S_{04},\ S_{12}S_{03},\ S_{12}S_{34},\\\ S_{34}S_{04}+S_{04}S_{34},\
S_{03}S_{04}+S_{04}S_{03},\ S_{43}S_{03}+S_{03}S_{43}.\end{array}$ (45)

The above presented list does not include the squared generator $S_{12}$ which
commutes with all the considered Hamiltonians by definition.

One more collection of scalars which includes the sums w.r.t. the repeating
indices $n=1,2$ and $m=1,2$ is presented in the following formulae:

$\displaystyle\begin{split}&S_{na}S_{nb}+S_{nb}S_{na}\end{split}$ (46)

and

$\displaystyle\begin{split}&\varepsilon_{3nm}(S_{na}S_{mb}+S_{mb}S_{na})\end{split}$
(47)

where $a,b=0,3,4$.

Equation (1) including arbitrary parameters (3) is transparently invariant
with respect the inversion of variables $x_{1}$ or $x_{2}$. The consequence of
this observation is that its scalar integrals of motion can be decoupled to
two linearly independent parts including proper scalars and pseudo scalars.
The proper scalars are collected in (45) and (46) while the pseudo scalars are
present in (47).

The generic scalar integral of motion is a linear combination of operators
(45), (46), and the term $P_{a}g(\tilde{r},x_{3})P_{a}$ (refer to (28)). The
pseudo scalar ones are the linear spans of operators (47) added by the same
term. Moreover, many terms of this linear combination should be omitted in
view of relations (23) and (29)-(32). As a result we will obtain the reduced
set of scalars which we write in terms of operators (21):

$\displaystyle P_{3}L_{3},\ DL_{3},\\{L_{3},K_{3}\\},$ (48) $\displaystyle
P_{3}^{2},\\{P_{3},K_{3}\\},\ D^{2},\ K_{3}^{2},\ \\{P_{3},D\\},\
\\{K_{3},D\\}.$ (49)

Notice that equation (46) collects the pseudo scalars while the proper scalars
are presented in (49).

One more way to specify the scalar integrals of motion is to use the formal
expressions for Killing tensors given by relations (12)-(41). Such integrals
are generated by the mentioned tensors in two cases: when the only nonzero
parameters $\lambda^{a}_{n}$ and $\lambda^{ab}_{n}$ correspond to $a=3$ and
$a=b=3$. One more possibility is $\lambda^{11}=\lambda^{22}\neq 0$, but it is
reduced to the previous one in view to the traceless condition for constant
tensors $\lambda^{ab}_{n}$. Going over the mentioned values of arbitrary
parameters we recover relations (48) and (49).

Let us specify the inequivalent linear combinations of pseudo scalars (48).
First we note that all of them are products of operator $L_{3}$ and linear
combinations of the following operators:

$\displaystyle S_{04}=D,\ S_{03}=\frac{1}{2}(K_{3}+P_{3}),\
S_{43}=\frac{1}{2}(K_{3}-P_{3})$ (50)

which satisfy the following relations

$\displaystyle[S_{04},S_{03}]=-{\mathrm{i}}S_{43},\
[S_{04},S_{43}]=-{\mathrm{i}}S_{03},\ [S_{04},S_{43}]=-{\mathrm{i}}S_{04}$
(51)

and so form a basis of Lie algebra so(1,2). Moreover, operators (50) commute
with $L_{3}$. It means that the number of inequivalent pseudoscalars which are
nothing but the products of linear combinations of operators (50) with
operator $L_{3}$ is equal to the number of inequivalent subalgebras of algebra
so(1,2), since just these subalgebras generate our integrals of motion.

The inequivalent subalgebras of algebra $so(1,2)$ are one dimensional and
include the following basis elements:

$\displaystyle S_{03},\ S_{43},\ S_{03}\pm S_{43}.$ (52)

Moreover, up to discrete transformation (22) which changes the sign of
$S_{03}$ we can restrict ourselves to the positive sign in the last term in
(52). As a result we obtain the following inequivalent symmetries (28) which
we present together with the related Killing tensors (33)-(41):

$\displaystyle\hat{Q}_{1}=P_{3}L_{3}+P_{a}\tilde{q}_{1}(\tilde{r},x_{3})P_{a},\
\mu^{ab}=\mu_{3}^{ab}$ (53)

and

$\displaystyle\hat{Q}_{2}=(K_{3}+P_{3})L_{3}+P_{a}\tilde{q}_{2}(\tilde{r},x_{3})P_{a},\
\mu^{ab}=\mu_{3}^{ab}+\mu_{8}^{ab},$ (54)
$\displaystyle\hat{Q}_{3}=(K_{3}-P_{3})L_{3}+P_{a}\tilde{q}_{3}(\tilde{r},x_{3})P_{a},\
\mu^{ab}=\mu_{3}^{ab}-\mu_{8}^{ab}$ (55)

we remind that the only nonzero parameters in tensors $\mu_{3}^{ab}$ and
$\mu_{8}^{ab}$ are $\lambda^{33}_{3}$ and $\lambda^{33}_{8}$.

Let us specify one more pseudo scalar operator

$\displaystyle\hat{Q}_{4}=\\{D,L_{3}\\}+P_{a}\tilde{q}_{4}(\tilde{r},x_{3})P_{a},\
\mu^{ab}=\mu_{4}^{ab},\lambda_{4}^{1}\neq 0.$ (56)

Operator (56) is equivalent to (54) and can be ignored in the analysis of
integrable PDM systems. However, we cannot ignore it in the case of
superintegrable systems.

The specific arguments of functions $q(.)$ in (72) are caused by the requested
symmetry of (72) with respect to rotations around the third coordinate axis.

The next task is to specify all inequivalent proper scalars. In accordance
with (48 ) and (51 they belong to the enveloping algebra of algebra so(1,2).
Moreover, it follows from the first of equations (30) that the related Casimir
operator of algebra so(1,2) takes the following form:

$\displaystyle
S_{03}^{2}+S_{04}^{2}-S_{43}^{2}=P_{a}\tilde{r}^{2}P_{a}-L_{3}^{2}.$ (57)

Thus the considered integrals of motion are linear combination of basic
elements of the mentioned enveloping algebra:

$\displaystyle Q=\sum c_{\alpha}Q_{\alpha}$ (58)

where $\alpha=1,...,6$,

$\displaystyle\begin{split}&Q_{1}=(S_{03}^{2}-S_{43}^{2}),\
Q_{2}=(S_{03}^{2}+S_{43}^{2}),\ Q_{3}=\\{S_{03},S_{43}\\},\\\
&Q_{4}=\\{S_{03},S_{04}\\},\ Q_{5}=\\{S_{04},S_{43}\\},\
Q_{6}=S_{04}^{2}-S_{03}^{2}\end{split}$ (59)

and $c_{1},...,c_{6}$ are real constants.

Notice that relation (57) can be rewritten in terms of operators $Q_{1},Q_{2}$
and $Q_{6}$ in the following way:

$\displaystyle Q_{1}+3Q_{2}+2Q_{6}=2(P_{a}\tilde{r}^{2}P_{a}-L_{3}^{2}).$ (60)

The l.h.s. of relation (60) includes the symmetry operator $L_{3}^{2}$ which
commutes with Hamiltonian (2) by definition and the term
$P_{a}\tilde{r}^{2}P_{a}$ which can be included to the last term of the
generic integral of motion (28). It means that operators $Q_{1},Q_{2}$ and
$Q_{6}$ can be treated as linearly dependent and so one of the coefficients
$c_{1},c_{2}$ and $c_{3}$ can be nullified without loss of generality.

The next step in the simplification of expression (58) can be made using the
equivalence transformations generated by the Lie group whose Lie algebra is a
linear span of base elements (50). The transformations generated by $S_{04}$
are the Lorentz transformations which look as:

$\displaystyle\begin{split}&Q_{2}\to
Q_{2}\cosh(2\lambda)+Q_{3}\sinh(2\lambda),\\\ &Q_{3}\to
Q_{3}\cosh(2\lambda)+Q_{2}\sinh(2\lambda),\end{split}$ (61)
$\displaystyle\begin{split}&Q_{4}\to
Q_{4}\cosh(\lambda)+Q_{5}\sinh(\lambda)\\\ &Q_{5}\to
Q_{4}\cosh(\lambda)+Q_{5}\sinh(\lambda).\end{split}$ (62)

In addition, $Q_{6}$ and $Q_{2}$ are transformed to their linear combinations,
but such transformations can be ignored in view of relation (60).

Let us use transformations (61) to specify the inequivalent versions of
coefficients $c_{2}$ and $c_{3}$ present in formula (58). For
$c_{2}^{2}>c_{3}^{2}$ or $c_{2}^{2}<c_{3}^{2}$ we can nullify $c_{3}$ or
$c_{2}$ correspondingly, while for the special case $c_{4}=\pm c_{5}$ we can
restrict ourselves to the case $c_{4}=-c_{5}$ up to the discrete equivalence
transformations (22). So the inequivalent versions of the pairs
$(c_{4},c_{5})$ are:

$\displaystyle c_{2}\neq 0,c_{3}=0;\ c_{2}=0,c_{3}\neq 0;\ c_{2}=-c_{3}.$ (63)

Whenever the last version presented in (63) is true we can use transformations
(62) to specify the inequivalent versions of constants $c_{4}$ and $c_{5}$
which are analogous to the ones presented in (63):

$\displaystyle c_{4}\neq 0,c_{5}=0;\ c_{4}=0,c_{5}\neq 0;\ c_{4}=-c_{5}$ (64)

provided $c_{2}=-c_{3}.$

The next step is to use the rotation like transformation generated by $S_{43}$
which look as follows:

$\displaystyle\begin{split}&Q_{3}\to Q_{3}\cos(\omega)-Q_{5}\sin(\omega),\\\
&Q_{5}\to Q_{5}\cos(\omega)+Q_{3}\sin(\omega),\\\ \end{split}$ (65)
$\displaystyle\begin{split}&Q_{6}\to Q_{6}\cos(2\omega)-Q_{4}\sin(2\omega),\\\
&Q_{4}\to Q_{4}\cos(2\omega)+Q_{6}\sin(2\omega)\end{split}$ (66)

Let at least one of the constants $c_{3},c_{3}$ is nontrivial than
transformations (65) can be used to transform them to the following form:

$\displaystyle c_{5}=0\ \ c_{3}\neq 0.$ (67)

Alternatively, both $c_{5}$ and $c_{3}$ can be equal to zero.

Summarising the results presented in (63), (64) and (67) and taking into
account that in view of relation (60) coefficients $c_{1},c_{2},c_{3}$ can be
treated as linearly dependent we come to the following inequivalent versions
of operators (58) which we present together with the related Killing tensors:

$\displaystyle Q_{1}=P_{3}^{2}+P_{a}q_{1}(\tilde{r},x_{3})P_{a}+\eta^{1},\
\mu^{ab}=\mu_{1}^{ab},$ (68) $\displaystyle
Q_{2}=\\{P_{3},D\\}+P_{a}q_{2}(\tilde{r},x_{3})P_{a}+\eta^{2},\
\mu^{ab}=\mu_{2}^{ab},$ (69) $\displaystyle
Q_{3}=\\{P_{3},K_{3}\\}+P_{a}q_{3}(\tilde{r},x_{3})P_{a},\
\mu^{ab}=\mu_{6}^{ab}+\eta^{3},$ (70) $\displaystyle Q_{4}=P_{3}^{2}\pm
K_{3}^{2}+P_{a}q_{4}(\tilde{r},x_{3})P_{a}+\eta^{4},\
\mu^{ab}=\mu_{1}^{ab}\pm\mu_{9}^{ab},$ (71)
$\displaystyle\begin{split}&Q_{5}=(K_{3}\pm
P_{3})^{2}+P_{a}q_{5}(\tilde{r},x_{3})P_{a}+\eta^{5},\
\mu^{ab}=\mu_{1}^{ab}\pm 2\mu_{6}^{ab}+\mu_{9}^{ab},\end{split}$ (72)
$\displaystyle Q_{6}=\\{P_{3},K_{3}\pm
P_{3}\\}+P_{a}q_{6}(\tilde{r},x_{3})P_{a}+\eta^{6},\mu^{ab}=\mu_{6}^{ab}\pm\mu_{1}^{ab},$
(73) $\displaystyle Q_{7}=K_{3}^{2}\pm
P_{3}^{2}+2n\\{K_{3},P_{3}\\}+P_{a}q_{7}(\tilde{r},x_{3})P_{a}+\eta^{7},\
\mu^{ab}=\mu^{ab}_{9}-\mu^{ab}_{1}+2n\mu^{ab}_{6},$ (74) $\displaystyle
Q_{8}=\\{D,(K_{3}\pm P_{3})\\}+P_{a}q_{8}(\tilde{r},x_{3})P_{a}+\eta^{8},\
\mu^{ab}=\mu_{7}^{ab}\pm\mu_{1}^{ab}$ (75)

and nonzero parameters in tensors
$\mu^{ab}_{1},\mu^{ab}_{2},\mu^{ab}_{6},\mu_{7}^{ab}$ and $\mu^{ab}_{9}$ are
$\lambda_{1}^{33},\lambda_{2}^{33},\lambda_{6}^{33},\lambda_{7}^{3}$ and
$\lambda_{9}^{33}$.

Thus to classify the scalar integrals of motion it is sufficient to solve the
determining equations (10) and (11) where $\mu^{ab}$ are the Killing tensors
fixed in (68) - (75).

###  6.2 Pseudo scalar integrals of motion

Let us search for solutions of the above defined determining equations for
pseudo scalar integrals of motion.

Let us start which the pseudo scalar integrals of motion fixed in (53). The
related determining equations (14) are reduced to the following form:

$\displaystyle x_{2}\partial_{3}M=\partial_{1}N,\
x_{1}\partial_{3}M=-\partial_{2}N,\ \partial_{3}N=0$ (76)

and are solved by the following functions:

$\displaystyle M=\frac{\nu x_{3}+F(\tilde{r})}{\tilde{r}^{2}},\
N=\frac{\nu}{2}\varphi$ (77)

where $\varphi=\arctan\left(\frac{x_{2}}{x_{1}}\right)$ is the Euler angle,
$\nu$ is the integration constant, and $F(\tilde{r})$ is an arbitrary function
of $\tilde{r}=\sqrt{x_{1}^{2}+x_{2}^{2}}.$ The corresponding functions
$\tilde{M}$ and $\tilde{N}$ which generate potential $V$ and function $\eta$
are solutions of the same equation and so have the same generic forms as given
in (77), but in general with different parameter $\nu$ and different function
$F(\tilde{r})$. Substituting the obtained results into equations (18) we come
to functions $f$ and $V$ presented in Item 1 of Table 1.

The next pseudoscalar operator which we consider is fixed by relations (54)
and (55). The corresponding determining equations (14) for
$M=M(\tilde{r}^{2},x_{3})$ and $N=N(\tilde{r}^{2},x_{3})$ are reduced to the
form:

$\displaystyle\partial_{3}N=0,\ \partial_{4}N=0,$ (78)
$\displaystyle\partial_{\varphi}N=\tilde{r}^{2}(2x_{3}\tilde{r}^{2}\partial_{4}M+2x_{3}M)-\frac{1}{2}(\tilde{r}^{2}-x_{3}^{2}\pm
1)\partial_{3}M$ (79)

where
$\partial_{3}=\frac{\partial}{\partial_{x_{3}}},\partial_{4}=\frac{\partial}{\partial_{x_{4}}},x_{4}=\tilde{r}^{2}$.

In accordance with (78) function $N$ depends on $\varphi$ only. Since $M$ by
definition does not depend on this variable, to solve equation (79) $N$ should
be a linear function, i.e., $N=c\varphi.$ Then equation (79) is reduced to the
following one:

$\displaystyle
2x_{3}\tilde{r}^{2}\partial_{4}M+2x_{3}M-\frac{1}{2}((\tilde{r}^{2}-x_{3}^{2}\pm
1)\partial_{3}M=c.$ (80)

Equation (80) is easily integrable. Its solutions are presented in Items 2 and
3 of Table 1. Notice that these solutions are qualitative different for
different sighs before 1 in the formulae presented above. And these signs are
the same as signs for $P_{3}$ in equations (54) and (55).

Let us fix also the PDM system admitting symmetry $DL_{3}+...$ which is
generated by the Killing tensor (37) with the only nonzero parameter
$\lambda_{4}^{3}$.. This symmetry is equivalent to the symmetry
$L_{3}(P_{3}+K_{3})+...$ considered above, see equation (53). However, its
presentation will be useful in searching for the systems admitting more than
one second order integrals of motion. We will not present the calculation
details but give the corresponding equations (14):

$\displaystyle\begin{array}[]{l}x_{4}x_{3}\partial_{3}M+2x_{4}\partial_{4}M+2M+\partial_{\varphi}N=0,\\\
\partial_{3}N=0,\ \partial_{4}N=0,\ \partial_{\varphi\varphi}N=0\end{array}$
(83)

and their solutions:

$\displaystyle M=\frac{F(\theta)-2\nu\ln(\tilde{r})}{\tilde{r}^{2}}$ (84)

where $F(\theta)$ and $G(\theta)$ are arbitrary functions, $\varphi$ and
$\theta$ are the Euler angle, $\mu$ and $\nu$ are arbitrary parameters.

Just these functions together with the corresponding integrals of motion are
presented in Item 4 of Table 1.

Thus we have found all inequivalent systems with position dependent mass which
admit pseudo scalar integrals of motion. They are defined up to pairs of
arbitrary functions. The related Hamiltonians commute also with the third
component of angular momentum and so the found systems are integrable. We will
se that for some particular arbitrary functions these systems are
superintegrable.

###  6.3 True scalar integrals of motion

Consider now the integrals of motion which are invariant with respect to the
space reflections. Their generic form is given by formulae (68)–(74).

Let us represent the corresponding determining equations. They include
functions $M$ and $N$ which depend on $\tilde{r}$ and $x_{3}$, but in some
cases it is reasonable to treat them as functions of $r$ and $x_{3}$. To unify
the representation we will use the notations $\tilde{r}^{2}=x_{4}$ and
$r^{2}=x_{5}$ .

Substituting the Killing tensors fixed in (68)-(74) into (14) we come to tjhe
following equations:

$\displaystyle\text{ For }\ Q_{1}:\ \partial_{\tilde{r}}N=0,\
\partial_{3}(M+N)=0;$ (85) $\displaystyle\text{ For }\ Q_{2}:\
\begin{array}[]{l}\partial_{3}M+2\partial_{5}N=0,\\\
2\partial_{5}(x_{5}M+x_{3}N)+\partial_{3}(x_{3}M+N)-2M=0;\end{array}$ (88)
$\displaystyle\text{ For }\ Q_{3}:\
\begin{array}[]{l}2\partial_{4}(x_{3}^{2}M+N)-x_{3}p_{3}M=0,\\\
\partial_{3}(\tilde{r}^{2}M+N)-2x_{3}\tilde{r}^{2}\partial_{4}M=0;\end{array}$
(91) $\displaystyle\text{ For }\ Q_{4}:\
\begin{array}[]{l}\partial_{4}(N+4x_{3}^{2}\tilde{r}^{2}M)-W_{\pm}\partial_{3}M=0,\\\
\partial_{3}(N+W_{\pm}^{2}M)-4x_{3}\tilde{r}^{2}W_{\pm}\partial_{4}M=0,\\\
W_{\pm}=\tilde{r}^{2}\pm 1-x_{3}^{2}=x_{4}-x_{3}^{2}\pm 1;\end{array}$ (95)
$\displaystyle\text{ For }\ Q_{5}:\
\begin{array}[]{l}\partial_{4}(N+4(x_{3}^{2}\pm
1)\tilde{r}^{2}M)-x_{3}(W_{\pm})\partial_{3}M=0,\\\
\partial_{3}N+((W\pm)^{2}+4\tilde{r}^{2})\partial_{3}M-4\partial_{4}(x_{3}\tilde{r}^{2}M(W\mp
2));\end{array}$ (98) $\displaystyle\text{ For }\ Q_{6}:\
\begin{array}[]{l}(2x_{3}^{2}\partial_{4}M-x_{3}\partial_{3})M=\partial_{4}N,\\\
(x_{4}^{2}-1)(\partial_{3}-2x_{3}x_{4}\partial_{4})M=\partial_{3}N;\end{array}$
(101) $\displaystyle\text{ For }\ Q_{7}:\
\begin{array}[]{l}x_{3}(x_{4}-x_{3}^{2}+n)\partial_{3}M+2(N+2-2x_{3}^{2})\partial_{4}(x_{4}M)=\partial_{4}N,\\\
(x_{3}^{4}-2(x_{4}+n)x_{3}^{2}+x_{4}^{2}-4x_{4}\pm 1)\partial_{3}M\\\
-4x_{3}(x_{4}-x_{3}^{2}+n)\partial_{4}x_{4}M=-\partial_{3}N,\end{array}$ (105)
$\displaystyle\text{ For }\ Q_{8}:\ \begin{array}[]{l}(3x_{3}^{2}-x_{4}\pm
1)\partial_{3}M+8x_{3}\partial_{4}(x_{4}M)=2\partial_{4}N,\\\
2(x_{4}-3x_{3}^{2}\pm
1)(x_{4}\partial_{4}M+\partial_{3}(x_{3}M))=-\partial_{3}N\end{array}$ (108)

The systems (85) and (91) are easy solvable. They are solved by the following
functions:

$\displaystyle M=F(\tilde{r})+G(x_{3}),\ N=-G(x_{3})$

and

$\displaystyle M=\frac{F(\theta)+G(R)}{r^{2}},\ N=-F(\theta)$

where $F(.)$ and $\theta$ are arbitrary functions,
$\theta=\arctan\left(\frac{\tilde{r}}{x_{3}}\right)$ is the Euler angle. These
solutions generate the inverse masses and potentials represented in Items 5
and 6 of Table 1

Equations (95) are a bit more complicated. Excluding unknown variable $N$ we
obtain the following second order equation for $M$:

$\displaystyle\partial_{33}M-4\partial_{55}(x_{5}M)+2\partial_{5}M=0.$ (109)

The obtained partial differential equation with variable coefficients appears
to be exactly solvable. To discover its exact solutions it is reasonable to
represent them as

$\displaystyle M=\frac{P}{\sqrt{x_{5}}}=\frac{P}{r}$ (110)

where $P$ are polynomials in $r$ and $x_{3}$.

The only first order polynomial satisfying (109), (109) is $c_{1}x_{3}$, the
second order polynomial is $c_{2}rx_{3}+c_{3}(r^{2}+x_{3}^{2}),$ the third
order polynomial is $c_{4}(x_{3}^{3}+3x_{3}r^{2})+c_{5}(r^{3}+3rx_{3}^{2})$,
etc. In this way we can find the infinite (but countable) set of exact
solutions for $M$. A more difficult step is to recognize the fact that such
polynomials are only particular case of the generic solution
$P=F(r+x_{3})+G(r-x_{3})$ were $F(.)$ and $G(.)$ are arbitrary functions of
the arguments fixed in the brackets. It means that rather complicated second
order equation with variable coefficients presented in (109) can be reduced to
the D’Alambert equation if we represent the dependent variable in form (110)
and change the independent variables $(x_{3},x_{5})$ to
$(x_{3},\sqrt{x_{5}})$.

The corresponding inverse masses and potentials are presented in Item 7 of
Table 1.

The remaining tasks, i.e., the constructions of exact solutions for the
determining equations (95)-(108) appear to be much more difficult problems
which, however, are solvable.

Consider equations (95). Excluding variable $N$ we obtain the following
compatibility conditions for this system:

$\displaystyle\begin{array}[]{l}-(W_{\pm}^{2}\mp\tilde{r}^{3}x_{3}^{2})\partial_{34}M+x_{3}W_{\pm}(4x_{4}\partial_{4}M-\partial_{33}M)\\\
+(20x_{4}-8x_{3}(x_{3}\mp
1))\partial_{4}M+(3W_{\pm}+6x_{3}^{2})\partial_{3}M+12x_{3}M.\end{array}$
(113)

This rather complicated partial differential equation of second order can be
solved by the following trick. Let us choose new independent variables. To
find the first of them we solve equations (95) for $N=0$ and obtain the
following subclass of solutions for our problem:

$\displaystyle M=\frac{G\left(\frac{r^{2}\pm 1}{\tilde{r}}\right)}{(r^{2}\pm
1)^{2}\mp 4\tilde{r}^{2}}.$ (114)

This solution is valid for any arbitrary function $G\left(\frac{r^{2}\pm
1}{\tilde{r}}\right)$, in particular, for $G\left(\frac{r^{2}\pm
1}{\tilde{r}}\right)=1$, when

$\displaystyle M=\frac{1}{(r^{2}\pm 1)^{2}\mp 4\tilde{r}^{2}}.$ (115)

Thus it is reasonable to search for solutions with non-trivial $N$ in the
form:

$\displaystyle M=\frac{F(r,x_{3})}{(r^{2}\pm 1)^{2}\mp 4\tilde{r}^{2}}.$ (116)

Substituting this form into (95) we immediately recognize that the latter
equations turn to the identities provided $N=-F(r,x_{3})$ and

$\displaystyle F=F\left(\frac{r^{2}\mp 1}{x_{3}}\right).$ (117)

In this way we come to the solutions presented in Item 8 of Table 1.

In complete analogy with the above one can solve equations (98). The inverse
masses and potentials generated by solutions of these equations are presented
in Item 9 of Table 1.

Let us consider equations (101). Excluding $N$ we obtain the following
compatibility condition for this system:

$\displaystyle((x_{4}-x_{3}^{2}-1)\partial_{3}\partial_{4}+4x_{3}x_{4}\partial_{4}\partial_{4}-x_{3}\partial_{3}\partial_{3}-3\partial_{3}+8x_{3}\partial_{4})M(x_{4},x_{3})=0.$
(118)

Like (109) equation (118) can be reduced to the D’Alambert equation if we
chose the following new independent and dependent variables:

$\displaystyle x=r^{2}\pm 1,\ y=x_{3}^{2},\ \tilde{M}(x,y)=\sqrt{x^{2}\mp
4y}M(\tilde{r},x_{3}).$ (119)

As a result we come to the following generic solution for (118):

$\displaystyle M=\frac{F(x_{-})+G(\tilde{x}_{-})}{\sqrt{(r^{2}\pm 1)^{2}\mp
4x_{3}^{2}}}$ (120)

where $x_{-}=\sqrt{x^{2}\mp 4y}+x,\tilde{x}_{-}=\sqrt{x^{2}\pm 4y}-x$, $F(.)$
and $G(.)$ are arbitrary functions. The corresponding potential $V$ (17) looks
as:

$\displaystyle
V=\frac{\tilde{F}(x_{-})+\tilde{G}(\tilde{x}_{-})}{F(x_{-})+G(\tilde{x}_{-})}.$
(121)

For any fixed $F(x_{-})$ and $G(\tilde{x}_{-})$ we can solve equations (101)
and find functions $N=N_{F,G}$ corresponding to $M$ defined in (120).
Unfortunately, it is seemed be impossible to represent functions $N_{F,G}$ in
closed form for $F$ and $G(\tilde{x}_{-})$ arbitrary. However, we can do it at
least for some rather extended classes of $F(.)$ and $G$. In particular it is
the case if $F(x_{-})+G(\tilde{x}_{-})$ is a homogeneous function of
$x_{3}^{2}$ and $y$, say, for

$\displaystyle
F(x_{-})+G(\tilde{x}_{-})={x_{-}^{n}+(-1)^{n+1}\tilde{x}_{-}^{n}}=\Phi_{n}.$
(122)

Substituting (122) into (119) and integrating the corresponding equations (14)
we find the related functions $M$ and $N$ in closed form for $n$ arbitrary:

$\displaystyle\begin{split}&M=M_{n}=\frac{\Phi_{n}}{\sqrt{(r^{2}-1)^{2}+4x_{3}^{2}}},\\\
&N=N_{n}=2^{n}x_{3}^{2}M_{n-1}.\end{split}$ (123)

Thus we found the infinite (but countable) set of solutions of the determining
equations (14) which generate integrals of motion (73). Let us represent
explicitly some of them:

$\displaystyle\begin{split}&M_{1}=r^{2}-1,\ M_{-1}=\frac{1}{x_{3}^{2}}\\\
&M_{2}=(r^{2}-1)^{2}+x_{3}^{2},\\\
&M_{3}=((r^{2}-1)((r^{2}-1)^{2}+2x_{3}^{2})\\\
&M_{4}=(r^{2}-1)^{4}+3(r^{2}-1)x_{3}^{4}+x_{3}^{6},\\\
&M_{5}=(r^{2}-1)^{5}+4(r^{2}-1)^{3}x_{3}^{2}+(r^{2}-1)x_{3}^{4},\\\ &...\\\
&M_{n}=\sum_{m\leq
n}\frac{(n-m)!x^{n-2m}x_{3}^{2m}}{(n-2m)!2^{m}m!}.\end{split}$ (124)

Two more particular solutions are:

$\displaystyle M=\frac{1}{\tilde{r}^{2}}\ \text{ and }\
M=\frac{x_{3}^{2}+1}{\tilde{r}^{4}}.$ (125)

Notice that linear combinations of solutions (123), (124) and (125) also solve
equations (14), see Item 10 of Table 1.

The last class of symmetries which we are supposed to study is represented in
equation (74). The corresponding determining equations (105) are the most
complicated, but in many aspects they are analogous to ones requested for the
systems admitting integrals of motion of type (73).

Excluding $N$ we come to the following compatibility condition for the system
(105):

$\displaystyle\begin{array}[]{l}((x_{4}-x_{3}^{2}+n)^{2}-n^{2}\pm
1-4x_{3}^{2}x_{4})\partial_{3}\partial_{4}-(n+x_{4}-x_{3}^{2})(x_{3}\partial_{3}\partial_{3}-4x_{3}x_{4}\partial_{4}\partial_{4})\\\
+3(n+x_{4}-3x_{3}^{2}r)\partial_{3}-(8n-8x_{3}^{2}+20x_{4})x_{3}\partial_{4}-12x_{3})M=0.\end{array}$
(128)

This frighteningly looking equation appears to be equivalent to the d’Alambert
one if we choose the new independent and dependent variables

$\displaystyle\begin{split}&y=\frac{1}{x_{4}}(\sqrt{z_{\pm}^{2}+2x_{3}^{2}x_{4}}+z_{\pm}),\
\tilde{y}=\frac{1}{x_{4}}(\sqrt{z_{\pm}^{2}+2x_{3}^{2}x_{4}}-z_{\pm}),\\\
&z_{\pm}=\frac{x_{4}^{2}\pm 1+2n(x_{4}-x_{3}^{2})}{2\sqrt{2(n^{2}\mp 1)}},\ \
1-n^{2}\leq
0,\hat{M}(y,\tilde{y})=\sqrt{z_{\pm}^{2}+2x_{3}^{2}x_{4}}M(\tilde{r},x_{3}).\end{split}$

As a result equation (128) takes the canonical form

$\displaystyle\partial_{y}\partial_{\tilde{y}}\hat{M}=0$ (129)

and so the related mass function has the following form:

$\displaystyle
M=\frac{F(y)+G(\tilde{y})}{x_{4}\sqrt{z_{\pm}^{2}+2x_{3}^{2}x_{4}}}$

with arbitrary functions $F(.)$ and $G(.)$.

For any particular functions $F(.)$ and $F_{2}(.)$ we can find the
corresponding functions $N$ by direct integration of equations (105).
Moreover, for some rather extended classes of the arbitrary functions it is
possible to present the related functions $N$ in closed form like it was done
in equations (122)-(125). In particular we can set:

$\displaystyle
M=M_{m}=\frac{1}{x_{4}\sqrt{z_{\pm}^{2}+2x_{3}^{2}x_{4}}}(y^{m}+(-1)^{m+1}\tilde{y}^{n})$
(130)

where $n$ are natural numbers. Then the corresponding functions $N=N_{m}$ can
be found in the closed form for arbitrary $m$ by directs integration of
equations (105):

$\displaystyle N_{m}=2^{m-1}x_{3}^{2}M_{m-1}.$ (131)

We again have a countable set of exact solutions whose linear combinations
$\sum_{m}c_{m}M_{m}$ are solutions also. The corresponding potentials looks as
follows:

$\displaystyle
V=\frac{\sum_{m}\tilde{c}_{m}(y^{m}+(-1)^{m+1}\tilde{y}^{m})}{\sum_{m}c_{m}(y^{m}+(-1)^{m+1}\tilde{y}^{m})}.$
(132)

The obtained in this way solutions are represented in Items 11 and 12 of Table
1. Solutions corresponding to $n^{2}-1<0$ which are represented in Items 13
and 14 of the same table can be obtained in analogous way.

Thus we find all inequivalent PDM systems which, in addition to the cylindric
symmetry, admit at least one second order scalar integral of motion and so are
integrable. The inverse masses and potentials of these systems are defined up
to two arbitrary functions for pseudo scalar integrals of motion and up to
four arbitrary functions if integrals of motion are true scalars. For some
particular classes of the mentioned arbitrary functions the number of the
admitted scalar integrals of motion can be extended as it is indicated in
Tables 3-5.

Table 1. Inverse masses, potentials and scalar integrals of motion for
integrable systems

No | $f$ | $V$ | Integrals of motion
---|---|---|---
1 | $\frac{\tilde{r}^{2}}{c_{1}x_{3}+F(\tilde{r})}$ | $\frac{c_{2}x_{3}+G(\tilde{r})}{c_{1}x_{3}+F(\tilde{r})}$ | $P_{3}L_{3}+\frac{c_{1}}{2}(\varphi\cdot H)-\frac{c_{2}}{2}\varphi$
2 | $\frac{\tilde{r}^{2}}{F\left(\frac{r^{2}-1}{\tilde{r}}\right)+c_{1}\text{arctanh}\left(\frac{r^{2}+1}{2x_{3}}\right)}$ | $\frac{G\left(\frac{r^{2}-1}{\tilde{r}}\right)+c_{2}\text{arctanh}\left(\frac{r^{2}+1}{2x_{3}}\right)}{F\left(\frac{r^{2}-1}{\tilde{r}}\right)+c_{1}\text{arctanh}\left(\frac{r^{2}+1}{2x_{3}}\right)}$ | $\begin{array}[]{c}\\{L_{3},(K_{3}+P_{3})\\}\\\ -2c_{1}\left(\varphi\cdot H\right)+2c_{2}\varphi\end{array}$
3 | $\frac{\tilde{r}^{2}}{c_{1}\arctan\left(\frac{r^{2}-1}{2x_{3}}\right)+F\left(\frac{r^{2}+1}{\tilde{r}}\right)}$ | $\frac{c_{2}\arctan\left(\frac{r^{2}-1}{2x_{3}}\right)+G\left(\frac{r^{2}+1}{\tilde{r}}\right)}{c_{1}\arctan\left(\frac{r^{2}-1}{2x_{3}}\right)+F\left(\frac{r^{2}+1}{\tilde{r}}\right)}$ | $\begin{array}[]{c}\\{L_{3},(K_{3}-P_{3})\\}\\\ -2c_{1}\left(\varphi\cdot H\right)+2c_{2}\varphi\end{array}$
4 | $\frac{\tilde{r}^{2}}{F(\theta)-2c_{1}\ln(\tilde{r})}$ | $\frac{G(\theta)-2c_{2}\ln(\tilde{r})}{F(\theta)-2c_{1}\ln(\tilde{r})}$ | $\begin{array}[]{c}DL_{3}+c_{1}(\varphi\cdot H)-c_{2}\varphi\end{array}$
5 | $\frac{1}{F(\tilde{r})+G(x_{3})}$ | $\frac{\tilde{F}(\tilde{r})+\tilde{G}(x_{3})}{F(\tilde{r})+G(x_{3})}$ | $\begin{array}[]{c}P_{3}^{2}-\left(G(x_{3})\cdot H\right)+\tilde{G}(x_{3})\end{array}$
6 | $\frac{r^{2}}{F(\theta)+G(r)}$ | $\frac{\tilde{F}(\theta)+\tilde{G}(r)}{F(\theta)+G(r)}$ | $L_{1}^{2}+L_{2}^{2}-(F(\theta)\cdot H)+\tilde{F}(\theta)$
7 | $\begin{array}[]{c}\frac{r^{2}}{F(r+x_{3})+G(r-x_{3})}\end{array}$ | $\frac{\tilde{F}(r+x_{3})+\tilde{G}(r-x_{3})}{F(r+x_{3})+G(r-x_{3})}$ | $\begin{array}[]{c}\\{P_{3},D\\}-P_{n}x_{3}P_{n}\\\ -((F(r+x_{3})-G(r-x_{3}))\cdot H)\\\ +\tilde{F}(r+x_{3})-\tilde{G}(r-x_{3})\end{array}$
8 | $\frac{F_{\pm}}{F\left(\frac{r^{2}\pm 1}{{\tilde{r}}}\right)-G\left(\frac{r^{2}\mp 1}{x_{3}}\right)}$ | $\frac{\tilde{F}\left(\frac{r^{2}\pm 1}{\tilde{r}}\right)+\tilde{G}\left(\frac{r^{2}\mp 1}{x_{3}}\right)}{F\left(\frac{r^{2}\pm 1}{{\tilde{r}}}\right)-G\left(\frac{r^{2}\mp 1}{x_{3}}\right)}$ | $\begin{array}[]{c}(K_{3}\mp P_{3})^{2}+\tilde{G}\left(\frac{r^{2}\mp 1}{x_{3}}\right)\\\ +\left(G\left(\frac{r^{2}\mp 1}{x_{3}}\right)\cdot H\right)-\hat{\eta}\end{array}$
9 | $\begin{array}[]{c}\vspace{1mm}\frac{\sqrt{F_{\pm}}}{F(x_{\pm})+G(\tilde{x}_{\pm})}\end{array}$ | $\frac{\tilde{F}(x_{\pm})+\tilde{G}(\tilde{x}_{\pm})}{F_{1}(x_{\pm})+F_{2}(\tilde{x}_{\pm})}$ | $\begin{array}[]{c}\\{P_{3},(K_{3}\pm P_{3})\\}+P_{a}x_{3}^{2}P_{a}\\\ +\left((g(F,G)\cdot H\right)-\eta((\tilde{F},\tilde{G})\end{array}$
10 | $\begin{array}[]{c}\sum_{n}{c_{n}}F_{n},\ F_{n}=\frac{x_{\pm}^{n}+(-1)^{n+1}\tilde{x}_{\pm}^{n}}{F_{+}}\end{array}$ | $\begin{array}[]{c}\frac{\sum_{n}\tilde{c}_{n}F_{n}}{\sum_{n}c_{n}F_{n}}\end{array}$ | $\begin{array}[]{c}\\{P_{3},(K_{3}\pm P_{3})\\}+P_{a}x_{3}^{2}P_{a}\\\ -\left(x_{3}^{2}\sum_{n}c_{n}2^{n}F_{n-1}\cdot H\right)\\\ +x_{3}^{2}\sum_{n}\tilde{c}_{n}2^{n}F_{n-1}\end{array}$
11 | $\frac{\sqrt{z_{\pm}^{2}+2x_{3}^{2}\tilde{r}^{2}}}{F({y_{\pm}})+G({\tilde{y}_{\pm}})}$ | $\frac{\tilde{F}({y_{\pm}})+\tilde{G}({\tilde{y}_{\pm}})}{F({y_{\pm}})+G({\tilde{y}_{\pm}})}$ | $\begin{array}[]{c}K_{3}^{2}\pm P_{3}^{2}+2n\\{K_{3},P_{3}\\}+4P_{a}\tilde{r}^{2}P_{a}\\\ +\left(g(F,G)\cdot H\right)-\tilde{\eta}(\tilde{F},\tilde{G})-\hat{\eta},\\\ 1-n^{2}\leq 0\end{array}$
12 | $\begin{array}[]{c}\sum_{m}c_{m}\tilde{F}_{m},\ \tilde{F}_{m}=\frac{y_{\pm}^{m}+(-1)^{m+1}\tilde{y}_{\pm}^{m}}{z_{\pm}}\end{array}$ | $\begin{array}[]{c}\frac{\sum_{m}\tilde{c}_{m}\tilde{F}_{m}}{\sum_{m}c_{m}\tilde{F}_{m}}\end{array}$ | $\begin{array}[]{c}K_{3}^{2}\pm P_{3}^{2}+2n\\{K_{3},P_{3}\\}+4P_{a}\tilde{r}^{2}P_{a}\\\ -\left(x_{3}^{2}\sum_{m}c_{m}2^{m-1}\tilde{F}_{m-1}\cdot H\right)\\\ +x_{3}^{2}\sum_{m}\tilde{c}_{m}2^{m-1}\tilde{F}_{m-1}-\hat{\eta},\\\ 1-n^{2}\leq 0\end{array}$
13 | $\frac{\sqrt{\tilde{z}_{\pm}^{2}\mp 2x_{3}^{2}\tilde{r}^{2}}}{F({s_{\pm}})+G({\tilde{s}_{\pm}})}$ | $\frac{\tilde{F}({s_{\pm}})+\tilde{G}({\tilde{s}_{\pm}})}{F({s_{\pm}})+G({\tilde{s}_{\pm}})}$ | $\begin{array}[]{c}K_{3}^{2}\pm P_{3}^{2}+2n\\{K_{3},P_{3}\\}+4P_{a}\tilde{r}^{2}P_{a}\\\ +\left(g(F,G)\cdot H\right)-\tilde{\eta}(\tilde{F},\tilde{G})-\hat{\eta},\\\ n^{2}-1<0\end{array}$
14 | $\begin{array}[]{c}\sum_{m}c_{m}\Phi_{m},\ \Phi_{m}=\frac{s_{\pm}^{m}+(-1)^{m+1}\tilde{s}_{\pm}^{m}}{z_{\pm}}\end{array}$ | $\begin{array}[]{c}\frac{\sum_{m}\tilde{c}_{m}\Phi_{m}}{\sum_{m}c_{m}\Phi_{m}}\end{array}$ | $\begin{array}[]{c}K_{3}^{2}\pm P_{3}^{2}+2n\\{K_{3},P_{3}\\}+4P_{a}\tilde{r}^{2}P_{a}\\\ -\left(x_{3}^{2}\sum_{m}c_{m}2^{m-1}\Phi_{m-1}\cdot H\right)\\\ +x_{3}^{2}\sum_{m}\tilde{c}_{m}2^{m-1}\Phi_{m-1}-\hat{\eta},\\\ n^{2}-1<0\end{array}$

In Table 1 and the following Tables 2, 3, 4 $F(.),G(.),\tilde{F}(.)$ and
$\tilde{G}(.)$ are arbitrary functions, $g(F,G)$ and
$\eta(\tilde{F},\tilde{G})$ are solutions of equation (101), (105) with given
$F_{1}(x),F_{2}(\tilde{x})$ and $G_{1},G_{2}$, $c_{1},c_{2},..$ are arbitrary
real parameters. In addition we use the following notations which also will be
used in all tables:

$\displaystyle\begin{split}&F_{\pm}=(r^{2}\pm 1)^{2}\mp 4x_{3}^{2},\
\tilde{F}_{\pm}=(r^{2}\pm 1)^{2}\pm 4x_{3}^{2},\\\
&\hat{\eta}=3(r^{2}-5x_{3}^{2}),\ G_{\pm}=(\tilde{r}^{2}-2x_{3}^{2})(r^{2}\pm
1)^{2}+4x_{3}^{2}\tilde{r}^{2},\\\ &z_{\pm}=\frac{r^{4}\pm
1+2n(\tilde{r}^{2}-x_{3}^{2})}{2\sqrt{2(n^{2}\mp 1)}},\ \ -n^{2}\leq 1,\\\
&y_{\pm}=\frac{1}{\tilde{r}^{2}}\left(\sqrt{z_{\pm}^{2}+2x_{3}^{2}\tilde{r}^{2}}+z_{\pm}\right),\
\tilde{y}_{\pm}=\frac{1}{\tilde{r}^{2}}\left(\sqrt{z_{\pm}^{2}+2x_{3}^{2}\tilde{r}^{2}}-z_{\pm}\right),\\\
&\tilde{z}_{\pm}=\frac{r^{4}\pm 1+2n(\tilde{r}^{2}-x_{3}^{2})}{2\sqrt{2(1\mp
n^{2})}},\ n^{2}<1\\\
&s_{\pm}=\frac{1}{\tilde{r}^{2}}\left(\sqrt{\tilde{z}_{\pm}^{2}\mp
2x_{3}^{2}\tilde{r}^{2}}+\tilde{z}_{\pm}\right),\
\tilde{s}_{\pm}=\frac{1}{\tilde{r}^{2}}\left(\sqrt{\tilde{z}^{2}\mp
2x_{3}^{2}\tilde{r}^{2}}-\tilde{z}_{\pm}\right),\\\ &x_{\pm}=\sqrt{(r^{2}\pm
1)^{2}\mp x_{3}^{2}}+r^{2}\pm 1,\tilde{x}_{\pm}=\sqrt{(r^{2}\pm 1)^{2}\mp
x_{3}^{2}}-r^{2}\mp 1.\end{split}$ (133)

The mysterious term $\hat{\eta}$ appears as a result of the following
identity:

$\displaystyle K_{3}^{2}=p_{a}\mu_{8}^{ab}p_{b}+\hat{\eta}$

where $\mu_{9}^{ab}$ is the Killing tensor (41) with $\lambda_{9}^{3}\neq 0$ ,
and it is the cost paid for our desire to represent the integrals of motion
via the anticommutators of operators (21).

##  7 Vector integrals of motion

The systems admitting vector integrals of motion can be classified in the way
analogous to one used for the case of scalar integrals. First, using our
knowledge of the requested Killing tensors we will find the generic form of
second order integrals of motion which transform as a vectors under rotations
with respect to the third coordinate axis. Then, using the discrete
equivalence transformations we will select the subsets of such integrals with
fixed parities. Then, using the continuous equivalence transformations we will
specify all inequivalent linear combinations of the Killing tensors which
should be considered. Finally we will solve the determining equations
generated by such inequivalent Killing tensors.

###  7.1 Generic vector integrals of motion and their parity properties

Vector integrals of motion are generated by linear combinations of Killing
tensors (33)-(41) with nonzero parameters $\lambda^{\alpha}_{n}$ and
$\lambda^{3\alpha}_{n}$, where $\alpha=1,2,n=1,...,9$. The corresponding
integrals of motion (28) are linear combinations of operator $Q^{(0)}$
generated by $\mu^{ab}_{0}$, i.e.,

$\displaystyle Q^{(0)}=P_{a}g({\bf x})P_{a}$ (134)

and the following operators:

$\displaystyle\\{S_{12},S_{3\alpha}\\},\ \\{S_{12},S_{4\alpha}\\},\
\\{S_{12},S_{0\alpha}\\},$ (135) $\displaystyle\\{S_{3\alpha},S_{03}\\},\
\\{S_{3\alpha},S_{43}\\},\ \\{S_{3\alpha},S_{04}\\},\
\\{S_{0\alpha},S_{43}\\},\ \\{S_{4\alpha},S_{03}\\}.$ (136)

Notice that operators (135) are invariant with respect to the reflections
$x_{\alpha}\to-x_{\alpha}$ while operators (136) change their size under this
operation. Thus operators (135) and (136) are true and pseudo vectors
respectively.

The generic vector integral of motion includes a linear combination of
operators (135) and (136). Since Hamiltonians (42) are invariant with respect
to the space reflection, such combinations can include either operators (135)
or (136), but never both of them. Thus it is necessary to consider the true
vector and pseudo vector integrals of motion separately, while their linear
combinations can be rearrangered in accordance with the following comment.

Let us note that pseudo vectors (135) have the following property: whenever we
change the sign of $x_{1}$, the first components of these vectors are not
changed but the second components change their sign. Analogously, if we change
the sign of $x_{2}$: the second components are kept unchanged while the first
components change their signs also. In other words, if $V_{\alpha}$ is one of
pseudo vectors (135), then $\tilde{V}_{\alpha}=\varepsilon_{\alpha\nu}V_{\nu}$
(were $\varepsilon_{\alpha\nu}$ is the antisymmetric unit tensor, $\alpha$ and
$\nu$ take the values 1 or 2) transforms as a true vector under the space
reflections.

That is why the true vectors and pseudo vectors the generic forms of the
related operator $Q^{(0)}$ (134) are

$\displaystyle Q^{(0)}_{\alpha}=P_{a}x_{\alpha}g({\tilde{r},x_{3}})P_{a}$
(137)

and

$\displaystyle\tilde{Q}^{(0)}_{\alpha}=P_{a}\varepsilon_{\alpha\nu}x_{\nu}g({\tilde{r},x_{3}})P_{a}$
(138)

correspondingly, were $\varepsilon_{\alpha\nu}$ is the antisymmetric unit
tensor, $\alpha$ and $\nu$ take the values 1 or 2.

Like in Section 5, we can use the equivalence transformations from group
SO(1,2) whose generators are $S_{03},S_{04}$ and $S_{43}$ to simplify the
integrals o motion. Let us discuss these possible simplifications.

In accordance with (33)-(41) it is possible to specify five classes of vector
integrals of motion, whose coefficients are zero, first, second, third or four
order polynomials in $x_{a}$ which are fixed in the following equations (139),
(140), (141), (142) and (143) correspondingly:

$\displaystyle
Q^{(1)}_{\alpha}=P_{3}P_{\alpha}+Q^{(0)}_{\alpha}+x_{\alpha}\tilde{\eta}_{1},\
\mu^{ab}=\mu^{ab}_{1},$ (139)
$\displaystyle\begin{split}&Q^{(2)}_{\alpha}=\\{P_{\alpha},D\\}+Q^{(0)}_{\alpha}+x_{\alpha}\tilde{\eta}_{2},\
\mu^{ab}=\mu^{ab}_{2},\\\
&Q^{(3)}_{\alpha}=\varepsilon_{\alpha\nu}\\{P_{3},L_{\nu}\\}+Q^{(0)}_{\alpha}+x_{\alpha}\tilde{\eta}_{3},\
\mu^{ab}=\mu^{ab}_{2}+\mu^{ab}_{3},\end{split}$ (140)
$\displaystyle\begin{split}&Q^{(4)}_{\alpha}=\\{K_{3},P_{\alpha}\\}+Q^{(0)}_{\alpha}+x_{\alpha}\tilde{\eta}_{4},\
\mu^{ab}=\mu^{ab}_{4},\\\
&Q^{(5)}_{\alpha}=\varepsilon_{\alpha\nu}\\{D,L_{\nu}\\}+Q^{(0)}_{\alpha}+x_{\alpha}\tilde{\eta}_{5},\
\mu^{ab}=\mu^{ab}_{4}+\mu^{ab}_{6},\end{split}$ (141)
$\displaystyle\begin{split}&Q^{(6)}_{\alpha}=\\{D,K_{\alpha}\\}+Q^{(0)}_{\alpha}+x_{\alpha}\tilde{\eta}_{6},\
\mu^{ab}=\mu^{ab}_{7}\\\
&Q^{(7)}_{\alpha}=\varepsilon_{\alpha\nu}\\{K_{3},L_{\nu}\\}+Q^{(0)}_{\alpha}+x_{\alpha}\tilde{\eta}_{7},\
\mu^{ab}=\mu^{ab}_{7}+\mu^{ab}_{8},\end{split}$ (142) $\displaystyle
Q^{(8)}_{\alpha}=\\{K_{3},K_{\alpha}\\}+Q^{(0)}_{\alpha}+x_{\alpha}\tilde{\eta}_{8},\
\mu^{ab}=\mu^{ab}_{9}$ (143)

where $\tilde{\eta}_{n}=\tilde{\eta}_{n}(\tilde{r},x_{3})$ are an unknown
functions, $\mu^{ab}_{1},..,\mu^{ab}_{8}$ are tensors (33)-(41) with the
following nonzero coefficients:

$\displaystyle\lambda^{\alpha}_{n},\ \lambda^{3\alpha}_{m},\ \alpha=1,2,\
n=2,4,6,m=1,3,7,8,9.$ (144)

In accordance with (139)- (143), to generate vector integrals of motion we
choose the special Killing tensors whose nonzero parameters are fixed in
(144). In addition, functions $\eta$ and $g$ in (7) and (134) and functions
$N,K$ in (19) should be changed in the following manner:

$\displaystyle\eta\to x_{\alpha}\tilde{\eta}(\tilde{r},x_{3}),\ N\to
x_{\alpha}\tilde{N}(\tilde{r},x_{3}),\ K\to
x_{\alpha}\tilde{K}(\tilde{r},x_{3}).$ (145)

Let us note that operators $Q^{(8)},Q^{(7)}$ and $Q^{(6)}$ are equivalent to
$Q^{(1)},Q^{(3)}$ and $Q^{(2)}$ respectively up to transformation (24) and so
it is sufficient to consider only the potential integrals of motion
represented in (139)-(141). The generic vector integral of motion is a linear
combination of operators (139)-(143). However, such linear combination can be
a priori simplified.

Operators (139)-(143) have different invariance properties with respect to
transformation $x_{3}\to-x_{3}$. Half part of them including
$Q^{(1)},Q^{(4)},Q^{(5)}$ and $Q^{(8)}$ change their signs together with
$x_{3}$ while the remaining operators keep their form.

Let us start with operators $Q^{(1)},Q^{(4)},Q^{(5)},Q^{(8)}$ and consider
their linear combination:

$\displaystyle
Q=c_{1}Q^{(1)}+c_{4}Q^{(4)}+c_{5}Q^{(5)}+c_{8}Q^{(8)}+Q^{(0)}_{\alpha}+x_{\alpha}\tilde{\eta}.$
(146)

Let the coefficients $c_{1}$ and $c_{8}$ are nontrivial, then up to the
dilatation transformation $c_{1}=\pm c_{8}$. Setting $c_{8}=2$ we come to the
following version of operator (146):

$\displaystyle Q=\\{K_{3},K_{\alpha}\\}\pm
2P_{3}P_{\alpha}+c_{4}\\{K_{3},P_{\alpha}\\}+c_{5}\varepsilon_{\alpha\nu}DL_{\nu}+Q^{(0)}_{\alpha}+x_{\alpha}\tilde{\eta}.$
(147)

If only one of the coefficients $c_{1},c_{8}$ is trivial and $c_{4}$ is
nontrivial then up to transformation (24) and dilatation we can reduce (146)
to the following form:

$\displaystyle Q=\\{(K_{3}\pm
P_{3}),P_{\alpha}\\}+c_{5}\varepsilon_{\alpha\nu}DL_{\nu}+Q^{(0)}_{\alpha}+x_{\alpha}\tilde{\eta},$
(148)

while the remaining inequivalent versions of operator (146) are

$\displaystyle
Q=P_{3}P_{\alpha}+c_{5}\varepsilon_{\alpha\nu}DL_{\nu}+Q^{(0)}_{\alpha}+x_{\alpha}\tilde{\eta}$
(149)

and

$\displaystyle
c_{4}\\{K_{3},P_{a}\\}+c_{5}\varepsilon_{\alpha\nu}DL_{\nu}+Q^{(0)}_{\alpha}+x_{\alpha}\tilde{\eta}$
(150)

Analogous speculations with the linear combinations of operators
$Q_{2},Q_{3},Q_{6}$ and $Q_{7}$ make it possible to specify the following
their inequivalent versions:

$\displaystyle\begin{split}&Q=c_{1}\\{D,(K_{\alpha}\pm
P_{\alpha})\\}+c_{2}\varepsilon_{\alpha\nu}\\{L_{\nu},(K_{3}\pm
P_{3})\\}+Q^{(0)}_{\alpha}+c_{3}\varepsilon_{\alpha\nu}\\{L_{\nu},P_{3}\\}+x_{\alpha}\tilde{\eta},\\\
&Q=c_{1}\\{D,P_{\alpha}\\}+c_{2}\varepsilon_{\alpha\nu}\\{L_{\nu},(K_{3}\pm
P_{3})\\}+Q^{(0)}_{\alpha}+x_{\alpha}\tilde{\eta},\\\
&Q=c_{1}\\{D,(K_{\alpha}\pm
P_{\alpha})\\}+c_{2}\varepsilon_{\alpha\nu}\\{L_{\nu},P_{3}\\}+Q^{(0)}_{\alpha}+x_{\alpha}\tilde{\eta},\\\
&Q=c_{1}\\{D,P_{\alpha}\\}+c_{2}\varepsilon_{\alpha\nu}\\{L_{\nu},P_{\alpha}\\}+Q^{(0)}_{\alpha}+x_{\alpha}\tilde{\eta}.\end{split}$
(151)

Concerning the linear combinations of operators (139) - (143) with different
parities with respect to the reflection of $x_{3}$ we can note that such
combinations of operator $Q^{(1)}_{\alpha}$ with $Q^{(2)}_{\alpha}$ or
$Q^{(3)}_{\alpha}$ are equivalent to $Q^{(2)}_{\alpha}$ or $Q^{(3)}_{\alpha}$
up to the shift transformations. The same is true for the bilinear
combinations of $Q^{(2)}_{\alpha}$ with $Q^{(5)}_{\alpha}$ and
$Q^{(4)}_{\alpha}$ with $Q^{(7)}_{\alpha}$: they can be reduced to
$Q^{(5)}_{\alpha}$ and $Q^{(7)}_{\alpha}$ respectively.

We will not discuss the remaining linear combinations since the computing
experiments show that they do not correspond to nontrivial integrals of
motion.

Thus to classify the PDM systems with cylindric symmetry which admit vector
integrals of motion we are supposed to solve the determining equations (14),
(16) which correspond to the symmetries represented in equations (139) - (141)
and (147) - (151).

###  7.2 Selected calculations

The next step is to solve the the determining equations indicated in the
above.

For the symmetry specified in (139) the determining equations (14) have
especially simple form. Namely, substituting the Killing tensor
$\mu^{3\alpha}_{1}$ for $\alpha=1$ into (14) we come to the following system:

$\displaystyle\begin{split}&\partial_{3}M(\tilde{r},x_{3})=\partial_{1}(x_{1}\tilde{N}(\tilde{r},x_{3})),\\\
&\partial_{2}\tilde{N}(\tilde{r},x_{3})=0,\
\partial_{1}M(\tilde{r},x_{3})=\partial_{3}(x_{1}\tilde{N}(\tilde{r},x_{3}))\end{split}$
(152)

which is easy solvable. Its generic solution is:

$\displaystyle M=c_{1}r^{2}+c_{2}x_{3}+c_{3},\ \tilde{N}=-c_{1}x_{3}-c_{2}.$
(153)

In contrast with the case of scalar integrals of motion (refer to Table 1) the
obtained solution do not include arbitrary functions but only arbitrary
parameters $c_{1},c_{2}$ and $c_{3}$. Moreover, parameter $c_{2}$ can be
reduced to zero via the shift of $x_{3}$ provided parameter $c_{1}$ is
nontrivial. Thus taking into account that at least one of nonzero parameters
in parameter in (153) can be reduced to the unity by simultaneous scaling of
all independent variables $x_{1},x_{2}$ and $x_{3}$, without loss of
generality we can rewrite the latter equation in the following form:

$\displaystyle M=\cos(\lambda)r^{2}+\sin(\lambda)x_{3}+c_{3},\
\tilde{N}=-\cos(\lambda)x_{3}-\sin(\lambda)$ (154)

where parameter $\lambda$ takes two values: $\lambda=0$ or
$\lambda=\frac{\pi}{2}$

The reason of the reduced freedom in the arbitrary elements $M$ lies in the
fact that the considered systems by definition admit as minimum two integrals
of motion which are the components of the bivector (139). We have found the
first component, i.e, choose $\alpha=1$ in (139) and (145) but the second one
can be obtained simple by the changes $P_{1}\to P_{2},x_{1}\to x_{2}$.

Moreover, it happens that our system admits two more integrals of motion and
is maximally superintegrable. Its integrals of motion together with the
admissible potential are represented in Item 4 of Table 5.

Consider now integral of motion $Q^{(2)}_{a}$ presented in (140) where
$Q^{(0)}_{\alpha}=-P_{a}x_{\alpha}P_{a}$ and $\alpha=1$. The related
determining equations (14) have the following form:

$\displaystyle\begin{split}&2\tilde{r}\partial_{4}M(\tilde{r}^{2},x_{3})+x_{3}\partial_{3}M(\tilde{r}^{2},x_{3})+M(\tilde{r}^{2},x_{3})+N(\tilde{r}^{2},x_{3})=0,\\\
&2x_{3}\partial_{4}M(\tilde{r}^{2},x_{3})-\partial_{3}M(\tilde{r}^{2},x_{3})+\partial_{3}N(\tilde{r}^{2},x_{3})=0,\
\partial_{4}N(\tilde{r}^{2},x_{3})=0\end{split}$ (155)

and are easy solvable also. Their generic solution is:

$\displaystyle M=c_{1}+\frac{c_{2}}{x_{3}^{2}},+\frac{c_{3}}{\tilde{r}},\
N=-c_{1}+\frac{c_{2}}{x_{3}^{2}}.$ (156)

The corresponding inverse mass and potential are presented in Item 5 of Table
4 where one more integral of motion is indicated also. Thus the considered
system is superintegrable.

The next (and the last) example we consider is the system admitting the
symmetry presented in (150), i.e.,
$c_{1}Q^{(4)}_{\alpha}+c_{2}Q^{(5)}_{\alpha}.$ Setting $\alpha=1$ we come to
the following determining equations (14):

$\displaystyle(c_{2}-c_{1})(x_{4}-x_{3}^{2})\partial_{3}M+(2c_{1}-c_{2})\partial_{4}M+2(2c_{1}-c_{2})x_{3}M-2x_{4}\partial_{4}N-N=0,$
$\displaystyle
2(c_{2}-2c_{1})x_{3}x_{4}\partial_{4}M+((c2-c1)x_{3}^{2}+c_{1}x_{4})\partial_{3}M+2(c_{2}-2c_{1}))x_{3}\partial_{3}M+N=0,$
$\displaystyle(c_{1}-c_{2})(x_{3}^{2}-x_{4})-\frac{1}{2}\partial_{3}Nc_{2}x_{3}\partial_{3}M+c_{2}M=0.$

Excluding $N$ we come to the following compatibility conditions of the system
presented above:

$\displaystyle(c_{2}-2c_{1})x_{3}x_{4}\partial_{34}M+\frac{1}{2}(c_{2}-c_{1})(x_{3}^{2}+c_{1}x_{4})\partial_{34}M+\frac{3}{2}(c_{2}-2c_{1})\partial_{4}M+\frac{1}{4}(c_{2}+2c_{1})\partial_{3}M=0,$
$\displaystyle(c_{1}-c_{2})(x_{4}-x_{3}^{2})\partial_{44}M+\frac{3}{2}\left(c_{2}-\frac{2}{3}\right)x_{3}\partial_{34}M+\left(\frac{5}{2}c_{2}-2c_{1}\right)\partial_{4}M-\frac{1}{4}\partial_{33}M=0.$

This system has two special and one regular solutions. Namely, for
$c_{1}=2c_{2}$ and $c_{1}=c_{2}$ we obtain:

$\displaystyle M=c_{3}r$ (157)

and

$\displaystyle M=\frac{c_{3}x_{3}}{r^{2}}+F(r)$ (158)

respectively, where $F(r)$ is an arbitrary function. If parameters $c_{1}$ and
$c_{2}$ do not satisfy the conditions presented above, the solution is

$\displaystyle M=\frac{c_{3}}{r}.$ (159)

The system with whose mass is given in (158) is maximally superintegrable, see
Item 1 of Table 5. The solutions (157) and (159) correspond to the
rotationally invariant systems which admit as minimum the three parametric Lie
group. Such systems are completely classified in paper [9] and we will not
discuss them here.

Solving step by step the determining equations corresponding to the remaining
symmetries (141), (142), (143), (147) - (151) we find all inequivalent PDM
systems with cylindric symmetry, which admit vector integrals of motion. The
obtained results are presented in Tables 2-5.

##  8 Tensor integrals of motion

The last class of symmetries which is supposed to be considered are the tensor
integrals of motion. The corresponding Killing tensors are given by equations
(33), (35), (38), (40) and (41) where where the indices $a,b$ of the nonzero
parameters $\lambda_{k}^{ab}$ independently take the values $1$ and $2$. In
other words, we have to consider exactly five linearly independent tensors and
their linear combinations. The corresponding integrals of motion include
operators (134) and the following bilinear combinations:

$\displaystyle P_{\alpha}P_{\nu},\ \\{K_{\alpha}K_{\nu}\\},\
\\{P_{\alpha}K_{\nu}\\}$ (160)

and

$\displaystyle\\{P_{\alpha},L_{\nu}\\},\ \\{K_{\alpha}L_{\nu}\\}$ (161)

were $\alpha$ and $\nu$ take the values 1 and 2.

Notice that operators (160) are true tensors while operators (161) change
their signs under the reflection and so are pseudo tensors. It means that we
have to consider linear combinations of symmetries (160) and (161) separately.

One more note is that the traceless symmetric tensor $Q^{\alpha\nu}$ in two
dimensions has exactly two linearly independent components, i.e., $Q^{12}$ and
$Q^{11}-Q^{22}$. We will work with the first components, i.e., $Q^{12}$ and
will omit the top index $12$.

###  8.1 True tensor integrals of motion

Up to equivalence transformation (22) we can specify the following
inequivalent linear combinations of true scalars (157) and the universal block
(134):

$\displaystyle Q_{1}=Q^{(0)}_{1}+P_{1}P_{2}+\eta_{1},$ (162) $\displaystyle
Q_{2}=Q^{(0)}_{2}\\{K_{1},K_{2}\\}\pm 2P_{1}P_{2}+\eta_{2},$ (163)
$\displaystyle Q_{3}=Q^{(0)}+\\{P_{1},P_{2}\pm K_{2}\\}+\eta_{3},$ (164)
$\displaystyle Q_{4}=Q^{(0)}+\\{P_{1}\pm K_{1},P_{2}\pm K_{2}\\}+\eta_{4}.$
(165)

The determining equations (14) for symmetry (162) are simple and have the
following form:

$\displaystyle\partial_{2}M+\partial_{1}N=0,\ \partial_{1}M+\partial_{2}N=0,\
\partial_{3}N=0.$ (166)

Taking into account that the generic form of functions $M$ and $N$ is given by
the following equation

$\displaystyle M=M(\tilde{r},x_{3}),\ N=x_{1}x_{2}\tilde{N}(\tilde{r},x_{3})$

the system (166) is easy integrated and solved by the following functions:

$\displaystyle M=c_{1}r^{2}+F(x_{3}),\ N=c_{1}x_{1}x_{2}.$

The corresponding PDM system appears to be superintegrable, see Item 2 of
Table 4.

Let us consider one more integral of motion which is specified in (165).
Denoting $N=x_{1}x_{2}\tilde{N}(\tilde{r}^{2},x_{3})$ we come to the following
corresponding equations (14):

$\displaystyle\partial_{4}((x_{4}-x_{3}^{2}\pm
1)^{2}M)+x_{3}(x_{4}-x_{3}^{2}\pm
1)\partial_{3}M+\frac{1}{2}\partial_{4}(x_{4}\tilde{N})=0,$ $\displaystyle
x_{3}(x_{4}-x_{3}^{2}-1)\partial_{4}M+\partial_{3}(x_{3}^{2}M)+\frac{1}{8}\partial_{3}\tilde{N}$

which are perfectly solved by the following functions:

$\displaystyle M=\frac{(r^{2}\pm 1)F\left(\frac{r^{2}\mp
1}{x_{3}}\right)-c_{1}r^{2}}{(r^{4}-1)^{2}},\
\tilde{N}=\frac{c_{1}x_{1}x_{2}}{(r^{2}\pm 1)^{2}}.$

The obtained results are represented in Item 6 of Table 2.

We will not discuss the determining equations for symmetries (163) and (164)
but mention that both of them can be solved exactly also. However, for
symmetry (163) the related solutions are functions of $x_{3}$ only, but for
symmetry (164) we obtain the mass function which depends only on $r$. Such
systems are out of the scop of the present paper since they were discussed in
papers [35] and [39].

Table 2. Inverse masses, potentials, vector and tensor integrals of motion for
integrable PDM systems

No | $f$ | $V$ | Integrals of motion
---|---|---|---
1 | $\frac{x_{3}^{2}}{c_{1}+c_{2}(r^{2}\pm 1)F_{\pm}^{\frac{1}{2}}+c_{3}x_{3}^{2}(r^{2}\mp 1)^{-2}}$ | $\frac{c_{4}+c_{5}(r^{2}\pm 1)F_{\pm}^{\frac{1}{2}}+c_{6}x_{3}^{2}(r^{2}\mp 1)^{-2}}{c_{1}+c_{2}(r^{2}\pm 1)F_{\pm}^{\frac{1}{2}}+c_{3}x_{3}^{2}(r^{2}\mp 1)^{-2}}$ | $\begin{array}[]{c}\\{(K_{3}\pm P_{3}),L_{\alpha}\\}+3\varepsilon_{\alpha\nu}x_{\nu}\\\ +\frac{2\varepsilon_{\alpha\nu}x_{\nu}(c_{4}(r^{2}\pm 1)+c_{5}(r^{2}\mp 1)^{2})}{x_{3}^{2}}\\\ -\left(\frac{2\varepsilon_{\alpha\nu}x_{\nu}(c_{1}(r^{2}\pm 1)+c_{2}(r^{2}\mp 1)^{2})}{x_{3}^{2}}\cdot H\right)\end{array}$
2 | $\frac{x_{3}^{2}\tilde{F}_{\pm}}{c_{2}(r^{2}\mp 1)+\tilde{F}_{\pm}(c_{3}x_{3}^{2}(r^{2}\pm 1)^{-2}+c_{1})}$ | $\frac{c_{5}(r^{2}\mp 1)+\tilde{F}_{\pm}(c_{6}x_{3}^{2}(r^{2}\pm 1)^{-2}+c_{4})}{c_{2}(r^{2}\mp 1)+\tilde{F}_{\pm}(c_{3}x_{3}^{2}(r^{2}\pm 1)^{-2}+c_{1})}$ | $\begin{array}[]{c}\\{(K_{\alpha}\pm P_{\alpha}),D\\}-15x_{\alpha}\\\ \mp\varepsilon_{\alpha\nu}\\{L_{3},(K_{\nu}\mp P_{\nu})\\}\pm 3x_{\alpha}\\\ -2x_{\alpha}\left(c_{5}\frac{r^{2}\mp 1}{(r^{2}\pm 1)^{2}}\mp\frac{2c_{6}}{\tilde{F}_{\pm}}\right)\\\ +2\left(x_{\alpha}\left(c_{2}\frac{r^{2}\mp 1}{(r^{2}\pm 1)^{2}}\mp\frac{2c_{3}}{\tilde{F}_{\pm}}\right)\cdot H\right)\\\ \end{array}$
3 | $\frac{x_{3}^{6}}{c_{1}(x_{3}^{2}+4\tilde{r}^{2})+c_{2}x_{3}^{4}+c_{3}x_{3}^{6}}$ | $\frac{(c_{4}x_{3}^{2}+4\tilde{r}^{2})+c_{5}x_{3}^{4}+c_{6}x_{3}^{6}}{c_{1}(x_{3}^{2}+4\tilde{r}^{2})+c_{2}x_{3}^{4}+c_{3}x_{3}^{6}}$ | $\begin{array}[]{c}\\{P_{\alpha},D\\}+\varepsilon_{\alpha\nu}\\{L_{3},P_{\nu}\\}\\\ +2\left(\frac{x_{\alpha}(c_{1}-c_{3}x_{3}^{4})}{x_{3}^{4}}\cdot H\right)\\\ -2\frac{x_{\alpha}(c_{4}-c_{6}x_{3}^{4})}{x_{3}^{4}}\end{array}$
4 | $\frac{r(r^{2}\pm 1)^{2}x_{3}^{2}}{c_{1}r(r^{2}\pm 1)^{2}+c_{2}rx_{3}^{2}+c_{3}(r^{2}\mp 1)x_{3}^{2}}$ | $\frac{c_{4}r(r^{2}\pm 1)^{2}+c_{5}rx_{3}^{2}+c_{6}(r^{2}\mp 1)x_{3}^{2}}{c_{1}r(r^{2}\pm 1)^{2}+c_{2}rx_{3}^{2}+c_{3}(r^{2}\mp 1)x_{3}^{2}}$ | $\begin{array}[]{c}\\{D,(K_{\alpha}\pm P_{\alpha})\\}-\frac{2x_{\alpha}c_{5}(r^{2}\mp 1)}{(r^{2}\pm 1)^{2}}\\\ -\frac{x_{\alpha}c_{6}((r^{2}\mp 1)\mp 4r^{2})}{r(r^{2}\pm 1)^{2}}-15x_{\alpha}\\\ +\left(\frac{2x_{\alpha}c_{2}((r^{2}\mp 1)}{(r^{2}\pm 1)^{2}}\cdot H\right)\\\ +\left(\frac{x_{\alpha}c_{3}((r^{2}\mp 1)^{2}+4r^{2})}{r(r^{2}\pm 1)^{2}}\cdot H\right)\end{array}$
5 | $\frac{x_{3}^{2}}{x_{3}^{2}F(r)+c_{1}}$ | $\frac{x_{3}^{2}G(r)+c_{2}}{x_{3}^{2}F(r)+c_{1}}$ | $\begin{array}[]{c}\\{L_{1},L_{2}\\}+\left(\frac{c_{1}x_{1}x_{2}}{x_{3}^{2}}\cdot H\right)-\frac{c_{2}x_{1}x_{2}}{x_{3}^{2}},\\\ L_{1}^{2}-L_{2}^{2}+\left(\frac{c_{1}(x_{1}^{2}-x_{2}^{2})}{2x_{3}^{2}}\cdot h\right)-\frac{c_{2}(x_{1}^{2}-x_{2}^{2})}{2x_{3}^{2}}\end{array}$
6 | $\frac{(r^{4}-1)^{2}}{(r^{2}\pm 1)^{2}F\left(\frac{r^{2}\mp 1}{x_{3}}\right)-c_{1}r^{2}}$ | $\frac{(r^{2}\pm 1)^{2}G\left(\frac{r^{2}\mp 1}{x_{3}}\right)+{c_{2}}r^{2}}{(r^{2}\pm 1)^{2}F\left(\frac{r^{2}\mp 1}{x_{3}}\right)-c_{1}r^{2}}$ | $\begin{array}[]{c}(K_{1}\pm P_{1})^{2}-(K_{2}\pm P_{2})^{2}\\\ +15(x_{1}^{2}-x_{2}^{2})\\\ +c_{1}\left(\frac{x_{1}^{2}-x_{2}^{2}}{(r^{2}\pm 1)^{2}}\cdot H\right)+{c_{2}}\frac{x_{1}^{2}-x_{2}^{2}}{(r^{2}\pm 1)^{2}},\\\ \\{(K_{1}\pm P_{1}),(K_{2}\pm P_{2})\\}\\\ +{c_{2}}\frac{x_{1}x_{2}}{(r^{2}\pm 1)^{2}}+15x_{1}x_{2}\end{array}$
7 | $\frac{(r^{2}\pm 1)^{2}F_{\pm}^{2}}{c_{1}G^{\pm}+2c_{2}x_{3}^{2}(r^{2}\pm 1)+c_{3}F_{\pm}^{2}}$ | $\frac{c_{4}G^{\pm}+2c_{5}x_{3}^{2}(r^{2}\pm 1)+2c_{6}F_{\pm}^{2}}{c_{1}G^{\pm}+2c_{2}x_{3}^{2}(r^{2}\pm 1)+c_{3}F_{\pm}^{2}}$ | $\begin{array}[]{c}\\{(K_{1}\pm P_{1}),L_{1}\\}\\\ -\\{(K_{2}\pm P_{2}),L_{2}\\}\\\ +2\left(x_{1}x_{2}(\frac{c_{1}x_{3}(r^{2}\pm 1)}{F_{\pm}^{2}}+c_{2})\cdot H\right)\\\ -2x_{1}x_{2}(\frac{c_{4}x_{3}(r^{2}\pm 1)}{F_{\pm}^{2}}+c_{5}),\\\ \\{(K_{1}\pm P_{1}),L_{2}\\}\\\ +\\{(K_{2}\pm P_{2}),L_{1}\\}\\\ +\left((x_{1}^{2}-x_{2}^{2})(\frac{c_{1}x_{3}(r^{2}\pm 1)}{F_{\pm}^{2}}+c_{2})\cdot H\right)\\\ -(x_{1}^{2}-x_{2}^{2})(\frac{c_{4}x_{3}(r^{2}\pm 1)}{F_{\pm}^{2}}+c_{5})\end{array}$

###  8.2 Pseudo tensor integrals of motion

Such integrals should include linear combinations of the terms presented in
(161). Up to the equivalence we can restrict ourselves to two types of such
linear combinations:

$\displaystyle Q_{6}=Q^{(0)}+\\{P_{1},L_{2}\\}+\\{P_{2},L_{1}\\}+\eta_{6}$
(167)

and

$\displaystyle Q_{7}=Q^{(0)}+\\{(P_{1}\pm K_{1}),L_{2}\\}+\\{(P_{2}\pm
K_{2}),L_{1}\\}+\eta_{7}.$ (168)

Let us denote $N=(x_{1}^{2}-x_{2}^{2})\tilde{N}(r,x_{3})$ then symmetry (167)
generates the following version of the determining equations (14):

$\displaystyle
2(x_{1}^{2}-x_{2}^{2})\partial_{4}\tilde{N}+4x_{3}\partial_{4}M-\partial_{3}M+2\tilde{N}=0,$
$\displaystyle 2\partial_{4}M-\partial_{3}N=0.$

which are solved by the following functions:

$\displaystyle
M=c_{1}+c_{2}(\tilde{r}^{2}+4x_{3}^{2}+2c_{3}x_{3}),\tilde{N}=-2c_{2}x_{3}-c_{3}.$

The corresponding PDM system appears to be maximally superintegrable, see Item
2 of Table 5.

The next (and the last in this section) symmetry which we consider is
specified in (168). The corresponding determining equations (14) are reduced
to the following system:

$\displaystyle
x_{3}(x_{4}-x_{3}^{2}-1)\partial_{4}M+\frac{1}{4}(1+3x_{3}^{2}-x_{4})\partial_{3}M-\frac{1}{2}(x_{4}\tilde{N})+x_{3}M=0,$
$\displaystyle(x_{4}+3x_{3}^{2}+1)\partial_{3}M-4x_{3}(x_{3}^{2}+1)\partial_{4}M+4x_{3}M-2\tilde{N},$
$\displaystyle(x_{4}-3x_{3}^{2}-1)\partial_{4}M+2\partial_{3}(x_{3}M)+\frac{1}{2}\tilde{N}=0$

which is solved by the following functions:

$\displaystyle M=c_{1}\frac{(\tilde{r}^{2}-2x_{3}^{2})(r^{2}\pm
1)^{2}+4x_{3}^{2}\tilde{r}^{2}}{(r^{2}\pm 1)^{2}((r^{2}\pm 1)^{2}\pm
4x_{3}^{2})}+2c_{2}\frac{x_{3}^{2}}{(r^{2}\pm 1)((r^{2}\pm 1)^{2}\pm
4x_{3}^{2})}+\frac{c_{3}}{(r^{2}\pm 1)^{2}},$
$\displaystyle\tilde{N}=c_{1}\frac{x_{3}}{((r^{2}\pm 1)^{2}\pm
4x_{3}^{2})}+c_{2}.$

The related PDM system is integrable, see Item 7 of Table 2.

##  9 Superintegrable systems

Thus we have specified all nonequivalent cylindrically invariant PDM systems
which admit a second order integral of motion. By definition such systems
admit two integrals of motion one of which is the generator of rotations
around the third coordinate axis, and so they are integrable. All the systems
admitting two (but no more) integrals of motion are presented in Tables 1 and
2.

However, some of the found systems automatically admit more than two integrals
of motion. In addition, for some particular versions of the arbitrary
elements, i.e., arbitrary functions and integration constants we also can find
additional integrals of motion. The related systems are superintegrable or
even maximally superintegrable.

The next (and final) step of our analysis is the specifications of just these
systems. It can be done in two ways.

The first way presupposes solving of the functional equations. Let you have
two systems any of which admit one second order integral of motion, and the
related mass functions are $M^{(1)}$ and $M^{(2)}$. Then we have to equate to
zero their linear combination:

$\displaystyle k_{1}M^{(1)}+k_{2}M^{(2)}=0.$ (169)

and this condition generates a functional equation for the arbitrary elements
present in $M^{(1)}$ and $M^{(2)}$.

The second way is to solve the extended systems of the determining equations
(14) corresponding to the inequivalent pairs of integrals of motion. In all
cases when such extended systems have nontrivial solutions we can obtain the
mass functions of the related superintegrable systems.

The same speculations with the inequivalent triplets of integrals of motion
and the related mass functions make it possible to specify the inequivalent
maximally superintegrable sytems. We will not present the complete set of the
corresponding calculations which do not include new elements in comparison
with the given above, but restrict ourselves to considering few examples.

###  9.1 Superintegrable systems admitting scalar integrals of motion

The systems admitting scalar integrals of motion are defined up to arbitrary
functions, and to specify such of them which admit as minimum two such
integrals it is preferable to use the second approach, i.e. to search for such
cases when arbitrary elements $M$ and $N$ satisfy more than one system of the
determining equations specified in (76), (78), (79), (84) and (85)-(108).

Let the PDM system includes arbitrary elements satisfying determining
equations (85), see Item 5 of Table 1. The corresponding position dependent
mass include two arbitrary functions. Let us search for such the arbitrary
functions that the related PDM system admits one more second order integral of
motion. It is the case when functions $M$ and $N$ satisfies one more system of
the determining equations mentioned in the above.

Solving step by step all pairs of equations including (85) and one of of
equations (76), (78), (79), (84) and (88)-(108) we obtain the following
results. Equation (85) is incompatible with equations (78), i.e, the
corresponding systems of equations have no solutions. The remaining pairs of
equations are consistent, but only three of them, namely, (85) and (88), (85)
and (91), (85) and (83) have good solutions, presented in Items 1, 3 and 3 of
Table 3. Solutions for the other pairs of equations are invariant with respect
to two parametric Lie groups and so can be ignored as being classified in our
previous papers.

Considering the doublets of equations which include (88) we recognize two more
inequivalent systems with good solutions. Those are the doublets including
systems (95) or (108). One of the mentioned solutions is represented in Item 4
of Table 3. The other solutions are too cumbersome to be placed in the table,
and we fix them in equations (171).

Table 3. Inverse masses, potentials and scalar integrals of motion for
superintegrable systems

No | $f$ | $V$ | Integrals of motion
---|---|---|---
1 | $\frac{\tilde{r}^{2}}{\tilde{r}^{2}(c_{1}+2c_{2}x_{3}+c_{3}(\tilde{r}^{2}+4x_{3}^{2}))+c_{4}}$ | $\frac{\tilde{r}^{2}(c_{5}+2c_{6}x_{3}+c_{7}(\tilde{r}^{2}+4x_{3}^{2}))+c_{8}}{\tilde{r}^{2}(c_{1}+2c_{2}x_{3}+c_{3}(\tilde{r}^{2}+4x_{3}^{2}))+c_{4}}$ | $\begin{array}[]{c}P_{3}^{2}+2c_{6}x_{3}+4c_{7}x_{3}^{2}+c_{4}\frac{x_{3}}{\tilde{r}^{2}}\\\ -\left((2c_{2}x_{3}+4c_{3}x_{3}^{2})\cdot H\right),\\\ \\{P_{3},D\\}-P_{n}x_{3}P_{n}\\\ +\left((-c_{1}x_{3}-c_{2}(2x_{3}^{2}+\tilde{r}^{2})\right.\\\ \left.-c_{3}x_{3}(4x_{3}^{2}+3\tilde{r}^{2}))\cdot H\right)\\\ +c_{5}x_{3}+c_{6}(2x_{3}^{2}+\tilde{r}^{2})\\\ -c_{8}\frac{x_{3}}{\tilde{r}^{2}}+c_{7}x_{3}(4x_{3}^{2}+3\tilde{r}^{2})\end{array}$
2 | $\frac{r\tilde{r}^{2}}{c_{1}r\tilde{r}^{2}+c_{2}r+c_{3}x_{3}+c_{7}\tilde{r}^{2}}$ | $\frac{c_{4}r\tilde{r}^{2}+c_{5}r+c_{6}x_{3}+c_{8}\tilde{r}^{2}}{c_{1}r\tilde{r}^{2}+c_{2}r+c_{3}x_{3}+c_{7}\tilde{r}^{2}}$ | $\begin{array}[]{c}\\{P_{3},D\\}-P_{n}x_{3}P_{n}\\\ +c_{2}\left(\frac{c_{2}x_{3}+c_{3}r}{\tilde{r}^{2}}\cdot H\right)\\\ -\frac{c_{5}x_{3}+c_{6}r}{\tilde{r}^{2}},\\\ L_{1}^{2}+L_{2}^{2}\\\ -(\frac{r^{2}}{\tilde{r}^{2}}(c_{2}+c_{3}\frac{x_{3}}{r})\cdot H)\\\ +(c_{5}+c_{6}\frac{x_{3}}{r})\frac{r^{2}}{\tilde{r}^{2}}\end{array}$
3 | $\frac{\tilde{r}^{2}x_{3}^{2}}{c_{1}x_{3}^{2}+{c_{2}}\tilde{r}^{2}+{c_{3}}x_{3}^{2}\tilde{r}^{2}+{c_{4}}r^{2}\tilde{r}^{2}x_{3}^{2}}$ | $\frac{{c_{5}}x_{3}^{2}+{c_{6}}\tilde{r}^{2}+{c_{7}}x_{3}^{2}\tilde{r}^{2}+{c_{8}}r^{2}\tilde{r}^{2}x_{3}^{2}}{c_{1}x_{3}^{2}+{c_{2}}\tilde{r}^{2}+{c_{3}}x_{3}^{2}\tilde{r}^{2}+{c_{4}}r^{2}\tilde{r}^{2}x_{3}^{2}}$ | $\begin{array}[]{c}\\{P_{3},K_{3}\\}+\frac{2{c_{6}}\tilde{r}^{2}}{x_{3}^{2}}-2\omega x_{3}^{2}\\\ +2\left(({c_{3}}x_{3}^{2}-{c_{2}}\frac{\tilde{r}^{2}}{x_{3}^{2}})\cdot H\right),\\\ P_{3}^{2}-\left((\frac{c_{2}}{x_{3}^{2}}+{c_{4}}x_{3}^{2})\cdot H\right)\\\ +\frac{c_{6}}{x_{3}^{2}}+{c_{8}}x_{3}^{2}\end{array}$
4 | $\frac{\tilde{r}^{2}x_{3}^{2}}{c_{1}\tilde{r}^{2}+x_{3}^{2}(c_{3}-2{c_{2}}\ln(\tilde{r}))}$ | $\frac{x_{3}^{2}(2c_{5}\ln(\tilde{r})+c_{6})+c_{4}\tilde{r}^{2}}{c_{1}\tilde{r}^{2}+x_{3}^{2}(c_{3}-2{c_{2}}\ln(\tilde{r}))}$ | $\begin{array}[]{c}P_{3}^{2}-\left((\frac{c_{1}}{x_{3}^{2}}\cdot H\right)+\frac{c_{4}}{x_{3}^{2}},\\\ DL_{3}+{c_{2}}(\varphi\cdot H)-c_{5}\varphi\end{array}$
5 | $\frac{\tilde{r}^{2}x_{3}^{2}\sqrt{F_{\pm}}}{c_{1}x_{3}^{2}\sqrt{F_{\pm}}+c_{2}r^{2}\sqrt{F_{\pm}}+c_{3}\tilde{r}^{2}(r^{2}\pm 1)}$ | $\frac{c_{4}x_{3}^{2}\sqrt{F_{\pm}}+c_{5}r^{2}\sqrt{F_{\pm}}+c_{6}\tilde{r}^{2}(r^{2}\pm 1)}{c_{1}x_{3}^{2}\sqrt{F_{\pm}}+c_{2}r^{2}\sqrt{F_{\pm}}+c_{3}\tilde{r}^{2}(r^{2}\pm 1)}$ | $\begin{array}[]{l}(K_{3}\pm P_{3})^{2}+\frac{F_{\pm}(c_{5}+c_{6}(r^{2}\pm 1))}{x_{3}^{2}}\\\ -\left(\frac{F_{\pm}(c_{2}+c_{3}(r^{2}\pm 1))}{x_{3}^{2}}\cdot H\right)-\hat{\eta},\\\ \\{P_{3},(K_{3}\pm P_{3})\\}+c_{5}\frac{(\tilde{r}^{2}\pm 1)}{x_{3}^{2}}\\\ +c_{6}\frac{x_{3}^{2}(\tilde{r}^{2}\mp 1)+(\tilde{r}^{2}\pm 1)^{2}}{x_{3}^{2}\sqrt{F_{\pm}}}-\left(c_{2}\frac{(\tilde{r}^{2}\pm 1)}{x_{3}^{2}}\cdot H\right)\\\ {}^{\prime}-\left(c_{3}\frac{x_{3}^{2}(\tilde{r}^{2}\mp 1)+(\tilde{r}^{2}\pm 1)^{2}}{x_{3}^{2}\sqrt{F_{\pm}}}\cdot H\right)\end{array}$
6 | $\frac{\tilde{r}^{2}x_{3}^{2}\sqrt{F_{\pm}}}{c_{1}x_{3}^{2}\sqrt{F_{\pm}}+c_{3}r^{2}\sqrt{F_{\pm}}+c_{3}(r^{2}\pm 1)\tilde{r}^{2}}$ | $\frac{c_{4}x_{3}^{2}\sqrt{F_{\pm}}+c_{5}r^{2}\sqrt{F_{\pm}}+c_{6}(r^{2}\pm 1)\tilde{r}^{2}}{c_{1}x_{3}^{2}\sqrt{F_{\pm}}+c_{2}r^{2}\sqrt{F_{\pm}}+c_{3}(r^{2}\pm 1)\tilde{r}^{2}}$ | $\begin{array}[]{c}K_{3}^{2}-P_{3}^{2}-\hat{\eta}+\frac{c_{5}(r^{4}-1)}{x_{3}^{2}}\\\ +\frac{c_{6}(r^{2}\mp 1)(r^{4}\pm 1+2\tilde{r}^{2})}{x_{3}^{2}\sqrt{F_{\pm}}}\\\ -\left(\frac{c_{2}(r^{4}-1)}{x_{3}^{2}}\cdot H\right)\\\ -\left(\frac{c_{3}(r^{2}\mp 1)(r^{4}\pm 1+2\tilde{r}^{2})}{x_{3}^{2}\sqrt{F_{\pm}}}\cdot H\right),\\\ (K_{3}\pm P_{3})^{2}-\hat{\eta}+\frac{c_{5}(r^{2}\pm 1)}{x_{3}^{2}}\\\ +\frac{c_{6}(r^{2}\pm 1)^{2}\sqrt{F_{\pm}}}{x_{3}^{2}}-\left(\frac{c_{2}(r^{2}\pm 1)}{x_{3}^{2}}\cdot H\right)\\\ -\left(\frac{c_{3}(r^{2}\pm 1)^{2}\sqrt{F_{\pm}}}{x_{3}^{2}}\cdot H\right)\end{array}$

We will not represent the following steps and the calculation details since
any of the related systems includes four or more equations for two unknowns
depending on two variables which make them rather simple solvable. The
obtained results are represented in Table 3 and equations (172)-(175). Let us
however represent the examples of solutions which are missing in the table and
mentioned formulae.

An example of solutions which are not ”good” since the related PDM system is
invariant with respect to two parametric group is:

$\displaystyle
M=\frac{c_{1}}{\tilde{r}^{2}}+\frac{c_{2}x_{3}}{r\tilde{r}^{2}}.$ (170)

Function (170) solves the system of equations (88), (108). However, the
related PDM Scrödinger equation is invariant with respect to two parametric
Lie group including rotations around the third coordinate axis and
dilatations. Such systems were discussed in [39].

One more nice solution looks as:

$\displaystyle
M=\frac{c_{1}}{x_{3}^{2}}+\frac{1}{\tilde{r}^{2}}\left(c_{2}+c_{3}\ln\left(\frac{\tilde{r}}{r^{2}}\right)\right).$

The corresponding integrals of motion are $\hat{Q}_{4}$ (56) and
$K_{3}^{2}+P_{a}q(\tilde{r},x_{3})P_{a}+\eta$. We will not represent the
explicit form of the related functions $q(\tilde{r},x_{3}),\eta$ and
$q_{4}(\tilde{r},x_{3}),\tilde{\eta}_{4}$ in (56) since the latter integral of
motion can be reduced to $Q_{1}$ via transformation (22). Notice that this
transformation keeps operator $\hat{Q}_{4}$ invariant.

$\displaystyle\begin{split}&f=\frac{\tilde{F}_{\pm}^{2}(r^{2}\mp
1)^{2}\tilde{r}^{2}}{c_{1}(r^{2}\mp
1)^{2}\tilde{F}_{\pm}^{2}+c_{2}x_{3}\tilde{r}^{2}(r^{4}\pm 1)(r^{2}\mp
1)+c_{3}\Phi_{\pm}+c_{4}\Lambda_{\pm}},\\\ &V=\frac{c_{5}(r^{2}\mp
1)\tilde{F}_{\pm}^{2}+c_{6}x_{3}\tilde{r}^{2}(r^{4}\pm 1)(r^{2}\mp
1)+c_{7}\Phi_{\pm}+c_{8}\Lambda_{\pm}}{c_{1}(r^{2}\mp
1)\tilde{F}_{\pm}^{2}+c_{2}x_{3}\tilde{r}^{2}(r^{4}\pm 1)(r^{2}\mp
1)+c_{3}\Phi_{\pm}+c_{4}\Lambda_{\pm}},\\\ &\Phi_{\pm}=(r^{2}\pm
1)^{2}((r^{2}\mp 1)^{2}\pm r)+4x_{3}(4r^{2}-\tilde{r}^{2}),\\\
&\Lambda_{\pm}=(\tilde{r}^{2}+4x_{3}^{2})(r^{2}\pm
1)^{2}-4x_{3}^{2}\tilde{r}^{2},\\\ &(K_{3}\pm
P_{3})^{2}-\hat{\eta}-\frac{1}{\tilde{F}_{\pm}^{2}}(c_{6}x_{3}(r^{2}\pm
1)F_{-}+4(c_{8}-3c_{7})x_{3}^{2}(r^{2}\pm 1)^{2})\\\
&+\left(\frac{1}{\tilde{F}_{\pm}^{2}}(c_{2}x_{3}(r^{2}\pm
1)F_{\pm}+4(c_{4}-3c_{3})x_{3}^{2}(r^{2}\pm 1)^{2}\cdot H\right),\\\
&\\{D,(K_{3}\mp
P_{3})\\}-15x_{3}-\frac{c_{6}F_{\pm}}{8\tilde{F}_{\pm}^{2}}+\left(-\frac{c_{2}F_{\pm}}{8\tilde{F}_{\pm}^{2}}\cdot
H\right)\\\ &-\frac{c_{7}x_{3}(r^{2}\mp 1)^{4}+(c_{8}-2c_{7}x_{3})((r^{2}\pm
1)^{2}+2\tilde{r}^{2})(2r^{2}-\tilde{r}^{2})(r^{2}\pm
1)}{\tilde{F}_{\pm}^{2}(r^{2}\mp 1)^{2}}\\\ &+\left(\frac{c_{3}x_{3}(r^{2}\mp
1)^{4}+(c_{4}-2c_{3}x_{3})((r^{2}\pm
1)^{2}+2\tilde{r}^{2})(2r^{2}-\tilde{r}^{2})(r^{2}\pm
1)}{\tilde{F}_{\pm}^{2}(r^{2}\mp 1)^{2}}\cdot H\right);\end{split}$ (171)
$\displaystyle\begin{split}&f=\frac{(r^{4}\pm
1)^{2}x_{3}^{2}\tilde{r}^{2}}{x_{3}^{2}(c_{1}(r^{4}\mp
1)^{2}\tilde{r}^{2}+c_{3}\tilde{r}^{2}(r^{4}\pm
1)+c_{4}r^{2}\tilde{r}^{2})+c_{2}(r^{4}\pm 1\mp
2r^{2}\tilde{r}^{2})\tilde{r}^{2}},\\\ &V=\frac{x_{3}^{2}(c_{5}(r^{4}\mp
1)^{2}\tilde{r}^{2}+c_{7}\tilde{r}^{2}(r^{4}\pm
1)+c_{8}r^{2}\tilde{r}^{2})+c_{6}(r^{4}\pm 1\mp
2r^{2}\tilde{r}^{2})\tilde{r}^{2}}{x_{3}^{2}(c_{1}(r^{4}\mp
1)^{2}\tilde{r}^{2}+c_{3}\tilde{r}^{2}(r^{4}\pm
1)+c_{4}r^{2}\tilde{r}^{2})+c_{2}(r^{4}\pm 1\mp
2r^{2}\tilde{r}^{2})\tilde{r}^{2}},\\\ &Q_{1}:=K_{3}^{2}\pm
P_{3}^{2}-\hat{\eta}-\frac{(1\pm r^{4})(c_{8}x_{3}^{4}+c_{6}((r^{4}\mp
1)^{2}\pm 4x_{3}^{4})+2c_{7}r^{2}x_{3}^{4})}{x_{3}^{2}(r^{4}\mp 1)^{2}}\\\
&+\left(\frac{(1\pm r^{4})(c_{4}x_{3}^{4}+c_{2}((r^{4}\mp 1)^{2}\pm
4x_{3}^{4})+2c_{3}r^{2}x_{3}^{4})}{x_{3}^{2}(r^{4}\mp 1)^{2}}\cdot
H\right),\end{split}$ (172)
$\displaystyle\begin{split}&Q_{2}=\\{K_{3},P_{3}\\}+\frac{c_{6}(\tilde{r}^{2}(r^{4}\mp
1)^{2}\pm 2x_{3}^{4}r^{2})-c_{7}(r^{4}\pm
1)x_{3}^{4}-2c_{8}r^{2}x_{3}^{4}}{2x_{3}^{2}(r^{4}\mp 1)^{2}}\\\
&-\left(\frac{c_{2}(\tilde{r}^{2}(r^{4}\mp
1)^{2}-2x_{3}^{4}r^{2})-c_{3}(r^{4}\pm
1)x_{3}^{4}-2c_{4}r^{2}x_{3}^{4}}{2x_{3}^{2}(r^{4}\mp 1)^{2}}\cdot
H\right);\end{split}$
$\displaystyle\begin{split}&f=\frac{\tilde{r}^{2}(r^{2}\pm
1)^{2}\sqrt{F_{\pm}}}{\sqrt{F_{\pm}}(c_{1}(r^{2}\pm 1)^{2}+c_{2}(r^{2}\pm
1)^{2}-2\tilde{r}^{2}))+c_{3}(r^{4}-1)(r^{2}\pm
1)+c_{4}x_{3}\tilde{r}^{2}},\\\
&V=\frac{c_{3}\tilde{F}_{\pm}^{2}+c_{4}x_{3}\tilde{r}^{2}(r^{2}\mp
1)}{c_{1}\tilde{F}_{\pm}^{2}+c_{2}x_{3}\tilde{r}^{2}(r^{2}\mp 1)},\\\
&Q_{1}=(K_{3}\pm
P_{3})^{2}-\hat{\eta}+\frac{8c_{6}x_{3}^{2}+c_{8}x_{3}\sqrt{F_{\pm}}}{(r^{2}\pm
1)^{2}}-\left(\frac{8c_{2}x_{3}^{2}+c_{4}x_{3}\sqrt{F_{\pm}}}{(r^{2}\pm
1)^{2}}\cdot H\right),\\\ &Q_{2}=\\{D,K_{3}\pm
P_{3}\\}-15x_{3}-4\left(\left(\frac{c_{3}x_{3}}{\sqrt{F_{\pm}}}+\frac{c_{1}x_{3}(r^{2}\mp
1)}{\sqrt{F_{\pm}}(r^{2}\pm 1)^{2}}\right)\cdot H\right)\\\
&-\frac{1}{4}\left(\frac{c_{4}(r^{2}\mp 1)((r^{2}\pm 1)^{2}\mp
8x_{3}^{2})}{(r^{2}\pm 1)^{2}\sqrt{F_{\pm}}}\cdot H\right)\\\
&+4\left(\frac{c_{7}x_{3}}{\sqrt{F_{\pm}}}+\frac{c_{5}x_{3}(r^{2}\mp
1)}{(r^{2}\pm 1)^{2}}\right)+\frac{c_{8}(r^{2}\mp 1)((r^{2}\pm 1)^{2}\mp
8x_{3}^{2})}{4(r^{2}\pm 1)^{2}\sqrt{F_{\pm}}};\end{split}$ (173)
$\displaystyle\begin{split}&f=\frac{\tilde{r}^{2}\sqrt{F_{\pm}}}{c_{1}\sqrt{F_{\pm}}+c_{2}(x_{3}-1)\sqrt{r^{2}\pm
1+2x_{3}}+c_{3}(x_{3}+1)\sqrt{r^{2}\pm 1-2x_{3}}+c_{4}(r^{2}\mp 1)},\\\
&V=\frac{c_{5}\sqrt{F_{\pm}}+c_{6}(x_{3}-1)\sqrt{r^{2}\pm
1+2x_{3}}+c_{7}(x_{3}+1)\sqrt{r^{2}\pm 1-2x_{3}}+c_{8}(r^{2}\mp
1)}{c_{1}\sqrt{F_{\pm}}+c_{2}(x_{3}-1)\sqrt{r^{2}\pm
1+2x_{3}}+c_{3}(x_{3}+1)\sqrt{r^{2}\pm 1-2x_{3}}+c_{4}(r^{2}\mp 1)},\\\
&Q_{1}=\\{D,K_{3}\pm P_{3}\\}-15x_{3}+\frac{c_{6}(r^{2}\pm 1)}{\sqrt{r^{2}\pm
1-2x_{3}}}+\frac{c_{7}(r^{2}\pm 1)}{\sqrt{r^{2}\pm
1+2x_{3}}}+\frac{4c_{8}}{\sqrt{F\pm}}\\\ &-\left(\left(\frac{c_{2}(r^{2}\pm
1)}{\sqrt{r^{2}\pm 1-2x_{3}}}+\frac{4c_{4}}{sqrt{F_{\pm}}}\right)\cdot
H\right)-\left(\frac{c_{3}(r^{2}\pm 1)}{\sqrt{r^{2}\pm 1+2x_{3}}}\cdot
H\right),\\\ &Q_{2}=\\{P_{3},K_{3}\pm
P_{3}\\}+\frac{c_{6}x_{3}}{\sqrt{r^{2}\pm
1-2x_{3}}}+\frac{c_{7}x_{3}}{\sqrt{r^{2}\pm 1+2x_{3}}}+\frac{c_{8}(r^{2}\pm
1)}{\sqrt{F_{\pm}}}\\\ &-\left(\frac{c_{3}x_{3}}{\sqrt{r^{2}\pm
1+2x_{3}}}\cdot H\right)-\left(\left(\frac{c_{2}x_{3}}{\sqrt{r^{2}\pm
1-2x_{3}}}+\frac{c_{4}(r^{2}\pm 1)}{\sqrt{F_{\pm}}}\right)\cdot
H\right);\end{split}$ (174) $\displaystyle\begin{split}&f=\frac{(r^{2}\pm
1)^{2}r\tilde{r}^{2}}{c_{1}(r^{2}\pm
1)^{2}+c_{2}r\tilde{r}^{2}+c_{3}x_{3}(r^{2}\pm
1)^{2}+c_{4}\tilde{r}^{2}(r^{2}\mp 1)},\\\ &V=\frac{c_{5}(r^{2}\pm
1)^{2}+c_{6}r\tilde{r}^{2}+c_{7}x_{3}(r^{2}\pm
1)^{2}+c_{8}\tilde{r}^{2}(r^{2}\mp 1)}{c_{1}(r^{2}\pm
1)^{2}+c_{2}r\tilde{r}^{2}+c_{3}x_{3}(r^{2}\pm
1)^{2}+c_{4}\tilde{r}^{2}(r^{2}\mp 1)},\\\
&Q_{1}=\\{K_{3},P_{3}\\}+\frac{c_{7}x_{3}(r^{2}\pm
1)^{2}-x_{3}^{2}(c_{6}r+c_{8}(r^{2}\mp 1))}{r(r^{2}\pm 1)^{2}}\\\
&-\left(\frac{c_{3}x_{3}(r^{2}\pm 1)^{2}-x_{3}^{2}(c_{2}r+c_{4}(r^{2}\mp
1))}{r(r^{2}\pm 1)^{2}}\cdot H\right),\\\ &Q_{2}=\\{D,(K_{3}\pm
P_{3})\\}-\hat{\eta}+\frac{2c_{6}x_{3}(r^{2}\mp 1)}{(r^{2}\pm
1)^{2}}+\frac{c_{8}(r^{4}\mp r^{2}+1)-c_{7}(r^{4}-1)(r^{2}\pm 1)}{r(r^{2}\pm
1)^{2}}\\\ &-\left(\frac{2c_{2}x_{3}r(r^{2}\mp 1)-c_{3}(r^{4}-1)(r^{2}\pm
1)+c_{4}(r^{4}\mp r^{2}+1)}{r(r^{2}\pm 1)^{2}}\cdot H\right).\end{split}$
(175) $\displaystyle\begin{split}&f=\frac{\tilde{F}_{\pm}^{2}(r^{2}\mp
1)^{2}\tilde{r}^{2}}{c_{1}(r^{2}\mp
1)^{2}\tilde{F}_{\pm}^{2}+c_{2}x_{3}\tilde{r}^{2}(r^{4}\pm 1)(r^{2}\mp
1)+c_{3}\Phi_{\pm}+c_{4}\Lambda_{\pm}},\\\ &V=\frac{c_{5}(r^{2}\mp
1)\tilde{F}_{\pm}^{2}+c_{6}x_{3}\tilde{r}^{2}(r^{4}\pm 1)(r^{2}\mp
1)+c_{7}\Phi_{\pm}+c_{8}\Lambda_{\pm}}{c_{1}(r^{2}\mp
1)\tilde{F}_{\pm}^{2}+c_{2}x_{3}\tilde{r}^{2}(r^{4}\pm 1)(r^{2}\mp
1)+c_{3}\Phi_{\pm}+c_{4}\Lambda_{\pm}},\\\ &\Phi_{\pm}=(r^{2}\pm
1)^{2}((r^{2}\mp 1)^{2}\pm r)+4x_{3}(4r^{2}-\tilde{r}^{2}),\\\
&\Lambda_{\pm}=(\tilde{r}^{2}+4x_{3}^{2})(r^{2}\pm
1)^{2}-4x_{3}^{2}\tilde{r}^{2},\\\ &Q_{1}=(K_{3}\pm
P_{3})^{2}-\hat{\eta}-\frac{1}{\tilde{F}_{\pm}^{2}}(c_{6}x_{3}(r^{2}\pm
1)F_{-}+4(c_{8}-3c_{7})x_{3}^{2}(r^{2}\pm 1)^{2})\\\
&+\left(\frac{1}{\tilde{F}_{\pm}^{2}}(c_{2}x_{3}(r^{2}\pm
1)F_{\pm}+4(c_{4}-3c_{3})x_{3}^{2}(r^{2}\pm 1)^{2}\cdot H\right),\\\
&Q_{2}=\\{D,(K_{3}\mp
P_{3})\\}-15x_{3}-\frac{c_{6}F_{\pm}}{8\tilde{F}_{\pm}^{2}}+\left(-\frac{c_{2}F_{\pm}}{8\tilde{F}_{\pm}^{2}}\cdot
H\right)\\\ &-\frac{c_{7}x_{3}(r^{2}\mp 1)^{4}+(c_{8}-2c_{7}x_{3})((r^{2}\pm
1)^{2}+2\tilde{r}^{2})(2r^{2}-\tilde{r}^{2})(r^{2}\pm
1)}{\tilde{F}_{\pm}^{2}(r^{2}\mp 1)^{2}}\\\ &+\left(\frac{c_{3}x_{3}(r^{2}\mp
1)^{4}+(c_{4}-2c_{3}x_{3})((r^{2}\pm
1)^{2}+2\tilde{r}^{2})(2r^{2}-\tilde{r}^{2})(r^{2}\pm
1)}{\tilde{F}_{\pm}^{2}(r^{2}\mp 1)^{2}}\cdot H\right).\end{split}$ (176)

Summarizing, we find all inequivalent superintegrable PDM systems admitting
scalar integrals of motion. The number of such systems appears to be rather
extended, see Table 3 and the following formulae (171)-(175). The latter
formulae collect such expressions which are too cumbersome to be placed in a
table.

###  9.2 Superintegrable systems with combined symmetries

The next class of superintegrable PDM systems which we consider are those ones
which admit both the scalar and vector or tensor integrals of motion. Since
the systems admitting the non-scalar integral of motion up to the only
exception has a rather simple shape defined up to arbitrary parameters, it is
relatively easy to solve the functional equations (169) (where $M^{(1)}$ and
and $M^{(2)}$ are position dependent masses of the systems with different
types of symmetries) in order to select the cases when there exist additional
integrals of motion and so the related systems are superintegrable. To verify
the completeness of the obtained in this way list of superintegrable systems
it is reasonable to verify the consistence of the systems of the determining
equations corresponding to the combinations of symmetries missing in this
list.

Table 4. Inverse masses, potentials and combined integrals of motion for
superintegrable systems

No | $f$ | $V$ | Integrals of motion
---|---|---|---
1 | $\frac{\tilde{r}x_{3}^{2}}{\tilde{r}F\left(\frac{r^{2}\pm 1}{x_{3}}\right)+c_{1}x_{3}^{2}(r^{2}\mp 1)F_{\pm}}$ | $\frac{\tilde{r}G\left(\frac{r^{2}\pm 1}{x_{3}}\right)+{c_{2}}x_{3}^{2}(r^{2}\mp 1)F_{\pm}}{\tilde{r}F\left(\frac{r^{2}\pm 1}{x_{3}}\right)+c_{1}x_{3}^{2}(r^{2}\mp 1)F_{\pm}}$ | $\begin{array}[]{c}\\{L_{3},(K_{\alpha}\pm P_{\alpha})\\}+3\varepsilon_{\alpha\nu}x_{\nu}\\\ +2c_{1}\left(\frac{\varepsilon_{\alpha\nu}x_{\nu}}{\tilde{r}}\cdot H\right)-2{c_{2}}\frac{\varepsilon_{\alpha\nu}x_{\nu}}{\tilde{r}},\\\ (K_{3}\pm p_{3})^{2}+\frac{G\left(\frac{r^{2}\pm 1}{x_{3}}\right)F_{\pm}}{x_{3}^{2}}\\\ -\left(\frac{F\left(\frac{r^{2}\pm 1}{x_{3}}\right)F_{\pm}}{x_{3}^{2}}\cdot H\right)\end{array}$
2 | $\frac{(r^{4}-1)^{2}x_{3}^{2}}{(r^{2}\pm 1)^{2}({c_{1}}(r^{2}\mp 1)^{2}+{c_{2}}x_{3}^{2})-c_{3}x_{3}^{2}r^{2}}$ | $\frac{(r^{2}\pm 1)^{2}({c_{4}}(r^{2}\mp 1)^{2}+{c_{5}}x_{3}^{2})-{c_{6}}x_{3}^{2}r^{2}}{(r^{2}\pm 1)^{2}({c_{1}}(r^{2}\mp 1)^{2}+{c_{2}}x_{3}^{2})-c_{1}x_{3}^{2}r^{2}}$ | $\begin{array}[]{c}(K_{1}\pm P_{1})^{2}-(K_{2}\pm P_{2})^{2}\\\ +15(x_{1}^{2}-x_{2}^{2})\\\ +c_{1}\left(\frac{x_{1}^{2}-x_{2}^{2}}{(r^{2}\pm 1)^{2}}\cdot H\right)+{c_{2}}\frac{x_{1}^{2}-x_{2}^{2}}{(r^{2}\pm 1)^{2}},\\\ \\{(K_{1}\pm P_{1}),(K_{2}\pm P_{2})\\}+15x_{1}x_{2}\\\ +c_{1}\left(\frac{x_{1}x_{2}}{(r^{2}\pm 1)^{2}}\cdot H\right)+{c_{2}}\frac{2x_{1}x_{2}}{(r^{2}\pm 1)^{2}},\\\ \\{L_{1},L_{2}\\}-\frac{{c_{4}}x_{1}x_{2}}{x_{3}^{2}}+\left(\frac{{c_{1}}x_{1}x_{2}}{x_{3}^{2}}\cdot H\right),\\\ L_{1}^{2}-L_{2}^{2}-\frac{c_{4}(x_{1}^{2}-x_{2}^{2})}{2x_{3}^{2}}\\\ +\left(\frac{c_{1}(x_{1}^{2}-x_{2}^{2})}{2x_{3}^{2}}\cdot H\right)\end{array}$
3 | $\frac{1}{c_{1}\tilde{r}^{2}+F(x_{3})}$ | $\frac{{c_{2}}\tilde{r}^{2}+G(x_{3})}{c_{1}\tilde{r}^{2}+F(x_{3})}$ | $\begin{array}[]{c}P_{1}P_{2}-(c_{1}x_{1}x_{2}\cdot H)+{c_{2}}x_{1}x_{2},\\\ P_{3}^{2}-\left(F(x_{3})\cdot H\right)+G(x_{3}),\\\ P_{1}^{2}-P_{2}^{2}+{c_{2}}(x_{1}^{2}-x_{2}^{2})-\\\ (c_{1}(x_{1}^{2}-x_{2}^{2})\cdot H)\end{array}$
4 | $\frac{x_{3}^{2}}{x_{3}^{2}(c_{1}r^{2}+{c_{2}})+{c_{1}}}$ | $\frac{x_{3}^{2}({c_{3}}r^{2}+{c_{2}})+c_{4}}{x_{3}^{2}(c_{1}r^{2}+{c_{2}})+{c_{1}}}$ | $\begin{array}[]{c}L_{1}^{2}-L_{2}^{2}-\frac{c_{4}(x_{1}^{2}-x_{2}^{2})}{2x_{3}^{2}}+\left(\frac{{c_{1}}(x_{1}^{2}-x_{2}^{2})}{2x_{3}^{2}}\cdot H\right),\\\ \\{L_{1},L_{2}\\}-\frac{c_{4}x_{1}x_{2}}{x_{3}^{2}}+\left(\frac{{c_{1}}x_{1}x_{2}}{x_{3}^{2}}\cdot H\right)\\\ P_{3}^{2}+{c_{3}}x_{3}^{2}+\frac{c_{4}}{x_{3}^{2}}-\left((c_{1}x_{3}^{2}+\frac{c_{2}}{x_{3}^{2}})\cdot H\right)\end{array}$
5 | $\frac{\tilde{r}}{F(x_{3})\tilde{r}+c_{1}}$ | $\frac{G(x_{3})\tilde{r}+{c_{2}}}{F(x_{3})\tilde{r}+c_{1}}$ | $\begin{array}[]{c}\\{L_{3},P_{\alpha}\\}-{c_{2}}\frac{\varepsilon_{\alpha\nu}x_{\nu}}{\tilde{r}}+c_{1}\varepsilon_{\alpha\nu}\left(\frac{x_{\nu}}{\tilde{r}}\cdot H\right)\\\ P_{3}^{2}-\left(F(x_{3})\cdot H\right)+G(x_{3})\end{array}$
6 | $\frac{x_{3}^{2}r}{(c_{1}x_{3}^{2}+c_{3})r+2c_{2}x_{3}^{2}}$ | $\frac{(c_{5}x_{3}^{2}+c_{4})r+c_{6}x_{3}^{2}}{(c_{1}x_{3}^{2}+c_{3})r+2c_{2}x_{3}^{2}}$ | $\begin{array}[]{c}\\{P_{\alpha},D\\}-2\left(x_{\alpha}(c_{1}+\frac{c_{2}}{r})\cdot H\right)\\\ +2x_{\alpha}(2c_{5}+\frac{c_{6}}{r}),\\\ P_{3}^{2}-\left(\frac{c_{3}^{2}}{x_{3}^{2}}\cdot H\right)+\frac{c_{4}}{x_{3}^{2}}\end{array}$
7 | $((r^{2}\pm 1)^{2}\mp 4x_{3}^{2})^{\frac{1}{2}}$ | $c_{1}$ | $\begin{array}[]{c}\\{(K_{3}\pm P_{3},P_{\alpha}\\}+P_{n}x_{3}x_{\alpha}P_{n},\\\ \\{P_{3},(K_{3}\pm P_{3})\\}+P_{a}x_{3}^{2}P_{a}\end{array}$

In this way we find the superintegrable systems presented in Table 4 which
collects the systems admitting integrals of motion which transform in
different way under rotation transformations. More exactly, there are such
systems which admit both scalar and vector or tensor integrals of motion.
However, the system represented in Item 2 admits only tensorial integrals of
motion.

The final step is to classify the maximally superintegrable system. It can be
done in the same way as in the superintegrable case. The classification
results are presented in Table 5 which includes two parts.

Table 5. Inverse masses, potentials and integrals of motion for maximally
superintegrable systems

No | $f$ | $V$ | Integrals of motion
---|---|---|---
1 | $\frac{r^{2}{\tilde{r}}}{\tilde{r}F(r)+c_{1}x_{3}}$ | $\frac{\tilde{r}(G(r)+{c_{2}}x_{3})}{\tilde{r}F(r)+c_{1}x_{3}}$ | $\begin{array}[]{c}\\{L_{3},L_{\alpha}\\}-2\left(\frac{c_{1}x_{\alpha}}{\tilde{r}}\cdot H\right)+\frac{2{c_{2}}x_{\alpha}}{\tilde{r}},\\\ L_{1}^{2}+L_{2}^{2}-c_{1}(\frac{x_{3}}{\tilde{r}}\cdot H)+\frac{{c_{2}}x_{3}}{\tilde{r}},\\\ \\{K_{3},P_{\alpha}\\}+2\varepsilon_{\alpha\nu}DL_{\nu}\\\ +\left(\frac{x_{\alpha}}{r^{2}}(x_{3}F(r)-c_{1}\tilde{r})\cdot H\right)\\\ -\frac{x_{\alpha}}{r^{2}}(G(r)-c_{2}\tilde{r})+P_{a}x_{3}x_{\alpha}P_{a}\end{array}$
2 | $\frac{1}{c_{1}+2c_{2}x_{3}+c_{3}(\tilde{r}^{2}+4x_{3}^{2})}$ | $\frac{c_{4}+2c_{5}x_{3}+c_{6}(\tilde{r}^{2}+4x_{3}^{2})}{c_{1}+2c_{2}x_{3}+c_{3}(\tilde{r}^{2}+4x_{3}^{2})}$ | $\begin{array}[]{c}P_{3}^{2}+2c_{5}x_{3}+4c_{6}x_{3}^{2}\\\ -\left((2c_{2}x_{3}+4c_{3}x_{3}^{2})\cdot H\right),\\\ \\{P_{3},D\\}-P_{n}x_{3}P_{n}+c_{4}x_{3}+c_{5}(2x_{3}^{2}+\tilde{r}^{2})\\\ +c_{6}x_{3}(4x_{3}^{2}+3\tilde{r}^{2})-\left((c_{1}x_{3}+c_{2}(2x_{3}^{2}+\tilde{r}^{2})\right.\\\ \left.+c_{3}x_{3}(4x_{3}^{2}+3\tilde{r}^{2}))\cdot H\right),\\\ P_{1}L_{1}-P_{2}L_{2}\\\ -(x_{1}x_{2}(2c_{3}x_{3}+c_{2})\cdot H)+x_{1}x_{2}(2c_{6}x_{3}+c_{5}),\\\ \\{P_{1},L_{2}\\}+\\{P_{2},L_{1}\\}\\\ -((x_{1}^{2}-x_{2}^{2})(2c_{3}x_{3}+c_{2})\cdot H)\\\ +(x_{1}^{2}-x_{2}^{2})(2c_{6}x_{3}+c_{5})\end{array}$
3 | $\frac{\tilde{r}^{2}}{c_{1}+c_{2}x_{3}\tilde{r}}$ | $\frac{c_{3}+c_{4}x_{3}\tilde{r}}{c_{1}+c_{2}x_{3}\tilde{r}}$ | $\begin{array}[]{c}\\{L_{3},P_{\alpha}\\}+c_{1}\varepsilon_{\alpha\nu}\left(\frac{x_{\nu}}{\tilde{r}}\cdot H\right)-{c_{3}}\frac{x_{2}}{\tilde{r}},\\\ P_{3}^{2}-\left(c_{2}x_{3}\cdot H\right)+c_{4}x_{3},\\\ \\{P_{3},D\\}+\left(\left(\frac{c_{1}}{x_{3}^{2}}-c_{2}(x_{3}^{2}+\frac{\tilde{r}^{2}}{2}\right)\cdot H\right)\\\ -\frac{c_{3}}{x_{3}^{2}}-c_{4}(x_{3}^{2}+\frac{\tilde{r}^{2}}{2})\end{array}$
4 | $\frac{1}{c_{1}r^{2}+2c_{2}x_{3}+c_{3}}$ | $\frac{c_{4}r^{2}+2c_{5}x_{3}+c_{6}}{c_{1}r^{2}+2c_{2}x_{3}+c_{3}}$ | $\begin{array}[]{c}P_{3}P_{\alpha}-\left(x_{\alpha}(c_{1}x_{3}+c_{2})\cdot H\right)\\\ +x_{\alpha}(c_{4}x_{3}+c_{5}),\\\ P_{3}^{2}-\left((c_{1}x_{3}^{2}+2c_{2}x_{3})\cdot H\right)\\\ +c_{4}x_{3}^{2}+2c_{5}x_{3},\\\ P_{1}P_{2}-c_{1}(x_{1}x_{2}\cdot H)+c_{4}x_{1}x_{2},\\\ (P_{1}^{2}-P_{2}^{2})+c_{4}(x_{1}^{2}-x_{2}^{2})\\\ -c_{1}((x_{1}^{2}-x_{2}^{2})\cdot H)\\\ \end{array}$
5 | $\frac{x_{3}^{2}}{c_{1}x_{3}^{2}r^{2}+{c_{2}}}$ | $\frac{{c_{3}}x_{3}^{2}r^{2}+{c_{4}}}{c_{1}x_{3}^{2}r^{2}+{c_{2}}}$ | $\begin{array}[]{c}P_{1}P_{2}-c_{1}(x_{1}x_{2}\cdot H)+{c_{3}}x_{1}x_{2},\\\ P_{1}^{2}-P_{2}^{2}+{c_{3}}(x_{1}^{2}-x_{2}^{2})-{c_{1}}((x_{1}^{2}-x_{2}^{2})\cdot H),\\\ \\{L_{1},L_{2}\\}+\left(\frac{2{c_{2}}x_{1}x_{2}}{x_{3}^{2}}\cdot H\right)-\frac{2{c_{4}}x_{1}x_{2}}{x_{3}^{2}},\\\ L_{1}^{2}-L_{2}^{2}+\left(\frac{{c_{2}}(x_{1}^{2}-x_{2}^{2})}{x_{3}^{2}}\cdot H\right)-\frac{{c_{4}}(x_{1}^{2}-x_{2}^{2})}{x_{3}^{2}},\\\ P_{3}^{2}-\left((c_{1}x_{3}^{2}+\frac{c_{2}}{x_{3}^{2}})\cdot H\right)+{c_{3}}x_{3}^{2}+\frac{c_{4}}{x_{3}^{2}}\end{array}$

Table 5 (continued). Inverse masses, potentials and integrals of motion for
maximally superintegrable systems

No | $f$ | $V$ | Integrals of motion
---|---|---|---
6 | $\frac{\tilde{r}^{2}\sqrt{F_{\pm}}}{c_{1}\sqrt{F_{\pm}}+c_{2}(r^{2}\mp 1)}$ | $\frac{c_{3}\sqrt{F_{\pm}}+c_{4}(r^{2}\mp 1)}{c_{1}\sqrt{F_{\pm}}+c_{2}(r^{2}\mp 1)}$ | $\begin{array}[]{l}K_{3}^{2}-P_{3}^{2}-\hat{\eta}-\frac{2c_{4}(r^{2}\mp 1)}{\sqrt{F_{\pm}}}+\left(\frac{2c_{2}(r^{2}\mp 1)}{\sqrt{F_{\pm}}}\cdot H\right),\\\ \\{D,(K_{3}\pm P_{3})\\}-15x_{3}\mp\frac{4c_{4}x_{3}}{\sqrt{F_{\pm}}}\pm\left(\frac{4c_{2}x_{3}}{\sqrt{F_{\pm}}}\cdot H\right),\\\ (K_{3}\pm P_{3})^{2}-\hat{\eta}\end{array}$
7 | $\frac{r^{2}\tilde{r}^{2}}{c_{1}{r^{2}}+c_{2}r^{4}\tilde{r}^{2}+c_{3}r^{2}\tilde{r}^{2}}$ | $\frac{c_{4}{r^{2}}+c_{5}r^{4}\tilde{r}^{2}+c_{6}r^{2}\tilde{r}^{2}}{c_{1}{r^{2}}+c_{2}r^{4}\tilde{r}^{2}+c_{3}r^{2}\tilde{r}^{2}}$ | $\begin{array}[]{c}\\{P_{3},K_{3}\\}-2P_{a}x_{3}P_{a}\\\ +\left((\frac{c_{2}x_{3}^{2}}{\tilde{r}^{2}}-c_{3}x_{3}^{2})\cdot H\right)\\\ -\frac{c_{5}x_{3}^{2}}{\tilde{r}^{2}}+c_{6}x_{3}^{2},\\\ L_{1}^{2}+L_{2}^{2}-(c_{1}\frac{r^{2}}{\tilde{r}^{2}}\cdot H)+c_{4}\frac{r^{2}}{\tilde{r}^{2}},\\\ P_{3}^{2}-\left(c_{2}x_{3}^{2}\cdot H\right)+c_{5}x_{3}^{2}\end{array}$
8 | $\frac{1}{{c_{2}}(\tilde{r}^{2}+4x_{3}^{2})+2{c_{3}}x_{3}+c_{1}}$ | $\frac{{c_{6}}x_{3}+{c_{5}}(\tilde{r}^{2}+4x_{3}^{2})+c_{4}}{{c_{2}}(\tilde{r}^{2}+4x_{3}^{2})+2{c_{3}}x_{3}+c_{1}}$ | $\begin{array}[]{c}P_{1}L_{1}-P_{2}L_{2}+x_{1}x_{2}(2{c_{5}}x_{3}+{c_{6}})\\\ -(x_{1}x_{2}(2{c_{2}}x_{3}+{c_{3}})\cdot H),\\\ \\{P_{1},L_{2}\\}+\\{P_{2},L_{1}\\}+(x_{1}^{2}-x_{2}^{2})(2{c_{5}}x_{3}+{c_{6}})\\\ -((x_{1}^{2}-x_{2}^{2})(2{c_{2}}x_{3}+{c_{3}})\cdot H)\\\ P_{3}^{2}-2\left((2{c_{2}}x_{3}^{2}+{c_{3}}x_{3})\cdot H\right)\\\ +4{c_{5}}x_{3}^{2}+{c_{6}}x_{3},\\\ P_{1}P_{2}-{c_{2}}(x_{1}x_{2}\cdot H)+{c_{5}}x_{1}x_{2},\\\ P_{1}^{2}-P_{2}^{2}-c_{1}((x_{1}^{2}-x_{2}^{2})\cdot H)+c_{4}(x_{1}^{2}-x_{2}^{2})\end{array}$
9 | $\frac{x_{3}^{2}}{c_{1}x_{3}^{2}(x_{3}^{2}+4\tilde{r}^{2})+c_{2}+c_{3}x_{3}^{2}}$ | $\frac{c_{4}+c_{5}x_{3}^{2}(x_{3}^{2}+4\tilde{r}^{2})+c_{6}x_{3}^{2}}{c_{1}x_{3}^{2}(x_{3}^{2}+4\tilde{r}^{2})+c_{2}+c_{3}x_{3}^{2}}$ | $\begin{array}[]{c}\\{P_{3},L_{\alpha}\\}+2\frac{\varepsilon_{\alpha\nu}x_{\nu}(c_{4}-c_{5}x_{3}^{4})}{x_{3}^{2}}+2\left(\frac{\varepsilon_{\alpha\nu}x\nu(c_{1}x_{3}^{4}-c_{2})}{x_{3}^{2}}\cdot H\right),\\\ P_{3}^{2}-\left((\frac{c_{1}}{x_{3}^{2}}+\frac{c_{2}}{x_{3}^{2}})\cdot H\right)\ +\frac{c_{4}}{x_{3}^{2}}+\frac{c_{5}}{x_{3}^{2}},\\\ P_{1}P_{2}-c_{1}(x_{1}x_{2}\cdot H)+c_{4}x_{1}x_{2},\\\ P_{1}^{2}-P_{2}^{2}-c_{1}((x_{1}^{2}-x_{2}^{2})\cdot H)+c_{4}(x_{1}^{2}-x_{2}^{2})\end{array}$
10 | $\frac{x_{3}^{2}(r^{2}\pm 1)}{c_{1}(r^{2}\pm 1)+c_{2}x_{3}^{2}}$ | $\frac{c_{3}(r^{2}\pm 1)+c_{4}x_{3}^{2}}{c_{1}(r^{2}\pm 1)+c_{2}x_{3}^{2}}$ | $\begin{array}[]{c}\\{D,(K_{\alpha}\pm P_{\alpha})\\}-15x_{\alpha}-\frac{2x_{\alpha}c_{4}(r^{2}\mp 1)}{(r^{2}\pm 1)^{2}}\\\ +\left(\frac{2x_{\alpha}c_{1}((r^{2}\mp 1)}{(r^{2}\pm 1)^{2}}\cdot H\right),\\\ \\{(K_{3}\pm P_{3}),L_{\alpha}\\}+3\varepsilon_{\alpha\nu}x_{\nu}-\frac{2c_{4}\varepsilon_{\alpha\nu}x_{\nu}(r^{2}\pm 1)}{x_{3}^{2}}\\\ +2c_{1}\left(\frac{\varepsilon_{\alpha\nu}x_{\nu}(r^{2}\pm 1)}{x_{3}^{2}}\cdot H\right),\\\ \\{L_{1},L_{2}\\}+\left(\frac{2c_{1}x_{1}x_{2}}{x_{3}^{2}}\cdot H\right)-\frac{2c_{3}x_{1}x_{2}}{x_{3}^{2}},\\\ L_{1}^{2}-L_{2}^{2}+\left(\frac{2c_{1}(x_{1}^{2}-x_{2}^{2})}{x_{3}^{2}}\cdot H\right)-\frac{2c_{3}(x_{1}^{2}-x_{2}^{2})}{x_{3}^{2}}\par\end{array}$

At this point the classification of cylindrically invariant PDM system
admitting second order integrals of motion has been completed.

##  10 Discussion

In contrast with the cases of the standard Schrödinger equations and 2D
Schrödinger equations with position dependent mass we still do not have the
completed description of second order integrals of motion for 3D PDM
Schrödinger equations. However, some steps to such description have been
already made: the maximally superintegrable sytems and separable systems have
been already classified [25, 26], there are successes in the classification of
nondegenerate and semidegenerate systems [28, 30]. Moreover, the systems
invariant w.r.t. three- and two-parametric Lie groups are classified
completely [35, 39].

In the present paper we make the next stem to the complete classification of
the mentioned integrals of motion. Namely, we present all inequivalent quantum
mechanical PDM systems which, in addition to the second order integrals of
motion, admit the fixed one parameter Lie symmetry group.

As it was shown in [17] there are six inequivalent one parametric Lie groups
which can be possessed by the PDM systems. We start with one of them, namely,
with the group of rotations around the fixed axis. In other words, we deal
with the cylindrically symmetric PDM systems and classify such of them which
admit at least one second order integral of motion.

Let us mention that the PDM systems with cylindric symmetry possess a rather
extended collection of the mentioned integrals. Namely, we have fixed as much
as 66 inequivalent systems and presented their integrals of motion. In
particular we specify 18 superintrgrable and 10 maximally superintgrable
systems. They are collected in five tables. In addition, the most cumbersome
of them are presented separately in formulae (171)-(175). Notice that any item
including the terms ”$\pm$” in fact represents two systems one of which
corresponds to the sign ”$+$” and the other to the sign”$-$”.

To optimize calculations we separate the integrals of motion to three
qualitatively different subclasses in accordance with their transformation
properties with respect to the rotations which by definition leaves the PDM
Hamiltonians invariant. The mentioned subclasses include the scalar, vector
and tensor versions of the integrals. Moreover, the scalar and tensor
integrals of motion can be effective separated in accordance with their parity
properties.

The systems admitting one scalar integral of motion are defined up to
arbitrary functions which depend on specific variables. Such (integrable)
systems are presented in Table 1.

The majority of systems admitting vector or tensor integrals of motion are
defined more strictly and includes only arbitrary parameters. The reason of it
is that for these subclasses the related PDM system is supposed to admit as
minimum two linearly independent integrals of motion. The same is true for the
case of superintegrable systems admitting integrals of motion of arbitrary
type.

Thus we have made an essential step to the complete classification of the 3D
PDM systems admitting second order integrals of motion. In spite of that we
consider only one out of six inequivalent one parametric Lie groups which can
be accepted by such systems, this step is very important since the number of
found systems is very large, maybe more large then the total number of all
systems admitting the other inequivalent Lie symmetries. The latter statement
is supported by our computing experiments.

We believe that the presented classification is complete. However, the
determining equations which we solve to find the inequivalent systems and
their integrals of motion are rather complicated systems of partial
differential equations with variable coefficients, and there is a danger to
overlook some special solutions additional to the found generic ones. That is
why we present the mentioned determining equations explicitly. And the
absolutely rigorous statement (which, however, is tot constructive) is that
the discussed systems should include the arbitrary elements which are
functions solving these equations.

The next planned steps to the complete classification of the integrals of
motion admitted by the 3D PDM systems presuppose the classifications of the
systems which possess symmetries with respect to the remaining inequivalent
one parametric Lie groups for such systems specified in [17]. Finally, we
plane to classify such integrals of motion for the systems which have no Lie
symmetry. The latter problem appears not to be catastrophically complicated
thanks to the existence of rather strong equivalence group.

Acknowledgement. I am indebted with Universitá del Piemonte Orientale and and
Dipartimento di Scienze e Innovazione Tecnologica for the extended stay as
Visiting Professor.

## References

  * [1] W. Pauli, Jr. Über das Wasserstoffspektrum vom Standpunkt der neuen Quantenmechanik. Z. Physik 36, 336-363 (1926).
  * [2] P. Winternitz, J. Smorodinsky, M. Uhlir̆, and I. Fris̆. Symmetry groups in classical and quantum mechanics. Yad. Fiz. 4, 625-635 (1966) (English translation: Sov. J. Nucl. Phys. 4, 444-450 (1967)).
  * [3] A. Makarov, J. Smorodinsky, Kh. Valiev and P. Winternitz. A systematic search for non-relativistic systems with dynamical symmetries. Nuovo Cim. A 52, 1061-1084 (1967).
  * [4] N. Evans. Group theory of the Smorodinsky-Winternitz system. J. Math. Phys. 32, 3369-3375 (1991).
  * [5] N. W. Evans. Super-integrability of the Winternitz system. Phys. Lett. 147, 483-486 (1990).
  * [6] W. Miller, Jr, S. Post and P. Winternitz Classical and Quantum Superintegrability with Applications J. Phys. A: Math. Theor. 46 423001 (2013).
  * [7] P. Winternitz and I. Yurdusen. Integrable and superintegrable systems with spin. J. Math. Phys. 47 103509 (2006).
  * [8] J. -F. Désilets, P. Winternitz and I. Yurdusen. Superintegrable systems with spin and second-order integrals of motion. Phys. A: Math. Theor. 45 475201 (2012).
  * [9] A. G. Nikitin. New exactly solvable systems with Fock symmetry. J. Phys. A: Math. Theor. 45 485204 (2012)
  * [10] A. G. Nikitin. Matrix superpotentials and superintegrable systems for arbitrary spin J. Phys. A: Math. Theor. 45 225205 (2012)
  * [11] A. M. Escobar-Ruiz, R. Linares and P. Winternitz, New infinite families of Nth-order superintegrable systems separating in Cartesian coordinates, J. Phys. A: Math. Theor. 53 445203 (2020).
  * [12] A. G. Nikitin, Higher-order symmetry operators for Schrödinger equation. In CRM Proceedings and Lecture Notes (AMS), 37 , pp. 137–144 (2004).
  * [13] W. Miller, Symmetry and separation of variables (Cambridge University Press, 1984).
  * [14] Ian Marquette and Pavel Winternitz. Higher order quantum superintegrability: a new Painleve conjecture. Integrability, Supersymmetry and Coherent States. Springer, Cham, pp. 103-131. (2019).
  * [15] O. Rosas-Ortiz, Position-dependent mass systems: Classical and quantum pictures. In Geometric Methods in Physics XXXVIII (pp. 351-361). Birkhäuser, Cham (2020).
  * [16] von Roos O Position-dependent effective masses in semiconductor theory Phys. Rev. B 27, 7547–7552 ( 1983).
  * [17] A. G. Nikitin and T. M. Zasadko, Superintegrable systems with position dependent mass, J. Math. Phys. 56 042101 (2015)
  * [18] A.G Nikitin and T. M. Zasadko, Group classification of Schrödinger equations with position dependent mass, J. Phys. A: Math. Theor. 49, 365204 (2016).
  * [19] A. G. Nikitin, Kinematical invariance groups of the 3d Schrödinger equations with position dependent masses, J. Math. Phys. 58, 083508 (2017).
  * [20] Kalnins E.G., Kress J.M., Pogosyan G.S., Miller Jr. W., Completeness of superintegrability in two-dimensional constant-curvature spaces, J. Phys. A: Math. Gen. 34 (2001), 4705–4720,
  * [21] Kalnins E.G., Miller Jr. W., Post S., Contractions of 2D 2nd order quantum superintegrable systems and the Askey scheme for hypergeometric orthogonal polynomials, SIGMA 9 , 057 (2013).
  * [22] E. G. Kalnins, G. S. Pogosyan and W. Miller, Jr., Completeness of multiseparable superintegrability in two dimensions, Phys. Atomic Nuclei 65, 1033–1035 (2002).
  * [23] J. M. Kress, Equivalence of superintegrable systems in two dimensions, Phys. Atomic Nuclei 70, 560–566 (2007).
  * [24] J. Kress and K. Schöbel, An algebraic geometric classification of superintegrable systems in the Euclidean plane. Journal of Pure and Applied Algebra, 223, 1728-1752 (2019).
  * [25] E. G. Kalnins, J. M. Kress and W. Miller Jr . Classification of superintegrable systems in three dimensions. Bulg. J. Phys, 33(s1), 174-189 (2006).
  * [26] E. G. Kalnins, J. M. Kress and W. Miller, Jr., Second order superintegrable systems in conformally flat spaces. V. Two- and three-dimensional quantum systems, J. Math. Phys. 47 , 093501 (2006).
  * [27] B. K. Berndson, E. G. Kalnins and W. Miller, Jr. Toward Classification of 2nd Order Superintegrable Systems in 3-Dimensional Conformally Flat Spaces with Functionally Linearly Dependent Symmetry Operators, SIGMA 16 (2020).
  * [28] J. J. Capel and J. M. Kress, Invariant Classification of Second-order Conformally Flat Superintegrable Systems, J. Phys. A: Math. Theor. 47, 495202 (2014).
  * [29] A. Vollmer, Stäckel equivalence of non-degenerate superintegrable systems and invariant quadrics, SIGMA 17, 015 (2021).
  * [30] M. A. Escobar-Ruiz and W. Miller, Toward a classification of semidegenerate 3D superintegrable systems. J. Phys. A: Math. Theor., 50(9), 095203, (2017).
  * [31] J. F. Cariñena, M. F. Rañada and M. Santander, Superintegrability of three-dimensional Hamiltonian systems with conformally Euclidean metrics. Oscillator-related and Kepler-related systems, J. Phys. A: Math. Theor. 54, 105201 (2021).
  * [32] J. F. Cariñena, M. F. Rañada and M. Santander, Superintegrability on the three-dimensional spaces with curvature. Oscillator-related and Kepler-related systems on the sphere S3 and on the hyperbolic space H3, J. Phys. A: Math. Theor. 54, 365201 (2021).
  * [33] A. Ballesteros, A. Enciso, F. J. Herranz, O. Ragnisco and D. Riglioni, Superintegrable Oscillator and Kepler Systems on Spaces of Nonconstant Curvature via the Stäckel Transform, SIGMA 7, 048 (2011).
  * [34] O. Ragnisco and D. Riglioni, A Family of Exactly Solvable Radial Quantum Systems on Space of Non-Constant Curvature with Accidental Degeneracy in the Spectrum, SIGMA 6 097 (2010).
  * [35] A. G. Nikitin, Superintegrable quantum mechanical systems with position dependent masses invariant with respect to three parametric Lie groups. J. Math. Phys. DOI: 1063/5.0147792
  * [36] . A. G. Nikitin. Superintegrable and shape invariant systems with position dependent mass. J. Phys. A: Math. Theor. 48, 335201 (2015)
  * [37] E. Ferraro, A. Messina, A. G. Nikitin. Exactly solvable relativistic model with the anomalous interaction. Phys. Rev. A 81, 042108 (2010).
  * [38] J. Niederle, A. G. Nikitin. Relativistic wave equations for interacting, massive particles with arbitrary half-integer spins. Phys. Rev. D 64 (12), 125013 (2001).
  * [39] A. G. Nikitin, Superintegrable quantum mechanical systems with position dependent masses invariant with respect to two parametric Lie groups. J. Phys. A: Math. Theor. 56, 395203 (2023).
  * [40] J. Hietarinta, Pure quantum integrability Phys. Lett. A 246 97-104 (1998).
  * [41] A. G. Nikitin, Generalized Killing tensors of arbitrary valence and order, Ukrainian Mathematical Journal 43, 734–743 (1991).