# Cylindrical first order superintegrability with complex magnetic fields

O. Kubů<EMAIL_ADDRESS><EMAIL_ADDRESS>L Šnobl  Czech
Technical University in Prague, Faculty of Nuclear Sciences and Physical
Engineering, Department of Physics, Břehová 7, 115 19 Prague 1, Czech Republic

###### Abstract

This article is a contribution to the study of superintegrable Hamiltonian
systems with magnetic fields on the three-dimensional Euclidean space
$\mathbb{E}_{3}$ in quantum mechanics. In contrast to the growing interest in
complex electromagnetic fields in the mathematical community following the
experimental confirmation of its physical relevance [X. Peng et al., Phys.
Rev. Lett. 114 (2015)], they were so far not addressed in the growing
literature on superintegrability. Here we venture into this field by searching
for additional first order integrals of motion to the integrable systems of
cylindrical type. We find that already known systems can be extended into this
realm by admitting complex coupling constants. In addition to them, we find
one new system whose integrals of motion also feature complex constants. All
these systems are multiseparable. Rigorous mathematical analysis of these
systems is challenging due to the non-Hermitian setting and lost gauge
invariance. We proceed formally and pose the resolution of these problems as
an open challenge.

††preprint: AIP/123-QED

##  I Introduction

This article is a contribution to the study of integrable and superintegrable
Hamiltonian systems with magnetic fields on the three-dimensional (3D)
Euclidean space $\mathbb{E}_{3}$ in quantum mechanics, with focus on complex
valued magnetic fields. More specifically, we assume a Hamiltonian of the form
(using units where $e=-1,\ m=1$)

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

with implicit summation over repeated indices $j=1,2,3$ (in the whole paper),
$\vec{p}=-i\hbar\vec{\nabla}$ is the momentum operator and
$\vec{A}=(A_{1}(\vec{x}),A_{2}(\vec{x}),A_{3}(\vec{x}))$ and $W(\vec{x})$ are
the vector and electrostatic potentials of the electromagnetic field.

Integrability then entails the existence of two algebraically independent
integrals of motion $X_{1},X_{2}$ (further specified below) mutually in
involution, i.e.

$[H,X_{1}]=[H,X_{2}]=[X_{1},X_{2}]=0.$ (2)

Superintegrability assumes additional one or two integrals, algebraically
independent of each other and the integrals needed for integrability. They are
usually considered to be polynomials in the momenta $p_{j}$, for computational
feasibility usually of a low order (typically 2).

Integrable (and especially superintegrable) systems are rare and distinguished
by the possibility to obtain the solution to their equations of motion in a
closed form. They are subsequently invaluable for gaining physical intuition
and serve as a starting point for modeling more complicated systems. Finding
and classifying these systems is therefore of utmost importance.

Despite its physical relevance, (super)integrability with magnetic fields was
mostly ignored due to computational difficulty. The first systematic result
remedying this omission was the article by Shapovalov on separable systems
Shapovalov  _et al._ (1972), followed by the articles in $\mathbb{E}_{2}$
Bérubé and Winternitz (2004); McSween and Winternitz (2000). Subsequent
articles in $\mathbb{E}_{3}$ assumed first order integrals Marchesiello  _et
al._ (2015) or separation of variables Zhalij (2015); Marchesiello and Šnobl
(2018); Bertrand and Šnobl (2019); Fournier  _et al._ (2020), but some recent
articles go beyond separation Marchesiello and Šnobl (2022); Kubů  _et al._
(2022). The non-relativistic quantum case with spin was also investigated, see
e.g. the recent article İ.Yurduşen  _et al._ (2021) and references therein.

In all the cases above the magnetic field was by assumption real. On the other
hand, there has recently been a growing interest in imaginary or complex
magnetic fields following the experimental confirmation of their physical
relevance by observing the Yang-Lee zeros Peng  _et al._ (2015). A recent
paper Fernández and García-Muñoz (2022) investigates an exact solution for
electron in graphene interacting with a complex magnetic field. A paper by
Jaramillo Jaramillo (2015) shows a formal analogy between stable marginally
outer trapped surfaces (MOTS) of black holes and non-relativistic charged
particle in complex electromagnetic field on closed surfaces. In the
mathematical community these systems are considered as pseudo-Mostafazadeh
(2002) or quasi-Hermitian Scholtz  _et al._ (1992); Dieudonné (1961) systems,
see also Krejčiřk and Siegl (2015); Krejčiřík  _et al._ (2015); Krejčiřík
(2019) for some more recent work and additional references.

As we have demonstrated in the previous paragraph, complex magnetic fields
merit further research. However, the non-Hermitian setting of these systems
poses several problems that we do not resolve in this paper and pose them as
an open challenge.

The most fundamental one is the rigorous definition of the magnetic
Hamiltonians. This is caused by the fact that the time-independent gauge
transformation $A^{\prime}(\vec{x})=A(\vec{x})+\nabla\chi(\vec{x}),$
$W^{\prime}(\vec{x})=W(\vec{x})$, manifesting itself in the quantum context as
the position dependent change of phase of the wave function

$U\psi(\vec{x})=\exp\left(\frac{i}{\hbar}\chi(\vec{x})\right)\psi(\vec{x}),$
(3)

is unbounded once $\chi(\vec{x})$ has a nonvanishing imaginary part, which may
be necessary to fully fix the gauge. In other words, we can fix the real part
of the vector potential by a unitary gauge transformation, but the choice of
the imaginary part involves an unbounded transformation, which may affect the
spectral properties.

Even if this problem is addressed, the spectrum of the corresponding
Hamiltonian can be complex, which complicates its physical interpretation as
only real spectra are usually measurable. The mathematically correct approach
to the problem entails pseudo-Mostafazadeh (2002) or quasi-Hermiticity (or
self-adjointness) Scholtz  _et al._ (1992); Dieudonné (1961), i.e. the
Hamiltonian $H$ satisfies $H^{\dagger}=\Theta H\Theta^{-1}$ where the metric
$\Theta$ and its inverse $\Theta^{-1}$ are bounded operators, and is self-
adjoint with respect to the modified scalar product
$\langle\cdot,\Theta\cdot\rangle$. However, in the more general pseudo-
Hermitian case, where the metric $\Theta$ is indefinite, the modified “scalar
product” is indefinite as well and we deal with the so called Krein space, see
e.g. Bognár (2012). Only when the metric is positive, i.e. in the more strict
quasi-Hermitian case, we can talk about a nonstandard representation of usual
quantum mechanics. (For bounded operators $H$ these definitions coincide if
$0\notin\overline{\\{\langle\psi,\Theta\psi\rangle|\psi\in\mathcal{H},\lVert\psi\rVert=1\\}}$,
see Williams (1969) for a proof.) For more detail concerning these notions see
Krejčiřk and Siegl (2015); Krejčiřík  _et al._ (2015) and references therein.

Following on our research with real magnetic fields Kubů , we search for
superintegrable system of the cylindrical type with additional first order
integrals whose electromagnetic fields are complex functions. In Section II we
review cylindrical-type integrable systems, i.e. the corresponding integrals
of motion and the magnetic field, in cylindrical coordinates. In Section III
we specify the problem of first order superintegrability with complex fields
and briefly comment on the calculations leading to our classification.
Interested reader can find the details of the computation in the Appendix A.
The core of our article is the Results section IV. There we analyze the new
system with complex integrals of motion in Subsection IV.1 and the old systems
extended by choosing complex coupling constants in Subsection IV.2. Our
concluding remarks are in Section V.

##  II Cylindrical–type system

Before we specify the corresponding integrals $X_{1},X_{2}$, we have to
introduce the formalism used for magnetic field in curvilinear coordinates in
classical mechanics, cf. Marchesiello  _et al._ (2018); Fournier  _et al._
(2020).

Defining the cylindrical coordinates

$x^{1}=r\cos(\phi),\quad x^{2}=r\sin(\phi),\quad x^{3}=Z,$ (4)

we represent the vector potential $A$ as a 1-form

