Czech Technical University in Prague

Faculty of Nuclear Sciences and Physical Engineering

Department of Physics
---
Field of study: Mathematical Physics

Integrabiln\́mathrm{i} a superintegrabiln\́mathrm{i} systémy cylindrického
typu v magnetických pol\́mathrm{i}ch

Integrable and superintegrable systems of cylindrical type in magnetic fields

MASTER’S THESIS

Author: | Bc. Ondřej Kubů
---|---
Supervisor: | doc. Ing. Libor Šnobl, Ph.D.
Consultant: | Prof. Ing. Pavel Winternitz, Ph. D.
Year: | 2020

See pages 1,2 of zadani.pdf

Prohlášen\́mathrm{i}

Prohlašuji, že jsem svoj\́mathrm{i} diplomovou práci vypracoval samostatně a
použil jsem pouze podklady (literaturu, projekty, SW atd.) uvedené v
přiloženém seznamu.

Nemám závažný důvod proti použit\́mathrm{i} tohoto školn\́mathrm{i}ho
d\́mathrm{i}la ve smyslu $\S$ 60 Zákona č. 121/2000 Sb., o právu autorském, o
právech souvisej\́mathrm{i}c\́mathrm{i}ch s právem autorským a o změně
některých zákonů (autorský zákon).

Declaration

I declare that I wrote my master’s thesis independently and that I have not
made use of any aid (literature, projects, SW etc.) other than those
acknowledged.

I agree with the usage of this thesis in the sense of the $\S$ 60 Act 121/2000
(Copyright Act of the Czech Republic).

V Praze dne ……………….. | ………………………………….
---|---
| Bc. Ondřej Kubů

Acknowledgement

First, I would like to thank my supervisor doc. Ing. Libor Šnobl, Ph.D. for
perfect collaboration and support in writing the thesis. The weekly
consultations were immensely helpful for overcoming smaller or bigger problems
encountered during the work on the thesis. I am also very grateful for support
and help during the not-so-smooth administration of my stay in Montréal.

Second, I would like to thank my consultant and supervisor during the stay at
Université de Montréal, Prof. Ing. Pavel Winternitz, Ph.D. for warm
hospitality during the stay. I am very grateful for the opportunity to study
abroad, which is an invaluable experience, and it would not be possible if
Prof. Ing. Pavel Winternitz, Ph. D. did not accept to supervise my thesis.

And third I would like to thank my parents whose support allows me to focus on
my studies.

This work was supported by the Grant Agency of the Czech Technical University
in Prague, grant No. SGS19/183/OHK4/3T/14, and by the Czech Science Foundation
(Grant Agency of the Czech Republic), project 17-11805S.

Bc. Ondřej Kubů

Název práce: |
---|---
Integrabiln\́mathrm{i} a superintegrabiln\́mathrm{i} systémy cylindrického
typu v magnetických pol\́mathrm{i}ch
Autor: | Bc. Ondřej Kubů
Studijn\́mathrm{i} program: | Aplikace př\́mathrm{i}rodn\́mathrm{i}ch věd
Obor: | Matematická fyzika
Druh práce: | Diplomová práce
Vedouc\́mathrm{i} práce: |  doc. Ing. Libor Šnobl, Ph.D. Katedra fyziky, Fakulta jaderná a fyzikálně inženýrská, České vysoké učen\́mathrm{i} technické v Praze
Konzultant: |  Prof. Ing. Pavel Winternitz, Ph. D. Département de mathématiques et de statistique & Centre de recherches mathématique, Université de Montréal
Abstrakt: C\́mathrm{i}lem této práce je hledán\́mathrm{i}
integrabiln\́mathrm{i}ch a superintegrabiln\́mathrm{i}ch systémů cylindrického
typu s magnetickým polem. Po zformulován\́mathrm{i} kvantově mechanických
určuj\́mathrm{i}c\́mathrm{i}ch rovnic pro integrály pohybu druhého řádu v
cylindrických souřadnic\́mathrm{i}ch jsou nalezeny všechny kvadraticky
integrabiln\́mathrm{i} systémy cylindrického typu. Mezi nimi jsou hledány
systémy, které připouštěj\́mathrm{i} dodatečné integrály pohybu. Nejprve jsou
př\́mathrm{i}mým řešen\́mathrm{i}m určuj\́mathrm{i}c\́mathrm{i}ch rovnic
nalezeny všechny systémy s dodatečným integrálem prvn\́mathrm{i}ho řádu v
klasické i kvantové mechanice. Ukazuje se, že všechny tyto systémy již byly
známé a žádné dalš\́mathrm{i} neexistuj\́mathrm{i}. Nalezeny jsou také všechny
klasické systémy s dodatečným integrálem typu $L^{2}+\ldots$, respektive
$L_{x}p_{y}-L_{y}p_{x}+\ldots$, z nichž většina zat\́mathrm{i}m nebyla
publikována. Všechny nalezené superintegrabiln\́mathrm{i} systémy
připouštěj\́mathrm{i} integrál prvn\́mathrm{i}ho řádu $L_{z}$ a jejich
Hamiltonovy-Jacobiho, respektive Schrödingerovy, rovnice jsou vyřešeny
separac\́mathrm{i} proměnných v cylindrických souřadnic\́mathrm{i}ch, u
systémů prvn\́mathrm{i}ho řádu i v kartézských.
Kl\́mathrm{i}čová slova: |  integrabilita, superintegrabilita, kvantová korekce, magnetické pole, cylindrické souřadnice
Title: |
Integrable and superintegrable systems of cylindrical type in magnetic fields
Author: | Bc. Ondřej Kubů
Abstract: The goal of this thesis is the search for integrable and
superintegrable systems with magnetic field. We formulate the quantum
mechanical determining equations for second order integrals of motion in the
cylindrical coordinates and we find all quadratically integrable systems of
the cylindrical type. Among them we search for systems admitting additional
integrals of motion. We find all systems with an additional first order
integral both in classical and quantum mechanics. It turns out that all these
systems have already been known and no other exist. We also find all systems
with an additional integral of type $L^{2}+\ldots$, respectively
$L_{y}p_{y}-L_{x}p_{y}+\ldots$, of which the majority is new to the
literature. All found superintegrable systems admit the first order integral
$L_{z}$ and we solve their Hamilton-Jacobi and Schrödinger equations by
separation of variables in the cylindrical coordinates, for the first order
systems in the Cartesian coordinates as well.
Key words: |  integrability, superintegrability, quantum correction, magnetic field, cylindrical coordinates

###### Contents

  1. Introduction
  2. 1 Quantum integrable systems in cylindrical coordinates
    1. 1.1 Integrability and superintegrability
    2. 1.2 Classical determining equations
    3. 1.3 Determining equations in quantum mechanics
    4. 1.4 Quantum integrable systems
      1. 1.4.1 Case b) $\mathrm{rank}\,(M)=2$: $\psi^{\prime}(\phi)\neq 0$ and $\mu(Z)=0$
      2. 1.4.2 Case b) $\mathrm{rank}\,(M)=1$: $\psi^{\prime}(\phi)\neq 0$ and $\mu(Z)=0$
  3. 2 Found superintegrable systems
    1. 2.1 First order integrals
      1. 2.1.1 General considerations
      2. 2.1.2 Case $\mathrm{rank}(M)=2$ subcase 2a)
      3. 2.1.3 Case $\mathrm{rank}(M)=2$ subcase 2b 1.1)
      4. 2.1.4 Case $\mathrm{rank}(M)=2$ subcase 2b 1.2)
      5. 2.1.5 Case $\mathrm{rank}(M)=1$ subcase 3a 1a)
      6. 2.1.6 Case $\mathrm{rank}(M)=1$ subcase 3a 1b)
      7. 2.1.7 Case $\mathrm{rank}(M)=1$ subcase 3a 2)
      8. 2.1.8 Case $\mathrm{rank}(M)=1$ subcase 3a 3)
      9. 2.1.9 Case $\mathrm{rank}(M)=1$ subcase 3a 4)
      10. 2.1.10 Case $\mathrm{rank}(M)=1$ subcase 3b)
    2. 2.2 Second order integrals
      1. 2.2.1 System with bi-quadratic potential and non-zero $B^{\phi}$
      2. 2.2.2 Systems with additional integral $L^{2}+\ldots$
      3. 2.2.3 Systems with additional integral $L_{x}p_{y}-L_{y}p_{x}+\ldots$
  4. Conclusions

## Introduction

This thesis is a contribution to the study of integrable and superintegrable
systems with electric and magnetic fields. These are Hamiltonian systems
distinguished by the existence of as many independent integrals of motion in
involution as degrees of freedom (integrability) or even some additional not
necessarily in involution with the rest (superintegrability) which give them
extraordinary properties: in classical mechanics the equations of motion can
be solved by quadratures (integrable) or even algebraically in closed form
(maximally superintegrable), in quantum mechanics the energy levels are
degenerate and it has been conjectured that all maximally superintegrable
systems are exactly solvable [26]. These properties make them invaluable as
physical models which allow us to develop insight into the principles
governing physical laws. They are also the basis for constructing more
complicated models, often by the method of perturbations. The prime example in
this regard is the periodic table [23], which is obtained as a perturbation to
the Coulomb model.

The most well-known superintegrable systems are the Kepler-Coulomb system and
the harmonic oscillator. As stated by Bertrand’s theorem [3], see also [13],
these are the only spherically symmetrical maximally superintegrable systems
without magnetic field. For analysis of these systems in quantum mechanics see
e.g. [14].

Most of the work in the field of superintegrability was done for the systems
with so-called natural Hamiltonian

$H=\frac{1}{2}\left(\vec{p}\right)^{2}+W(\vec{x}),$ (1)

mainly with the assumption that the integrals of motion are polynomial in
momenta. The best studied cases are those on the Euclidean spaces
$\mathbb{E}_{2}$ and $\mathbb{E}_{3}$, for which all second order
superintegrable systems were found [28, 12, 10]. The subsequent developments
include studies in higher dimensional Euclidean spaces, more general spaces,
e.g. Riemannian or pseudo-Riemannian, and of course higher order integrals,
see the review article [23] and references therein.

Another natural generalization, which we consider in this thesis, is the
Hamiltonian for the electromagnetic field, that is with the vector potential
$\vec{A}(\vec{x})$ in addition to the scalar potential $W(\vec{x})$, namely

$H=\frac{1}{2}\left(\vec{p}+\vec{A}(\vec{x})\right)^{2}+W(\vec{x}).$ (2)

(We consider an electron in electromagnetic field, so we choose the units of
measurement so that $e=-1$ and $m=1$.) The earlier work, which focused mainly
on the $\mathbb{E}_{2}$ case, see e.g. [5, 22, 8], was followed by the
$\mathbb{E}_{3}$ case, see [30, 20, 17, 18, 4] and references therein.

The case with the Hamiltonian (1) is related to the separation of Hamilton-
Jacobi equation in classical mechanics and Schrödinger equation in quantum
mechanics. More specifically, in [28] it was shown that in $\mathbb{E}_{3}$
there are 11 pairs of commuting integrals corresponding to 11 coordinate
systems in which Hamilton-Jacobi or Schrödinger equations separate determined
in [9]. Even though the correspondence between second order integrals and
separation no longer holds when magnetic field is present [2, 20], we still
talk about these 11 classes because the highest order terms have the same
structure, namely lie in the enveloping algebra of the Euclidean Lie algebra.

The focus of our thesis is one of the 11 coordinate systems, namely the
cylindrical with the defining transformations $x=r\cos(\phi),\ y=r\sin(\phi),\
z=Z$. If we have vanishing magnetic field and scalar potential of the form

$W(\vec{x})=W_{1}(r)+\frac{1}{r^{2}}W_{2}(\phi)+W_{3}(Z),$ (3)

the Hamilton-Jacobi equation separates in these coordinates. In the article
[11] the authors obtained all classical integrable cases with magnetic field
in this class. In this thesis we will continue this research program by
studying the quantum integrable cases and determine which of the integrable
cases admit additional first order independent integrals of motion, hence
being superintegrable. We use Maple™ [16] for the calculations and
occasionally verify some results in Mathematica® [29].

The structure of the thesis is as follows: After reviewing definitions of
integrability and superintegrability in Section 1.1, the first chapter is
focused on integrable systems in quantum mechanics. In Section 1.2 we review
classical determining equations for integrals of second order and in Section
1.3 we extend them to quantum mechanics by calculating the quantum correction
in the cylindrical coordinates with the focus on the first order and
cylindrical integrals. In Section 1.4 we find all quadratic cylindrical
integrable systems in quantum mechanics and analyse those differing from the
classical case considered in [11].

In Chapter 2 we search for additional integrals of motion to make the
integrable systems superintegrable. In Section 2.1 we find all cylindrical
systems with at least one additional first order integral and in Section 2.2
we consider second order integrals. Due to the computational complexity of the
problem we were not able to solve the second order case in general. After
analysing one system we found by chance in Subsection 2.2.1, we restrict to a
physically motivated _ansatz_ and find all cylindrical systems admitting
integrals of the form $L^{2}+\ldots$ in Subsection 2.2.2 and
$L_{x}p_{y}-L_{y}p_{x}+\ldots$ in Subsection 2.2.3. We conclude by a summary
of results.

##  Chapter 1 Quantum integrable systems in cylindrical coordinates

In this chapter, we consider quantum integrable systems in cylindrical
coordinates. In particular, we derive the determining equations in cylindrical
coordinates (Section 1.3), compare them with the classical ones from [11],
which we include in Section 1.2 for reference, and solve them for the so-
called cylindrical case, i.e. integrals with the highest order term
$L_{z}^{2}$ or $P_{z}^{2}$ (Section 1.4).

But first we review the key notions for our thesis: integrability and
superintegrability.

###  1.1 Integrability and superintegrability

We will use the standard definitions in the field, see e.g. [20], which are,
however, slightly different from those in the review article [23].

We start in classical mechanics: A finite-dimensional classical Hamiltonian
system in a $2n$-dimensional phase space is called integrable (or Liouville
integrable) if it allows $n$ integrals of motion $\\{X_{0}=H,...,X_{n-1}\\}$
(including the Hamiltonian $H$), which are

  1. 1.

well-defined analytic functions on the phase space (possibly with exceptions
of lower dimensional manifolds, such as the $z$-axis or $xy$-plane),

  2. 2.

in involution, that is pairwise Poisson commute,
$\\{X_{i},X_{j}\\}_{\text{P.B.}}=0$, where the Poisson bracket is defined in
any canonical coordinates $(q_{i},p_{i})$ as

$\\{X_{1},X_{2}\\}_{\text{P.B.}}=\sum_{i=1}^{n}\left(\frac{\partial
X_{1}}{\partial q_{i}}\frac{\partial X_{2}}{\partial p_{i}}-\frac{\partial
X_{2}}{\partial q_{i}}\frac{\partial X_{1}}{\partial p_{i}}\right),$ (1.1)

  3. 3.

and are functionally independent, i.e. their Jacobian matrix
$\left(\frac{\partial X_{i}}{\partial q_{j}},\frac{\partial X_{i}}{\partial
p_{k}}\right)$ has the maximal possible rank (here $n$) on the region where
the integrals are well defined and locally analytic.

A finite dimensional Hamiltonian system is superintegrable if it is integrable
with integrals $\\{X_{0}=H,X_{1},...,X_{n-1}\\}$ and admits additional
integrals of motion $Y_{1},\ldots,Y_{k},$ $1\leq k\leq n-1$ so that the set
$\\{X_{0}=H,X_{1},\ldots,X_{n-1},Y_{1},\ldots,Y_{k}\\}$ is functionally
independent. (Integrals $Y_{i}$ need not Poisson commute with any other
integral except the Hamiltonian.) We call the system minimally superintegrable
if $k=1$ and maximally superintegrable if $k=2n-1$. (In our case we will have
$n=3$ and maximal superintegrability means $k=2$, so there is no other
option.)

In quantum mechanics the notions in previous definitions must be slightly
modified:

  1. 1.

The integrals of motion must be well-defined operators in the enveloping
algebra of the Heisenberg Lie algebra, i.e. polynomials in phase space
coordinate operators $\hat{Q}_{i}$ and $\hat{P}_{i}$, or convergent series in
them. (Here we must work in the Cartesian coordinates because there are some
problems with quantizing the momenta in other coordinate systems.)

  2. 2.

The Poisson bracket is replaced by the commutator of operators. (We proceed
formally, ignoring complications arising from unbounded operators.)

  3. 3.

Functional independence is replaced by the so-called algebraic (or polynomial)
independence which means that no non-trivial fully symmetrized (Jordan)
polynomial of the integrals vanishes.

###  1.2 Classical determining equations

We are working with magnetic fields so let us review a few basic notions: The
magnetic field $B$, which is defined in the Cartesian coordinates as

$\vec{B}=\nabla\times\vec{A},\quad\text{i.e.}\quad
B_{j}=\varepsilon_{jkl}\frac{\partial A_{l}}{\partial x_{k}},$ (1.2)

where the Levi-Civita symbol $\varepsilon_{jkl}$ is totally antisymmetric with
$\varepsilon_{123}=1$, is gauge invariant, i.e. does not change under
transformation of the potentials

$\vec{A}^{\prime}(\vec{x})=\vec{A}(\vec{x})+\nabla\chi(\vec{x}),\quad
V^{\prime}(\vec{x})=V(\vec{x})$ (1.3)

with an arbitrary choice of the scalar function $\chi(\vec{x})$. (We consider
the time independent case only.) For the effect of gauge transformation in
quantum mechanics see eq. 1.38.

Because the systems we consider are gauge invariant, it is useful to write
their integrals of motion in terms of the so-called covariant momenta (using
units $e=-1$, $m=1$)

$p_{j}^{A}=p_{j}+A_{j}.$ (1.4)

Using these momenta, the commutation relations in the Cartesian coordinates
read

$\\{p_{j}^{A},p_{k}^{A}\\}_{\text{P.B.}}=\varepsilon_{jkl}B_{l},\quad\\{p_{j}^{A},q_{k}\\}_{\text{P.B.}}=-\delta_{jk},$
(1.5)

where $\delta_{jk}$ is the Kronecker delta.

To work with momenta and magnetic fields in the cylindrical coordinates
$(r,\phi,Z)$ defined by the transformation

$x=r\cos(\phi),\quad y=r\sin(\phi),\quad z=Z,$ (1.6)

it is convenient to use the formalism of differential forms introduced in
[21]: Given the structure of the canonical 1-form

$\lambda=p_{x}\mathrm{d}x+p_{y}\mathrm{d}y+p_{z}\mathrm{d}z=p_{r}\mathrm{d}r+p_{\phi}\mathrm{d}\phi+p_{Z}\mathrm{d}Z,$
(1.7)

we obtain the following transformation for the linear momenta

$\displaystyle p_{x}=\cos(\phi)p_{r}-\frac{\sin(\phi)}{r}p_{\phi},\quad
p_{y}=\sin(\phi)p_{r}+\frac{\cos(\phi)}{r}p_{\phi},\quad p_{z}=p_{Z},$ (1.8)

and the components of the vector potential $A$ transform the same way. (We
directly associate the Cartesian vector components to the 1-form components.)
This enables us to use the gauge invariant momenta from eq. 1.4 in the
cylindrical coordinates as well.

On the other hand, components of the magnetic field 2-form $B=\mathrm{d}A$ are

$\displaystyle\begin{aligned}
B&=B^{x}(\vec{x})\,\mathrm{d}y\wedge\mathrm{d}z+B^{y}(\vec{x})\,\mathrm{d}z\wedge\mathrm{d}x+B^{z}(\vec{x})\,\mathrm{d}x\wedge\mathrm{d}y\\\
&=B^{r}(r,\phi,Z)\,\mathrm{d}\phi\wedge\mathrm{d}Z+B^{\phi}(r,\phi,Z)\,\mathrm{d}Z\wedge\mathrm{d}r+B^{Z}(r,\phi,Z)\,\mathrm{d}r\wedge\mathrm{d}\phi.\end{aligned}$
(1.9)

This leads to the following transformation

$\displaystyle B^{x}(\vec{x})$
$\displaystyle=\frac{\cos(\phi)}{r}B^{r}(r,\phi,Z)-\sin(\phi)B^{\phi}(r,\phi,Z),$
$\displaystyle B^{y}(\vec{x})$
$\displaystyle=\frac{\sin(\phi)}{r}B^{r}(r,\phi,Z)+\cos(\phi)B^{\phi}(r,\phi,Z),$
(1.10) $\displaystyle B^{z}(\vec{x})$
$\displaystyle=\frac{1}{r}B^{Z}(r,\phi,Z),$

so that the components of the magnetic field are computed in the same way as
in the Cartesian coordinates, namely

$B^{r}=\frac{\partial A_{Z}}{\partial\phi}-\frac{\partial A_{\phi}}{\partial
Z},\quad B^{\phi}=\frac{\partial A_{r}}{\partial Z}-\frac{\partial
A_{Z}}{\partial r},\quad B^{Z}=\frac{\partial A_{\phi}}{\partial
r}-\frac{\partial A_{r}}{\partial\phi}.$ (1.11)

Now we search for integrals of motion. In classical mechanics, we use the
determining equations from [11], which were obtained by setting the Poisson
bracket $\\{X,H\\}_{\text{P.B.}}=0$ in the cylindrical coordinates (for
definition see eq. 1.1) with the Hamiltonian $H$ written in terms of the
covariant momenta from eq. 1.4,

$H=\frac{1}{2}\left(\left(p_{r}^{A}\right)^{2}+\frac{\left(p_{\phi}^{A}\right)^{2}}{r^{2}}+\left(p_{Z}^{A}\right)^{2}\right)+W(r,\phi,Z),$
(1.12)

and with a general integral of motion of the second order

$\displaystyle X={}$ $\displaystyle
h^{r}\left(r,\phi,Z\right)\left(p_{r}^{A}\right)^{2}+h^{\phi}\left(r,\phi,Z\right)\left(p_{\phi}^{A}\right)^{2}+h^{Z}\left(r,\phi,Z\right)\left(p_{Z}^{A}\right)^{2}+$
$\displaystyle+n^{r}\left(r,\phi,Z\right)p_{\phi}^{A}p_{Z}^{A}+n^{\phi}\left(r,\phi,Z\right)p_{r}^{A}p_{Z}^{A}+n^{Z}\left(r,\phi,Z\right)p_{\phi}^{A}p_{r}^{A}+$
(1.13)
$\displaystyle+s^{r}\left(r,\phi,Z\right)p_{r}^{A}+s^{\phi}\left(r,\phi,Z\right)p_{\phi}^{A}+s^{Z}\left(r,\phi,Z\right)p_{Z}^{A}+m\left(r,\phi,Z\right),$

where the functions $h,$ $n,$ $s,$ $m$ are to be determined. We separate the
coefficients of obtained powers of momenta $p_{r},$ $p_{\phi},$ $p_{Z}$. The
determining equations obtained from the third order terms are

$\displaystyle\partial_{r}h^{r}$
$\displaystyle=0,\quad\partial_{\phi}h^{r}=-r^{2}\partial_{r}n^{Z},\quad\partial_{Z}h^{r}=-\partial_{r}n^{\phi},$
$\displaystyle\partial_{r}h^{\phi}$
$\displaystyle=-\frac{1}{r^{2}}\partial_{\phi}n^{Z}-\frac{2}{r^{3}}h^{r},\quad\partial_{\phi}h^{\phi}=-\frac{1}{r}n^{Z},\quad\partial_{Z}h^{\phi}=-\frac{1}{r^{2}}\partial_{\phi}n^{r}-\frac{1}{r^{3}}n^{\phi},$
$\displaystyle\partial_{r}h^{Z}$
$\displaystyle=-\partial_{Z}n^{\phi},\quad\partial_{\phi}h^{Z}=-r^{2}\partial_{Z}n^{r},\quad\partial_{Z}h^{Z}=0,$
(1.14) $\displaystyle\partial_{\phi}n^{\phi}$
$\displaystyle=-r^{2}(\partial_{Z}n^{Z}+\partial_{r}n^{r}).$

We simplify the second order term equations using the third order section 1.2
and rewrite derivatives of $\vec{A}$ in terms of $\vec{B}$ using eq. 1.11 to
obtain

$\displaystyle\partial_{r}s^{r}$ $\displaystyle=n^{\phi}B^{\phi}-n^{Z}B^{Z},$
$\displaystyle\partial_{\phi}s^{r}$
$\displaystyle=r^{2}(n^{r}B^{\phi}-2h^{\phi}B^{Z}-\partial_{r}s^{\phi})-n^{\phi}B^{r}+2h^{r}B^{Z},$
$\displaystyle\partial_{r}s^{Z}$
$\displaystyle=n^{Z}B^{r}-\partial_{Z}s^{r}-n^{r}B^{Z}+2h^{Z}B^{\phi}-2h^{r}B^{\phi},$
$\displaystyle\partial_{\phi}s^{\phi}$
$\displaystyle=-n^{r}B^{r}+n^{Z}B^{Z}-\frac{1}{r}s^{r},$ (1.15)
$\displaystyle\partial_{\phi}s^{Z}$
$\displaystyle=r^{2}(2h^{\phi}B^{r}-n^{Z}B^{\phi}-\partial_{Z}s^{\phi})-2h^{Z}B^{r}+n^{\phi}B^{Z},$
$\displaystyle\partial_{Z}s^{Z}$ $\displaystyle=n^{r}B^{r}-n^{\phi}B^{\phi}.$

After the same simplification the first and zeroth order terms imply

$\displaystyle\partial_{r}m$
$\displaystyle=s^{Z}B^{\phi}-s^{\phi}B^{Z}+n^{\phi}\partial_{Z}W+n^{Z}\partial_{\phi}W+2h^{r}\partial_{r}W,$
$\displaystyle\partial_{\phi}m$
$\displaystyle=s^{r}B^{Z}-s^{Z}B^{r}+r^{2}(n^{r}\partial_{Z}W+2h^{\phi}\partial_{\phi}W+n^{Z}\partial_{r}W),$
(1.16) $\displaystyle\partial_{Z}m$
$\displaystyle=s^{\phi}B^{r}-s^{r}B^{\phi}+2h^{Z}\partial_{Z}W+n^{r}\partial_{\phi}W+n^{\phi}\partial_{r}W,$

and

$s^{r}\partial_{r}W+s^{\phi}\partial_{\phi}W+s^{Z}\partial_{Z}W=0,$ (1.17)

respectively.

The third order equations (1.2) can be readily integrated, because they do not
depend on the magnetic field $\vec{B}$ nor the scalar potential $W$. The
solution is therefore the same as in the vanishing magnetic field: The highest
order terms lie in the enveloping algebra of the Euclidean Lie algebra [20]
generated by $p_{x},$ $p_{y},$ $p_{z},$ $L_{x},$ $L_{y},$ $L_{z}$. For
convenience we write it in terms of gauge covariant expressions and in the
Cartesian coordinates:

$X=\sum_{1\leq a\leq b\leq
6}\alpha_{ab}Y^{A}_{a}Y^{A}_{b}+\sum_{j=1}^{3}s_{j}(\vec{x})p_{j}^{A}+m(\vec{x}),$
(1.18)

where $Y^{A}=(p_{x}^{A},p_{y}^{A},p_{z}^{A},L_{x}^{A},L_{y}^{A},L_{z}^{A})$.
The coordinates of the gauge covariant angular momentum vector
$\vec{L}^{A}=(L_{1}^{A},L_{2}^{A},L_{3}^{A})\equiv(L_{x}^{A},L_{y}^{A},L_{z}^{A})$
are defined as

$L_{j}^{A}=\sum_{k,l}\varepsilon_{jkl}x_{k}p_{l}^{A}$ (1.19)

with the covariant momenta from eq. 1.4.

Transforming the second order terms to the cylindrical coordinates using eq.
1.6 and eq. 1.8 followed by collecting the momenta $p^{r}$, $p^{\phi}$,
$p^{Z}$, we obtain

$\displaystyle\begin{split}h^{r}={}&\left[\left(\alpha_{55}-\alpha_{44}\right){Z}^{2}+\left(\alpha_{24}+\alpha_{15}\right)Z+\alpha_{11}-\alpha_{22}\right]\left(\cos(\phi)\right)^{2}-\\\
&-\left[\alpha_{45}{Z}^{2}+\left(\alpha_{14}-\alpha_{25}\right)Z-\alpha_{12}\right]\sin(\phi)\cos(\phi)+\alpha_{44}{Z}^{2}-\alpha_{24}Z+\alpha_{22},\end{split}$
(1.20)
$\displaystyle\begin{split}{h^{\phi}}={}&\frac{\left[\left(\alpha_{44}-\alpha_{55}\right){Z}^{2}-\left(\alpha_{24}+\alpha_{15}\right)Z+\alpha_{22}-\alpha_{11}\right]\left(\cos(\phi)\right)^{2}}{r^{2}}+\\\
&+\frac{\left[\left(\alpha_{45}{Z}^{2}+\left(\alpha_{14}-\alpha_{25}\right)Z-\alpha_{12}\right)\sin(\phi)-\left(\alpha_{46}Z-\alpha_{26}\right)r\right]\cos(\phi)}{r^{2}}+\\\
&+\frac{-\left(\alpha_{56}Z+\alpha_{16}\right)r\sin(\phi)+\alpha_{55}{Z}^{2}+\alpha_{66}r^{2}+\alpha_{15}Z+\alpha_{11}}{r^{2}},\end{split}$
(1.21)
$\displaystyle\begin{split}{h^{Z}}={}&r^{2}\left[\alpha_{44}-\alpha_{45}\cos(\phi)\sin(\phi)-\left(\alpha_{44}-\alpha_{55}\right)\left(\cos(\phi)\right)^{2}\right]+\\\
&\alpha_{34}r\sin(\phi)-\alpha_{35}r\cos(\phi)+\alpha_{33},\end{split}$ (1.22)
$\displaystyle\begin{split}{n^{r}}={}&\left(2\alpha_{45}Z-\alpha_{25}+\alpha_{14}\right)\left(\cos(\phi)\right)^{2}+\\\
&+\left[2\left(\alpha_{55}-\alpha_{44}\right)Zr+\alpha_{24}+\alpha_{15}\right]\sin(\phi)\cos(\phi)+\\\
&+\frac{\left(-\alpha_{56}r^{2}-\alpha_{34}Z+\alpha_{23}\right)\cos(\phi)+\left(\alpha_{46}r^{2}-\alpha_{35}Z-\alpha_{13}\right)\sin(\phi)}{r}-\\\
&-\left(\alpha_{45}Z+\alpha_{14}-\alpha_{36}\right),\end{split}$ (1.23)
$\displaystyle\begin{split}{n^{\phi}}={}&\left[2\left(\alpha_{44}-\alpha_{55}\right)Z-\alpha_{24}-\alpha_{15}\right]r\left(\cos(\phi)\right)^{2}+\\\
&+\left(2\alpha_{45}Z-\alpha_{25}+\alpha_{14}\right)r\sin(\phi)\cos(\phi)+\\\
&+\left(\alpha_{35}Z+\alpha_{13}\right)\cos(\phi)-\left(\alpha_{34}Z-\alpha_{23}\right)r\sin(\phi)-\left(2\alpha_{44}Z-\alpha_{24}\right),\end{split}$
(1.24)
$\displaystyle\begin{split}{n^{Z}}={}&-\frac{2\left[\alpha_{45}{Z}^{2}+\left(\alpha_{14}-\alpha_{25}\right)Z-\alpha_{12}\right]\left(\cos(\phi)\right)^{2}}{r}-\\\
&-\frac{2\left[\left(\alpha_{55}-\alpha_{44}\right){Z}^{2}+\left(\alpha_{24}+\alpha_{15}\right)Z-\alpha_{22}+\alpha_{11}\right]\sin(\phi)\cos(\phi)}{r}-\\\
&-\left(\alpha_{46}Z-\alpha_{26}\right)\sin(\phi)-\left(\alpha_{56}Z+\alpha_{16}\right)\cos(\phi)+\\\
&+\frac{\alpha_{45}{Z}^{2}+\left(\alpha_{14}-\alpha_{25}\right)Z-\alpha_{12}}{r}.\end{split}$
(1.25)

If we suppose that the integrals of motion are of the first order,

$X=s^{r}\left(r,\phi,Z\right)p_{r}^{A}+s^{\phi}\left(r,\phi,Z\right)p_{\phi}^{A}+s^{Z}\left(r,\phi,Z\right)p_{Z}^{A}+m\left(r,\phi,Z\right),$
(1.26)

which corresponds to setting the functions $h,$ $n$ in section 1.2 to 0, the
determining equations simplify as follows: There are no equations of the third
order (they are satisfied identically). The second order equations read

$\displaystyle\begin{aligned}
&\partial_{r}{s^{r}}=0,\quad\partial_{\phi}{s^{r}}=-r^{2}\partial_{r}{s^{\phi}},\quad\partial_{r}{s^{Z}}=-\partial_{Z}{s^{r}},\\\
&\partial_{\phi}{s^{\phi}}=-{\frac{{s^{r}}}{r}},\quad\partial_{\phi}{s^{Z}}=-r^{2}\partial_{Z}{s^{\phi}},\quad\partial_{Z}{s^{Z}}=0.\end{aligned}$
(1.27)

The first order equations read

$\displaystyle\partial_{r}m$
$\displaystyle={s^{Z}}{B^{\phi}}-{s^{\phi}}{B^{Z}},$
$\displaystyle\partial_{\phi}m$ $\displaystyle={s^{r}}{B^{Z}}-{s^{Z}}{B^{r}},$
(1.28) $\displaystyle\partial_{Z}m$
$\displaystyle={s^{\phi}}{B^{r}}-{s^{r}}{B^{\phi}},$

and the zeroth order equation is

${s^{r}}\partial_{r}W+{s^{\phi}}\partial_{\phi}W+{s^{Z}}\partial_{Z}W=0.$
(1.29)

The first thing we note is that the second order equations (1.27) do not
depend on either the magnetic field or the scalar potential in the same way as
the third order equations in the previous case. We can therefore solve them
for all the cases now. The solution reads

$\displaystyle s^{r}$
$\displaystyle=(k_{5}Z+k_{1})\cos(\phi)-(k_{4}Z-k_{2})\sin(\phi),$ (1.30)
$\displaystyle s^{\phi}$
$\displaystyle=-\frac{(k_{4}Z-k_{2})\cos(\phi)+(k_{5}Z+k_{1})\sin(\phi)-k_{6}r}{r},$
(1.31) $\displaystyle s^{Z}$
$\displaystyle=-k_{5}r\cos(\phi)+k_{4}r\sin(\phi)+k_{3},$ (1.32)

which corresponds to the integral with the first order term in the Euclidean
Lie algebra (in the Cartesian coordinates and gauge covariant form)

$Y=k_{1}p_{x}^{A}+k_{2}p_{y}^{A}+k_{3}p_{z}^{A}+k_{4}L_{x}^{A}+k_{5}L_{y}^{A}+k_{6}L_{z}^{A}+m(x,y,z).$
(1.33)

###  1.3 Determining equations in quantum mechanics

Let us start with some basic notions regarding the considered systems in
quantum mechanics.

We begin in the Cartesian coordinates. We work in the Schrödinger
representation, where the coordinates and linear momenta are

$\hat{Q}_{j}=x_{j},\quad\hat{P}_{j}=-\mathrm{i}\hbar\partial_{j}\equiv-\mathrm{i}\hbar\nabla_{j},$
(1.34)

respectively. (The operator $\hat{Q}$ is a multiplication by the corresponding
coordinate.)

The convention for operator ordering in the field of integrability is to
symmetrize properly [23], therefore the Hamiltonian for a system with magnetic
field in the Cartesian coordinates reads

$\displaystyle\hat{H}$
$\displaystyle=\frac{1}{2}\sum_{j}\left(\hat{P}_{j}\hat{P}_{j}+\hat{P}_{j}\hat{A}_{j}(\vec{x})+\hat{A}_{j}(\vec{x})\hat{P}_{j}+\hat{A}_{j}(\vec{x})\hat{A}_{j}(\vec{x})\right)+\hat{V}(\vec{x})=$
$\displaystyle=-\frac{\hbar^{2}}{2}\Delta-\mathrm{i}\hbar\hat{\vec{A}}\cdot\vec{\nabla}-\frac{\mathrm{i}\hbar}{2}\widehat{\textrm{div}\vec{A}}+\frac{1}{2}\hat{\vec{A}}\cdot\hat{\vec{A}}+\hat{V}(\vec{x}),$
(1.35)

where $\hat{A}_{j}(\vec{x})$, $\widehat{\textrm{div}\vec{A}}$ and
$\hat{V}(\vec{x})$ are the operators of multiplication by the $j$-th
coordinate of the vector potential $\vec{A}(\vec{x})$, the divergence of
$\vec{A}(\vec{x})$ and the scalar potential $V(\vec{x})$, respectively.

We follow the same convention for the integrals of motion as well, namely we
write the general second order integral as

$\hat{X}=\sum_{j=1}^{3}\\{h_{j}(\vec{x}),\hat{P}_{j}^{A}\hat{P}_{j}^{A}\\}+\sum_{j,k,l=1}^{3}\frac{|\varepsilon_{jkl}|}{2}\\{n_{j}(\vec{x}),\hat{P}_{k}^{A}\hat{P}_{l}^{A}\\}+\sum_{j=1}^{3}\\{s_{j}(\vec{x}),\hat{P}_{j}^{A}\\}+m(\vec{x}),$
(1.36)

where $\\{,\\}$ denotes the symmetrization

$\\{\hat{a},\hat{b}\\}=\frac{1}{2}\left(\hat{a}\hat{b}+\hat{b}\hat{a}\right).$
(1.37)

All choices of symmetrization are equivalent up to redefinition of the lower
order terms [20], so the choice is without loss of generality.

In quantum mechanics the gauge transformation (1.3) demonstrates itself as a
unitary transformation $\hat{U}$ of the underlying Hilbert space. Defining
$\hat{U}$ as

$\hat{U}\psi(\vec{x})=\exp\left(\frac{\mathrm{i}}{\hbar}\chi(\vec{x})\right)\cdot\psi(\vec{x}),$
(1.38)

the states and observables describing the same physical system are

$\psi\to\psi^{\prime}=\hat{U}\psi,\quad\hat{O}\to\hat{O}^{\prime}=\hat{U}\hat{O}\hat{U}^{\dagger}.$
(1.39)

The potentials and momenta transform the following way

$(\hat{P}_{j}+\hat{A}_{j})\to\hat{U}(\hat{P}_{j}+\hat{A}_{j})\hat{U}^{\dagger}=\hat{P}_{j}+\hat{A}^{\prime}_{j},\quad\hat{V}=\hat{U}\hat{V}\hat{U}^{\dagger}.$
(1.40)

We can therefore use the gauge covariant momenta

$\hat{P}_{j}^{A}=\hat{P}_{j}+\hat{A}_{j}$ (1.41)

in quantum mechanics as well. Their commutation relations in the Cartesian
coordinates read

$[\hat{P}_{j}^{A},\hat{P}_{k}^{A}]=-\mathrm{i}\hbar\varepsilon_{jkl}B_{l},\quad[\hat{P}_{j}^{A},\hat{Q}_{k}]=-\mathrm{i}\hbar\delta_{jk}.$
(1.42)

Now let us turn to the cylindrical coordinates $(r,\phi,Z)$, see eq. 1.6.
Quantization in non-Cartesian coordinates poses a problem, because the Dirac
quantization rule

$\\{A,B\\}_{\text{P.B}}\to\frac{1}{\mathrm{i}\hbar}[\hat{A},\hat{B}]$ (1.43)

is not invariant under canonical transformations (not even coordinate
transformations) in the sense that Schrödinger equations obtained from
different canonical coordinates are inequivalent [25, Exercise 7.4.10].
Experiments show that the correct Schrödinger equation is obtained by
quantization of the classical Hamiltonian in the Cartesian coordinates and
subsequent transformation into the other coordinate systems. Succinctly, the
rule is to quantize and transform (in that order).

Let us therefore transform the quantized Hamiltonian (1.35). The gradient
operator in the cylindrical coordinates reads

$\nabla f=\left(\frac{\partial f}{\partial r},\frac{1}{r}\frac{\partial
f}{\partial\phi},\frac{\partial f}{\partial Z}\right)$ (1.44)

and the Laplace operator is

$\Delta=\frac{1}{r}\frac{\partial}{\partial r}\left(r\frac{\partial}{\partial
r}\right)+\frac{1}{r^{2}}\frac{\partial^{2}}{\partial\phi^{2}}+\frac{\partial^{2}}{\partial
Z^{2}}.$ (1.45)

The divergence is

$\vec{\nabla}\cdot\vec{A}=\frac{1}{r}\frac{\partial(rA_{r})}{\partial
r}+\frac{\partial A_{\phi}}{\partial\phi}+\frac{\partial A_{Z}}{\partial Z},$
(1.46)

(the right hand side in term of 1-form components) so the final form of the
Hamiltonian in the cylindrical coordinates reads ($A_{i}$ and $V$ functions of
$r,$ $\phi,$ $Z$)

$\displaystyle\hat{H}={}$
$\displaystyle-\frac{\hbar^{2}}{2}\left(\frac{1}{r}\frac{\partial}{\partial
r}\left(r\frac{\partial}{\partial
r}\right)+\frac{1}{r^{2}}\frac{\partial^{2}}{\partial\phi^{2}}+\frac{\partial^{2}}{\partial
Z^{2}}\right)-\mathrm{i}\hbar\left(\hat{A}_{r}\frac{\partial}{\partial
r}+\frac{\hat{A}_{\phi}}{r^{2}}\frac{\partial}{\partial\phi}+\hat{A}_{Z}\frac{\partial}{\partial
Z}\right)-$
$\displaystyle-\frac{\mathrm{i}\hbar}{2}\left(\frac{1}{r}\frac{\partial(rA_{r})}{\partial
r}+\frac{\partial A_{\phi}}{\partial\phi}+\frac{\partial A_{Z}}{\partial
Z}\right)+\frac{1}{2}\left((\hat{A}_{r})^{2}+\frac{(\hat{A}_{\phi})^{2}}{r^{2}}+(\hat{A}_{Z})^{2}\right)+\hat{V}.$
(1.47)

The same procedure must be applied to the integrals of motion as well, we thus
quantize in the Cartesian coordinates, symmetrize and transform into the
cylindrical coordinates. For this we need the transformation rules for momenta
form eq. 1.8, which we readily invert to get

$p_{r}=\cos(\phi)p_{x}+\sin(\phi)p_{y},\quad
p_{\phi}=-r\sin(\phi)p_{x}+r\cos(\phi)p_{y},\quad p_{Z}=p_{z}$ (1.48)

and the momentum operators $\hat{P}_{i}$ in the cylindrical coordinates read

$\hat{P}_{x}=-\mathrm{i}\hbar\left(\cos(\phi)\partial_{r}-\frac{\sin(\phi)}{r}\partial_{\phi}\right),\
\hat{P}_{y}=-\mathrm{i}\hbar\left(\sin(\phi)\partial_{r}+\frac{\cos(\phi)}{r}\partial_{\phi}\right),\
\hat{P}_{z}=-\mathrm{i}\hbar\partial_{Z}.$ (1.49)

The correct way to obtain the determining equations in the cylindrical
coordinates is to transform the Cartesian equations from [20]. Let us list
them for reference:

The third order equations read

$\displaystyle\partial_{x}h_{x}=0,\quad\partial_{y}h_{x}=-\partial_{x}n_{x},\quad\partial_{z}h_{x}=-\partial_{x}n_{y},$
$\displaystyle\partial_{x}h_{y}=-\partial_{y}n_{z},\quad\partial_{y}h_{y}=0,\quad\partial_{z}h_{y}=-\partial_{y}n_{x},$
(1.50)
$\displaystyle\partial_{x}h_{z}=-\partial_{z}n_{z},\quad\partial_{y}h_{z}=-\partial_{z}n_{x},\quad\partial_{z}h_{z}=0,$
$\displaystyle\nabla\cdot\vec{n}=0.$

The second order equations are

$\displaystyle\partial_{x}s_{x}=n_{y}B^{y}-n_{z}B^{z},$
$\displaystyle\partial_{y}s_{y}=n_{z}B^{z}-n_{x}B^{x},$
$\displaystyle\partial_{z}s_{z}=n_{x}B^{x}-n_{y}B^{y},$
$\displaystyle\partial_{y}s_{x}+\partial_{x}s_{y}=n_{x}B^{y}-n_{y}B^{y}+2(h_{x}-h_{y})B^{z},$
(1.51)
$\displaystyle\partial_{z}s_{x}+\partial_{x}s_{z}=n_{z}B^{x}-n_{x}B^{z}+2(h_{z}-h_{x})B^{y},$
$\displaystyle\partial_{y}s_{z}+\partial_{z}s_{y}=n_{y}B^{z}-n_{z}B^{y}+2(h_{y}-h_{z})B^{x}$

with the following consequence