$A=A_{1}\mathrm{d}x^{1}+A_{2}\mathrm{d}x^{2}+A_{3}\mathrm{d}x^{3}=A_{r}\mathrm{d}r+A_{\phi}\mathrm{d}\phi+A_{Z}\mathrm{d}Z.$
(5)

Hence we obtain the following transformations

$\displaystyle A_{1}=\cos(\phi)A_{r}-\frac{\sin(\phi)}{r}A_{\phi},\quad
A_{2}=\sin(\phi)A_{r}+\frac{\cos(\phi)}{r}A_{\phi},\quad A_{3}=A_{Z}.$ (6)

As components of the canonical 1-form $\lambda=p_{j}{\operatorname{d}}x^{j}$,
the momenta $p_{j}$ transform in the same way and we can define the magnetic.
momenta by $p^{A}_{j}=p_{j}+A_{j}$ in both Cartesian and cylindrical
coordinates. (We used to call $p^{A}_{j}$ the gauge covariant momenta, the
terminology coined probably in Marchesiello  _et al._ (2015) but implicitly
used as velocities $p_{j}^{A}=\dot{x}^{j}$ since the first papers on
superintegrability with magnetic field Dorizzi  _et al._ (1985). With complex
magnetic field it is no longer appropriate terminology because, as we
mentioned in the introduction, the gauge transformation (3) may be unbounded.)

Components of the magnetic field 2-form $B=\mathrm{d}A$ are

$\displaystyle\begin{aligned}
B&=B^{1}(\vec{x})\mathrm{d}x^{2}\wedge\mathrm{d}z+B^{2}(\vec{x})\mathrm{d}z\wedge\mathrm{d}x^{1}+B^{3}(\vec{x})\mathrm{d}x^{1}\wedge\mathrm{d}x^{2}\\\
&=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}$
(7)

which leads to the following transformation

$\displaystyle B^{1}(\vec{x})$
$\displaystyle=\frac{\cos(\phi)}{r}B^{r}(r,\phi,Z)-\sin(\phi)B^{\phi}(r,\phi,Z),$
$\displaystyle B^{2}(\vec{x})$
$\displaystyle=\frac{\sin(\phi)}{r}B^{r}(r,\phi,Z)+\cos(\phi)B^{\phi}(r,\phi,Z),$
(8) $\displaystyle B^{3}(\vec{x})$ $\displaystyle=\frac{1}{r}B^{Z}(r,\phi,Z).$

We can use the same formalism and notation in quantum mechanics as well, we
just have to quantize the equations and the (properly symmetrized) integrals
in Cartesian coordinates and subsequently transform our equations into
cylindrical ones. For example, the transformed momenta read

$p_{1}=-i\hbar\left(\cos(\phi)\partial_{r}-\frac{\sin(\phi)}{r}\partial_{\phi}\right),\
p_{2}=-i\hbar\left(\sin(\phi)\partial_{r}+\frac{\cos(\phi)}{r}\partial_{\phi}\right),\
p_{3}=-i\hbar\partial_{Z},$ (9)

i.e. transformation is the same as (6) upon defining
$p_{r,\phi,Z}=-i\hbar\partial_{r,\phi,Z}$.

We can now introduce integrals of motion of the cylindrical type, i.e.
integrals that imply separation of Schrödinger (or, classically, Hamilton-
Jacobi) equation in the cylindrical coordinates in the limit of vanishing
magnetic field $\vec{B}$. Expressed in the cylindrical coordinates they read

$\begin{split}&X_{1}=(p_{\phi}^{A})^{2}+\frac{1}{2}\sum_{\alpha=r,\phi,Z}\left(s_{1}^{\alpha}(r,\phi,Z)p_{\alpha}^{A}+p_{\alpha}^{A}s_{1}^{\alpha}(r,\phi,Z)\right)+m_{1}(r,\phi,Z),\\\
&X_{2}=(p_{Z}^{A})^{2}+\frac{1}{2}\sum_{\alpha=r,\phi,Z}\left(s_{2}^{\alpha}(r,\phi,Z)p_{\alpha}^{A}+p_{\alpha}^{A}s_{2}^{\alpha}(r,\phi,Z)\right)+m_{2}(r,\phi,Z).\end{split}$
(10)

The functions $s^{r,\phi,Z}_{1,2},m_{1,2}$ together with the electromagnetic
fields $B,W$ are to be determined from the integrability conditions (2).

Since the integrals of motion as well as the Hamiltonian and their commutators
are differential operators, the integrability conditions (2) can be separated
into the coefficients of derivatives which must vanish independently, yielding
the so-called determining equations. The second order ones can be solved in
terms of 5 auxiliary functions
$\rho(r),\sigma(r),\psi(\phi),\tau(\phi),\mu(Z)$ of one variable each that we
use in the subsequent calculations, namely

$\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}$ (11)
$\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}$
(12)

cf. Fournier  _et al._ (2020); Kubů ; Kubů and Šnobl (2022). We further use
primes for derivatives of these functions with respect to their respective
independent variable.

##  III First order superintegrability in the cylindrical case with complex
fields

We are looking for first order superintegrability in the cylindrical case.
This means that we impose the integrability conditions (2) with the
cylindrical integrals (LABEL:cyl_integrals). First order superintegrability
then amounts to the existence of an additional integral of the form (in
Cartesian coordinates)

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

Here we have the magnetic momenta $p^{A}_{j}=p_{j}+A_{j}(x^{1},x^{2},x^{3})$
and angular momenta $L^{A}_{j}=\epsilon_{jkl}x^{k}p^{A}_{l}.$ (We sum over
both indices $k,l$ and $\epsilon_{jkl}$ is the totally antisymmetric Levi-
Civita symbol.) The constants $k_{j}$ and the function $m$ may be complex but
the coordinates $x^{1},x^{2},x^{3}$ and the corresponding momenta are assumed
to be real.

To summarize our problem, we assume that

$[X_{1},H]=[X_{2},H]=[X_{1},X_{2}]=[Y,H]=0$ (14)

with the Hamiltonian (1), and integrals (LABEL:cyl_integrals) and (13). As
mentioned above, we can separate these equations into coefficients of
derivatives such as $\partial_{x_{1}x_{3}}$ (or $\partial_{rZ}$ in cylindrical
coordinates) that must vanish separately. We thus obtain the determining
equations which we solve for the fields $\vec{B},W$ together with functions
$s^{r,\phi,Z}_{1,2},m_{1,2},m$ and constants $k_{j}$ in the integrals. (We
note that if the constants $k_{6}$ and $k_{3}$ are not intertwined with any
other constant $k_{j}$ in (13), i.e. all other constants $k_{j}$ in (13) can
be set to 0, then the corresponding integrals,
$\tilde{X}_{1}=L_{3}^{A}+\tilde{m}(x^{1},x^{2},x^{3})=p_{\phi}^{A}+\tilde{m}(r,\phi,Z)$
and
$\tilde{X}_{2}=p_{3}^{A}+\tilde{m}(x^{1},x^{2},x^{3})=p_{Z}^{A}+\tilde{m}(r,\phi,Z)$,
respectively, can be seen as a “square root” of the cylindrical integrals
(LABEL:cyl_integrals). In this case the cylindrical integrals
(LABEL:cyl_integrals) are dependent on these first order integrals and we do
not have the required superintegrability.)

We proceed in cylindrical coordinates. We substitute the partially reduced
form of magnetic field $\vec{B}$ from (12) into the equations for $Y$, which
we use to constrain $\vec{B}$ further and determine $W$. In the process we
split into several cases depending on whether the constants $k_{j}$ vanish or
not, and subsequently use conditions on $X_{1},X_{2}$ to constrain the fields
further and determine their final form. We include the details of our
calculation in the Appendix A.

Taking into account that the cases with real electromagnetic fields and
integrals were considered in O. Kubů’s Master thesis Kubů , we are mainly
interested in solutions not found there, especially those with complex
constants $k_{j}$ in the integral $Y$.

We find one such system, which we analyze in Subsection IV.1. We also confirm
that the systems found earlier can be extended into complex systems by simply
allowing complex coupling constant. We analyze those in Subsection IV.2.

##  IV Results

###  IV.1 The new system

The corresponding electromagnetic field in Cartesian coordinates is given by