$\nabla\cdot\vec{s}=0.$ (1.52)

The first order

$\displaystyle\partial_{x}m=2h_{x}\partial_{x}W+n_{z}\partial_{y}W+n_{y}\partial_{z}W+s_{z}B^{y}-s_{y}B^{z},$
$\displaystyle\partial_{y}m=2h_{y}\partial_{y}W+n_{z}\partial_{x}W+n_{x}\partial_{z}W+s_{x}B^{z}-s_{z}B^{x},$
(1.53)
$\displaystyle\partial_{z}m=2h_{z}\partial_{z}W+n_{y}\partial_{x}W+n_{x}\partial_{y}W+s_{y}B^{x}-s_{x}B^{y}.$

The zeroth order equation containing the $\hbar^{2}$–proportional quantum
correction reads

$\begin{split}\vec{s}\cdot\nabla
W+\frac{\hbar^{2}}{4}\left(\partial_{z}n_{x}\partial_{z}B^{x}-\partial_{y}n_{x}\partial_{y}B^{x}+\partial_{x}n_{y}\partial_{x}B^{y}-\partial_{z}n_{y}\partial_{z}B^{y}\right.+\\\
+\left.\partial_{y}n_{z}\partial_{y}B^{z}-\partial_{x}n_{z}\partial_{x}B^{z}+\partial_{x}n_{x}\partial_{y}B^{y}-\partial_{y}n_{y}\partial_{x}B^{x}\right)=0.\end{split}$
(1.54)

Although the last two terms suggest that the quantum correction is not
invariant under Euclidean transformations, which are symmetries the system
possesses clearly, it is not the case due to the identity [20]

$\partial_{x}n_{x}\partial_{y}B^{y}-\partial_{y}n_{y}\partial_{x}B^{x}=\partial_{y}n_{y}\partial_{z}B^{z}-\partial_{z}n_{z}\partial_{y}B^{y}=\partial_{z}n_{z}\partial_{x}B^{x}-\partial_{x}n_{x}\partial_{z}B^{z},$
(1.55)

which follows from

$\nabla\cdot\vec{B}=0,\quad\nabla\cdot\vec{n}=0.$ (1.56)

The last missing piece needed to perform the transformation is

$\begin{split}h_{x}&=h^{r}(\cos(\phi))^{2}+h^{\phi}r^{2}(\sin(\phi))^{2}-n^{Z}r\cos(\phi)\sin(\phi),\\\
h_{y}&=h^{r}(\sin(\phi))^{2}+h^{\phi}r^{2}(\cos(\phi))^{2}+n^{Z}r\cos(\phi)\sin(\phi),\\\
h_{z}&=h^{Z},\\\ n_{x}&=n^{r}r\cos(\phi)+n^{\phi}\sin(\phi),\\\
n_{y}&=-n^{r}r\sin(\phi)+n^{\phi}\cos(\phi),\\\
n_{z}&=2h^{r}\sin(\phi)\cos(\phi)-2h^{\phi}r^{2}\sin(\phi)\cos(\phi)+n^{Z}r[(\cos(\phi))^{2}-(\sin(\phi))^{2}],\\\
s_{x}&=s^{r}\cos(\phi)-s^{\phi}r\sin(\phi),\\\
s_{x}&=s^{r}\sin(\phi)+s^{\phi}r\cos(\phi),\\\ s_{z}&=s^{Z}.\end{split}$
(1.57)

It is obtained (in classical mechanics) by transforming the linear momenta
using eq. 1.8 in the second order integral of motion $X$, see section 1.2, and
the definition of the functions $n,$ $h,$ $s$ as the coefficients of the
corresponding momenta ($h_{x}$ is the coefficient of $p_{x}^{2}$, $n_{x}$ of
$p_{y}p_{z}$, $s_{x}$ of $p_{x}$ etc. and analogously in the cylindrical
coordinates).

Due to the transformation rules (1.57), the third, second and first order
equations (1.3), (1.3) and (1.3), which do not have any quantum correction,
get the same form as in classical mechanics, namely the form of equations
(1.2), (1.2) and (1.2), respectively. The transformed zeroth order equation
with derivatives of the functions $h$ eliminated using the third order
equations (1.2) reads

$\displaystyle\begin{split}&s^{r}\partial_{r}W+s^{\phi}\partial_{\phi}W+s^{Z}\partial_{Z}W+\\\
&+\frac{\hbar^{2}}{4}\left[-{\frac{\partial_{Z}n^{Z}B^{r}}{r}}+{\frac{\left({r}\partial_{\phi}n^{r}+n^{\phi}\right)B^{\phi}}{r^{2}}}+{\frac{\left(r\partial_{r}n^{Z}+n^{Z}\right)B^{Z}}{r^{2}}}-\right.\\\
&-{\frac{\left(r\partial_{\phi}n^{r}+n^{\phi}\right)\partial_{\phi}B^{r}}{r^{3}}}+\frac{\left(r\partial_{Z}n^{Z}+r\partial_{r}n^{r}+n^{r}\right)\partial_{r}B^{r}}{r}+\left(\partial_{Z}n^{r}\right)\partial_{Z}B^{r}-\\\
&-\left(\partial_{Z}n^{\phi}\right)\partial_{Z}B^{\phi}+\frac{\left(r\partial_{r}n^{r}+n^{r}\right)\partial_{\phi}B^{\phi}}{r}+\left(\partial_{r}n^{\phi}\right)\partial_{r}B^{\phi}+\\\
&+\left.{\frac{\left(r\partial_{\phi}n^{Z}+2h^{r}-2r^{2}h^{\phi}\right)\partial_{\phi}B^{Z}}{r^{3}}}-{\frac{\left(r\partial_{r}n^{Z}+n^{Z}\right)\partial_{r}B^{Z}}{r}}\right]=0.\end{split}$
(1.58)

If we use the identity eq. 1.55 to change the last 2 terms in the quantum
correction in eq. 1.54, we find that only the terms with diagonal derivatives
such as $\partial_{r}B^{r}$ change. By taking $\tfrac{1}{3}$ of each term in
eq. 1.55 instead of the last 2 terms in eq. 1.54 and using the identity

$\partial_{r}B^{r}+\partial_{\phi}B^{\phi}+\partial_{Z}B^{Z}=0,$ (1.59)

which follows from exactness of the 2-form $B={\operatorname{d}}A$, we
simplify the diagonal terms

$\frac{\left(r\partial_{Z}n^{Z}+r\partial_{r}n^{r}+n^{r}\right)\partial_{r}B^{r}}{r}+\frac{\left(r\partial_{r}n^{r}+n^{r}\right)\partial_{\phi}B^{\phi}}{r}=(\partial_{Z}n^{Z})\partial_{r}B^{r}-\frac{(r\partial_{r}n^{r}+n^{r})\partial_{Z}B^{Z}}{r}.$
(1.60)

We note that brackets of the type $\left(r\partial_{r}n^{Z}+n^{Z}\right)$
vanish if $n^{Z}\sim\frac{1}{r}$. This means vanishing quantum correction for
the integrals of motion with the highest order term $P_{i}P_{j},\
i,j\in\\{x,y,z\\}$, taking into account eq. 1.57. (The result can be clearly
seen from eq. 1.54, so this merely checks the transformation.)

For comparison, we include the quantum correction with collected functions $n$
and their derivatives. We simplified it using the identities (1.59) and
(1.60).

$\displaystyle\begin{split}&\frac{\hbar^{2}}{4}\left[{\frac{\left(B^{Z}-r\partial_{r}B^{Z}\right)(n^{Z}+r\partial_{r}n_{Z})}{r^{2}}}+{\frac{\left(rB^{\phi}-\partial_{\phi}B^{r}\right)(n^{\phi}+r\partial_{\phi}n^{r})}{r^{3}}}-\right.\\\
&-\frac{\left(\partial_{Z}B^{Z}\right)(n^{r}+r\partial_{r}n^{r})}{r}-\frac{\left(B^{r}-r\partial_{r}B^{r}\right)\partial_{Z}n^{Z}}{r}+{\frac{\left(\partial_{\phi}B^{Z}\right)\partial_{\phi}n^{Z}}{r^{2}}}-\\\
&-(\partial_{Z}B^{\phi})\partial_{Z}n^{\phi}+\left(\partial_{r}B^{\phi}\right)\partial_{r}n^{\phi}+(\partial_{Z}B^{r})\partial_{Z}n^{r}\left.+\frac{2\partial_{\phi}B^{Z}\left(h^{r}-r^{2}h^{\phi}\right)}{r^{3}}\right].\end{split}$
(1.61)

Let us have a look at the effects of the correction (LABEL:corr_cyl) on the
cylindrical integrable systems. From eq. 1.54 follows that the determining
equations for the first order integrals are the same as in the classical case
($n^{i}=0$, $i\in\\{x,y,z\\}$). As we have already mentioned, the correction
vanishes for integrals of motion with the highest order term
$\hat{P}_{i}\hat{P}_{j},$ $i,j\in\\{x,y,z\\},$ whose functions $n^{i}$ are
constant. This includes our integral $\hat{P}_{Z}^{2}+\ldots$ On the other
hand, the integral $\hat{P}_{\phi}^{2}+\ldots$ with $h^{\phi}=1$ corresponds
to

$h^{x}=y^{2},\quad h^{y}=x^{2},\quad h^{z}=0,\quad n^{z}=-{2xy},\quad
n^{x}=0,\quad n^{y}=0$ (1.62)

and the quantum correction reads

$-\frac{\hbar^{2}}{2r}\partial_{\phi}B^{Z}.$ (1.63)

(It can be seen from LABEL:corr_cyl directly.) The correction is non-zero for
the general form of $B^{Z}$ obtained from the highest order determining
equations [11],

$B^{Z}=-\frac{1}{2r^{2}}(\psi(\phi)+\psi^{\prime\prime}(\phi))+r\mu(Z)-\frac{1}{2}\rho^{\prime}(r),$
(1.64)

see eq. 1.88 below.

There could be another quantum correction arising from
$[\hat{X}_{1},\hat{X}_{2}]=0$, which is the quantum analogue of involutivity
condition. Let us show that is not the case:

We first assume that both integrals of motion are of the first order, that is

$\displaystyle\hat{X}_{1}=\sum_{j=1}^{3}\frac{1}{2}\left(s_{j}(x,y,z)\hat{P}_{j}^{A}+\hat{P}_{j}^{A}s_{j}(x,y,z)\right)+m(x,y,z),$
(1.65)
$\displaystyle\hat{X}_{2}=\sum_{k=1}^{3}\frac{1}{2}\left(\tilde{s}_{k}(x,y,z)\hat{P}_{k}^{A}+\hat{P}_{k}^{A}\tilde{s}_{k}(x,y,z)\right)+\tilde{m}(x,y,z).$
(1.66)

After some calculations using the commutation relations from eq. (1.42), we
obtain

$\displaystyle\begin{split}[X_{1},X_{2}]={}&\sum_{j,k}\mathrm{i}\hbar\left(\tilde{s}_{k}\partial_{k}s_{j}\hat{P}_{j}^{A}-s_{j}\partial_{j}\tilde{s}_{k}\hat{P}_{k}^{A}-\varepsilon_{jkl}s_{j}\tilde{s}_{k}B^{l}-s_{j}\partial_{j}\tilde{m}+\tilde{s}_{k}\partial_{k}m\right)+\\\
&+\sum_{j,k}\frac{\hbar^{2}}{2}(\tilde{s}_{k}\partial_{jk}s_{j}-s_{j}\partial_{jk}\tilde{s}_{k}).\end{split}$
(1.67)

The first line corresponds to the Poisson bracket, the second line is an
apparent quantum correction. Collecting the coefficients of $\hat{P}_{j}^{A}$
(renaming the dummy indices in the second term), we obtain the first order
equations

$\sum_{k=1}^{3}(\tilde{s}_{k}\partial_{k}s_{j}-s_{k}\partial_{k}\tilde{s}_{j})=0.$
(1.68)

The quantum correction can be rewritten as follows

$\sum_{j,k}\frac{\hbar^{2}}{2}\partial_{j}(\tilde{s}_{k}\partial_{k}s_{j}-s_{k}\partial_{k}\tilde{s}_{j}),$
(1.69)

because the terms with the first order derivatives cancel if we rename the
indices. The bracket contains the first order equations; therefore, the
quantum correction vanishes.

Now let us turn to the cylindrical case, i.e. the integrals of motion of the
form

$\displaystyle\hat{X}_{1}=\left(\hat{L}_{z}^{A}\right)^{2}+\sum_{j=1}^{3}\frac{1}{2}\left(s_{j}(x,y,z)\hat{P}_{j}^{A}+\hat{P}_{j}^{A}s_{j}(x,y,z)\right)+m(x,y,z),$
(1.70)
$\displaystyle\hat{X}_{2}=\left(\hat{P}_{z}^{A}\right)^{2}+\sum_{k=1}^{3}\frac{1}{2}\left(\tilde{s}_{k}(x,y,z)\hat{P}_{k}^{A}+\hat{P}_{k}^{A}\tilde{s}_{k}(x,y,z)\right)+\tilde{m}(x,y,z).$
(1.71)

We impose the relation $[\hat{X}_{1},\hat{X}_{2}]=0$ and separate the equation
according to the powers of $P$. The second order terms

$\displaystyle\hat{P}_{z}\hat{P}_{z}$ $\displaystyle:\partial_{z}s_{z}$
$\displaystyle=0,$ $\displaystyle\hat{P}_{z}\hat{P}_{x}$
$\displaystyle:y(x\partial_{y}-y\partial_{x})\tilde{s}_{z}+2y(yB_{x}+yB_{y})+\partial_{z}s_{x}$
$\displaystyle=0,$ $\displaystyle\hat{P}_{z}\hat{P}_{y}$
$\displaystyle:x(x\partial_{y}-y\partial_{x})\tilde{s}_{z}+2x(yB_{x}+yB_{y})-\partial_{z}s_{y}$
$\displaystyle=0,$ (1.72) $\displaystyle\hat{P}_{x}\hat{P}_{x}$
$\displaystyle:y\tilde{s}_{y}+y(x\partial_{y}-y\partial_{x})\tilde{s}_{x}$
$\displaystyle=0,$ $\displaystyle\hat{P}_{y}\hat{P}_{y}$
$\displaystyle:x\tilde{s}_{x}-y(x\partial_{y}-y\partial_{x})\tilde{s}_{y}$
$\displaystyle=0,$

give the same equations as would be obtained in the classical case. (The
equation for $\hat{P}_{x}\hat{P}_{y}$ vanishes as a consequence of those for
$\hat{P}_{x}\hat{P}_{x}$ and $\hat{P}_{y}\hat{P}_{y}$.)

The first order equations obtain some rather complicated apparent corrections,
which nevertheless vanish when we use the second order equations (1.3). The
equations (without the apparent corrections) read

$\displaystyle
2y[(x\partial_{y}-y\partial_{x})\tilde{m}+(xB_{x}+yB_{y})\tilde{s}_{z}-(x\tilde{s}_{x}+y\tilde{s}_{y})B_{z}]+\sum_{j=1}^{3}(\tilde{s}_{j}\partial_{j}s_{x}-\tilde{s}_{j}\partial_{j}s_{x})$
$\displaystyle=0,$ $\displaystyle
2x[(x\partial_{y}-y\partial_{x})\tilde{m}+(xB_{x}+yB_{y})\tilde{s}_{z}-(x\tilde{s}_{x}+y\tilde{s}_{y})B_{z}]+\sum_{j=1}^{3}(\tilde{s}_{j}\partial_{j}s_{y}-\tilde{s}_{j}\partial_{j}s_{y})$
$\displaystyle=0,$ (1.73) $\displaystyle
2\partial_{z}m-2(s_{y}B_{x}-s_{x}B_{y})+\sum_{j=1}^{3}(\tilde{s}_{j}\partial_{j}s_{z}-s_{j}\partial_{j}\tilde{s}_{z})$
$\displaystyle=0.$

The same happens in the zeroth order: All apparent corrections vanish when we
use the higher order conditions (1.3) and (1.3), so the equation has the
classical form

$\tilde{s}_{x}(s_{y}B_{z}-s_{z}B_{y})+\tilde{s}_{y}(s_{z}B_{x}-s_{x}B_{z})+\tilde{s}_{z}(s_{y}B_{x}-s_{x}B_{y})+\sum_{j=1}^{3}(\tilde{s}_{j}\partial_{j}m-s_{j}\partial_{j}\tilde{m})=0.$
(1.74)

To sum up: The apparent quantum corrections to the involutivity condition
$[\hat{X}_{1},\hat{X}_{2}]=0$ vanish once we impose the higher order
conditions from the same commutator for both the first order and the
cylindrical integrals. We note that we do not need to use the determining
equations for the integrals to prove the assertion.

In the cylindrical case there is, therefore, only one equation with quantum
correction, namely the zeroth order equation for $X_{1}$, which now reads

$s_{1}^{r}\partial_{r}W+s_{1}^{\phi}\partial_{\phi}W+s_{1}^{Z}\partial_{Z}W-\frac{\hbar^{2}}{2r}\partial_{\phi}B^{Z}=0.$
(1.75)

In order to consider the quantum correction for integrals containing angular
momenta in the second order term, we suggest using the Cartesian form of the
correction from eq. 1.54 or transforming it into more suitable coordinates
because using the correction from LABEL:corr_cyl seems not tractable.

###  1.4 Quantum integrable systems

Here we solve the quantum determining equations for the case of cylindrical
integrals in the cylindrical coordinates. As we have seen in the last section,
the quantum correction has arisen in one equation only, namely the zeroth
order equation for the integral $\hat{X}_{1}=(\hat{P}_{\phi}^{A})^{2}+\ldots$,
see eq. (1.75). Because there are so little changes, we will closely follow
the analysis from [11] and focus mainly on the cases which are affected by the
correction.

We start with reducing the determining equations (1.2)–(1.17) to the case of
cylindrical integrals, which in classical mechanics read

$\begin{split}&X_{1}=(p_{\phi}^{A})^{2}+s_{1}^{r}(r,\phi,Z)p_{r}^{A}+s_{1}^{\phi}(r,\phi,Z)p_{\phi}^{A}+s_{1}^{Z}(r,\phi,Z)p_{Z}^{A}+m_{1}(r,\phi,Z),\\\
&X_{2}=(p_{Z}^{A})^{2}+s_{2}^{r}(r,\phi,Z)p_{r}^{A}+s_{2}^{\phi}(r,\phi,Z)p_{\phi}^{A}+s_{2}^{Z}(r,\phi,Z)p_{Z}^{A}+m_{2}(r,\phi,Z).\end{split}$
(1.76)

To obtain their quantum mechanical form, we would have to transform the
integrals into the Cartesian coordinates, quantize with proper symmetrization
and transform back into the cylindrical coordinates. Because we have
transformed the determining equations from the Cartesian form, we do not need
the explicit form of the integrals here.

The reduction is done by substituting the appropriate values for the functions
$h,$ $n$, namely the only non-zero functions

$h^{\phi}_{1}=1,\quad h_{2}^{Z}=1,$ (1.77)

into the general determining equations (1.2)–(1.17) and taking into account
the quantum correction (1.63) in the zeroth order equation. We obtain the
following:

The third order equations (1.2) are trivially satisfied for both integrals.
The second order equations for the integrals $\hat{X}_{1}$ and $\hat{X}_{2}$
read

$\begin{split}\partial_{r}s_{1}^{r}&=0,\quad\partial_{\phi}s_{1}^{\phi}=-\frac{s_{1}^{r}}{r},\\\
\partial_{\phi}s_{1}^{r}&=-r^{2}(\partial_{r}s_{1}^{\phi}+2B^{Z}),\quad\partial_{\phi}s_{1}^{Z}=-r^{2}(\partial_{Z}s_{1}^{\phi}-2B^{r}),\\\
\partial_{r}s_{1}^{Z}&=-\partial_{Z}s_{1}^{r},\quad\partial_{Z}s_{1}^{Z}=0,\end{split}$
(1.78)

and

$\begin{split}\partial_{r}s_{2}^{r}&=0,\quad\partial_{\phi}s_{2}^{\phi}=-\frac{s_{2}^{r}}{r},\\\
\partial_{\phi}s_{2}^{r}&=-r^{2}\partial_{r}s_{2}^{\phi},\quad\partial_{\phi}s_{2}^{Z}=-r^{2}\partial_{Z}s_{2}^{\phi}-2B^{r},\\\
\partial_{r}s_{2}^{Z}&=-\partial_{Z}s_{2}^{r}+2B^{\phi},\quad\partial_{Z}s_{2}^{Z}=0,\end{split}$
(1.79)

respectively.

The first order equations reduce to those for $\hat{X}_{1}$

$\displaystyle\partial_{r}m_{1}$
$\displaystyle=s_{1}^{Z}B^{\phi}-s_{1}^{\phi}B^{Z},$
$\displaystyle\partial_{\phi}m_{1}$
$\displaystyle=s_{1}^{r}B^{Z}-s_{1}^{Z}B^{r}+2r^{2}\partial_{\phi}W,$ (1.80)
$\displaystyle\partial_{Z}m_{1}$
$\displaystyle=s_{1}^{\phi}B^{r}-s_{1}^{r}B^{\phi},$

and those for $\hat{X}_{2}$

$\displaystyle\partial_{r}m_{2}$
$\displaystyle=s_{2}^{Z}B^{\phi}-s_{2}^{\phi}B^{Z},$
$\displaystyle\partial_{\phi}m_{2}$
$\displaystyle=s_{2}^{r}B^{Z}-s_{2}^{Z}B^{r},$ (1.81)
$\displaystyle\partial_{Z}m_{2}$
$\displaystyle=s_{2}^{\phi}B^{r}-s_{2}^{r}B^{\phi}+2\partial_{Z}W.$

One of the zeroth order equations is the only one to change with respect to
the classical case:

$\displaystyle
s_{1}^{r}\partial_{r}W+s_{1}^{\phi}\partial_{\phi}W+s_{1}^{Z}\partial_{Z}W-\frac{\hbar^{2}}{2r}\partial_{\phi}B^{Z}=0,$
(1.82) $\displaystyle
s_{2}^{r}\partial_{r}W+s_{2}^{\phi}\partial_{\phi}W+s_{2}^{Z}\partial_{Z}W=0.$
(1.83)

In addition to the commutation with the Hamiltonian, i.e. the condition on
integrals of motion, we impose the commutation of the integrals
$[\hat{X}_{1},\hat{X}_{2}]=0$, which corresponds to the classical notion of
involution. We obtain additional equations for each order in momenta, namely
those of the second order