$\vec{B}(\vec{x})=\left(-\frac{ib}{(x^{1}-ix^{2})^{3}},-\frac{b}{(x^{1}-ix^{2})^{3}},0\right),\quad
W(\vec{x})=\frac{w_{1}}{2(x^{1}-ix^{2})^{2}}-\frac{b^{2}}{8(x^{1}-ix^{2})^{4}},$
(15)

where both coupling constants $b$ and $w_{1}$ may be complex.

With complex magnetic field, we no longer have the gauge freedom. Let us write
all the results in terms of magnetic momenta $p_{j}^{A}$.

This system is maximally superintegrable. We have 3 first order integrals,

$Y_{1}=ip^{A}_{1}+p^{A}_{2},\quad
Y_{2}=L^{A}_{2}+iL^{A}_{1}-\frac{b}{(x^{1}-ix^{2})},\quad\tilde{X}_{2}=p^{A}_{3}+\frac{b}{2(x^{1}-ix^{2})^{2}},$
(16)

i.e. one of the cylindrical integrals reduces to a first order one,
$X_{2}=\tilde{X}_{2}^{2}$, followed by the Hamiltonian and the second
cylindrical integral

$\displaystyle H={}$
$\displaystyle\frac{1}{2}((p^{A}_{1})^{2}+(p^{A}_{2})^{2}+(p^{A}_{3})^{2})+\frac{w_{1}}{2(x^{1}-ix^{2})^{2}}-\frac{b^{2}}{8(x^{1}-ix^{2})^{4}},$
(17) $\displaystyle X_{1}={}$
$\displaystyle(L^{A}_{3})^{2}-\frac{b(x^{1}+ix^{2})}{x^{1}-ix^{2}}p^{A}_{3}-b^{2}\frac{(x^{1}+ix^{2})}{2(x^{1}-ix^{2})^{3}}+w_{1}\frac{(x^{1}+ix^{2})}{x^{1}-ix^{2}}.$
(18)

The algebraic independence of all 5 integrals is obvious once we notice that
$L_{3}^{2}$ contains $(x^{3})^{2}p_{1}^{2}$ and $x^{3}$ is absent from the
other integrals.

The determining equations for the classical and quantum (super)integrability
are very similar if we choose a suitable form of symmetrization. For first
order integrals there is no change at all, for second order integrals only the
zeroth order equation differs, namely there is an a priori nonvanishing
additional term proportional to $\hbar^{2}$ Marchesiello  _et al._ (2015). In
Kubů  we showed that the correction for the integral $X_{1}=L_{3}^{2}+\ldots$
has the following form
$-\frac{\hbar^{2}}{2}(x^{1}\partial_{x^{2}}-x^{2}\partial_{x^{1}})B^{3}$ and
vanishes for the Cartesian integrals $P_{i}^{2}+\ldots$, including
$X_{1}=P_{3}^{2}+\ldots$ In our case $B^{3}=0$, therefore the integrals have
the same form in classical mechanics as well. However, interpretation of our
system in that context is unclear due to complex coefficients in (15).

Comparing the form of the Hamiltonian (17) with Case IIc in Shapovalov  _et
al._ (1972), we see that the system separates in cylindrical coordinates, with
the corresponding commuting integrals $\tilde{X}_{1}$ and $X_{2}$.

Let us introduce “complex” coordinates for the $x^{1}x^{2}$ plane,
$z=x^{1}+ix^{2}$, $\bar{z}=x^{1}-ix^{2}$. The corresponding magnetic momenta
read

$p^{A}_{z}=\frac{1}{2}(p^{A}_{1}-ip^{A}_{2}),\quad
p^{A}_{\bar{z}}=\frac{1}{2}(p^{A}_{1}+ip^{A}_{2}).$ (19)

We note that the vanishing $B^{z}$, see (15), implies that $p_{1}^{A}$ and
$p_{2}^{A}$ commute. The same clearly holds also for $p_{z}^{A}$ and
$p_{\bar{z}}^{A}$.

The transformed electromagnetic field reads

$B=\frac{b}{\bar{z}^{3}}{\operatorname{d}}\bar{z}\wedge{\operatorname{d}}x^{3},\quad
W(z,\bar{z},x^{3})=\frac{w_{1}}{2\bar{z}^{2}}-\frac{b^{2}}{8\bar{z}^{4}},$
(20)

and the integrals are as follows

$\displaystyle Y_{1}={}$ $\displaystyle 2ip^{A}_{z},\quad
Y_{2}=2x^{3}p^{A}_{z}-\bar{z}p^{A}_{3}-\frac{b}{\bar{z}},\quad\tilde{X}_{2}=p^{A}_{3}+\frac{b}{2\bar{z}^{2}},$
(21) $\displaystyle H={}$
$\displaystyle\frac{1}{2}(4p^{A}_{z}p^{A}_{\bar{z}}+(p^{A}_{3})^{2})+\frac{w_{1}}{2\bar{z}^{2}}-\frac{b^{2}}{8\bar{z}^{4}},$
(22) $\displaystyle X_{1}={}$
$\displaystyle-(zp^{A}_{z}-\bar{z}p^{A}_{\bar{z}})^{2}-\frac{bz}{\bar{z}}p^{A}_{3}-b^{2}\frac{z}{2\bar{z}^{3}}+\frac{w_{1}z}{\bar{z}}.$
(23)

The minus sign in $X_{1}$ is obtained from the complex unit $i$, as the
angular momentum $L^{A}_{3}=x^{1}p^{A}_{2}-x^{2}p^{A}_{1}$ transforms into
$i(zp^{A}_{z}-\bar{z}p^{A}_{\bar{z}})$.

The forms of the Hamiltonian and the integral $X_{1}$ imply that this system
can be seen as an extension of a 2D system (at least in the classical
context), a method to obtain 3D system with magnetic field from 2D systems
discussed in classical setting without complex coordinates in Marchesiello and
Šnobl (2020); Kubů  _et al._ (2021). There we can set the integral of motion
$\tilde{X}_{2}$ to a constant $\kappa$ and obtain the corresponding 2D system
by merging the first and zeroth order terms. The corresponding 2D system lacks
the linear term and thus has no magnetic field. We see that the use of complex
coordinates $z=x^{1}+ix^{2}$, $\bar{z}=x^{1}-ix^{2}$ in the 2D system does not
interfere with the method.

However, this method is formal in the quantum context because the purely
continuous spectrum of $p_{3}$ precludes a reasonable choice of the constant
$\kappa$.

Let us analyze the algebra of integrals. The commutators are as follows, with
only the leading order terms for the dependent integrals $R_{1}$ and $R_{2}$
written for brevity:

$\displaystyle[Y_{1},Y_{2}]={}$ $\displaystyle
0,\quad[Y_{1},\tilde{X}_{2}]=0,\quad[Y_{2},\tilde{X}_{2}]=-2i\hbar p_{z}=\hbar
Y_{1},\quad[X_{1},\tilde{X}_{2}]=0$ (24) $\displaystyle[Y_{1},X_{1}]={}$
$\displaystyle
2i\hbar(i\\{p_{z}^{A},zp_{z}^{A}-\bar{z}p_{\bar{z}}^{A}\\})+\ldots=:i\hbar
R_{1},$ (25) $\displaystyle[Y_{2},X_{1}]={}$ $\displaystyle
i\hbar\\{2x^{3}p_{z}^{A}-\bar{z}p_{3}^{A},zp_{z}^{A}-\bar{z}p_{\bar{z}}^{A}\\}+\ldots=:i\hbar
R_{2},$ (26)

where $\\{\cdot,\cdot\\}$ denotes the anticommutator, $\\{a,b\\}=ab+ba.$ The
higher order commutators are as follows

$\displaystyle[Y_{1},R_{1}]={}$ $\displaystyle-2i\hbar
Y_{1}^{2},\quad[Y_{2},R_{1}]=-i\hbar\\{Y_{1},Y_{2}\\},\quad[\tilde{X}_{2},R_{1}]=0,\quad$
(27) $\displaystyle[Y_{1},R_{2}]={}$
$\displaystyle-i\hbar\\{Y_{1},Y_{2}\\},\quad[Y_{2},R_{2}]=-i\hbar(2Y_{2}^{2}-6b\tilde{X}_{2}+2w_{1}),\quad[\tilde{X}_{2},R_{2}]=-\hbar
R_{1},$ (28) $\displaystyle[X_{1},R_{1}]={}$ $\displaystyle
i\hbar(2\\{X_{1},Y_{1}\\}-\hbar^{2}Y_{1}),\quad[X_{1},R_{2}]=i\hbar(2\\{Y_{2},X_{1}\\}-\hbar^{2}Y_{2}),$
(29) $\displaystyle[R_{1},R_{2}]={}$ $\displaystyle
i\hbar(4ibH-6ib\tilde{X}_{2}^{2}+4iw_{1}\tilde{X}_{2}),$ (30)

so the algebra closes. We notice that the terms with $\hbar^{3}$ in (29) do
not arise in the classical mechanics, i.e. they modify the Poisson brackets.

A similar 2D system arises as $E_{8}$ with $a_{2}=0$ in the classification of
2D superintegrable systems (see e.g. the review article Miller  _et al._
(2013))

$H=4p_{z}p_{\bar{z}}+a_{1}z\bar{z}-\frac{a_{3}}{\bar{z}^{2}},$ (31)

(notice the choice $m=\frac{1}{2}$ instead of $m=1$ in Miller  _et al._
(2013)), was recently analyzed from the point of view of dynamical symmetry
algebras Marquette and Quesne (2022). We did not find this more general 2D
system, which can be extended into 3D as well, due to our first order
additional integral ansatz. In contrast to the system (31), which exhibits an
equidistant, oscillator-like spectrum (determined using ladder operators
Marquette and Quesne (2022)), for our system (20) the quadratic integral
$p_{z}^{2}+\frac{a_{1}}{4}\bar{z}^{2}$ reduces to a first order one
$p_{z}=-i\hbar\partial_{z}$ and there are no ladder operators. The free motion
in the $z$ direction of our system suggests a purely essential spectrum.

####  IV.1.1 Physical interpretation of the system

There are several problems with the rigorous definition of the Hamiltonian
(17). First the choice of the vector potential $A$: we no longer have the
freedom to choose the imaginary part as the gauge transformation is unbounded.
Second, if we cannot choose $A$, we must retain the $(x^{1}-ix^{2})^{-4}$ term
in $W$ (15). Due to its strong singularity, it cannot serve as a perturbation
with respect to magnetic Laplacian
$(p^{A}_{1})^{2}+(p^{A}_{2})^{2}+(p^{A}_{3})^{2}$ in the standard definition.

We leave these problems as an open challenge. In this article, we proceed
formally by choosing the simplest gauge for (15) satisfying
$\operatorname{div}\vec{A}=0$, i.e.

$\vec{A}(\vec{x})=\left(0,0,-\frac{b}{2(x^{1}-ix^{2})^{2}}\right),$ (32)

which also eliminates the most singular term in $W$ (15).

The Hamiltonian (17) now reads

$H=\frac{1}{2}(p_{1}^{2}+p_{2}^{2}+p_{3}^{2})-\frac{b}{2(x^{1}-ix^{2})^{2}}p_{3}+\frac{w_{1}}{2(x^{1}-ix^{2})^{2}}.$
(33)

This Hamiltonian is pseudo-Hermitian Mostafazadeh (2002), i.e. $\Theta
H\Theta^{-1}=H^{\dagger}$, with $\Theta=\mathcal{P}_{2}$ a Hermitian (i.e.
bounded and self-adjoint) invertible operator that changes the sign of
$x^{2}\to-x^{2}$ and $p_{2}\to-p_{2}$. However, the partial parity
$\mathcal{P}_{2}$ is not a positive operator, therefore our Hamiltonian is not
quasi-Hermitian in the sense of Scholtz  _et al._ (1992); Dieudonné (1961).

The Hamiltonian is neither invariant under the parity $\mathcal{P}$ nor the
time reversal $\mathcal{T}$. It is $\mathcal{PT}$-self-adjoin without being
self-adjoint with respect to $\mathcal{P}$ or $\mathcal{T}.$

We therefore conclude that the Hamiltonian (17) is pseudo-Hermitian but not
quasi-Hermitian and the physical interpretation of this system is not clear.

The Schrödinger equation $H\Psi=E\Psi$ for the system (33) can be solved in
the complex coordinates $z=x^{1}+ix^{2}$, $\bar{z}=x^{1}-ix^{2}$ by separation
of variables with the following smooth non-normalizable solution

$\Psi\left(z,\bar{z},x^{3}\right)\sim
N\exp\left(Cz+\frac{\lambda_{3}^{2}-2E}{4C\hbar^{2}}\bar{z}+\frac{w_{1}-\lambda_{3}b}{4C\hbar^{2}}\frac{1}{\bar{z}}\right)\exp\left(\frac{i}{\hbar}\lambda_{3}x^{3}\right),$
(34)

where $C$ is the separation constant and $\lambda_{3}\in\mathbb{R}$
corresponds to the momentum $p_{3}$ by the Fourier transform.

We are not able to prove that the spectrum is purely essential. Both
exponentials in (34) need to be purely imaginary for standard construction of
the Weyl sequence, but we do not know whether $\lambda$ is real or not. In the
oscillating case (31) the eigenfunctions $\psi$ are bounded by the Gaussian
term
$\exp(-\frac{\lambda}{2}z\bar{z})=\exp[-\frac{\lambda}{2}((x^{1})^{2}+(x^{2})^{2})]$
so the previously discussed condition is not necessary.

###  IV.2 Systems found earlier admit complex coupling constants

All other systems have already been found in Marchesiello  _et al._ (2015) or
Marchesiello and Šnobl (2020) (and considered also in Kubů ), i.e. the form of
their integrals is a standard one. We have the following 3 systems, where the
coupling constants $b_{j},w_{j}$ may be complex.

  1. 1.

The constant magnetic field and vanishing scalar potential

$\vec{B}(\vec{x})=(0,0,b),\quad W(\vec{x})=0,$ (35)

i.e. our Hamiltonian is

$H=\frac{1}{2}((p^{A}_{1})^{2}+(p^{A}_{2})^{2}+(p^{A}_{3})^{2}).$ (36)

Proceeding formally, we present 2 choices of gauge, the first one suitable for
separation in Cartesian coordinates, the second one for the cylindrical ones,

$\vec{A}=(0,bx^{1},0),\quad\vec{A}=\left(-\tfrac{b}{2}x^{2},\tfrac{b}{2}x^{1},0\right).$
(37)

For a general complex $b=\alpha+i\beta$ the Hamiltonian is symmetric under the
parity $\mathcal{P}$ but not $\mathcal{P}$-self-adjoint,
$\mathcal{P}H\mathcal{P}=H\neq H^{\dagger}$, and neither $\mathcal{T}$\- nor
$\mathcal{PT}$-symmetric.

Excluding the self-adjoint real case, the purely imaginary $b=i\beta$ yields
an interesting example of $\mathcal{PT}$-symmetric (also $\mathcal{P}$\- and
$\mathcal{T}$-symmetric) operator which is neither $\mathcal{T}$\- nor
$\mathcal{P}$-self-adjoint, i.e.

$\mathcal{PT}H\mathcal{PT}=\mathcal{P}H\mathcal{P}=\mathcal{T}H\mathcal{T}=H\neq
H^{\dagger}.$ (38)

The Hamiltonian is therefore not pseudo-Hermitian with respect to a metric
$\mathcal{P}$, $\mathcal{T}$ nor $\mathcal{PT}$, but we cannot exclude a more
complicated metric. This is another example showing that
$\mathcal{PT}$-symmetry and pseudo- or quasi-Hermiticity are not necessarily
equivalent. (In fact, the inequivalence was demonstrated in Siegl , see also
Krejčiřk and Siegl (2015); Krejčiřík  _et al._ (2015).)

This system is maximally superintegrable with 4 first order integrals