$\displaystyle\partial_{\phi}s_{2}^{\phi}$
$\displaystyle=0,\quad\partial_{\phi}s_{2}^{r}=0,\quad\partial_{Z}s_{1}^{r}=0,\quad\partial_{\phi}s_{2}^{Z}=\partial_{Z}s_{1}^{\phi}-2B^{r},$
(1.84)

of the first order

$\displaystyle
s_{2}^{Z}\partial_{Z}s_{1}^{r}+s_{2}^{\phi}\partial_{\phi}s_{1}^{r}$
$\displaystyle=0,$ $\displaystyle-
s_{1}^{\phi}(2B^{r}+\partial_{\phi}s_{2}^{Z})+s_{2}^{Z}\partial_{Z}s_{1}^{Z}-s_{1}^{Z}\partial_{Z}s_{2}^{Z}+$
$\displaystyle+s_{2}^{\phi}\partial_{\phi}2_{1}^{Z}+s_{1}^{r}(2B^{\phi}-\partial_{r}s_{2}^{Z})+2\partial_{Z}m_{1}$
$\displaystyle=0,$ (1.85) $\displaystyle-
s_{2}^{Z}(2B^{r}-\partial_{Z}s_{1}^{\phi})+s_{2}^{\phi}\partial_{\phi}s_{1}^{\phi}-$
$\displaystyle-
s_{1}^{Z}\partial_{Z}s_{2}^{\phi}-s_{1}^{r}\partial_{r}s_{2}^{\phi}-2\partial_{\phi}m_{2}$
$\displaystyle=0,$

and of the zeroth order

$\displaystyle-
s_{1}^{r}\partial_{r}m_{2}+s_{2}^{\phi}\partial_{\phi}m_{1}-s_{1}^{\phi}\partial_{\phi}m_{2}+s_{2}^{Z}\partial_{Z}m_{1}-s_{1}^{Z}\partial_{Z}m_{2}+$
$\displaystyle+B^{r}(s_{2}^{\phi}s_{1}^{Z}-s_{1}^{\phi}s_{2}^{Z})+B^{\phi}s_{1}^{r}s_{2}^{Z}-B^{Z}s_{1}^{r}s_{2}^{\phi}$
$\displaystyle=0.$ (1.86)

None of them has any quantum correction, as was considered at the end of
Section 1.3 (in the Cartesian coordinates).

The second order equations (1.78), (1.79) and (1.84) can be solved for the
functions $s_{j}$ and the magnetic field $B$ in terms of 5 functions of one
variable each, which we call the auxiliary functions:

$\displaystyle\begin{split}s_{1}^{r}&=\frac{\mathrm{d}}{\mathrm{d}\phi}\psi(\phi),\quad
s_{1}^{\phi}=-\frac{\psi(\phi)}{r}-r^{2}\mu(Z)+\rho(r),\quad
s_{1}^{Z}=\tau(\phi),\\\ s_{2}^{r}&=0,\quad s_{2}^{\phi}=\mu(Z),\quad
s_{2}^{Z}=-\frac{\tau(\phi)}{r^{2}}+\sigma(r),\end{split}$ (1.87)
$\displaystyle\begin{split}B^{r}&=-\frac{r^{2}}{2}\frac{\mathrm{d}}{\mathrm{d}Z}\mu(Z)+\frac{1}{2r^{2}}\frac{\mathrm{d}}{\mathrm{d}\phi}\tau(\phi),\quad
B^{\phi}=\frac{\tau(\phi)}{r^{3}}+\frac{1}{2}\frac{\mathrm{d}}{\mathrm{d}r}\sigma(r),\\\
B^{Z}&=\frac{-\psi(\phi)}{2r^{2}}+r\mu(Z)-\frac{1}{2}\frac{\mathrm{d}}{\mathrm{d}r}\rho(r)-\frac{1}{2r^{2}}\frac{\mathrm{d}^{2}}{\mathrm{d}\phi^{2}}\psi(\phi).\end{split}$
(1.88)

From now on we use primes for derivatives of functions of one variable, but we
use dot for derivatives with respect to time.

We substitute the result into the remaining determining equations to replace
the functions $s$ and the magnetic field $B$ by the auxiliary functions
$\rho(r),$ $\sigma(r),$ $\tau(\phi),$ $\psi(\phi),$ $\mu(Z)$. The first order
equations (1.4), (1.4) and (1.4) give us one direct condition on the auxiliary
functions $\psi(\phi)$ and $\mu(Z)$ and several conditions on the derivatives
of $m_{1}$ and $m_{2}$:

$\displaystyle\left(r\rho(r)-\psi(\phi)-r^{3}\mu(Z)\right)\left(\psi^{\prime\prime}(\phi)+r^{2}\rho^{\prime}(r)\right)+$
$\displaystyle+\left(r^{3}\mu(Z)+r\rho(r)\right)\psi(\phi)-\psi(\phi)^{2}+r^{3}\tau(\phi)\sigma^{\prime}(r)+$
$\displaystyle+2r^{6}\mu(Z)^{2}-2r^{4}\rho(r)\mu(Z)+2\tau(\phi)^{2}-2r^{3}\partial_{r}m_{1}$
$\displaystyle{}={}$ $\displaystyle 0,$
$\displaystyle\psi^{\prime}(\phi)\left(2r^{3}\mu(Z)-r^{2}\rho^{\prime}(r)-\psi(\phi)-\psi^{\prime\prime}(\phi)\right)+$
$\displaystyle+\tau(\phi)\left(r^{4}\mu^{\prime}(Z)-\tau^{\prime}(\phi)\right)+4r^{4}W_{\phi}-2r^{2}\partial_{\phi}m_{1}$
$\displaystyle{}={}$ $\displaystyle 0,$
$\displaystyle\left(\tau^{\prime}(\phi)-r^{4}\mu^{\prime}(Z)\right)\left(r\rho(r)-\psi(\phi)-r^{3}\mu(Z)\right)-$
$\displaystyle-\psi^{\prime}(\phi)\left(r^{3}\sigma^{\prime}(r)+2\tau(\phi)\right)-2r^{3}\partial_{Z}m_{1}$
$\displaystyle{}={}$ $\displaystyle 0,$ $\displaystyle
r^{3}\mu(Z)(\psi^{\prime\prime}(\phi)+\psi(\phi))+r^{3}\sigma^{\prime}(r)\left(r^{2}\sigma(r)-\tau(\phi)\right)-2r^{6}\mu(Z)^{2}+$
(1.89)
$\displaystyle+r^{5}\mu(Z)\rho^{\prime}(r)+2r^{2}\sigma(r)\tau(\phi)-2\tau(\phi)^{2}-2r^{5}\partial_{r}m_{2}$
$\displaystyle{}={}$ $\displaystyle 0,$
$\displaystyle\left(r^{4}\mu^{\prime}(Z)-\tau^{\prime}(\phi)\right)\left(r^{2}\sigma(r)-\tau(\phi)\right)-2r^{4}\partial_{\phi}m_{2}$
$\displaystyle{}={}$ $\displaystyle 0,$
$\displaystyle-r^{4}\mu(Z)\mu^{\prime}(Z)+\mu(Z)\tau^{\prime}(\phi)+4r^{2}W_{Z}-2r^{2}\partial_{Z}m_{2}$
$\displaystyle{}={}$ $\displaystyle 0,$
$\displaystyle\mu(Z)\psi^{\prime\prime}(\phi)$ $\displaystyle{}={}$
$\displaystyle 0,$
$\displaystyle\left(r^{2}\rho(r)-r\psi(\phi)-r^{4}\mu(Z)\right)\mu^{\prime}(Z)+\mu(Z)\tau^{\prime}(\phi)+2\partial_{Z}m_{1}$
$\displaystyle{}={}$ $\displaystyle 0,$
$\displaystyle\tau^{\prime}(\phi)\left(r^{2}\sigma(r)-\tau(\phi)\right)+r^{4}\tau(\phi)\mu^{\prime}(Z)+r^{3}\mu(Z)\psi^{\prime}(\phi)+2r^{4}\partial_{\phi}m_{2}$
$\displaystyle{}={}$ $\displaystyle 0.$

We see that we have two equations for each of the derivatives
$\partial_{Z}m_{1}$ and $\partial_{\phi}m_{2}$, so we obtain constraints for
the auxiliary functions (the two values of $\partial_{Z}m_{1}$ and
$\partial_{\phi}m_{2}$ must coincide). Assuming $m_{1}$ and $m_{2}$ to be
sufficiently smooth, we impose the Clairaut compatibility conditions
$\partial_{ba}m_{i}=\partial_{ab}m_{i}$ on the second derivatives and obtain
equations for the mixed second derivatives of the scalar potential $W$, which
we list shortly, after we analyse the zeroth order equations.

Only one of the zeroth order equations (1.82) and (1.4) obtains the quantum
correction:

$\displaystyle\begin{aligned}
&\left(r^{2}\sigma(r)-\tau(\phi)\right)W_{Z}+r^{2}\mu(Z)W_{\phi}&{}={}&0,&\\\
&\left(r\rho(r)-\psi(\phi)-r^{3}\mu(Z)\right)W_{\phi}+r\left(\psi^{\prime}(\phi)W_{r}+\tau(\phi)W_{Z}\right)+&\\\
&+\hbar^{2}\frac{\psi^{\prime\prime\prime}(\phi)+\psi^{\prime}(\phi)}{4r^{2}}&{}={}&0,&\\\
&2r^{4}\left(r^{3}\mu(Z)-r\rho(r)+\psi(\phi)\right)\partial_{\phi}m_{2}+2r^{3}\left(r^{2}\sigma(r)-\tau(\phi)\right)\partial_{Z}m_{1}+&\\\
&+r^{3}\mu(Z)\psi^{\prime}(\phi)\psi^{\prime\prime}(\phi)+2r^{5}\mu(Z)\partial_{\phi}m_{1}-2r^{5}\tau(\phi)\partial_{Z}m_{2}+&\\\
&+\left[r^{3}\left(r^{2}\sigma(r)-\tau(\phi)\right)\sigma^{\prime}(r)-2r^{6}\mu(Z)^{2}+r^{5}\mu(Z)\rho^{\prime}(r)+\right.&\\\
&\left.+r^{3}\mu(Z)\psi(\phi)+2r^{2}\sigma(r)\tau(\phi)-2\tau(\phi)^{2}-2r^{5}\partial_{r}m_{2}\right]\psi^{\prime}(\phi)-&\\\
&-\left[r^{2}\left(r^{3}\mu(Z)-r\rho(r)+\psi(\phi)\right)\sigma(r)\right.+&\\\
&+\left.\tau(\phi)\left(r\rho(r)-\psi(\phi)\right)\right]\left(r^{4}\mu^{\prime}(Z)-\tau^{\prime}(\phi)\right)&{}={}&0.&\end{aligned}$
(1.90)

Substituting for the derivatives of $m$ from section 1.4 into eq. (1.90), we
obtain a system of linear inhomogeneous algebraic equations for the first
derivatives of the scalar potential $W$.

The final form of the reduced equations, which we analyse hereafter, is as
follows. (Indexes of $W$ mean partial derivatives.)

$\displaystyle\psi^{\prime}(\phi)\left(r^{3}\sigma^{\prime}(r)+2\tau(\phi)\right)-\tau^{\prime}(\phi)\left(r\rho(r)-\psi(\phi)\right)={}$
$\displaystyle 0,$ (1.91)
$\displaystyle\mu(Z)\psi^{\prime}(\phi)+r^{3}\sigma(r)\mu^{\prime}(Z)={}$
$\displaystyle 0,$ (1.92) $\displaystyle W_{r\phi}={}$
$\displaystyle-\frac{2}{r}W_{\phi}+\frac{1}{4r^{5}}\left[\psi^{\prime}(\phi)\left(r^{3}(\rho^{\prime\prime}(r)-\mu(Z))-r^{2}\rho^{\prime}(r)+r\rho(r)-3\psi^{\prime\prime}(\phi)-4\psi(\phi)\right)+\right.$
$\displaystyle+\left.\tau^{\prime}(\phi)\left(r^{3}\sigma^{\prime}(r)+2\tau(\phi)\right)-2r^{4}\tau(\phi)\mu^{\prime}(Z)-\psi^{\prime\prime\prime}(\phi)\left(\psi(\phi)-r\rho(r)\right)\right],$
$\displaystyle W_{\phi Z}={}$
$\displaystyle-\frac{1}{4r^{2}}\left[r^{2}\mu^{\prime\prime}(Z)\left(\tau(\phi)-r^{2}\sigma(r)\right)+\tau^{\prime\prime}(\phi)\mu(Z)\right],$
(1.93) $\displaystyle W_{rZ}={}$
$\displaystyle\frac{1}{4r^{3}}\left[r\mu^{\prime}(Z)\left(r^{2}\rho^{\prime}(r)+\psi(\phi)-2r^{3}\mu(Z)\right)+2\mu(Z)\tau^{\prime}(\phi)\right],$

$\displaystyle\begin{pmatrix}0&r^{2}\mu(Z)&r^{2}\sigma(r)-\tau(\phi)\\\
\psi^{\prime}(\phi)&\rho(r)-r^{2}\mu(Z)-\frac{\psi(\phi)}{r}&\tau(\phi)\\\
0&4r^{7}\mu(Z)&-4r^{5}\tau(\phi)\end{pmatrix}\cdot\begin{pmatrix}W_{r}\\\
W_{\phi}\\\ W_{Z}\end{pmatrix}=\begin{pmatrix}0\\\
-\frac{\hbar^{2}(\psi^{\prime\prime\prime}(\phi)+\psi^{\prime}(\phi))}{4r^{3}}\\\
\alpha(r,\phi,Z)\end{pmatrix}$ (1.94)

with

$\displaystyle\alpha(r,\phi,Z)={}$
$\displaystyle\psi^{\prime}(\phi)\left[\left(r^{5}\sigma(r)-r^{3}\tau(\phi)\right)\sigma^{\prime}(r)+r^{5}\mu(Z)\rho^{\prime}(r)-2\tau(\phi)^{2}\right.+$
$\displaystyle\left.+2r^{2}\sigma(r)\tau(\phi)-r^{3}\mu(Z)\left(r^{3}\mu(Z)+r\rho(r)-2\psi(\phi)\right)\right]+$
$\displaystyle+\tau^{\prime}(\phi)\left[\left(r\rho(r)-\psi(\phi)\right)\tau(\phi)+r^{2}\sigma(r)\left(r^{3}\mu(Z)-r\rho(r)+\psi(\phi)\right)\right]-$
$\displaystyle-r^{4}\mu^{\prime}(Z)\tau(\phi)\left(r\rho(r)-\psi(\phi)\right).$
(1.95)

We denote the matrix in eq. 1.94 by $M$.

The only change with respect to the classical case from [11] is the non-zero
RHS in the system of linear inhomogeneous algebraic equations (1.94),
corresponding to eq. (39) in the original article.

We separate the analysis of the reduced system into several cases with respect
to the rank of the matrix $M$, following [11]. It can be either $3$, $2$ or
$1$. Rank $0$ would mean that all the auxiliary functions vanish and with them
the magnetic field as well, see eq. 1.88, so we rule this case out.

If the rank is $3$, then the determinant of $M$,

$\det(M)=4r^{9}\psi^{\prime}(\phi)\mu(Z)\sigma(r),$ (1.96)

is not zero and it implies a unique solution for each first derivative of $W$.
The analysis of this case from [11] shows that this assumption leads to a
contradiction and the reduced system is inconsistent. The analysis remains
valid because it uses the matrix $M$ and equations (1.91) and (1.92) only,
which are unaffected by the quantum correction.

If the rank is either $2$ or $1$ instead, then $\det(M)=0$, and from eq. 1.96
we see that there are a priori three possible cases:

  1. a)

$\psi^{\prime}(\phi)=0$,

  2. b)

$\psi^{\prime}(\phi)\neq 0$ and $\mu(Z)=0$,

  3. c)

$\psi^{\prime}(\phi)\neq 0$, $\mu(Z)\neq 0$ and $\sigma(r)=0$. We rule this
case out due to inconsistency with eq. 1.92.

The assumption of case a) implies that the quantum correction vanishes; thus,
we obtain the classical systems from [11]. Therefore, we continue with the
case b) only, i.e. $\psi^{\prime}(\phi)\neq 0$ and $\mu(Z)=0$. (The key
results for case a) are summarised in the corresponding subsections of Section
2.1, where we search for additional integrals of motion.)

Before we go to the specific subcases, we show that $\alpha=0$ for $\alpha$
from section 1.4. In the considered case it follows from eq. 1.91 and eq. 1.92
only, so the considerations are the same as in [11]:

Assuming $\psi^{\prime}(\phi)\neq 0$ and $\mu(Z)=0$, eq. 1.92 is satisfied
trivially. By differentiating eq. 1.91 with respect to $r$, we obtain

$\left(r^{3}\sigma^{\prime}(r)\right)^{\prime}=\frac{\tau^{\prime}(\phi)}{\psi^{\prime}(\phi)}\left(r\rho(r)\right)^{\prime}.$
(1.97)

This leads to 3 possibilities:

  * •

$\left(r^{3}\sigma^{\prime}(r)\right)^{\prime}=\left(r\rho(r)\right)^{\prime}=0$,
i.e.

$\sigma(r)=\frac{C_{\sigma}}{r^{2}}+\tilde{C}_{\sigma},\qquad\rho(r)=\frac{C_{\rho}}{r}.$
(1.98)

Substituting eq. 1.98 into eq. 1.91 we find

$2\left(\tau(\phi)-C_{\sigma}\right)\psi^{\prime}(\phi)+\left(\psi(\phi)-C_{\rho}\right)\tau^{\prime}(\phi)=0$
(1.99)

which directly implies that $\alpha$ defined in section 1.4 vanishes.

  * •

$\left(r^{3}\sigma^{\prime}(r)\right)^{\prime}=\tau^{\prime}(\phi)=0$, i.e.

$\sigma(r)=\frac{C_{\sigma}}{r^{2}}+\tilde{C}_{\sigma},\qquad\tau(\phi)=C_{\tau}.$
(1.100)

Substituting eq. 1.100 into eq. 1.91 we find $C_{\sigma}=C_{\tau}$ and that
together with eq. 1.100 implies again that $\alpha=0$ in section 1.4.

  * •

$\frac{\left(r^{3}\sigma^{\prime}(r)\right)^{\prime}}{\left(r\rho(r)\right)^{\prime}}=\frac{\tau^{\prime}(\phi)}{\psi^{\prime}(\phi)}=\lambda\neq
0,$ implying that

$\rho(r)=\frac{1}{\lambda}r^{2}\sigma^{\prime}(r)+\frac{C_{\rho}}{r},\qquad\tau(\phi)=\lambda\psi(\phi)+C_{\tau}.$
(1.101)

Substituting eq. 1.101 into eq. 1.91 and differentiating it with respect to
$\phi$ we arrive at a contradiction with our assumptions, namely
$\lambda\psi^{\prime}(\phi)=0$.

####  1.4.1 Case b) $\mathrm{rank}\,(M)=2$: $\psi^{\prime}(\phi)\neq 0$ and
$\mu(Z)=0$

Here the second and third equations in (1.4) reduce to

$\partial_{rZ}W=0,\quad\partial_{\phi Z}W=0,$ (1.102)

implying separation of the $Z$ coordinate from $r$ and $\phi$ in the scalar
potential $W$, i.e.

$W(r,\phi,Z)=W_{12}(r,\phi)+W_{3}(Z).$ (1.103)

The reduced row echelon form of the extended matrix of the system of equations
(1.94) reads

$\begin{pmatrix}r\psi^{\prime}(\phi)&r\rho(r)-\psi(\phi)&0&-\frac{\hbar^{2}(\psi^{\prime\prime\prime}(\phi)+\psi^{\prime}(\phi))}{4r^{2}}\\\
0&0&\sigma(r)&0\\\ 0&0&\tau(\phi)&0\end{pmatrix}$ (1.104)

Our assumption $\psi^{\prime}(\phi)\neq 0$ implies that we have two
possibilities to have $\mathrm{rank}\,(M)=2$, namely either $\sigma(r)\neq 0$
or $\tau(\phi)\neq 0$. Both of them imply $W_{3}(Z)^{\prime}=0$ and without
loss of generality we can absorb the constant into redefinition of $W_{12}$,
so we have $W(r,\phi,Z)=W_{12}(r,\phi)\equiv W(r,\phi)$.

  1. 1)

$\sigma(r)\neq 0$: We follow the further splitting into the subcases from
[11]. We use our assumption $\psi^{\prime}(\phi)\neq 0$ to rewrite eq. 1.91 in
the following way:

$r^{3}\sigma^{\prime}(r)+2\tau(\phi)-\frac{\tau^{\prime}(\phi)}{\psi^{\prime}(\phi)}(r\rho(r)-\psi(\phi))=0.$
(1.105)

Differentiation with respect to $\phi$ leads to the equation:

$3\tau^{\prime}(\phi)+\psi(\phi)\frac{\tau^{\prime\prime}(\phi)\psi^{\prime}(\phi)-\tau^{\prime}(\phi)\psi^{\prime\prime}(\phi)}{\psi^{\prime}(\phi)^{2}}=r\rho(r)\frac{\tau^{\prime\prime}(\phi)\psi^{\prime}(\phi)-\tau^{\prime}(\phi)\psi^{\prime\prime}(\phi)}{\psi^{\prime}(\phi)^{2}}.$
(1.106)

If
$\tau^{\prime\prime}(\phi)\psi^{\prime}(\phi)-\tau^{\prime}(\phi)\psi^{\prime\prime}(\phi)\neq
0$, we can separate the variables $r$ and $\phi$, else the expression vanishes
and we conclude from eq. 1.106 that $\tau^{\prime}(\phi)=0$, thus
$\tau(\phi)=\tau_{0}$ is a constant. We treat the subcases separately.

    1. 1.1)

$\tau^{\prime\prime}(\phi)\psi^{\prime}(\phi)-\tau^{\prime}(\phi)\psi^{\prime\prime}(\phi)=0$,
so $\tau(\phi)=\tau_{0}$: From eq. 1.91, which now reads
$r^{3}\sigma^{\prime}(r)=-2\tau_{0}$, follows

$\sigma(r)=\frac{\tau_{0}}{r^{2}}+\sigma_{0}.$ (1.107)

The yet unsolved equations from eq. (1.91)–(1.94) constraining the scalar
potential $W$ read

$\displaystyle
r\psi^{\prime}(\phi)W_{r}+\left(r\rho(r)-\psi(\phi)\right)W_{\phi}+\frac{\hbar^{2}(\psi^{\prime\prime\prime}(\phi)+\psi^{\prime}(\phi))}{4r^{2}}$
$\displaystyle=0,$ (1.108)
$\displaystyle\psi^{\prime}(\phi)\left(r^{3}\rho^{\prime\prime}(r)-r^{2}\rho^{\prime}(r)+r\rho(r)-3\psi^{\prime\prime}(\phi)-4\psi(\phi)\right)+$
$\displaystyle+\psi^{\prime\prime\prime}(\phi)\left(r\rho(r)-\psi(\phi)\right)-4r^{5}W_{r\phi}-8r^{4}W_{\phi}$
$\displaystyle=0.$ (1.109)

The magnetic field takes the form

$B^{\phi}=0,\quad B^{r}=0,\quad
B^{Z}=-\frac{\rho^{\prime}(r)}{2}-\frac{\psi^{\prime\prime}(\phi)+\psi(\phi)}{2r^{2}}.$
(1.110)

Thus, we have a free motion in the $z$-direction (in the direction of the
$z$-axis) and motion in the $xy$-plane under the influence of the scalar
potential $W(r,\phi)$ constrained by eq. 1.108 and eq. 1.109 and the
perpendicular magnetic field $B^{Z}(r,\phi)$. Such 2D problem was discussed by
McSween and Winternitz [22] in classical mechanics and by Bérubé and
Winternitz [5] in quantum mechanics.

For both motions we have one integral of motion in addition to the separable
Hamiltonian, namely the integral $X_{1}$ for the motion in the $xy$-plane,
which we list shortly together with corresponding magnetic field $B$ and
potential $W$, and $X_{2}$ for the $z$-direction. The integral $X_{2}$ reduces
to the first order one

$\tilde{X}_{2}=p_{Z}^{A}+\frac{\sigma_{0}}{2},$ (1.111)

which in a suitably chosen gauge reads $p_{Z}$.

In [5] it was shown that there are 2 solutions to equations (1.108) and
(1.109) for $\psi^{\prime}(\phi)\neq 0$, in both cases the magnetic field
$B^{Z}$ does not depend on $Z$. Omitting the lengthy details, we present the
results only.

      1. i.

The magnetic field independent of $\phi$ and $Z$ corresponds to

$\begin{split}\psi(\phi)&=\psi_{0}+\psi_{1}\cos(\phi)+\psi_{2}\sin(\phi),\\\
\rho(r)&=3\rho_{2}r^{4}-\rho_{1}r^{2}+\rho_{0}\end{split}$ (1.112)

which means that this quantum system is the same as the classical (the
correction vanishes). The resulting magnetic field $B^{Z}$ and potential $W$
are

$\begin{split}B^{Z}&=-6\rho_{2}r^{3}+\rho_{1}r,\\\
W&=-2\rho_{2}r(\psi_{1}\cos(\phi)+\psi_{2}\sin(\phi))-\rho_{2}^{2}r^{6}+\frac{\rho_{2}\rho_{1}}{2}r^{4}-\rho_{2}W_{0}r^{2}.\end{split}$
(1.113)

The corresponding integral of motion in the Cartesian coordinates is

$\displaystyle\begin{aligned}
X_{1}={}&(L_{z}^{A})^{2}+(3\rho_{2}r^{4}-\rho_{1}r^{2}+W_{0})L_{z}^{A}-\psi_{2}p_{x}^{A}+\psi_{1}p_{y}^{A}+\\\
&+\psi_{1}(2\rho_{2}r^{2}-\rho_{1})x+\psi_{2}(2\rho_{2}r^{2}-\rho_{1})y+\\\
&+\frac{(3\rho_{2}r^{4}-\rho_{1}r^{2}+2W_{0})(3\rho_{2}r^{2}-\rho_{1})r^{2}}{4},\end{aligned}$
(1.114)

where we write $r=\sqrt{x^{2}+y^{2}}$ for brevity.

      2. ii.

The magnetic field depends on $\phi$ as well as $r$. We have

$\rho(r)=\frac{\rho_{0}}{r}.$ (1.115)

This leads to the special case of the following subcase 1)1.2) with
$\tau_{1}=0$. We therefore refer the reader there for details and list the
results only.

The magnetic field $B$ together with the potential $W$ read

$\displaystyle B^{r}={}$ $\displaystyle 0,\quad B^{\phi}={}0,\quad
B^{Z}=\frac{2\beta_{1}\beta(\phi)^{2}+\beta_{2}}{4r^{2}\beta(\phi)^{5}},$
(1.116) $\displaystyle W={}$
$\displaystyle\frac{W_{0}}{r^{2}\beta(\phi)^{2}}-\frac{\beta_{2}}{32\beta(\phi)^{4}r^{4}}+\hbar^{2}\frac{f_{1}\beta(\phi)^{4}-12\beta_{1}\beta(\phi)^{2}-5\beta_{2}}{32r^{2}\beta(\phi)^{6}},$
(1.117)

where $\beta_{1}$, $\beta_{2}$ are constants and
$\beta(\phi)=\psi(\phi)-\rho_{0}$ must satisfy

$4\beta(\phi)^{4}\beta^{\prime}(\phi)^{2}+4\beta(\phi)^{6}-4\beta_{1}\beta(\phi)^{2}+f_{1}\beta(\phi)^{4}=\beta_{2}.$
(1.118)

The corresponding integral $X_{1}$ is defined by

$\displaystyle s_{1}^{r}={}$
$\displaystyle\frac{\sqrt{4\beta_{1}\beta(\phi)^{2}+\beta_{2}-4\beta(\phi)^{6}-\beta(\phi)^{4}f_{1}}}{2\beta(\phi)^{2}},$
$\displaystyle s_{1}^{\phi}={}$ $\displaystyle-\frac{\beta(\phi)}{r},\qquad
s_{1}^{Z}=0,$ (1.119) $\displaystyle m_{1}={}$
$\displaystyle\frac{2W_{0}}{\beta(\phi)^{2}}-\frac{2\beta_{1}\beta(\phi)^{2}+\beta_{2}}{8\beta(\phi)^{4}r^{2}}+\hbar^{2}\frac{f_{1}\beta(\phi)^{4}-12\beta_{1}\beta(\phi)^{2}-5\beta_{2}}{16\beta(\phi)^{6}}.$

    2. 1.2)

$\tau^{\prime\prime}(\phi)\psi^{\prime}(\phi)-\tau^{\prime}(\phi)\psi^{\prime\prime}(\phi)\neq
0$. Separating the variables $r$ and $\phi$ in eq. 1.106,

$\frac{3\tau^{\prime}(\phi)\psi^{\prime}(\phi)^{2}}{\tau^{\prime\prime}(\phi)\psi^{\prime}(\phi)-\tau^{\prime}(\phi)\psi^{\prime\prime}(\phi)}+\psi(\phi)=\rho_{0}=r\rho(r)$
(1.120)

with separation constant $\rho_{0}$, we obtain

$\rho(r)=\frac{\rho_{0}}{r},\qquad\tau(\phi)=\tau_{0}+\frac{\tau_{1}}{(\psi(\phi)-\rho_{0})^{2}}.$
(1.121)

Using these results in eq. 1.91, we get
$\sigma(r)=\frac{\tau_{0}}{r^{2}}+\sigma_{0}$.

There are two unsolved equations among eq. (1.4)–(1.94). Simplifying them
using the previous results, we get

$\displaystyle
r\psi^{\prime}(\phi)W_{r}+\left(\rho_{0}-\psi(\phi)\right)W_{\phi}+\hbar^{2}\frac{\psi^{\prime\prime\prime}(\phi)+\psi^{\prime}(\phi)}{4r^{2}}$
$\displaystyle=0,$ $\displaystyle
3\psi^{\prime}(\phi)\psi^{\prime\prime}(\phi)+4\psi^{\prime}(\phi)\left(\psi(\phi)-\rho_{0}\right)+\psi^{\prime\prime\prime}(\phi)(\psi(\phi)-\rho_{0})+$
(1.122)
$\displaystyle+\frac{4\tau_{1}^{2}}{\left(\psi(\phi)-\rho_{0}\right)^{5}}\psi^{\prime}(\phi)+4r^{5}W_{r\phi}+8r^{4}W_{\phi}$
$\displaystyle=0.$

We substitute $\beta(\phi)=\psi(\phi)-\rho_{0}$ and integrate the second
equation with respect to $\phi$ to obtain

$\displaystyle
r\beta^{\prime}(\phi)W_{r}-\beta(\phi)W_{\phi}+\hbar^{2}\frac{\beta^{\prime\prime\prime}(\phi)+\beta^{\prime}(\phi)}{4r^{2}}={}$
$\displaystyle 0,$ (1.123)
$\displaystyle\beta(\phi)\beta^{\prime\prime}(\phi)+\beta^{\prime}(\phi)^{2}+2\beta(\phi)^{2}-\frac{\tau_{1}^{2}}{\beta(\phi)^{4}}+$
$\displaystyle+4r^{5}W_{r}+8r^{4}W(r,\phi)+f(r)={}$ $\displaystyle 0,$ (1.124)

where $f(r)$ is an arbitrary function arising in the integration.

Substituting for $W_{r}$ from eq. 1.124 into eq. 1.123 we find expressions for
both $W_{r}$ and $W_{\phi}$, where only the expression for $W_{\phi}$ has the
correction,

$\partial_{\phi}W=-\frac{2\beta^{\prime}W}{\beta}-\frac{\beta^{\prime}(2\beta^{6}+\beta^{5}\beta^{\prime\prime}+\beta^{4}\beta^{\prime
2}+\beta^{4}f\left(r\right)-\tau_{1}^{2})}{4r^{4}\beta^{5}}+\hbar^{2}\frac{\beta^{\prime\prime\prime}+\beta^{\prime}}{4r^{2}\beta}.$
(1.125)

Substituting them into item 1)1.2) we obtain the following equation

$\beta^{\prime}(\phi)\left(12\beta(\phi)^{2}+7\beta\beta^{\prime\prime}(\phi)+4\beta^{\prime}(\phi)^{2}-f^{\prime}(r)r+4f(r)\right)+\beta(\phi)^{2}\beta^{\prime\prime\prime}(\phi)=0,$
(1.126)

which we differentiate by $r$ and get

$f(r)=\frac{f_{1}}{4}+f_{2}r^{4},$ (1.127)

where $f_{1}$, $f_{2}$ are integration constants. With the form of $f(r)$ from
eq. 1.127, we can determine the dependence of the scalar potential $W$ on $r$
from eq. 1.124, treating the yet unknown function $\beta(\phi)$ as a
parameter:

$W=-\frac{f_{2}}{8}+\frac{\tilde{W}(\phi)}{r^{2}}+\frac{\beta(\phi)\beta^{\prime\prime}(\phi)+\beta^{\prime}(\phi)^{2}+\frac{f_{1}}{4}-\frac{\tau_{1}^{2}}{\beta(\phi)^{4}}+2\beta(\phi)^{2}}{8r^{4}}.$
(1.128)

We determine the function $\tilde{W}(\phi)$ from an equation obtained by
inserting $W$ from eq. (1.128) into eq. 1.125 and subtracting eq. 1.126 with
$f(r)$ from eq. 1.127, i.e. we subtract

$\displaystyle\beta^{\prime}(\phi)\left(7\beta(\phi)\beta^{\prime\prime}(\phi)+4\beta^{\prime}(\phi)^{2}+12\beta(\phi)^{2}+f_{1}\right)+\beta(\phi)^{2}\beta^{\prime\prime\prime}(\phi)=0.$
(1.129)

The resulting equation

$\frac{\hbar^{2}\beta^{\prime\prime\prime}(\phi)+(\hbar^{2}-8\tilde{W}(\phi))\beta^{\prime}(\phi)-4(\tilde{W}^{\prime}(\phi))\beta(\phi)}{4r^{2}\beta(\phi)}=0$
(1.130)

yields

$\tilde{W}(\phi)=\frac{W_{0}}{\beta(\phi)^{2}}+\hbar^{2}\frac{2\beta(\phi)\beta^{\prime\prime}(\phi)+\beta(\phi)^{2}-\beta^{\prime}(\phi)^{2}}{8\beta(\phi)^{2}},$
(1.131)

where $W_{0}$ is an arbitrary constant. The final form of the scalar potential
after a constant shift eliminating $-\frac{f_{1}}{8}$ is

$\begin{split}W={}&\frac{W_{0}}{r^{2}\beta(\phi)^{2}}+\frac{\beta(\phi)\beta^{\prime\prime}(\phi)+\beta^{\prime}(\phi)^{2}+\frac{f_{1}}{4}-\frac{\tau_{1}^{2}}{\beta(\phi)^{4}}+2\beta(\phi)^{2}}{8r^{4}}+\\\
&+\hbar^{2}\frac{2\beta(\phi)\beta^{\prime\prime}(\phi)+\beta(\phi)^{2}-\beta^{\prime}(\phi)^{2}}{8r^{2}\beta(\phi)^{2}}.\end{split}$
(1.132)

We have therefore almost explicit form of the potential, the only undetermined
function is $\beta(\phi)$, which must satisfy the only remaining equation
(1.129). We reduce its order as follows [11]: we multiply by $\beta(\phi)$ and
integrate, then multiply by $\beta^{\prime}(\phi)\beta(\phi)$ and integrate
again. The result is

$4\beta(\phi)^{4}\beta^{\prime}(\phi)^{2}+4\beta(\phi)^{6}-4\beta_{1}\beta(\phi)^{2}+f_{1}\beta(\phi)^{4}=\beta_{2}$
(1.133)

where $\beta_{1},$ $\beta_{2}$ are the constants of integration. According to
[5], eq. 1.133 can be written as a quadrature expressing the independent
variable $\phi$ as a function of $\beta$ in terms of elliptic integrals but
the authors analyse further some special cases only. For more details see the
cited article, we will later consider only the special case $\beta_{2}=0$.

We can use eq. 1.133 to further simplify the scalar potential (1.132)

$W=\frac{W_{0}}{r^{2}\beta(\phi)^{2}}-\frac{(4\tau_{1}^{2}+\beta_{2})}{32r^{4}\beta(\phi)^{4}}+\hbar^{2}\frac{f_{1}\beta(\phi)^{4}-12\beta_{1}\beta(\phi)^{2}-5\beta_{2}}{32r^{2}\beta(\phi)^{6}}.$
(1.134)

The magnetic field is also expressed in terms of the function $\beta(\phi)$
and its derivatives

$\displaystyle
B^{r}=-\frac{\tau_{1}\beta^{\prime}(\phi)}{r^{2}\beta(\phi)^{3}},\quad
B^{\phi}=\frac{\tau_{1}}{r^{3}\beta(\phi)^{2}},\quad
B^{Z}=-\frac{\beta(\phi)+\beta^{\prime\prime}(\phi)}{2r^{2}},$ (1.135)

or substituting for the derivatives from eq. (1.129) and its integrated form
eq. (1.133)

$\displaystyle B^{r}={}$
$\displaystyle-\tau_{1}\frac{\sqrt{4\beta_{1}\beta(\phi)^{2}+\beta_{2}-4\beta(\phi)^{6}-f_{1}\beta(\phi)^{4}}}{2r^{2}\beta(\phi)^{5}},$
$\displaystyle B^{\phi}={}$
$\displaystyle\frac{\tau_{1}}{r^{3}\beta(\phi)^{2}},\quad
B^{Z}=\frac{2\beta_{1}\beta(\phi)^{2}+\beta_{2}}{4r^{2}\beta(\phi)^{5}}.$
(1.136)

(The square root appears due to substitution for $\beta^{\prime}(\phi)$ from
eq. (1.133), so the sign depends on the choice of the branch in that
equation.)

We note that eq. (1.133) is the same in classical and quantum cases. This
implies that the magnetic field is the same in both cases, only the scalar
potential $W$ obtains an $\hbar^{2}$-proportional correction depending on
$\beta$.

The integrals from LABEL:cyl_integrals are determined by

$\displaystyle s_{1}^{r}={}$
$\displaystyle\frac{\sqrt{4\beta_{1}\beta(\phi)^{2}+\beta_{2}-4\beta(\phi)^{6}-\beta(\phi)^{4}f_{1}}}{2\beta(\phi)^{2}},$
$\displaystyle s_{1}^{\phi}={}$ $\displaystyle-\frac{\beta(\phi)}{r},\qquad
s_{1}^{Z}=\tau_{0}+\frac{\tau_{1}}{\beta(\phi)^{2}},$ $\displaystyle m_{1}={}$
$\displaystyle\frac{2W_{0}}{\beta(\phi)^{2}}-\frac{4\beta(\phi)^{2}\tau_{0}\tau_{1}+2\beta_{1}\beta(\phi)^{2}+4\tau_{1}^{2}+\beta_{2}}{8\beta(\phi)^{4}r^{2}}+$
$\displaystyle+\hbar^{2}\frac{f_{1}\beta(\phi)^{4}-12\beta_{1}\beta(\phi)^{2}-5\beta_{2}}{16r^{2}\beta(\phi)^{6}},$
(1.137) $\displaystyle s_{2}^{r}={}$ $\displaystyle 0,\qquad
s_{2}^{\phi}=0,\qquad
s_{2}^{Z}=\sigma_{0}-\frac{\tau_{1}}{r^{2}\beta(\phi)^{2}},$ $\displaystyle
m_{2}={}$
$\displaystyle\frac{\tau_{1}}{\beta(\phi)^{2}r^{2}}\left(\frac{\tau_{1}}{4\beta(\phi)^{2}r^{2}}-\frac{\sigma_{0}}{2}\right).$

Writing the integral $X_{2}$ in full, we get

$X_{2}=\left(p_{Z}^{A}-\frac{\tau_{1}}{2r^{2}\beta(\phi)^{2}}\right)^{2}+\sigma_{0}\left(p_{Z}^{A}-\frac{\tau_{1}}{2r^{2}\beta(\phi)^{2}}\right).$
(1.138)

so this integral of motion reduces to the first order one

$\tilde{X}_{2}=p_{Z}^{A}-\frac{\tau_{1}}{2r^{2}\beta(\phi)^{2}}$ (1.139)

and is in suitable gauge $X_{2}=p_{Z}$, as can be checked in eq. (1.135).

We can rewrite eq. 1.133 using a substitution $\gamma(\phi)=\beta(\phi)^{2}$
to obtain

$\gamma(\phi)\gamma^{\prime}(\phi)^{2}+4\gamma(\phi)^{3}-4\beta_{1}\gamma(\phi)+f_{1}\gamma(\phi)^{2}=\beta_{2}.$
(1.140)

Let us analyse the special case $\beta_{2}=0$. We obtain the following
solution:

$\gamma(\phi)=\frac{\sqrt{64\beta_{1}+f_{1}^{2}}\sin\left(2(\phi-\phi_{0})\right)-f_{1}}{8}.$
(1.141)

Under the assumptions that $f_{1}<0$, $-\frac{f_{1}^{2}}{64}<\beta_{1}<0$, the
corresponding $\beta(\phi)$ is well defined, bounded and positive,

$\beta(\phi)=\sqrt{\frac{\sqrt{64\beta_{1}+f_{1}^{2}}\sin\left(2(\phi-\phi_{0})\right)-f_{1}}{8}}.$
(1.142)

The magnetic field now reads:

$\displaystyle B^{r}={}$
$\displaystyle-\frac{8\tau_{1}\sqrt{64\beta_{1}+f_{1}^{2}}\cos\left(2(\phi-\phi_{0})\right)}{r^{2}\left(\sqrt{64\beta_{1}+f_{1}^{2}}\sin\left(2(\phi-\phi_{0})\right)-f_{1}\right)^{2}},$
$\displaystyle B^{\phi}={}$
$\displaystyle\frac{8\tau_{1}}{r^{3}\left(\sqrt{64\beta_{1}+f_{1}^{2}}\sin\left(2(\phi-\phi_{0})\right)-f_{1}\right)},$
(1.143) $\displaystyle B^{Z}={}$
$\displaystyle\frac{\beta_{1}}{2r^{2}}\left(\frac{\sqrt{64\beta_{1}+f_{1}^{2}}\sin\left(2(\phi-\phi_{0})\right)-f_{1}}{8}\right)^{-\frac{3}{2}}$

For solution from eq. 1.141 we find

$\displaystyle W$
$\displaystyle=\frac{8W_{0}}{r^{2}\left(a\sin\left(2(\phi-\phi_{0})\right)-f_{1}\right)}-\frac{8\tau_{1}^{2}}{r^{4}\left(a\sin\left(2(\phi-\phi_{0})\right)-f_{1}\right)^{2}}+$
$\displaystyle+\hbar^{2}\frac{a^{2}f_{1}(\sin(2(\phi-\phi_{0}))^{2}+1)-2a^{3}\sin(2(\phi-\phi_{0}))+32\beta_{1}f_{1}-320\beta_{2}}{4r^{2}\left(a\sin(2(\phi-\phi_{0}))-f_{1}\right)^{3}},$
(1.144)

where $a=\sqrt{f_{1}^{2}+64\beta_{1}}$.

Choosing the coordinates so that $\phi_{0}=0$ and transforming the magnetic
field into the Cartesian coordinates using section 1.2, we get

$\displaystyle B^{x}(\vec{x})$
$\displaystyle=-\frac{8\tau_{1}\left(\sqrt{f_{1}^{2}+64\beta_{1}}x-f_{1}y\right)}{\left(f_{1}(x^{2}+y^{2})-2\sqrt{f_{1}^{2}+64\beta_{1}}xy\right)^{2}},$
(1.145) $\displaystyle B^{y}(\vec{x})$
$\displaystyle=\frac{8\tau_{1}\left(\sqrt{f_{1}^{2}+64\beta_{1}}y-f_{1}x\right)}{\left(f_{1}(x^{2}+y^{2})-2\sqrt{f_{1}^{2}+64\beta_{1}}xy\right)^{2}},$
(1.146) $\displaystyle B^{z}(\vec{x})$
$\displaystyle=-\frac{8\sqrt{2}\beta_{1}}{\left(f_{1}(x^{2}+y^{2})-2\sqrt{f_{1}^{2}+64\beta_{1}}xy\right)^{\frac{3}{2}}},$
(1.147) $\displaystyle W(\vec{x})$
$\displaystyle=-\frac{8W_{0}\left(f_{1}(x^{2}+y^{2})-2\sqrt{f_{1}^{2}+64\beta_{1}}xy\right)+8\tau_{1}^{2}}{\left(f_{1}(x^{2}+y^{2})-2\sqrt{f_{1}^{2}+64\beta_{1}}xy\right)^{2}}+$
(1.148)
$\displaystyle+\hbar^{2}\frac{(r^{4}+4x^{2}y^{2})f_{1}^{3}-4a^{3}r^{2}xy+\beta_{1}f_{1}(96r^{4}+256y^{2}x^{2})-320r^{4}\beta_{2}}{4(f_{1}r^{2}-2ayx)^{3}}.$

  2. 2)