$Y_{1}=p_{1}^{A}-bx^{2},\quad
Y_{2}=p_{2}^{A}+bx^{1},\quad\tilde{X}_{1}=L_{3}^{A}+\frac{b((x^{1})^{2}+(x^{2})^{2})}{2},\quad\tilde{X}_{2}=p_{3}^{A},$
(39)

and a nonpolynomial one Marchesiello  _et al._ (2015)

$X_{5}=(p_{2}^{A}+bx^{1})\sin\left(\frac{bx^{3}}{p_{3}^{A}}\right)-(p_{1}^{A}-bx^{2})\cos\left(\frac{bx^{3}}{p_{3}^{A}}\right),$
(40)

where we assume that $p_{3}^{A}\neq 0$. Interpretation of the last integral in
the quantum mechanical context remains unclear due to the
$\frac{x^{3}}{p_{3}^{A}}$ term in the trigonometric functions and purely
continuous spectrum of $p_{3}^{A}=p_{3}$.

As in the real case analyzed in Kubů , the stationary Schrödinger equation
$\hat{H}\psi=E\psi$ separates in the Cartesian coordinates in the first
(Landau) gauge (37)

$\displaystyle\Psi(\vec{x})=f(x)\exp\left(\frac{i}{\hbar}\lambda_{2}y\right)\exp\left(\frac{i}{\hbar}\lambda_{3}z\right),$
(41)
$\displaystyle\hbar^{2}f^{\prime\prime}(x)=\left((bx-\lambda_{2})^{2}+\lambda_{3}^{2}-2E\right)f(x),$
(42) $\displaystyle
Y_{2}\psi(\vec{x})=\lambda_{2}\Psi(\vec{x}),\quad\tilde{X}_{2}\Psi(\vec{x})=\lambda_{3}\Psi(\vec{x}).$
(43)

(The separation can be done rigorously using the (unitary) Fourier transform.)

The reduced Schrödinger equation for $f(x)$ can be solved in terms of
confluent hypergeometric functions Kubů , which under some condition reduce to
Hermite polynomials (see e.g. (DLMF, , Chapter 13)), i.e. the form known from
(shifted) 1D harmonic oscillator. This happens if the energy $E$ satisfies

$E=\frac{\lambda_{3}^{2}}{2}+\hbar b\left(n+\frac{1}{2}\right),\quad
n=0,1,\ldots,$ (44)

and the corresponding non-normalizable eigenfunctions are

$\Psi_{n,\lambda_{2},\lambda_{3}}(\vec{x})=K_{n}\exp\left(\frac{i}{\hbar}(\lambda_{2}y+\lambda_{3}z)\right)H_{n}\left(\sqrt{\frac{b}{\hbar}}\left(x-\frac{\lambda_{2}}{b}\right)\right)\exp\left(-\frac{b}{2\hbar}\left(x-\frac{\lambda_{2}}{b}\right)^{2}\right).$
(45)

The result so far suggests that all the energies are in the essential
spectrum,

$[\hbar b/2,+\infty)\subset\sigma_{ess}(H).$ (46)

However, the standard argument, i.e. construction of Weyl sequences, does not
work as the last exponential in (45) is not bounded due to the nonvanishing
imaginary part ($b=\alpha+i\beta$)

$-\frac{\beta\left(\alpha^{2}x^{2}+\beta^{2}x^{2}-\lambda_{2}^{2}\right)}{2\hbar\left(\alpha^{2}+\beta^{2}\right)}.$
(47)

We are therefore not able to prove the result above.

  2. 2.

Constant magnetic field with a nontrivial scalar potential governing the
motion in direction of the $x^{3}$ axis,

$\vec{B}(\vec{x})=(0,0,b),\quad W(\vec{x})=W_{3}(x^{3}),$ (48)

i.e. our Hamiltonian is

$H=\frac{1}{2}((p^{A}_{1})^{2}+(p^{A}_{2})^{2}+(p^{A}_{3})^{2})+W_{3}(x^{3}).$
(49)

We again continue formally, i.e. we choose the gauge (37). The properties of
the Hamiltonian (49) with respect the parity $\mathcal{P}$ and/or time-
reversal $\mathcal{T}$ carry over from the previous case if and only if
$W_{3}$ has the same property.

In this case the second cylindrical integral remains quadratic,

$X_{2}=\left(p_{3}^{A}\right)^{2}+W_{3}(x^{3}),$ (50)

the other integrals in (39) remain unchanged. The nonpolynomial integral from
(40) cannot be used.

Some special choices of the potential $W_{3}(x^{3})$ admit additional
integrals, e.g.

$W_{3}(x^{3})=\frac{c}{(x^{3})^{2}}+\frac{b^{2}(x^{3})^{2}}{8},\quad
W_{3}(x^{3})=\frac{b^{2}}{2}(x^{3})^{2},$ (51)

are known to admit the fifth quadratic integral Marchesiello and Šnobl (2017).
These algebraic results carry over to the case with complex coupling
constants, but we again have the problem with rigorous definition of the
system.

For (49) we also obtain purely essential spectrum consisting of infinitely
degenerate energies (44) where we must replace $\lambda_{3}^{2}$ by the
eigenvalue $\xi$ corresponding to the spectral problem
$(\left(p_{3}^{A}\right)^{2}+W_{3}(x^{3}))\psi=\xi\psi$.

  3. 3.

The system with non-constant magnetic field,

$\displaystyle\vec{B}=\left({\frac{4b}{(x^{2})^{3}}},0,0\right),\quad
W=-4\left(\frac{b^{2}}{2(x^{2})^{4}}+\frac{w_{0}}{(x^{2})^{2}}\right),$ (52)

i.e. our Hamiltonian is

$H=\frac{1}{2}((p^{A}_{1})^{2}+(p^{A}_{2})^{2}+(p^{A}_{3})^{2})-4\left(\frac{b^{2}}{2(x^{2})^{4}}+\frac{w_{0}}{(x^{2})^{2}}\right).$
(53)

Formally choosing the gauge satisfying $\operatorname{div}\vec{A}=0$,

$\vec{A}=(0,0,-2b(x^{2})^{-2}),$ (54)

the Hamiltonian is even in $x^{2}$ and therefore invariant under the parity
$\mathcal{P}$. The $\mathcal{PT}$-symmetry then depends on the invariance of
constants under time-reversal $\mathcal{T}$, which means that they have to be
real. On the other hand, the Hamiltonian is $\mathcal{P}$-pseudo-Hermitian if
we choose them purely imaginary. As the parity $\mathcal{P}$ is not positive
definite, the Hamiltonian is not quasi-Hermitian with the metric
$\mathcal{P},\mathcal{T}$ or $\mathcal{PT}$, though we cannot exclude a more
complicated metric.

The corresponding first order integrals of motion are

$\tilde{X}_{2}=p_{3}^{A}+{\frac{2b}{(x^{2})^{2}}},\quad Y_{1}=p_{1}^{A}.$ (55)

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

$X_{1}=\left(L_{3}^{A}\right)^{2}-{\frac{4b\left((x^{1})^{2}+(x^{2})^{2}\right)}{(x^{2})^{2}}\left(p_{3}^{A}+\frac{2b}{(x^{2})^{2}}\right)}-{\frac{8w_{0}\left((x^{1})^{2}+(x^{2})^{2}\right)}{(x^{2})^{2}}},\\\
$ (56)

This case in the chosen coordinates corresponds to the Case Ib) or Ic) in
Marchesiello and Šnobl (2020) with constants

$a_{1}=a_{2}=0,\quad a_{3}=2b,\quad b_{1}=b_{2}=0,\quad b_{3}=-4w_{0}.$ (57)

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

$\displaystyle X_{3}$
$\displaystyle=L_{3}^{A}p_{2}^{A}-\frac{4bx^{1}}{(x^{2})^{2}}\left(p_{3}^{A}+{\frac{2b}{(x^{2})^{2}}}\right)-\frac{8w_{0}x^{1}}{(x^{2})^{2}},$
(58)

but it is not independent due to the relation

$(\tilde{X}_{2}^{2}+Y_{1}^{2}-2H)X_{1}+X_{3}^{2}=-4(\tilde{X}_{2}b+2w_{0})(2H-\tilde{X}_{2}^{2}).$
(59)

The solution of the stationary Schrödinger equation in terms of Bessel
functions (Kubů, , equation (2.147)) does not yield a dense set of orthogonal
polynomials. We have not obtained the full spectrum; we expect a purely
essential one but are not able to prove this result.

##  V Conclusion

Motivated by the recent observation that the complex magnetic field is
physically relevant at least in statistical mechanics Peng  _et al._ (2015),
in this article we looked for the first order superintegrable systems of the
cylindrical type that admit complex-valued magnetic field, i.e.
$B^{j}:\mathbb{R}^{3}\to\mathbb{C}$.

We have obtained four such systems, three of them known from Kubů  but
allowing complex coupling constants and one new system (15), where also the
integrals of motion feature complex coefficients.

Drawing on results from Kubů , all these systems separate in cylindrical
coordinates and at least one further coordinate system, namely the Cartesian
or the extension of “complex” coordinate system
$z=x^{1}+ix^{2},\bar{z}=x^{1}-ix^{2},x^{3}$ for the new system. This is a
typical behavior for first order superintegrable systems with magnetic field
if the integrals have the standard form. On the other hand, the
superintegrable system from Kubů  _et al._ (2022) has a nonstandard set of
first order integrals and does not separate in any coordinates system.

The new system is maximally superintegrable with first order integrals, the
Hamiltonian and one quadratic cylindrical integral. The algebra of integrals
closes with two additional quadratic integrals, dependent on the previously
mentioned.

The rigorous definition of these systems poses a formidable challenge. In the
standard Hermitian setting, the calculations are performed in the gauge-
covariant form. With complex magnetic fields we are in the non-Hermitian
setting where the gauge freedom no longer holds, and it is not clear how we
should choose the corresponding vector potential. Leaving the resolution of
this problem as an open challenge, we proceeded by choosing the simplest gauge
that also eliminates the most singular terms.

Even in this case we face another problem: The spectra of the found systems
are not guaranteed to be real (and measurable) as the resulting Hamiltonians
are usually only pseudo-Hermitian and therefore self-adjoint only with respect
to an indefinite scalar product. Standard methods used to infer the spectrum
do not work in this setting. Unlike in the quasi-Hermitian case, where the
modified scalar product is positive definite, we thus cannot interpret these
systems as nonstandard representations of the usual quantum mechanics.

## Acknowledgements

OK thanks David Krejčiřk for encouragement to investigate the topic of complex
magnetic fields and discussions about the spectral analysis of the found
systems.

Computations in this paper were performed using MapleTM 2022.1 by Maplesoft, a
division of Waterloo Maple Inc., Waterloo, Ontario.

OK’s research is supported by the Grant Agency of the Czech Technical
University in Prague, grant No. SGS22/178/OHK4/3T/14.

##  Appendix A Detailed solution to the determining equations

The condition $[Y,H]=0$ separates into coefficients of the second, first and
zeroth order derivatives. The first order ones can be solved for the partial
derivatives $\partial_{\alpha}m,\alpha\in\\{r,\phi,Z\\}$. Assuming sufficient
smoothness, we use the compatibility conditions

$\partial_{\alpha\beta}m=\partial_{\beta\alpha}m,\quad\alpha,\beta\in\\{r,\phi,Z\\}$
(60)

to arrive at the division into subcases. Namely, differentiating
$\partial_{r\phi}m=\partial_{\phi r}m$ twice with respect to $Z$, we get the
simple condition

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

Together with the condition derived by differentiating
$\partial_{rZ}m=\partial_{Zr}m$ with respect to $Z$ twice, namely

$\displaystyle\left(\left({k_{4}}Z-{k_{2}}\right)\mu^{\prime\prime\prime}\left(Z\right)+4{k_{4}}\mu^{\prime\prime}\left(Z\right)\right)\cos\left(\phi\right)+\left(\left({k_{5}}Z+{k_{1}}\right)\mu^{\prime\prime\prime}\left(Z\right)+4k_{5}\mu^{\prime\prime}\left(Z\right)\right)\sin\left(\phi\right)=0,$
(62)

we find that the assumption $k_{4}\neq 0$ leads to $\mu(Z)=\mu_{1}Z+\mu_{2}$.
The case $k_{5}\neq 0$ is equivalent to $k_{4}\neq 0$ by a rotation of the
coordinate axes.

If both $k_{4}=k_{5}=0$, we similarly consider
$\partial_{Z}(\partial_{r\phi}m=\partial_{\phi r}m)$ and
$\partial_{Z}(\partial_{rZ}m=\partial_{Zr}m)$, which read

$-2k_{3}\mu^{\prime\prime}\\!\left(Z\right)r^{4}=0,\quad(k_{1}\sin\\!\left(\phi\right)-k_{2}\cos\\!\left(\phi\right))\mu^{\prime\prime}\\!\left(Z\right)r^{4}=0.$
(63)

Eliminating all three constants here would leave only $k_{6}$, which implies
that the integral $Y$ is a reduced version of $X_{1}$, not an additional
independent integral. We can therefore conclude that

$\mu\left(Z\right)=\mu_{1}Z+\mu_{2}$ (64)

holds in both cases.

Differentiating $\partial_{r\phi}m=\partial_{\phi r}m$ with respect to $r$ and
$Z$ now yields

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

We therefore split our considerations into two cases:

  1. 1.

$k_{4}=k_{5}=0$,

  2. 2.

$k_{4}\neq 0$ without loss of generality, i.e.
$\rho(r)=\rho_{3}r^{2}+\rho_{1}+\rho_{2}r^{-1}$.

###  A.1 Case 1 \- $k_{4}=k_{5}=0$

In this case we shall consider two complementary possibilities: $k_{1}\neq 0$
or $k_{2}\neq 0$, and $k_{1}=k_{2}=0$ with $k_{3}=ck_{6}$ a nontrivial
combination, i.e. $c\neq 0$ and $k_{6}\neq 0$.

Let us treat the latter case first. It does not yield anything new in the end.
Substituting the assumptions $k_{1}=k_{2}=0$ and $k_{3}=ck_{6}$ into the
compatibility conditions (60), we immediately see that the magnetic field $B$
depends only on $r$. A straightforward calculation using the equations for $Y$
and the first order equations from $[X_{1},H]=0,[X_{2},H]=0$ and
$[X_{1},X_{2}]=0$ yield that also the scalar potential $W$ depends on $r$
only, which means the coordinates $\phi$ and $Z$ are cyclic, e.g. we get no
additional integral, only the reduced cylindrical ones needed for
integrability.

We can therefore proceed with $k_{1}\neq 0$ without loss of generality. (This
is equivalent to $k_{2}\neq 0$ by a rotation of the coordinate system.) The
compatibility conditions (60) differentiated with respect to $r$ twice yield

$\rho\left(r\right)=\frac{{\rho_{1}}r+{\rho_{2}}\ln\left(r\right)r+{\rho_{3}}r^{3}+{\rho_{4}}}{r},\quad\sigma\left(r\right)=\frac{{k_{2}}\mu_{1}r^{2}}{2{k_{1}}}+\frac{{\sigma_{1}}}{2r}+\frac{{\sigma_{2}}}{6r^{2}}+{\sigma_{3}}r+{\sigma_{4}}.$
(66)

Plugging this result back into the compatibility conditions, we immediately
find that $\sigma_{3}=0$ and $k_{3}\mu_{1}=0$.

When we choose $\mu_{1}\neq 0$, the compatibility conditions for $m$ can be
solved for the remaining auxiliary functions $\tau(\phi)$ and $\psi(\phi)$,
but the remaining equations for integrals $X_{1}$ and $X_{2}$ lead to a
contradiction.

We therefore continue with $\mu_{1}=0.$ The terms without $r$ in the
compatibility conditions for $m$ can be solved for $\psi(\phi)$ and
$\tau(\phi)$ in general as they do not contain $k_{3}$ nor $k_{6}$. The last
remaining terms in the same equations yield further constraints on the
constants in these functions, leading to 2 possible magnetic fields (absorbing
$\rho_{3}$ into $\mu_{2}$ by a redefinition),