$\tau(\phi)\neq 0$: We can use considerations from the previous case as in the
classical case, because we have never divided by $\sigma(r)$, which may now
vanish, and never assumed $\tau(\phi)=0$. We only need some constraints on the
constants to assure $\tau(\phi)\neq 0$. We therefore obtain the system which
splits into independent motions in $xy$-plane and $z$-direction (the motion in
2D was considered in [5]) with magnetic field from eq. 1.113 and the system
with the magnetic field from item 1)1.2) and scalar potential from eq. 1.134.
(If $\tau_{1}=0$ in the second system, we get the separating system from eq.
1.116 in subcase 1)(1.1))i.)

####  1.4.2 Case b) $\mathrm{rank}\,(M)=1$: $\psi^{\prime}(\phi)\neq 0$ and
$\mu(Z)=0$

From section 1.4 follows the separation of the scalar potential

$W(r,\phi,Z)=W_{12}(r,\phi)+W_{3}(Z).$ (1.149)

Let us have a look at extended matrix eq. 1.104 for the system of equations
(1.94). For the (non-extended) matrix of the system $M$ to have rank 1, the
auxiliary functions $\tau(\phi)$ and $\sigma(r)$ must be zero, therefore
$W_{3}(Z)$ remains unconstrained, which is in contrast with respect to rank
$(M)=2$. The remaining equations of the system (1.91)–(1.94) are the same as
eq. 1.108 and eq. 1.109, which read

$\displaystyle
r\psi^{\prime}(\phi)W_{r}+\left(r\rho(r)-\psi(\phi)\right)W_{\phi}+\hbar^{2}\frac{\psi^{\prime\prime\prime}(\phi)+\psi^{\prime}(\phi)}{4r^{2}}=0,$
$\displaystyle\psi^{\prime}(\phi)\left(r^{3}\rho^{\prime\prime}(r)-r^{2}\rho^{\prime}(r)+r\rho(r)-3\psi^{\prime\prime}(\phi)-4\psi(\phi)\right)+$
(1.150)
$\displaystyle+\psi^{\prime\prime\prime}(\phi)\left(r\rho(r)-\psi(\phi)\right)-4r^{5}W_{r\phi}-8r^{4}W_{\phi}=0.$

The magnetic field reads

$B^{r}=0,\quad B^{\phi}=0,\quad
B^{Z}=-\frac{1}{2r^{2}}\left(\rho^{\prime}(r)r^{2}+\psi^{\prime\prime}(\phi)+\psi(\phi)\right).$
(1.151)

We have therefore obtained a separable motion – motion in the $xy$-plane under
the influence of the scalar potential $W(r,\phi)$, which has a quantum
correction, plus the perpendicular magnetic field $B^{Z}(r,\phi)$ and the
motion in the direction of the $z$-axis under the influence of the scalar
potential $W_{3}(Z)$, which is not constrained, in contrast with the case of
rank $(M)=2$. The $xy$-plane motion was analysed in [5] and [22]. For both
motions we have one integral of motion in addition to the separable
Hamiltonian, namely integral $X_{1}$ for the motion in the $xy$-plane and
$X_{2}$ for the $z$-direction.

There are 2 systems in 2D that are interesting for us in this case [5]. They
are those in case 1)1.1), so namely with the potential $W_{12}(r,\phi)$ and
magnetic field $B$ from eq. 1.113 and eq. 1.116, the form of the integral
$X_{1}$ can be seen under the cited equations. The $Z$ component of the
potential, $W_{3}(Z)$, remains unconstrained and the corresponding integral of
motion $X_{2}$ will not reduce to a linear one,

$X_{2}=(p_{Z}^{A})^{2}+2W_{3}(Z),$ (1.152)

as can be seen from section 1.4, where all $s^{i}$ vanish.

To sum up: In this section we obtained all quantum quadratically integrable
systems of cylindrical type by solving the corresponding determining
equations. We followed the analysis from the classical case [11] because only
one of the equations differs form the classical case, namely eq. (1.82). The
determining equations were reduced to equations (1.91)–(1.94) containing 5
auxiliary functions of one variable $\rho(r)$, $\sigma(r)$, $\psi(\phi)$,
$\tau(\phi)$ and $\mu(Z)$ which determine the functions $s^{i}$ and $B^{i}$,
see eq. (1.87) and (1.88).

Equation (1.94) contains a matrix which depends on the auxiliary functions
only and allows us to split the considerations according to its rank. In [11]
it was shown that rank 0 and 3 are either impossible (assuming a non-vanishing
magnetic field) or inconsistent and the arguments remain valid in quantum
mechanics. Ranks 1 and 2 split further into subcases a)
$\psi^{\prime}(\phi)=0$ and b) $\psi^{\prime}(\phi)\neq 0$, $\mu(Z)=0$. Case
a) implies vanishing of the quantum correction and the obtained systems are
therefore the same as in classical mechanics. Those were not analysed further
here, and we refer the reader to [11]. (The key results are cited in the
corresponding subsections of Section 2.1, where we search for the first order
superintegrable system among them.) In case b) the quantum correction is _a
priori_ non-trivial, but vanishes due to the consistency conditions on the
scalar potential $W$ in case 1)1.1) subcase 1)(1.1))i. and the corresponding
subcase in this subsection. The analysis shows that the remaining subcases
have the same magnetic field as their classical counterparts and only the
scalar potential $W$ is modified by a $\hbar^{2}$-proportional correction.

In all cases there is at least one free parameter, a constant or even a
function. In the next chapter we search for values of the constants and
functions so that the system admits additional integrals of motion, which
makes the system superintegrable.

##  Chapter 2 Found superintegrable systems

After the analysis of the cylindrical integrable systems both in classical
mechanics [11] and in quantum mechanics (Section 1.4) has been completed, we
go one step further, namely look for superintegrable systems among the
integrable ones.

The goal is to find additional integrals of motion. We use the usual _ansatz_
that the integrals of motion are first or second order polynomials in momenta.
We tried to solve the second order case in general, but we encountered
considerable computational difficulties. Although the third order equations
(1.2) do depend on neither magnetic nor electric field and can be solved, see
eq. (1.20)–(1.23) for the solution, we still have to solve 10 differential
equations (1.2)–(1.17) with 20 arbitrary constants $\alpha_{ij}$ from eq.
(1.20)–(1.23) and 4 unknown functions $s^{r},$ $s^{\phi},$ $s^{Z},$ $m$ of 3
variables $r,$ $\phi,$ $Z$ each in addition to the magnetic field $B$ and
scalar potential $W$ (simplified by assuming integrability).

Despite a lot of effort, we obtained a rather limited results by analysing the
second order case in general, which we present in Subsection 2.2.1. (We
obtained some first order systems here as well). Therefore, we have to be less
ambitious and consider the first order case (Section 2.1) and a physically
motivated second order _ansatz_ , namely $L^{2}+\ldots$ in Subsection 2.2.2
and $L_{x}p_{y}-L_{y}p_{x}+\ldots$ in Subsection 2.2.3. We solve these in a
systematic way, albeit the second order systems in classical mechanics only
because the quantum corrections for these integrals are non-trivial.

###  2.1 First order integrals

####  2.1.1 General considerations

Before we start analysing the cases, let us outline the general procedure and
prove that the first order integrals can contain neither $L_{x}$ nor $L_{y}$.

The task is to solve equations (1.27)–(1.29), where we substitute the magnetic
field from eq. (1.88) with the form of the auxiliary functions corresponding
to the considered case.

We start from the highest order equations (1.27), which do not depend on the
magnetic field nor the scalar potential and can be solved in general. As we
have already mentioned in Section 1.2, the solution (1.30) implies that the
first order terms are in the Euclidean Lie algebra generated by

$p_{x},\ p_{y},\ p_{z},\ L_{x},\ L_{y},\ L_{z},$ (2.1)

the solution therefore depends on $6$ constants $k_{1},\ldots,$ $k_{6}.$ We
associate the constants with the gauge invariant form of generators of the
Euclidean algebra, namely

$Y=k_{1}p_{x}^{A}+k_{2}p_{y}^{A}+k_{3}p_{z}^{A}+k_{4}L_{x}^{A}+k_{5}L_{y}^{A}+k_{6}L_{z}^{A}+m(x,y,z).$
(2.2)

We next solve the first order equations (1.2) with the solution (1.30)
substituted. The best way to proceed is to assume that the function $m$ is
smooth enough and impose the Clairaut compatibility conditions
$\partial_{ab}m=\partial_{ba}m$ on the second derivatives of $m(r,\phi,Z)$.
Taking into account the solution (1.30) to the second order equations (1.27)
and the form of magnetic field in terms of the auxiliary functions (1.88), the
compatibility equations read

$\displaystyle[(k_{4}Z-k_{2})\cos(\phi)+(k_{5}Z+k_{1})\sin(\phi)-k_{6}r](\psi^{\prime\prime\prime}(\phi)+\psi^{\prime}(\phi))+$
$\displaystyle+3[(k_{5}Z+k_{1})\cos(\phi)-(k_{4}Z-k_{2})\sin(\phi)]\psi^{\prime\prime}(\phi)+$
$\displaystyle+[(k_{5}Z+k_{1})\cos(\phi)-(k_{4}Z-k_{2})\sin(\phi)]r^{2}(\rho^{\prime}(r)-r\rho^{\prime\prime}(r))+$
(2.3)
$\displaystyle+[-3k_{5}r\cos(\phi)+3k_{4}r\sin(\phi)+2k_{3}]r^{4}\mu^{\prime}(Z)+$
$\displaystyle+r(k_{5}\cos(\phi)-k_{4}\sin(\phi))\tau^{\prime}(\phi)-r^{4}(k_{4}\cos(\phi)+k_{5}\sin(\phi))\sigma^{\prime}(r)+$
$\displaystyle+3[(k_{5}Z+k_{1})\cos(\phi)-(k_{4}Z-k_{2})\sin(\phi)]\psi(\phi)-$
$\displaystyle-2[k_{4}\cos(\phi)+k_{5}\sin(\phi)]r\tau(\phi)$
$\displaystyle=0,$
$\displaystyle[-(k_{5}Z+k_{1})\cos(\phi)+(k_{4}Z-k_{2})\sin(\phi)]r^{4}\sigma^{\prime\prime}(r)+$
$\displaystyle+(k_{4}\cos(\phi)+k_{5}\sin(\phi))r\psi^{\prime\prime}(\phi)+$
$\displaystyle+[3(k_{4}Z-k_{2})\cos(\phi)+3(k_{5}Z+k_{1})\sin(\phi)-2k_{6}r]\tau^{\prime}(\phi)-$
$\displaystyle-((k_{4}Z-k_{2})\cos(\phi)+\sin(\phi)(k_{5}Z+k_{1}))r^{4}\mu^{\prime}(Z)+$
(2.4) $\displaystyle+[k_{4}\cos(\phi)+k_{5}\sin(\phi)]r^{3}\rho^{\prime}(r)+$
$\displaystyle+6[(k_{5}Z+k_{1})\cos(\phi)-(k_{4}Z-k_{2})\sin(\phi)]\tau(\phi)-$
$\displaystyle-[k_{4}\cos(\phi)+k_{5}\sin(\phi)]r(2r^{3}\mu(Z)-\psi(\phi))$
$\displaystyle=0,$
$\displaystyle[(k_{4}Z-k_{2})\cos(\phi)+(k_{5}Z+k_{1})\sin(\phi)-k_{6}r]\tau^{\prime\prime}(\phi)+$
$\displaystyle+(-k_{5}r\cos(\phi)+k_{4}r\sin(\phi)+k_{3})r^{5}\mu^{\prime\prime}(Z)-$
$\displaystyle-(k_{5}\cos(\phi)-k_{4}\sin(\phi))r\psi^{\prime\prime}(\phi)+$
$\displaystyle+3[(k_{5}Z+k_{1})\cos(\phi)-\sin(\phi)(k_{4}Z-k_{2})]\tau^{\prime}(\phi)-$
$\displaystyle-[(k_{4}Z-k_{2})\cos(\phi)+(k_{5}Z+k_{1})\sin(\phi)]r^{3}\sigma^{\prime}(r)+$
(2.5)
$\displaystyle+[(k_{5}Z+k_{1})\cos(\phi)-(k_{4}Z-k_{2})\sin(\phi)]r^{4}\mu^{\prime}(Z)-$
$\displaystyle-4(k_{5}\cos(\phi)-k_{2}\sin(\phi))(r^{3}\rho^{\prime}(r))-$
$\displaystyle-2[(k_{4}Z+k_{2})\cos(\phi)+(k_{5}Z+k_{1})\sin(\phi)]\tau(\phi)-$
$\displaystyle-[k_{4}\sin(\phi)-k_{5}\cos(\phi)]r(2r^{3}\mu(Z)-\psi(\phi))$
$\displaystyle=0.$

(We have multiplied the equations by the power of $r$ which was originally in
the denominator. Because the equations must be satisfied for all $r$, $\phi$,
$Z$, we implicitly assume $r\neq 0$ and we do these simplifications
automatically, often without mentioning them.)

Although the equations are long, they are manageable because they can be
separated by taking derivatives with respect to $r$ and/or $Z$ (and
dividing/multiplying by $r$).

There is one thing that can be proved using the equations above only, i.e. for
all integrable cases: If we assume that the magnetic field is real and does
not vanish, these systems do not admit the first order integrals corresponding
to the constants $k_{4}$, $k_{5}$, namely $L_{x}$ and $L_{y}$. We prove this
fact by contradiction and at the same time illustrate the procedure used in
the following subsections, which consists in repeated solving of equations
with the highest powers of $Z$ and $r$ ($\mu(Z)$ reduces to a polynomial) and
substituting the results.

To prove the previous assertion, we assume that at least one of the constants
$k_{4}$ and $k_{5}$ does not vanish and we show it implies zero or complex
magnetic field. Taking the second derivative with respect to $Z$ of equation
(2.3), we obtain (dividing by $r^{4}$)

$[-3k_{5}r\cos(\phi)+3k_{4}r\sin(\phi)+2k_{3}]\mu^{\prime\prime\prime}(Z)=0.$
(2.6)

Differentiating it with respect to $r$, our assumption together with the fact
that the equation must be satisfied for all $r,$ $\phi,$ $Z$ implies

$\mu(Z)=\mu_{2}Z^{2}+\mu_{1}Z+\mu_{0}.$ (2.7)

Inserting it into eq. (2.4)–(2.5) and differentiating twice with respect to
$Z$, we obtain

$k_{5}\mu_{2}\sin(\phi)+k_{4}\mu_{2}\cos(\phi)=0,\quad
k_{5}\mu_{2}\sin(\phi)-k_{4}\mu_{2}\cos(\phi)=0,$ (2.8)

which imply linearity of $\mu(Z)$, $\mu_{2}=0$.

Using the result in eq. (2.3) differentiated with respect to $Z$ and $r$, we
obtain

$k_{5}(r^{3}\rho^{\prime\prime\prime}(r)+2r^{2}\rho^{\prime\prime}(r)-2r\rho^{\prime}(r))\cos(\phi)-k_{4}(r^{3}\rho^{\prime\prime\prime}(r)+2r^{2}\rho^{\prime\prime}(r)-2r\rho^{\prime}(r))\sin(\phi)=0,$
(2.9)

which has the following solution

$\rho(r)=\rho_{0}+\frac{\rho_{1}}{r}+\rho_{2}r^{2}.$ (2.10)

Differentiating eq. (2.4) with respect to $Z$ and $r$, we get

$(k_{5}r\sigma^{\prime\prime\prime}(r)+4k_{5}\sigma^{\prime\prime}(r)+12k_{4}\mu_{1})\cos(\phi)-(k_{4}r\sigma^{\prime\prime\prime}(r)+4k_{4}\sigma^{\prime\prime}(r)-12k_{5}\mu_{1})\sin(\phi)=0.$
(2.11)

Solving the brackets separately, the solutions are consistent only if
$\mu_{1}=0$ and read

$\sigma(r)=\frac{\sigma_{2}}{r^{2}}+\sigma_{1}r+\sigma_{0}.$ (2.12)

Inserting $\sigma(r)$ in into eq. (2.5), we get

$\sigma_{1}\left(k_{4}\cos(\phi)+k_{5}\sin(\phi)\right)=0,$ (2.13)

so we see that $\sigma_{1}=0$.

We now solve eq. (2.3)–(2.5) differentiated with respect to $Z$ only to obtain
$\psi(\phi)$ and $\tau(\phi)$ with the result

$\displaystyle\psi(\phi)$
$\displaystyle=\rho_{1}+\frac{c_{1}\cos(2\phi)+c_{2}\sin(2\phi)+c_{3}}{2(k_{4}\cos(\phi)+k_{5}\sin(\phi))},$
(2.14) $\displaystyle\tau(\phi)$
$\displaystyle=\sigma_{2}+\frac{2c_{4}-(k_{4}^{2}+k_{5}^{2})}{2(k_{4}\cos(\phi)+k_{5}\sin(\phi))^{2}}.$
(2.15)

With this solution, the coefficient of $Z$ in equations (2.3)–(2.5) vanishes,
but the equations are not satisfied yet. Differentiating eq. (2.4) with
respect to $r$ four times, we obtain

$(\mu_{0}-\rho_{2})(k_{4}\cos(\phi)+k_{5}\sin(\phi))=0,$ (2.16)

therefore $\rho_{2}=\mu_{0}$.

Inserting this result, we differentiate eq. (2.3)–(2.5) once with respect to
$r$, take the numerators only and collect the coefficients of the polynomial
in sines and cosines. Thus, we obtain a system of 7 algebraic equations:

(2.17) (2.18) (2.19) (2.20) (2.21) (2.22) (2.23)

It has 3 real solutions excluding $k_{4}=k_{5}=0$, namely

  1. 1.

$c_{1}=c_{3},\quad c_{4}=\frac{\sigma_{2}}{2}k_{4}^{2},\quad k_{5}=0,$

  2. 2.

$c_{1}=-c_{3},\quad c_{4}=\frac{\sigma_{2}}{2}k_{5}^{2},\quad k_{4}=0,$

  3. 3.

$c_{4}=\frac{\sigma_{2}}{2}(k_{4}^{2}+k_{5}^{2}),\quad
c_{2}=\frac{(-c_{1}+c_{3})k_{4}^{2}+(c_{1}+c_{3})k_{5}^{2}}{2k_{5}k_{4}}.$

(The remaining constants are unconstrained in all three cases.) All of them,
however, lead to vanishing of the magnetic field. Therefore, we have proved
that in order to have real non-vanishing magnetic field we need

$k_{4}=k_{5}=0.$ (2.24)

The remaining part of equations (2.3)–(2.5) can be solved for each integrable
system of cylindrical type in the same manner. In each step we have to
consider which of the constants $k_{i}$ from solution (1.30) should be set to
zero and which non-zero, because any vanishing constant means that the
corresponding integral is not allowed by the system, but non-vanishing
constants constrain auxiliary functions in such a way that the corresponding
systems have zero magnetic field.

Having solved the compatibility conditions (2.3)–(2.5), we solve eq. (1.2) and
eq. (1.29) to obtain $m(r,\phi,Z)$ and $W(r,\phi,Z)$, respectively. We have to
check compatibility of the obtained scalar potential $W(r,\phi,z)$ with the
restrictions imposed by the considered cylindrical system.

Let us now treat the specific cases from [11] and Section 1.4.

####  2.1.2 Case $\mathrm{rank}(M)=2$ subcase 2a)

The assumptions of this subcase read

$\psi(\phi)=0,\quad\mu(Z)=\mu_{0},\quad\tau(\phi)=\tau_{0},$ (2.25)

with the corresponding magnetic field and scalar potential

$W=W(r),\quad B^{r}=0,\quad
B^{\phi}=\frac{\tau_{0}}{r^{3}}+\frac{1}{2}\sigma^{\prime}(r),\quad
B^{Z}=\mu_{0}r-\frac{1}{2}\rho^{\prime}(r).$ (2.26)

In this case we assume $\psi(\phi)=0$, so the quantum correction vanishes and
therefore the system is the same in classical and quantum mechanics.

From [11] we know, that this system admits the integrals

$\tilde{X}_{1}=p_{\phi}^{A}+\frac{\rho(r)}{2}-\frac{\mu_{0}r^{2}}{2},\quad\tilde{X}_{2}=p_{Z}^{A}+\frac{\sigma(r)}{2}-\frac{\tau_{0}}{2r^{2}},$
(2.27)

so we set $k_{3}=k_{6}=0$.

The additional conditions which split the case into further subcases in [11]
all imply the scalar potential $W$ and magnetic field $B$ from eq. (2.26) with
different constraints on the functions and constants. However, we lose no
generality by redefining $\sigma(r)$ and $\rho(r)$ to set $\mu_{0}=\tau_{0}=0$
and proceeding with the assumption $\sigma(r)\rho(r)\neq 0$. (We can always
choose the constant term which does not appear in the magnetic field and can
be eliminated by the choice of gauge in the integrals.) Therefore, we do not
need to split into subcases.

Recalling eq. (2.24), the compatibility conditions (2.3)–(2.5) can be solved
at once with the result

$\rho(r)=\rho_{1}r^{2}+\rho_{0},\quad\sigma(r)=\sigma_{0},$ (2.28)

where $\rho_{1}\neq 0$ because otherwise the corresponding magnetic field

$B^{r}(r,\phi,Z)=0,\quad B^{\phi}(r,\phi,Z)=0,\quad
B^{Z}(r,\phi,Z)=-r\rho_{1}$ (2.29)

vanishes. We write the magnetic field also in the Cartesian coordinates
(transformed using eq. (1.2))

$B^{x}(x,y,z)=0,\quad B^{y}(x,y,z)=0,\quad B^{z}(x,y,z)=-\rho_{1}.$ (2.30)

We can now solve eq. (1.2) to get

$m=\rho_{1}r(-k_{1}\sin(\phi)+k_{2}\cos(\phi))-k_{3}\frac{\sigma_{0}}{2}-k_{6}\frac{\rho_{0}}{2}.$
(2.31)

(The last term is the integration constant, which we choose so that
$\tilde{X}_{1}$ and $\tilde{X}_{2}$ do not contain unnecessary constants.) The
zeroth order equation (1.29) becomes