$\begin{split}{B^{Z}}\left(r,\phi,Z\right)={}&\mu_{2}r-\frac{2k_{1}^{2}{\psi_{2}}+2k_{1}k_{2}{\psi_{4}}+2k_{2}^{2}{\psi_{3}}-k_{2}^{3}\rho_{4}}{2\left(k_{2}\cos\left(\phi\right)-k_{1}\sin\left(\phi\right)\right)^{3}r^{2}},\\\
{B^{\phi}}\left(r,\phi,Z\right)={}&-\frac{\left(k_{1}^{2}+k_{2}^{2}\right){\sigma_{2}}-12{\tau_{1}}}{12\left(k_{2}\cos\left(\phi\right)-k_{1}\sin\left(\phi\right)\right)^{2}r^{3}},\\\
{B^{r}}\left(r,\phi,Z\right)={}&-\frac{\left(\left(k_{1}^{2}+k_{2}^{2}\right){\sigma_{2}}-12{\tau_{1}}\right)\left(k_{1}\cos\left(\phi\right)+k_{2}\sin\left(\phi\right)\right)}{12\left(k_{2}\cos\left(\phi\right)-k_{1}\sin\left(\phi\right)\right)^{3}r^{2}},\end{split}$
(67)

where $\tau_{1},\psi_{i}$ are constants in $\tau(\phi)$ and $\psi(\phi),$ and

$\displaystyle B^{Z}(r,\phi,Z)={}$ $\displaystyle\mu_{2}r,\quad
B^{\phi}(r,\phi,Z)=0,\quad B^{r}(r,\phi,Z)=0.$ (68)

Let us treat these cases separately.

####  A.1.1 Nonconstant magnetic field

As we have noted, the compatibility conditions (60) have been solved, so we
can obtain the form of $m$ from the first order equations for integral $Y$.
The first order equations for $X_{2}$ now yield the function $m_{2}$, and
imply $\mu_{2}=0$ with $W=W(r,\phi)$, which the zeroth order equation for $Y$
reduces to a function of one variable,

$W=W(k_{1}r\sin(\phi)-k_{2}r\cos(\phi))=W(k_{1}y-k_{2}x).$ (69)

At this moment the equations for $Y$ and $X_{2}$ are solved. The commutation
$[X_{1},X_{2}]=0$ is solved once we set $m_{1}(r,\phi,Z)=m_{1}(r,\phi)$, so
only the equation $[X_{1},H]=0$ remains to be solved.

The coefficient of $\partial_{Z}$ in $[X_{1},H]=0$ is now only a constraint on
the constants. We obtain 2 solutions, both yield a magnetic field rational in
$y$ after transforming into Cartesian coordinates. The coefficient of
$\partial_{r}$ in $[X_{1},H]=0$ reduces $m_{1}(r,\phi)$ further. The last
first order equation, coefficient of $\partial_{\phi}$, and the zeroth order
equation both contain $W$. For them to be compatible with (69), the magnetic
field reduces to the form considered in Case 3 of Subsection IV.2 or to a
constant magnetic field, which we treat in the following subsection.

####  A.1.2 Constant magnetic field

Here the calculations proceed in the same way as for real magnetic field, see
e.g. Kubů . There is no problem to solve all first order equations, which
yield the functions $m,m_{1},m_{2}$ as well as the separation of $W$,

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

The zeroth order equation for $Y$ now reads

$\left({k_{1}}\cos\\!\left(\phi\right)+{k_{2}}\sin\\!\left(\phi\right)\right)W_{1}^{\prime}(r)+k_{3}W_{3}^{\prime}\\!\left(Z\right)+\frac{{k_{6}}W_{2}^{\prime}\\!\left(\phi\right)}{r^{2}}\\\
-\frac{\left({k_{1}}\sin\left(\phi\right)-{k_{2}}\cos\left(\phi\right)\right)W_{2}^{\prime}\left(\phi\right)+2\left({k_{1}}\cos\\!\left(\phi\right)+{k_{2}}\sin\\!\left(\phi\right)\right){W_{2}}\\!\left(\phi\right)}{r^{3}}=0.$
(71)

Zeroth order equation for $X_{1}$ reads

$W_{2}^{\prime}(\phi)\mu_{2}+\left({c_{1}}\cos\\!\left(\phi\right)+{c_{2}}\sin\\!\left(\phi\right)\right)W_{1}^{\prime}(r)+c_{3}W_{3}^{\prime}\\!\left(Z\right)+\frac{{c_{6}}W_{2}^{\prime}\\!\left(\phi\right)}{r^{2}}\\\
-\frac{\left({c_{1}}\sin\left(\phi\right)-{c_{2}}\cos\left(\phi\right)\right)W_{2}^{\prime}\left(\phi\right)+2\left({c_{1}}\cos\\!\left(\phi\right)+{c_{2}}\sin\\!\left(\phi\right)\right){W_{2}}\\!\left(\phi\right)}{r^{3}}=0,$
(72)

where $c_{i}$ are redefined constants from the auxiliary functions
$\rho(r),\sigma(r),\tau(\phi),\psi(\phi)$ not appearing in the magnetic field,
i.e. they appear in the first order terms of integrals $X_{1},X_{2}.$ The
assumption $k_{1}\neq 0$ allows us to solve for $W_{1}(r)$ by multiplying
equation (71) by $r^{3}$ and differentiating 3 times with respect to $r$.
Separating the equations (71) and (72) into coefficients of $r$ now eliminates
some constants from $W_{1}(r)$ and with $\mu_{2}\neq 0$ (nonvanishing magnetic
field) yields $W_{2}(\phi)$ equal to a constant that exactly cancels the
$r^{-2}$ term in $W_{1}(r)$. We conclude that

$W(r,\phi,Z)=w_{1}\text{ or }W(r,\phi,Z)=W_{3}(Z)$ (73)

depending on the values of constants $k_{3}$ and $c_{3}$. The constant in
$w_{1}$ in the first case is irrelevant, so we set it to zero.

We therefore arrived at systems analyzed in Kubů ; Marchesiello and Šnobl
(2020) and listed as Cases 1 and 2, respectively, in Subsection IV.2.

We conclude that the case $k_{4}=k_{5}=0$ yields only systems that were
already known from Kubů  but with complex coupling constants. We analyzed them
further in Subsection IV.2.

###  A.2 Case 2 \- $k_{4}\neq 0$

Differentiating all compatibility conditions (60) with respect to $r$ and $Z$,
the functions $\rho(r)$ and $\sigma(r)$ are determined to read

$\rho(r)=\rho_{3}r^{2}+\rho_{1}+\rho_{2}r^{-1},\quad\sigma(r)=\sigma_{3}+\sigma_{1}r^{-2},$
(74)

and we also find that $\mu_{1}=0$.

We plug these results into the compatibility conditions (60) again. The
highest order terms in $r$ imply $\mu_{2}=\rho_{3}$. Next, the terms with $Z$
can be solved for $\tau(\phi)$ and $\psi(\phi)$. The remaining terms in the
compatibility conditions now contain constants, some of which are therefore
eliminated, leading to only one solution with nontrivial magnetic field.
Namely, our auxiliary functions read

$\displaystyle\mu\left(Z\right)=\rho_{3},\quad\rho(r)=\rho_{3}r^{2}+\rho_{1}+\rho_{2}r^{-1},\quad\sigma(r)=\sigma_{3}+\sigma_{1}r^{-2},$
(75)
$\displaystyle\psi\left(\phi\right)=2\psi_{3}\cos(\phi)+i\psi_{2}e^{i\phi}+\rho_{2},\quad\tau(\phi)=\sigma_{1}-\tau_{1}e^{2\phi},$
(76)

and the corresponding magnetic field reads

$B^{r}=-i\frac{\tau_{2}e^{2i\phi}}{r^{2}},\quad
B^{\phi}=-\frac{\tau_{2}e^{2i\phi}}{r^{3}},\quad B^{Z}=0.$ (77)

We thus assume $\tau_{2}\neq 0$ in this subsection.

The compatibility conditions for $m$ are now solved. Substituting the
auxiliary functions into the first order equation for $Y$ yields