$\left(k_{1}\cos(\phi)+k_{2}\sin(\phi)\right)W^{\prime}(r)=0,$ (2.32)

we thus see that in order to have additional integrals (one of the constants
$k_{1},$ $k_{2}$ non-zero) we need $W$ constant, without loss of generality

$W=0.$ (2.33)

This is a well-known system, considered e.g. in [15]. It was considered in
[20, eq. (42)], with different choice of the frame of reference. The integrals
read (in the Cartesian coordinates)

$Y_{1}=p_{x}^{A}-\rho_{1}y,\quad
Y_{2}=p_{y}^{A}+\rho_{1}x,\quad\tilde{X}_{1}=L_{z}^{A}+\frac{\rho_{1}(x^{2}+y^{2})}{2},\quad\tilde{X}_{2}=p_{z}^{A},$
(2.34)

where $L_{z}^{A}=xp_{y}^{A}-yp_{x}^{A}$. The integrals are clearly mutually
independent but the Hamiltonian can be rewritten as

$H=\frac{1}{2}\left(Y_{1}^{2}+Y_{2}^{2}+\tilde{X}_{2}^{2}\right)-\rho_{1}X_{4}=\frac{1}{2}\left(\left(p_{x}^{A}\right)^{2}+\left(p_{y}^{A}\right)^{2}+\left(p_{z}^{A}\right)^{2}\right).$
(2.35)

The fifth independent integral, which makes the system maximally
superintegrable, is not polynomial in momenta [20], namely

$X_{5}=(p_{y}^{A}+\rho_{1}x)\sin\left(\frac{\rho_{1}z}{p_{z}^{A}}\right)-(p_{x}^{A}-\rho_{1}y)\cos\left(\frac{\rho_{1}z}{p_{z}^{A}}\right),$
(2.36)

where we assume that $p_{z}^{A}$. (Otherwise the system collapses to 2D.)

Let us follow [20] and solve the Hamiltonian equations of motion and the
Schrödinger equation. We choose the gauge so that two of the integrals
$Y_{1},$ $Y_{2},$ $\tilde{X}_{2}$ are simply momenta

$\vec{A}(\vec{x})=\left(0,-\rho_{1}x,0\right).$ (2.37)

The Hamiltonian equations of motion read

$\displaystyle\dot{x}=p_{x},\quad\dot{y}=p_{y}+\rho_{1}x,\quad\dot{z}=p_{z},$
(2.38)
$\displaystyle\dot{p}_{x}=-\rho_{1}(p_{y}+\rho_{1}x),\quad\dot{p}_{y}=0,\quad\dot{p}_{z}=0.$
(2.39)

The solution to these equations with the usual initial conditions
$x(0)=x_{0},$ $p_{x}(0)=p_{x}^{0}$ etc. is

$\displaystyle\begin{aligned}
x(t)&=\frac{p_{x}^{0}}{\rho_{1}}\sin(\rho_{1}t)+\left(x_{0}+\frac{p_{y}^{0}}{\rho_{1}}\right)\cos(\rho_{1}t)-\frac{p_{y}^{0}}{\rho_{1}},\\\
y(t)&=\left(x_{0}+\frac{p_{y}^{0}}{\rho_{1}}\right)\sin(\rho_{1}t)-\frac{p_{x}^{0}}{\rho_{1}}\cos(\rho_{1}t)+y_{0}+\frac{p_{x}^{0}}{\rho_{1}},\\\
z(t)&=p_{z}^{0}t+z_{0},\\\
p_{x}(t)&=-\left(\rho_{1}x_{0}+p_{y}^{0}\right)\sin(\rho_{1}t)+p_{x}^{0}\cos(\rho_{1}t),\quad
p_{y}(t)=p_{y}^{0},\quad p_{z}(t)=p_{z}^{0}.\end{aligned}$ (2.40)

If $p_{z}^{0}\neq 0$, the trajectory is a helix with its axis parallel to the
$z$-axis, otherwise it collapses to a circle in the plane $z=z_{0}$.

Let us now turn to the Schrödinger equation. Writing the Hamiltonian (2.35) in
full

$\hat{H}=-\frac{\hbar^{2}}{2}\left(\partial_{xx}+\partial_{yy}+\partial_{zz}\right)-\mathrm{i}\hbar\rho_{1}x\partial_{y}+\frac{\rho_{1}^{2}}{2}x^{2},$
(2.41)

we see that the stationary Schrödinger equation $\hat{H}\psi=E\psi$ separates
in the Cartesian coordinates

$\displaystyle\psi(\vec{x})=f(x)\exp\left(\frac{\mathrm{i}}{\hbar}\lambda_{2}y\right)\exp\left(\frac{\mathrm{i}}{\hbar}\lambda_{3}z\right),$
(2.42)
$\displaystyle\hbar^{2}f^{\prime\prime}(x)=\left((\rho_{1}x-\lambda_{2})^{2}+\lambda_{3}^{2}-2E\right)f(x),$
(2.43) $\displaystyle
Y_{2}\psi(\vec{x})=\lambda_{2}\psi(\vec{x}),\quad\tilde{X}_{2}\psi(\vec{x})=\lambda_{3}\psi(\vec{x}).$
(2.44)

The reduced Schrödinger equation corresponds to the 1D harmonic oscillator
with energy $E-\frac{\lambda^{2}_{3}}{2}$, angular frequency $\omega=\rho_{1}$
with the centre of force located at $x=\frac{\lambda_{2}}{\rho_{1}}$.
Therefore, the energy spectrum of the Hamiltonian (2.35)

$E=\frac{\lambda_{3}^{2}}{2}+\hbar\rho_{1}\left(n+\frac{1}{2}\right),$ (2.45)

is continuous due to the arbitrary momentum $\lambda_{3}$ and the generalized
eigenvectors are

$\psi_{n,\lambda_{2},\lambda_{3}}(\vec{x})=K_{n}\exp\left(\frac{\mathrm{i}}{\hbar}(\lambda_{2}y+\lambda_{3}z)\right)H_{n}\left(\sqrt{\frac{\rho_{1}}{\hbar}}\left(x-\frac{\lambda_{2}}{\rho_{1}}\right)\right)\exp\left(-\frac{\rho_{1}}{2\hbar}\left(x-\frac{\lambda_{2}}{\rho_{1}}\right)^{2}\right),$
(2.46)

where $H_{n}$ are the Hermite polynomials and $K_{n}$ is their normalization
constant.

Let us see if the Schrödinger equation separates in the cylindrical
coordinates as well. Using the gauge (2.37), the Hamiltonian reads

$\begin{split}\hat{H}={}&-\frac{\hbar^{2}}{2}\left(\frac{1}{r}\partial_{r}(r\partial_{r})+\frac{1}{r^{2}}\partial_{\phi\phi}+\partial_{ZZ}\right)-\\\
&-\mathrm{i}\hbar\rho_{1}r\cos(\phi)\left(\sin(\phi)\partial_{r}+\frac{\cos(\phi)}{r}\partial_{\phi}\right)+\frac{\rho_{1}}{2}r^{2}(\cos(\phi))^{2}.\end{split}$
(2.47)

It is clear that we can separate the $Z$ coordinate from the $r$ and $\phi$.
However, it is not possible to separate $r$ and $\phi$: The equation to
separate is (divided by $f(r)g(\phi)$)

$\begin{split}-\frac{\hbar^{2}}{2}\left(\frac{g^{\prime\prime}(\phi)}{g(\phi)}-\frac{f^{\prime\prime}(r)}{f(r)}\right)-\mathrm{i}\hbar\rho_{1}r\sin(\phi)\cos(\phi)\frac{f^{\prime}(r)}{f(r)}=\\\
=\mathrm{i}\hbar\rho_{1}\cos(\phi)^{2}\frac{g^{\prime}(\phi)}{g(\phi)}-\frac{\rho_{1}^{2}}{2}r^{2}\cos(\phi)^{2}+E,\end{split}$
(2.48)

where the last term on the left-hand side and the second term on the right-
hand side will always depend both on $r$ an $\phi$.

However, if we choose the gauge which makes $\tilde{X}_{1}$ simply $L_{z}$
[6], namely

$\vec{A}(\vec{x})=\left(\frac{\rho_{1}y}{2},-\frac{\rho_{1}x}{2},0\right),\qquad
A_{r}=0,\quad A_{\phi}=-\frac{\rho_{1}r^{2}}{2},\quad A_{Z}=0,$ (2.49)

the Hamiltonian reads

$\begin{split}\hat{H}={}&-\frac{\hbar^{2}}{2}\left(\frac{1}{r}\partial_{r}(r\partial_{r})+\frac{1}{r^{2}}\partial_{\phi\phi}+\partial_{ZZ}\right)+\mathrm{i}\hbar\frac{\rho_{1}}{2}\partial_{\phi}+\frac{\rho_{1}^{2}r^{2}}{8}.\end{split}$
(2.50)

This time $L_{z}=-\mathrm{i}\hbar\partial_{\phi}$ commutes with the
Hamiltonian and therefore the wave-function

$\psi(r,\phi,Z)=f(r)\exp\left(\frac{\mathrm{i}}{\hbar}\eta_{2}\phi\right)\exp\left(\frac{\mathrm{i}}{\hbar}\eta_{3}Z\right)$
(2.51)

leads to a separation of variables, with the $r$ coordinate equation reading

$-\hbar^{2}\left(f^{\prime\prime}(r)+\frac{1}{r}f^{\prime}(r)\right)+\left(\frac{\eta_{2}^{2}}{r^{2}}+\eta_{3}^{2}-\rho_{1}\eta_{2}\right)f(r)+\frac{\rho_{1}^{2}r^{2}f(r)}{4}=2Ef(r).$
(2.52)

Its solution in terms of Whittaker functions $M_{\mu,\nu}(z),$
$W_{\mu,\nu,}(z)$ [1, Chapter 13], [7, Section 13.14] reads

$f(r)=\frac{1}{r}\left(c_{1}M_{\mu,\nu}\left(\frac{\rho_{1}}{\sqrt{2}\hbar}r^{2}\right)+c_{2}M_{\mu,\nu}\left(\frac{\rho_{1}}{\sqrt{2}\hbar}r^{2}\right)\right),$
(2.53)

where $\mu$, $\nu$ are constants

$\mu=\frac{(\rho_{1}\eta_{2}-\eta_{3}^{2}+2E)}{\sqrt{2}\hbar\rho_{1}},\quad\nu=\frac{\eta_{2}}{\sqrt{2}\hbar}.$
(2.54)

We solve the stationary Hamilton-Jacobi equation as well. Since we work with
time-independent Hamiltonians, we can separate the time coordinate $t$ from
the Hamilton’s principal function $S$ with separation constant $-E$ and we use
Hamilton’s characteristic function $U$, namely $S=U-Et$. Using the Hamiltonian
(2.35) and the gauge (2.37) the equation reads

$\frac{1}{2}\left[\left(\frac{\partial U}{\partial
x}\right)^{2}+\left(\frac{\partial U}{\partial
y}-\rho_{1}x\right)^{2}+\left(\frac{\partial U}{\partial
z}\right)^{2}\right]=E.$ (2.55)

We have 2 cyclic coordinates $y$ and $z$; thus, we have the separation

$U(x,y,z)=u(x)+p_{y}y+p_{z}z,$ (2.56)

with $u(x)$ satisfying

$(u^{\prime}(x))^{2}+(p_{y}-\rho_{1}x)^{2}+p_{z}^{2}=2E.$ (2.57)

For later use let us consider more generally an equation of type

$(u^{\prime}(x))^{2}+\alpha_{1}x^{2}+\alpha_{2}x+\alpha_{3}=0$ (2.58)

whose solution is

$\begin{split}u(x)={}&\pm\left[\left(\frac{\alpha_{2}^{2}}{8\alpha_{1}^{3/2}}-\frac{\alpha_{3}}{2\sqrt{\alpha_{1}}}\right)\arctan\left(\frac{2\alpha_{1}x+\alpha_{2}}{2\sqrt{\alpha_{1}}\sqrt{-\alpha_{1}x^{2}-\alpha_{2}x-\alpha_{3}}}\right)\right.+\\\
&+\left.\left(\frac{x}{2}+\frac{\alpha_{2}}{4\alpha_{1}}\right)\sqrt{-\alpha_{1}x^{2}-\alpha_{2}x-\alpha_{3}}\right]+c_{1}.\end{split}$
(2.59)

In particular, the Hamilton’s characteristic function reads in our case

$\begin{split}U={}&\frac{(p_{z}^{2}-2E)}{2\rho_{1}}\arctan\left(\frac{(-\rho_{1}x+p_{y})}{\sqrt{-\rho_{1}^{2}x^{2}+2p_{y}\rho_{1}x+a}}\right)-\\\
&-\frac{(\rho_{1}x-p_{y})}{2\rho_{1}}\sqrt{-\rho_{1}^{2}x^{2}+2p_{y}\rho_{1}x+a}+p_{y}y+p_{z}z,\end{split}$
(2.60)

where $a=2E-p_{y}^{2}-p_{z}^{2}$.

In the cylindrical coordinates we use the gauge (2.49), in which the
Hamiltonian does not depend on $\phi$ and $Z$, therefore Hamilton’s
characteristic function

$U(r,\phi,Z)=u(r)+p_{\phi}\phi+p_{Z}Z$ (2.61)

reduces the Hamilton-Jacobi equation to

$\frac{1}{2}\left[u^{\prime}(r)^{2}+\frac{(2p_{\phi}-\rho_{1}r^{2})^{2}}{4r^{2}}+p_{Z}^{2}\right]=E.$
(2.62)

Let us again consider more general form of the equation for later use, namely

$(u^{\prime}(r))^{2}=\frac{ar^{4}+br^{2}+c}{r^{2}}.$ (2.63)

Assuming we take the square of a positive real number, the result to the last
equation is of the following type ($v=r^{2}$)

$\begin{split}u(v)=\int\\!{\frac{\sqrt{a{v}^{2}+bv+c}}{v}}{\rm
d}v=&\sqrt{a{v}^{2}+bv+c}+{\frac{b}{2\sqrt{a}}\ln\left({\frac{b+2av}{2\sqrt{a}}}+\sqrt{a{v}^{2}+bv+c}\right)}\\\
&-\sqrt{c}\arctan\left(\frac{bv+2c}{2\sqrt{c}\sqrt{a{v}^{2}+bv+c}}\right)+C.\end{split}$
(2.64)

The integration constant can absorb the purely imaginary constants arising
from the possibly negative value of $a$ in the second term, so we can use this
form of the solution in that case as well with $a$ replaced by $|a|$. If
$c<0$, we must change the third term to argtanh and $c$ to $|c|$. The explicit
expression for the solution of the Hamilton-Jacobi equation becomes too
complicated for practical use and we have already calculated the trajectories.
Thus, we do not present it here.

####  2.1.3 Case $\mathrm{rank}(M)=2$ subcase 2b 1.1)

This case corresponds to

$\tau(\phi)=\tau_{0},\quad\sigma(r)=\frac{\tau_{0}}{r^{2}}+\sigma_{0},\quad\mu(Z)=0,\quad\psi^{\prime}(\phi)\neq
0,$ (2.65)

which yields the following magnetic field

$B^{\phi}=0,\quad B^{r}=0,\quad
B^{Z}=-\rho^{\prime}(r)-\frac{\psi^{\prime\prime}(\phi)+\psi(\phi)}{2r^{2}}$
(2.66)

and the separable scalar potential with vanishing $Z$ component
$W(r,\phi,Z)=W_{12}(r,\phi)\equiv W(r,\phi)$ satisfying eq. (1.108) and eq.
(1.109).

Recalling eq. (2.24), we have only one non-trivial compatibility equation,
namely eq. (2.3). Differentiating it twice with respect to $r$ to eliminate
$\psi(\phi)$, we get an equation for $\rho(r)$

$(k_{1}\cos(\phi)+k_{2}\sin(\phi))[r^{3}\rho^{(4)}(r)+5r^{2}\rho^{\prime\prime\prime}(r)+2r\rho^{\prime\prime}(r)-2\rho^{\prime}(r)]=0.$
(2.67)

Because the vanishing of the constants $k_{1},$ $k_{2}$ would mean no
additional integrals of motion (we would the second order cylindrical
integrals or their reduced first order form only), we need $\rho(r)$ of the
form

$\rho(r)=\rho_{3}\ln(r)+\rho_{2}r^{2}+\frac{\rho_{1}}{r}+\rho_{0}.$ (2.68)

Using this in eq. (2.3) once differentiated with respect to $r$, we obtain

$k_{6}(\psi^{\prime\prime\prime}(\phi)+\psi^{\prime}(\phi))-2\rho_{3}(k_{2}\sin(\phi)+k_{1}\cos(\phi))=0,$
(2.69)

which splits the considerations into two subcases: $k_{6}\neq 0$ and
$k_{6}=0$.

  1. I.

$k_{6}\neq 0$: In this subcase we solve eq. (2.69) for $\psi(\phi)$ and get

$\psi(\phi)=-\frac{1}{k_{6}}[((k_{1}\phi+k_{2})\rho_{3}+\psi_{2}k_{6})\cos(\phi)+((k_{2}\phi-
k_{1})\rho_{3}-\psi_{1}k_{6})\sin(\phi)]+\psi_{3}.$ (2.70)

Using this in eq. (2.3) (without differentiation), we get an algebraic
equation with only one solution (excluding $k_{1}=k_{2}=0$):
$\psi_{3}=\rho_{1}$ and $\rho_{3}=0$, i.e.

$\psi(\phi)=-\psi_{2}\cos(\phi)+\psi_{1}\sin(\phi)+\rho_{1},\quad\rho(r)=\frac{\rho_{2}r^{3}+\rho_{1}}{r}+\rho_{0},$
(2.71)

so we obtained the constant magnetic field

$B^{r}=0,\quad B^{\phi}=0,\quad B^{Z}=-\rho_{2}r.$ (2.72)

Translating this into the Cartesian coordinates using eq. (1.2), it reads

$B^{x}=0,\quad B^{y}=0,\quad B^{z}=-\rho_{2}.$ (2.73)

We note that $\psi(\phi)$ from eq. (2.71) implies vanishing of the quantum
correction, so the quantum version of the system is identical to the classical
version in this subcase.

Having solved the compatibility conditions, we solve the original equations
(1.2) and (1.29) to obtain

$\displaystyle W(r,\phi)$
$\displaystyle=W\left(\frac{k_{6}r^{2}}{2}+k_{2}r\cos(\phi)-k_{1}r\sin(\phi)\right),$
(2.74) $\displaystyle m(r,\phi,Z)$
$\displaystyle=\rho_{2}\left(k_{2}r\cos(\phi)-k_{1}r\sin(\phi)+\frac{k_{6}r^{2}}{2}\right).$
(2.75)

(We set the integration constant to 0 in $m$ so that we get
$\tilde{X}_{1}=p_{Z}^{A}$ in accordance with eq. (2.34) in Subsection 2.1.2.)

The potential must, however, satisfy eq. (1.108) and eq. (1.109), where the
quantum correction vanishes. Inserting the obtained form of the scalar
potential into the first of them, we get

$W^{\prime}\cdot[(k_{2}\rho_{2}\sin(\phi)+k_{1}\rho_{2}\cos(\phi))r^{3}+((k_{2}\rho_{0}-k_{6}\psi_{2})\sin(\phi)+\\\
+(k_{1}\rho_{0}-k_{6}\psi_{1})\cos(\phi))r+k_{1}\psi_{2}-k_{2}\psi_{1}]=0.$
(2.76)

If we assume that the scalar potential $W$ is not constant, the terms in
square bracket must vanish separately for all powers of $r$. The $r^{3}$ term
in the bracket implies either $k_{1}=k_{2}=0$, which means only 2 reduced
cylindrical integrals, or $\rho_{2}=0$ and vanishing magnetic field. Thus, the
only interesting case for us is the constant scalar potential

$W(r,\phi)=W_{0}.$ (2.77)

Equation (1.109) is now also satisfied. Thus, we have arrived at the system
(2.34) and we refer the reader to Subsection 2.1.2.

  2. II.

$k_{6}=0$: With this assumption eq. (2.68) implies $\rho_{3}=0$ and together
with eq. (2.3) we get

$\displaystyle\rho(r)$
$\displaystyle=\frac{\rho_{2}r^{3}+\rho_{1}}{r}+\rho_{0},$ (2.78)
$\displaystyle\begin{split}\psi(\phi)&=\frac{c_{1}\cos(2\phi)+c_{2}\sin(2\phi)+c_{3}}{k_{2}\cos(\phi)-k_{1}\sin(\phi)}+\rho_{1}.\end{split}$
(2.79)

The corresponding magnetic field $B$ reads

$B^{r}=0,\quad B^{\phi}=0,\quad
B^{Z}=-\rho_{2}r+\frac{\xi}{r^{2}(k_{1}\sin(\phi)-k_{2}\cos(\phi))^{3}},$
(2.80)

where $\xi=-(c_{1}+c_{3})k_{1}^{2}-2c_{2}k_{1}k_{2}+k_{2}^{2}(c_{1}-c_{3})$.
Translating this into the Cartesian coordinates using eq. (1.2), it reads

$B^{x}=0,\quad B^{y}=0,\quad B^{z}=-\rho_{2}+\frac{\xi}{(k_{1}y-k_{2}x)^{3}}.$
(2.81)

In order to simplify the equations, we use the formula

$k_{1}\sin(\phi)-k_{2}\cos(\phi)=\sqrt{k_{1}^{2}+k_{2}^{2}}\cos(\phi+\phi_{0}),\quad\cos(\phi_{0})=-\frac{k_{2}}{\sqrt{k_{1}^{2}+k_{2}^{2}}},$
(2.82)

rotate the coordinate system so that $\phi_{0}=0$ and redefine the constants
$c_{i}$ in $\psi(\phi)$ to get

$\psi(\phi)=\frac{c_{1}\cos(2\phi)+c_{2}\sin(2\phi)+c_{3}}{\cos(\phi)}+\rho_{1}.$
(2.83)

The corresponding magnetic field $B$ reads in the cylindrical coordinates

$B^{r}=0,\quad B^{\phi}=0,\quad
B^{Z}=-\rho_{2}r+\frac{c_{1}-c_{3}}{r^{2}(\cos(\phi))^{3}}$ (2.84)

and in the Cartesian coordinates

$B^{x}=0,\quad B^{y}=0,\quad B^{z}=-\rho_{2}+\frac{c_{1}-c_{3}}{x^{3}}.$
(2.85)

Having solved the compatibility conditions, we solve the original equations
(1.2) and (1.29) to obtain