$m\left(r,\phi,Z\right)=\frac{\tau_{2}(2ik_{4}re^{i\phi}+ik_{3}e^{2i\phi})}{2r^{2}}.$
(78)

We proceed to solve the first order equations for $X_{1}$. Those corresponding
to $\partial_{\phi}$ and $\partial_{r}$ yield $m_{1}$ and $W$. We subsequently
substitute the results into the equation corresponding to $\partial_{Z}$,
differentiate it with respect to $r$ and obtain the terms multiplying
$e^{2i\phi},e^{3i\phi}$

$\frac{8\left(i\psi_{2}+\psi_{3}\right)\tau_{2}}{r},\quad 2i\rho_{1}\tau_{2},$
(79)

which due to our assumption $\tau_{2}\neq 0$ imply $\rho_{1}=0$ and
$\psi_{2}=i\psi_{3}$.

The zeroth order equation for $X_{1}$ and $Y$ yield the final form of the
scalar potential

$W(r,\phi,Z)=\frac{w_{1}e^{2i\phi}}{2r^{2}}-\frac{\tau_{2}^{2}e^{4i\phi}}{8r^{4}}.$
(80)

The equations for $X_{2}$ can be solved for $m_{2}$ and imply $\rho_{3}=0$.

Thus we obtain the new system written in Cartesian coordinates (15) and
analysed in Subsection (IV.1), with $\tau_{2}=b$.

## References

  * Shapovalov  _et al._ (1972) V. N. Shapovalov, V. G. Bagrov, and A. G. Meshkov, Sov. Phys. J. 15, 1115 (1972).
  * Bérubé and Winternitz (2004) J. Bérubé and P. Winternitz, J. Math. Phys. 45, 1959 (2004).
  * McSween and Winternitz (2000) E. McSween and P. Winternitz, J. Math. Phys. 41, 2957 (2000).
  * Marchesiello  _et al._ (2015) A. Marchesiello, L. Šnobl, and P. Winternitz, J. Phys. A 48, 395206, 24 (2015).
  * Zhalij (2015) A. Zhalij, J. Phys.: Conf. Ser. 621, 012019 (2015).
  * Marchesiello and Šnobl (2018) A. Marchesiello and L. Šnobl, SIGMA Symmetry Integrability Geom. Methods Appl. 14, 092, 11 (2018).
  * Bertrand and Šnobl (2019) S. Bertrand and L. Šnobl, J. Phys. A: Math. Theor. 52, 195201, 25 (2019).
  * Fournier  _et al._ (2020) F. Fournier, L. Šnobl, and P. Winternitz, J. Phys. A: Math. Theor. 53, 085203 (2020).
  * Marchesiello and Šnobl (2022) A. Marchesiello and L. Šnobl, J. Phys. A: Math. Theor. 55, 145203 (2022).
  * Kubů  _et al._ (2022) O. Kubů, A. Marchesiello, and L. Šnobl, “New classes of quadratically integrable systems in magnetic fields: the generalized cylindrical and spherical cases,” (2022), 2206.15305 .
  * İ.Yurduşen  _et al._ (2021) İ.Yurduşen, O. O. Tuncer, and P. Winternitz, J. Phys. A: Math. Theor. 54, 305201 (2021).
  * Peng  _et al._ (2015) X. Peng, H. Zhou, B. B. Wei, J. Cui, J. Du, and R. B. Liu, Phys. Rev. Lett. 114, 010601 (2015).
  * Fernández and García-Muñoz (2022) D. J. Fernández and J. D. García-Muñoz, Eur. Phys. J. Plus 137 (2022), 10.1140/epjp/s13360-022-03221-5.
  * Jaramillo (2015) J. L. Jaramillo, Classical Quantum Gravity 32, 132001 (2015).
  * Mostafazadeh (2002) A. Mostafazadeh, J. Math. Phys. 43, 205 (2002).
  * Scholtz  _et al._ (1992) F. G. Scholtz, H. B. Geyer, and F. J. W. Hahne, Ann. Physics 213, 74 (1992).
  * Dieudonné (1961) J. Dieudonné, in  _Proc. Internat. Sympos. Linear Spaces (Jerusalem, 1960)_ , Vol. 115122 (Pergamon, Oxford, 1961).
  * Krejčiřk and Siegl (2015) D. Krejčiřk and P. Siegl, in  _Non-Selfadjoint Operators in Quantum Physics : Mathematical Aspects_ , edited by F. Bagarello, J.-P. Gazeau, F. H. Szafraniec, and M. Znojil (Wiley, 2015) pp. 241–291.
  * Krejčiřík  _et al._ (2015) D. Krejčiřík, P. Siegl, M. Tater, and J. Viola, J. Math. Phys. 56, 103513 (2015).
  * Krejčiřík (2019) D. Krejčiřík, SIAM J. Math. Anal. 51, 790 (2019).
  * Bognár (2012) J. Bognár,  _Indefinite Inner Product Spaces_ (Springer Berlin, Heidelberg, 2012).
  * Williams (1969) J. P. Williams, Proc. Amer. Math. Soc. 20, 121 (1969).
  * (23) O. Kubů, Master Thesis, Czech Technical University in Prague, 2020, 2210.02393 .
  * Marchesiello  _et al._ (2018) A. Marchesiello, L. Šnobl, and P. Winternitz, J. Phys. A 51, 135205, 24 (2018).
  * Dorizzi  _et al._ (1985) B. Dorizzi, B. Grammaticos, A. Ramani, and P. Winternitz, J. Math. Phys. 26, 3070 (1985).
  * Kubů and Šnobl (2022) O. Kubů and L. Šnobl, “Quantum cylindrical integrability in magnetic fields,” (2022), to appear in Proceedings of 34th International Colloquium on Group Theoretical Methods in Physics, 2210.03468 .
  * Marchesiello and Šnobl (2020) A. Marchesiello and L. Šnobl, SIGMA Symmetry Integrability Geom. Methods Appl. 16, 015, 35 pages (2020).
  * Kubů  _et al._ (2021) O. Kubů, A. Marchesiello, and L. Šnobl, J. Phys. A: Math. Theor. 54, 425204 (2021).
  * Miller  _et al._ (2013) W. Miller, Jr., S. Post, and P. Winternitz, J. Phys. A 46, 423001, 97 (2013).
  * Marquette and Quesne (2022) I. Marquette and C. Quesne, J. Phys. A: Math. Theor. 55, 415203 (2022).
  * (31) P. Siegl, Master thesis, Czech Technical University in Prague, 2008\.
  * (32) DLMF, “NIST Digital Library of Mathematical Functions,” Release 1.0.25 of 2019-12-15, F. W. J. Olver et al., eds.
  * Marchesiello and Šnobl (2017) A. Marchesiello and L. Šnobl, J. Phys. A 50, 245202, 24 (2017).


* (100) L. Marsicano, M. Battaglieri, M. Bondi’, C. D. R. Carvajal, A. Celentano, M. De Napoli et al., “Dark photon production through positron annihilation in beam-dump experiments”, _Phys. Rev. D_ 98 (2018) 015031 [1802.03794].
  * (101) C. W. Fabjan and F. Gianotti, “Calorimetry for particle physics”, _Rev. Mod. Phys._ 75 (2003) 1243.
  * (102) A. Belyaev, N. D. Christensen and A. Pukhov, “CalcHEP 3.4 for collider physics within and beyond the Standard Model”, _Comput. Phys. Commun._ 184 (2013) 1729 [1207.6082].
  * (103) J. Alwall, R. Frederix, S. Frixione, V. Hirschi, F. Maltoni, O. Mattelaer et al., “The automated computation of tree-level and next-to-leading order differential cross sections, and their matching to parton shower simulations”, _JHEP_ 07 (2014) 079 [1405.0301].
  * (104) D. Forbes, C. Herwig, Y. Kahn, G. Krnjaic, C. Mantilla Suarez, N. Tran et al., “New Searches for Muonphilic Particles at Proton Beam Dump Spectrometers”, 2212.00033.
  * (105) O. Mattelaer and K. Ostrolenk, “Speeding up MadGraph5_aMC@NLO”, _Eur. Phys. J. C_ 81 (2021) 435 [2102.00773].