$\displaystyle W={}$ $\displaystyle W(r\cos(\phi))=W(x),$ (2.86)
$\displaystyle m={}$
$\displaystyle\rho_{2}r-\frac{c_{1}-c_{3}}{2r^{2}(\cos(\phi))^{2}}=\rho_{2}\sqrt{x^{2}+y^{2}}-\frac{c_{1}-c_{3}}{2x^{2}}.$
(2.87)

(We omit the constant of integration, which is only added to the integral of
motion.)

We need to satisfy 2 more equations to ensure integrability, namely eq.
(1.108) and eq. (1.109). The first of them is

$\begin{split}&\cos(\phi)^{3}r^{3}\left[((\rho_{2}r^{3}+\rho_{0}r)\cos(\phi)+2(c_{1}-c_{3}))\sin(\phi)-2c_{2}\cos(\phi)\right]W^{\prime}(r\cos(\phi))\\\
&+\tfrac{3}{2}\hbar^{2}(c_{1}-c_{3})\sin(\phi)=0.\end{split}$ (2.88)

(We differentiate $W$ with respect to its argument, i.e. $r\cos(\phi)$.)

In the classical case, $\hbar\to 0$, the conditions for non-vanishing magnetic
field $c_{1}\neq c_{3}$ and $\rho_{2}\neq 0$ imply constant scalar potential
$W$.

In the quantum case, there is another possibility in addition to $W=W_{0}$: If
$c_{1}\neq c_{3}$, we must have

$\rho_{2}=0,\quad\rho_{0}=0,\quad c_{2}=0,\quad
W=-\frac{3\hbar^{2}}{4(r\cos(\phi))^{2}}=-\frac{3\hbar^{2}}{4x^{2}}.$ (2.89)

Inserting this potential $W$ into eq. (1.109) with $\rho(r)$ from eq. (2.78)
and $\psi$ from eq. (2.83), the numerator becomes

$(c_{1}-c_{3})[(\rho_{2}r^{3}+\rho_{0}r)\sin(\phi)\cos(\phi)+2(c_{1}-c_{3})\sin(\phi)-2c_{2}\cos(\phi)]=0.$
(2.90)

Thus, we have a contradiction due to $c_{1}\neq c_{3}$ and the potential must
be constant in quantum mechanics as well.

If we insert constant potential $W=W_{0}$ into eq. (1.109), it is satisfied if
and only if $c_{1}=c_{3}$, because we obtain the same numerator as in eq.
(2.90). The result is, therefore, that both the magnetic field $B^{z}$ and
scalar potential $W$ must be constant. We have again the well-known maximally
superintegrable system from [15, 20], see Subsection 2.1.2.

####  2.1.4 Case $\mathrm{rank}(M)=2$ subcase 2b 1.2)

This case means

$\psi^{\prime}(\phi)\neq
0,\quad\mu(Z)=0,\quad\rho(r)=\frac{\rho_{0}}{r},\quad\tau(\phi)=\tau_{0}+\frac{\tau_{1}}{(\psi(\phi)-\rho_{0})^{2}},\quad\sigma(r)=\frac{\tau_{0}}{r^{2}}+\sigma_{0},$
(2.91)

thus, the following magnetic field

$B^{r}(r,\phi,Z)=-\frac{\tau_{1}\beta^{\prime}(\phi)}{r^{2}\beta(\phi)^{3}},\quad
B^{Z}(r,\phi,Z)=-\frac{\beta(\phi)+\beta^{\prime\prime}(\phi)}{2r^{3}},\quad
B^{\phi}(r,\phi,Z)=\frac{\tau_{1}}{r^{3}\beta(\phi)^{2}},$ (2.92)

and the scalar potential (1.132)

$\begin{split}W={}&\frac{W_{0}}{r^{2}\beta(\phi)^{2}}+\frac{\beta(\phi)\beta^{\prime\prime}(\phi)+\beta^{\prime}(\phi)^{2}+\frac{f_{1}}{4}-\frac{\tau_{1}^{2}}{\beta(\phi)^{4}}+2\beta(\phi)^{2}}{8r^{4}}+\\\
&+\frac{\hbar^{2}(2\beta(\phi)\beta^{\prime\prime}(\phi)+\beta(\phi)^{2}-\beta^{\prime}(\phi)^{2})}{8r^{2}\beta(\phi)^{2}},\end{split}$
(2.93)

where $\beta(\phi)$ is the auxiliary function $\psi(\phi)$ after a shift

$\beta(\phi)=\psi(\phi)-\rho_{0}.$ (2.94)

The function $\beta(\phi)$ must satisfy eq. (1.129) or equivalently its
reduced form (1.133). Using this equation, we can simplify the potential from
eq. (2.95) and obtain the potential from eq. (1.134), namely

$W=\frac{W_{0}}{r^{2}\beta(\phi)^{2}}-\frac{(4\tau_{1}^{2}+\beta_{2})}{32r^{4}\beta(\phi)^{4}}+\hbar^{2}\frac{f_{1}\beta(\phi)^{4}-12\beta_{1}\beta(\phi)^{2}-5\beta_{2}}{32r^{2}\beta(\phi)^{6}}.$
(2.95)

where $\beta_{1}$ and $\beta_{2}$ are integration constants in the reduced
equation (1.133).

Let us continue solving equations (2.3)–(2.5), taking eq. (2.24) into account.
Differentiating eq. (2.4) with respect to $r$, we obtain (multiplied by
$\beta(\phi)^{3}$, which must be non-zero)

$k_{6}\tau_{1}\beta^{\prime}(\phi)=0.$ (2.96)

Thus, we split our considerations into 4 subcases.

  1. I.

$\tau_{1}=0$, $k_{6}\neq 0$: Differentiating eq. (2.3) with respect to $r$, we
obtain

$k_{6}(\beta^{\prime}(\phi)+\beta^{\prime\prime\prime}(\phi))=0.$ (2.97)

Assumption $k_{6}\neq 0$ leads to

$\beta(\phi)=c_{1}\cos(\phi)+c_{2}\sin(\phi)+c_{3},$ (2.98)

but the constant $c_{3}$ must vanish to satisfy eq. (2.3), implying no
magnetic field.

  2. II.

$\tau_{1}=0$, $k_{6}=0$: From the last remaining compatibility equation (2.3)
we obtain $\beta(\phi)$ of the form

$\beta(\phi)=\frac{a_{1}+a_{2}\sin(2\phi)+a_{3}\cos(2\phi)}{k_{1}\sin(\phi)-k_{2}\cos(\phi)},$
(2.99)

leading to a system with non-vanishing magnetic field in the $z$-direction
only

$B^{Z}(r,\phi,Z)=-\frac{(a_{1}+a_{3})k_{1}^{2}+2a_{2}k_{2}k_{1}+k_{2}^{2}(a_{1}-a_{3})}{r^{2}(k_{1}\sin(\phi)-k_{2}\cos(\phi))^{3}},$
(2.100)

where $a_{i}$ are constants. Transforming the magnetic field into the
Cartesian coordinates using eq. (1.2), we get a somewhat simpler expression

$B^{z}(x,y,z)=-\frac{(a_{1}+a_{3})k_{1}^{2}+2a_{2}k_{2}k_{1}+k_{2}^{2}(a_{1}-a_{3})}{(k_{1}y-k_{2}x)^{3}}.$
(2.101)

Without loss of generality, we can simplify the previous expressions by using
eq. (2.82), choosing our coordinate system and redefining the constants
$a_{1},$ $a_{2}$ and $a_{3}$ so that $k_{1}=1$, $k_{2}=0$, i.e.

$\beta(\phi)=\frac{a_{1}+a_{2}\sin(2\phi)+a_{3}\cos(2\phi)}{\sin(\phi)}.$
(2.102)

The magnetic field simplifies to

$B^{Z}(r,\phi,Z)=-\frac{(a_{1}+a_{3})}{r^{2}(\sin(\phi))^{3}},\quad
B^{z}(x,y,z)=-\frac{(a_{1}+a_{3})}{y^{3}},$ (2.103)

with all other components vanishing both in the cylindrical and Cartesian
coordinates.

We have to check if $\beta(\phi)$ of the form eq. (2.102) satisfies eq.
(1.129). For that to be the case the constants must satisfy (omitting the
trivial $a_{1}=a_{2}=a_{3}=0$)

$a_{1}=-a_{3},\quad f_{1}=-16(a_{2}^{2}+a_{3}^{2}),$ (2.104)

which, however, implies vanishing of the magnetic field. This subcase,
therefore, does not admit any superintegrable system with additional first
order integrals and non-vanishing magnetic field.

  3. III.

$\tau_{1}\neq 0$, $k_{6}\neq 0$: Here equation (2.96) implies
$\beta(\phi)=\gamma=const$, which satisfies eq. (1.129) trivially. However,
eq. (2.4) in this case reads

$\tau_{1}(k_{1}\cos(\phi)+k_{2}\sin(\phi))=0,$ (2.105)

which (recall $\tau_{1}\neq 0$) implies no additional integral in this case.

  4. IV.

$\tau_{1}\neq 0$, $k_{6}=0$: From eq. (2.3) we obtain

$\beta(\phi)=c_{1}(k_{1}\sin(\phi)-k_{2}\cos(\phi))$ (2.106)

Substituting it into the differential equation (1.129), we obtain

$c_{1}\left(4(k_{1}^{2}+k_{2}^{2})c_{1}^{2}+f_{1}\right)(k_{1}\sin(\phi)+k_{2}\cos(\phi))=0,$
(2.107)

so we need

$c_{1}=\pm\sqrt{\frac{-f_{1}}{4(k_{1}^{2}+k_{2}^{2})}},\quad f_{1}<0,$ (2.108)

because vanishing $\beta(\phi)$ would imply division by 0.

We note that the quantum correction in eq. (2.95) vanishes for $\beta(\phi)$
of the form (2.106). The system defined by the following magnetic field $B$
and scalar potential $W$ is, therefore, the same in classical and quantum
mechanics.

$\displaystyle W$
$\displaystyle=-\left(\frac{2\tau_{1}^{2}\left(k_{2}^{2}+k_{1}^{2}\right)^{2}}{\left(k_{2}\cos(\phi)-k_{1}\sin(\phi)\right)^{4}f_{1}^{2}r^{4}}+\frac{4W_{0}\left(k_{2}^{2}+k_{1}^{2}\right)}{\left(k_{2}\cos(\phi)-k_{1}\sin(\phi)\right)^{2}f_{1}r^{2}}\right),$
(2.109) $\displaystyle B^{r}$
$\displaystyle=-{\frac{4\left(k_{2}\sin(\phi)+k_{1}\cos(\phi)\right)\left(k_{2}^{2}+k_{1}^{2}\right)\tau_{1}}{\left(k_{2}\cos(\phi)-k_{1}\sin(\phi)\right)^{3}r^{2}f_{1}}},$
(2.110) $\displaystyle{B^{\phi}}$
$\displaystyle=-{\frac{4\left(k_{2}^{2}+k_{1}^{2}\right)\tau_{1}}{r^{3}\left(k_{2}\cos(\phi)-k_{1}\sin\left(\phi\right)\right)^{2}f_{1}}},\quad{B^{Z}}=0.$
(2.111)

We note that we can write

$\frac{k_{2}\cos(\phi)-k_{1}\sin(\phi)}{\sqrt{k_{2}^{2}+k_{1}^{2}}}=\cos(\phi_{0})\cos(\phi)-\sin(\phi_{0})\sin(\phi)=\cos(\phi+\phi_{0}),$
(2.112)

which corresponds to the rotational symmetry of the system. Without loss of
generality, we can rotate the coordinate system so that
$\phi_{0}=\frac{3\pi}{2}$, which corresponds to $k_{2}=0,$ $k_{1}=1$. The
Cartesian form of the magnetic field and the scalar potential with the
adjusted coordinate axis reads

$\displaystyle
W=-4\left(\frac{\tau_{1}^{2}}{2f_{1}^{2}y^{4}}+\frac{W_{0}}{f_{1}y^{2}}\right),\quad
B^{x}={\frac{4\tau_{1}}{f_{1}y^{3}}},\quad{B^{y}}=0,\quad B^{z}=0$ (2.113)

The corresponding first order integrals of motion are

$\tilde{X}_{2}=p_{z}^{A}+{\frac{2\tau_{1}}{f_{1}y^{2}}},\quad
Y_{1}=p_{x}^{A}.$ (2.114)

The other cylindrical integral $X_{1}$ is of the second order and reads

$X_{1}=\left(xp_{y}^{A}-yp_{x}^{A}\right)^{2}-{\frac{4\tau_{1}\left({x}^{2}+y^{2}\right)}{f_{1}y^{2}}\left(p_{z}^{A}+\frac{2\tau_{1}}{f_{1}y^{2}}\right)}-{\frac{8W_{0}\left({x}^{2}+y^{2}\right)}{f_{1}y^{2}}},\\\
$ (2.115)

This case in the chosen coordinates corresponds to the Case Ib) or Ic) in [19]
with constants

$a_{1}=a_{2}=0,\quad a_{3}=\frac{2\tau_{1}}{f_{1}},\quad b_{1}=b_{2}=0,\quad
b_{3}=-\frac{4W_{0}}{f_{1}}.$ (2.116)

There is, therefore, another second order integral of motion, namely

$\displaystyle X_{3}$
$\displaystyle=\left(xp_{y}^{A}-yp_{x}^{A}\right)p_{y}^{A}-{\frac{4\tau_{1}x}{f_{1}y^{2}}\left(p_{z}^{A}+{\frac{2\tau_{1}}{f_{1}y^{2}}}\right)}-{\frac{8W_{0}x}{f_{1}y^{2}}},$
(2.117)

which in the natural choice of gauge

$A_{x}(x,y,z)=0,\quad A_{y}(x,y,z)=0,\quad
A_{z}(x,y,z)=-\frac{2\tau_{1}}{f_{1}y^{2}}$ (2.118)

simplifies to

$\displaystyle X_{3}$
$\displaystyle=\left(xp_{y}-yp_{x}\right)p_{y}-{\frac{4\tau_{1}x}{f_{1}y^{2}}p_{z}}-{\frac{8W_{0}x}{f_{1}y^{2}}}.$
(2.119)

There are, however, only 4 independent integrals of motion, because

$(\tilde{X}_{2}^{2}+Y_{1}^{2}-2H)X_{1}+X_{3}^{2}=-\frac{4}{f_{1}}(\tilde{X}_{2}\tau_{1}+2W_{0})(2H-\tilde{X}_{2}^{2}).$
(2.120)

The non-vanishing Poisson brackets read

$\displaystyle\\{X_{3},X_{1}\\}_{\text{P.B.}}$
$\displaystyle=2{p_{x}^{A}}\left(xp_{y}^{A}-yp_{x}^{A}\right)^{2}-{\frac{16W_{0}\left({x}^{2}+y^{2}\right){p_{x}^{A}}}{f_{1}y^{2}}}-$
(2.121)
$\displaystyle-{\frac{8\tau_{1}\left({x}^{2}+y^{2}\right){p_{x}^{A}}}{f_{1}y^{2}}\left({p_{z}^{A}}+{\frac{2\tau_{1}}{f_{1}y^{2}}}\right)}=2Y_{1}X_{1},$
$\displaystyle\\{X_{3},Y_{1}\\}_{\text{P.B.}}$
$\displaystyle=\left(p_{y}^{A}\right)^{2}-{\frac{4\tau_{1}}{f_{1}y^{2}}\left(p_{z}^{A}+{\frac{2\tau_{1}}{f_{1}y^{2}}}\right)}-{\frac{8W_{0}}{f_{1}y^{2}}}=2H-\tilde{X}_{2}^{2}-Y_{1}^{2},$
(2.122) $\displaystyle\\{X_{1},Y_{1}\\}_{\text{P.B.}}$
$\displaystyle=2p_{y}^{A}\left(xp_{y}^{A}-yp_{x}^{A}\right)-{\frac{8\tau_{1}x}{f_{1}y^{2}}\left(p_{z}^{A}+{\frac{2\tau_{1}}{f_{1}y^{2}}}\right)}-\frac{16W_{0}x}{f_{1}y^{2}}=2X_{3},$
(2.123)

so the Poisson algebra is already closed in the sense that the Poisson bracket
of any pair of integrals is a function (polynomial) in previously known
integrals, and thus generates no new independent integral. Therefore, the
system is not second order maximally superintegrable, as it was found in [19].

Let us now compute the trajectories and see if the bounded ones are closed,
which is typical for maximally superintegrable systems. The Hamilton’s
equations corresponding to the gauge-fixed Hamiltonian

$H=\frac{1}{2}\left((p_{x})^{2}+(p_{y})^{2}+\left(p_{z}-\frac{2\tau_{1}}{f_{1}y^{2}}\right)^{2}\right)-{\frac{4W_{0}f_{1}y^{2}+2\tau_{1}^{2}}{f_{1}^{2}y^{4}}}$
(2.124)

read

$\displaystyle\dot{x}=p_{x},\quad\dot{y}=p_{y},\quad\dot{z}=p_{z}-\frac{2\tau_{1}}{f_{1}y^{2}},$
(2.125)
$\displaystyle\dot{p}_{x}=0,\quad\dot{p}_{y}=-\frac{4\tau_{1}p_{z}+8W_{0}}{f_{1}y^{3}},\quad\dot{p}_{z}=0.$
(2.126)

The equations determining $x$ and $p_{x}$ correspond to the free motion of the
particle in the $x$-direction

$x(t)=p_{x0}t+x_{0}.$ (2.127)

Because the momentum $p_{z}$ is conserved, we find that $y(t)$ must satisfy
the following ODE

$\ddot{y}=-\frac{4\tau_{1}p_{z}+8W_{0}}{f_{1}y^{3}}.$ (2.128)

We note that the equation is not defined at $y(t)=0$. If
$\tau_{1}p_{z}+2W_{0}=0$, we have the free motion in this direction as well.
We continue with the assumption that it is not the case.

Equation (2.128) does not depend on the independent variable $t$, so it admits
$\dot{y}$ as an integrating factor and (assuming it is non-zero) we obtain

$\dot{y}(t)=\pm\frac{\sqrt{f_{1}(-C_{1}f_{1}y^{2}+4\tau_{1}p_{z}+8W_{0})}}{f_{1}y},$
(2.129)

where the integration constant $C_{1}$ is real and we assume that the square
root is well defined. Solving the separable ODE (2.129), we get

$y(t)=\pm\frac{\sqrt{C_{1}f_{1}\left(C_{1}^{2}f_{1}(t-t_{0})^{2}-4p_{z}\tau_{1}-8W_{0}\right)}}{C_{1}f_{1}}.$
(2.130)

We remind the reader that eq. (2.128) is not defined for $y(t)=0$, so we must
restrict the independent variable $t$ so that this does not occur. The
integration constant $t_{0}$ translates the origin of time and $C_{1}$ depends
on the initial values and constants of the system

$C_{1}=-\frac{f_{1}y_{0}^{2}\pm\sqrt{f_{1}^{2}y_{0}^{4}+16t_{0}^{2}f_{1}p_{z}\tau_{1}+32W_{0}t_{0}^{2}f_{1}}}{2t_{0}^{2}f_{1}}.$
(2.131)

The sign of the square root must be chosen so that the square root in eq.
(2.130) is well defined (positive argument) for $t$ in a suitably chosen
interval, which also determines the possible values of $t_{0}$.

The sign of $C_{1}$ determines the motion of the particle in the $y$
coordinate: If $C_{1}>0$, then the particle escapes to the infinity. If
$C_{1}<0$, then it falls on the origin $y=0$, where eq. (2.128) has a
singularity.

To solve the case $C_{1}=0$, we must return to the eq. (2.129), which leads to
the singular solution

$y(t)=\frac{2\sqrt{f_{1}(t-t_{0})\sqrt{f_{1}(p_{z}\tau_{1}+2W_{0})}}}{f_{1}}.$
(2.132)

This case is unbounded or goes to $y=0$ (where eq. (2.128) is ill defined) for
$f_{1}>0$ and $f_{1}<0$, respectively.

We can now solve the remaining equation in eq. (2.125), namely

$\dot{z}=p_{z}-\frac{2\tau_{1}}{f_{1}y^{2}}.$ (2.133)

The solution using $y(t)$ from eq. (2.130) is

$z(t)=-\frac{\tau_{1}C_{1}{\rm
argtanh}\left({\frac{C_{1}^{2}f_{1}\left(t-t_{0}\right)}{2\sqrt{C_{1}^{2}f_{1}\left(\tau_{1}p_{z}+2W_{0}\right)}}}\right)}{\sqrt{C_{1}^{2}f_{1}\left(\tau_{1}p_{z}+2W_{0}\right)}}+p_{z}t+C_{2}.$
(2.134)

The form of $z(t)$ above assumes that the argument of argtanh is real, i.e.

$f_{1}\left(\tau_{1}p_{z}+2W_{0}\right)>0,$ (2.135)

and not equal to $\pm 1$ (domain of argtanh). On the other hand, if

$f_{1}\left(\tau_{1}p_{z}+2W_{0}\right)<0,$ (2.136)

we can rewrite the argtanh using the identity

${\rm argtanh}(z)=\mathrm{i}\arctan(\mathrm{i}z),\quad\forall
z\in\mathbb{C}\setminus\\{\pm 1,\pm\mathrm{i}\\}$ (2.137)

and cancel the imaginary unit with the one from the denominator. Because

$\lim_{x\to\pm 1}{\rm argtanh}(x)=\pm\infty$ (2.138)

and it is dominant with respect to $p_{z}t$, $z(t)$ is unbounded in the first
case. In the second case, arctan is bounded and $z(t)$ is unbounded if and
only if $p_{z}\neq 0$.

Using the singular solution eq. (2.132) in eq. (2.133), we get

$z(t)=-\frac{\tau_{1}\ln(t-t_{0})}{2\sqrt{f_{1}(p_{z}\tau_{1}+2W_{0})}}+p_{z}t+C_{2},$
(2.139)

which is clearly unbounded as $t$ goes to $\infty$ or $t_{0}$.

The last possibility, $y(t)=p_{y}^{0}t+y_{0}$, leads to

$z(t)=p_{z}t+\frac{2\tau_{1}}{f_{1}p_{y}^{0}(p_{y}^{0}t+y_{0})}+C_{3},$
(2.140)

which is unbounded if $p_{z}\neq 0$ and singular at $t=0$ if $p_{y}\neq 0$.
(If $p_{y}^{0}=0$, the second term of $z(t)$ would be
$\frac{2\tau_{1}t}{f_{1}y_{0}^{2}}$, i.e. free motion with a modified
momentum.)