# Integrable systems in magnetic fields: the generalized parabolic cylindrical
case

O. Kubů1, A. Marchesiello2 and L. Šnobl1

###### Abstract

This article is a contribution to the classification of quadratically
integrable systems with vector potentials whose integrals are of the
nonstandard, nonseparable type. We focus on generalized parabolic cylindrical
case, related to non-subgroup-type coordinates. We find 3 new systems, two
with magnetic fields polynomial in Cartesian coordinates and one with
unbounded exponential terms. The limit in the parameters of the integrals
yields a new parabolic cylindrical system; the limit of vanishing magnetic
fields leads to the free motion. This confirms the conjecture that non-
subgroup type integrals can be related to separable systems only in a trivial
manner.

1Czech Technical University in Prague, Faculty of Nuclear Sciences and
Physical Engineering, Department of Physics, Břehová 7, 115 19 Prague 1, Czech
Republic

2Czech Technical University in Prague, Faculty of Information Technology,
Department of Applied Mathematics, Thákurova 9, 160 00 Prague 6, Czech
Republic

E-mail<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
Keywords: integrability, magnetic field, generalized parabolic cylindrical,
nonseparable

##  1 Introduction

This article is a step forward to the classification of quadratically
integrable and superintegrable systems with vector potentials that possess
integrals not connected to separability. We restrict our considerations to the
classical nonrelativistic particle on the three-dimensional (3D) Euclidean
space $\mathbb{E}_{3}$ and assume time-independent Hamiltonian.

Let us recall that an $n$-dimensional classical Hamiltonian system is called
integrable if it admits $n$ functionally independent first integrals
(including the Hamiltonian) that are in involution with respect to the Poisson
bracket. If there exist more than $n$ independent integrals, out of which $n$
are in involution, the system is said to be superintegrable. Maximally
superintegrable system has $2n-1$ independent integrals.

In 3D space, integrability implies the existence of two independent integrals
$X_{1}$, $X_{2}$ in involution in addition to the Hamiltonian, i.e.

$\\{X_{1},H\\}=0,\quad\\{X_{2},H\\}=0,\quad\\{X_{1},X_{2}\\}=0.$ (1.1)

The solvability of these systems (in quadratures for integrable, algebraically
for maximally superintegrable) together with their stability under
perturbations makes these systems a bedrock for the study of more complicated
and less regular systems. For these reasons, a thorough classification and
investigation of their properties merits substantial effort.

Most of the results in the field of superintegrability and the related
classification of integrable systems, starting with the seminal papers from
the 1960s by Smorodinsky, Winternitz et al. [9, 33, 18], have focused on the
so called “natural” Hamiltonians on $\mathbb{E}_{3}$, i.e. Hamiltonians with
quadratic but no linear terms in momenta in Euclidean 3D space
$\mathbb{E}_{3}$. Classification of such integrable and superintegrable
systems with integrals of motion at most quadratic in momenta, based on the
1:1 correspondence between (multi)separable and (super)integrable systems, has
been completed on 2D and 3D spaces of constant curvature including
$\mathbb{E}_{3}$, see the review [27] and the monograph [25]. Naturally, many
articles considering higher order integrals followed, see e.g. [28, 19, 6, 7]
and references therein. Although systems with only higher order integrals do
not separate on the coordinate space, separating coordinates on the phase
space for the Post-Winternitz system [28] and its Drach-Holt generalization
[4] were found using the Haantjes algebra approach in [32] and [30],
respectively. See also [15] for the first steps toward the classification of
“natural” quadratic maximally superintegrable systems in an arbitrary
dimension using methods of algebraic geometry.

However, the “natural” Hamiltonians certainly do not exhaust all physically
relevant systems. For example, the linear terms in the Hamiltonian

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

are necessary for models of a charged particle in an external electromagnetic
field or systems in uniform rotation [34]. Despite their relevance, such
systems have not been classified yet, probably due to the computational
difficulty stemming from the broken 1:1 correspondence between separability
and quadratic integrability [3, 24]. Namely, the separation now requires at
least one first-order integral [31, 2]. (We note that there is a rather vast
literature on magnetic monopole systems in many different spaces, see e.g.
[13] for references. However, these are exceptional, particular configurations
of the magnetic field, not a classification.)

As a consequence of the broken correspondence, there exist pairs of quadratic
integrals on $\mathbb{E}_{3}$ leading to integrability that do not correspond
to any orthogonal separable coordinates [21]. As the algebraic classification
of quadratic terms showed [20], there exist classes of integrals with
additional quadratic terms extending the original, separable ones, and also
with no connection to separability at all. We call such integrals nonstandard.
To confirm that these classes are not empty, we have to consider the lower
order terms in momenta as well in each class; we have done it for some of them
[16, 11, 12]. Quite a few nonstandard integrable systems were found, though a
complete study of all classes has not been achieved yet. The helical undulator
in infinite solenoid [16], a peculiar system with 3 first order integrals not
separable in orthogonal coordinates on the configuration space, remains the
sole superintegrable example (also in the relativistic setting [10]) and is
relevant e.g. for the study of free electron lasers [1].

In this paper, we consider one of the more complicated classes of [20], case
(j), which generalizes the parabolic cylindrical structure of integrals by
adding leading order terms that depend on three additional parameters. Our
goal is to classify the integrable systems in this class as well as to shed
more light on the conjecture in our previous papers [16, 12]: In the presence
of magnetic fields, there appears to be a discrepancy between the subgroup-
and non-subgroup-type coordinates111The subgroup type coordinates in 3D are
defined by the property that some subgroup chain
$E_{3}\supset\tilde{G}\supset G_{M}$ (1.3) exists where $E_{3}$ is the
Euclidean group and $G_{M}$ is its Abelian subgroup such that the leading
order terms of their corresponding quadratic integrals of motion are
second–order Casimir operators of subgroups [14, 26]. already at the level of
quadratic integrals. Namely, all known quadratically integrable systems with
vector potentials separate in Cartesian or cylindrical coordinates in the
limit of vanishing magnetic field, even when the integrals generalize the
leading order structure corresponding to separation in a different coordinate
system. In other words, the limit does not lead to separability in non-
subgroup coordinates (except the trivial free motion).

The structure of our article is as follows. In Section 2 we introduce our
formalism for magnetic fields and the generalized parabolic cylindrical case,
i.e. the explicit form of the integrals $X_{1}$ and $X_{2}$ depending on
parameters $a,b,c$. The following sections contain the detailed classification
of integrable systems in this class. Section 3 deals with the case $a\neq 0$
and contains one new family, (3.20)–(3.21), while Section 4 corresponding to
$a=0$ yields 2 more families, (4.6)–(4.7) and (4.15)–(4.16). In Section 5 we
summarize our results and consider the limits of vanishing magnetic field as
well as vanishing parameters. The latter yields a new standard, not
generalized parabolic cylindrical system. In Appendix A we present the scalar
term in integral $X_{1}$ for the system (3.20)–(3.21), obtained by solving
equations (3.24).

##  2 The generalized parabolic cylindrical case

Before we present the pair of integrals $X_{1}$ and $X_{2}$ defining our case,
we have to recall a few facts about magnetic fields and introduce our
notation.

As we have noted, the Hamiltonian of a system modeling a charged particle in
the three-dimensional Euclidean space $\mathbb{E}_{3}$ immersed in external
electromagnetic field is (1.2), where we have chosen the units $m=1,e=-1$,
considering the electron as our main example. We denote the Cartesian
coordinates by $\vec{x}=(x,y,z)$ and the corresponding linear momenta
$\vec{p}=(p_{x},p_{y},p_{z}).$

The key feature of such systems is the invariance of the dynamically relevant
magnetic field

$\vec{B}(\vec{x})=\nabla\times\vec{A}(\vec{x})$ (2.1)

with respect to the time-independent gauge transformation

$\vec{A}^{\prime}(\vec{x})=\vec{A}(\vec{x})+\mathrm{d}\chi(\vec{x}),\quad
W^{\prime}(\vec{x})=W(\vec{x}).$ (2.2)

On the other hand, the scalar potential $V(\vec{x})$ (i.e., the momentum-free
term in the Hamiltonian (1.2))

$V(\vec{x})=W(\vec{x})+\frac{1}{2}\left(A_{1}(\vec{x})^{2}+A_{2}(\vec{x})^{2}+A_{3}(\vec{x})^{2}\right)$
(2.3)

is affected by the transformation (2.2), as well as the first and zero order
terms in momenta of quadratic integrals.

In light of that we find it useful to write the integrals of motion in terms
of gauge covariant linear and angular momenta

$p_{i}^{A}=p_{i}+A_{i}(\vec{x}),\quad\ell^{A}_{i}=\epsilon_{ijk}x_{j}p^{A}_{k},$
(2.4)

with $\epsilon_{ijk}$ the totally anti-symmetric Levi-Civita symbol,
$\epsilon_{123}=1$. In this way, the conditions coming from (1.1) become
explicitly gauge invariant and depend only on the magnetic field
$\vec{B}(\vec{x})$, not on the vector potential $\vec{A}(\vec{x})$.

Let us now introduce case (j) from [20]: The initial form of the integrals as
derived in [20] is

$\displaystyle X_{1}={}$
$\displaystyle\ell^{A}_{x}p^{A}_{x}+a\ell^{A}_{y}p^{A}_{y}-(a+1)\ell^{A}_{z}p^{A}_{z}+\omega\left(\ell^{A}_{x}p^{A}_{z}-\ell^{A}_{z}p^{A}_{x}\right)+2bp^{A}_{x}p^{A}_{y}+c\left((p^{A}_{y})^{2}-(p^{A}_{z})^{2}\right)$
$\displaystyle+\sum_{j}s_{1}^{x_{j}}(\vec{x})p_{j}^{A}+m_{1}(\vec{x}),$ (2.5)
$\displaystyle X_{2}={}$
$\displaystyle(p^{A}_{x})^{2}+\frac{6\omega}{4a-1}p^{A}_{x}p^{A}_{z}+\frac{a+2}{4a-1}(p^{A}_{y})^{2}-\frac{5a+1}{4a-1}(p^{A}_{z})^{2}+\sum_{j}s_{2}^{x_{j}}(\vec{x})p_{j}^{A}+m_{2}(\vec{x}),$
(2.6)

with the parameters

$\omega=\sqrt{1+a-2a^{2}},\quad-\frac{1}{2}<a\leq 0,\quad 0\leq b,\quad
c\in\mathbb{R}.$ (2.7)

Our task is to find the functions $s_{i}^{x_{j}}$ and $m_{i}$ together with
the magnetic field $\vec{B}$ and the electrostatic potential $W$ such that
these integrals satisfy the integrability conditions (1.1).

The integrals contain 3 independent parameters $a,b,c$ and the quantity
$\omega$ is nonvanishing within the admissible range of $a$. The analysis in
[20] implies that the choice $a=0,\omega=1,c=b=0$ leads to the standard
parabolic cylindrical case with integrals $Y_{1}\equiv\ell_{z}p_{x}+\ldots$
and $Y_{2}\equiv p_{z}^{2}+\ldots$ upon a suitable rotation of the coordinate
system and a linear combination of integrals $X_{1},X_{2}$ and $H$. As we are
mainly interested in the integrals of the generalized form, we will hereafter
assume that at least one of the parameters $a,b,c$ is nonvanishing. As we
noted for the elliptic cylindrical case in [12], the determining equations for
the standard integrals (i.e. $a=b=c=0$) are less overdetermined and therefore
harder to solve. Their solutions include the partially separated (2+1)
dimensional system, but the parabolic 2D case has not been completely solved
yet [29], giving little hope for the full classification in 3D at this moment
as well.

The form of the integrals $X_{1}$ and $X_{2}$ above was suitable during the
algebraic study of the quadratic terms in [20], but for our purposes a
rotation around the $y$-axis dependent on the parameter $a$,

$\displaystyle X={}$ $\displaystyle\frac{\omega x}{\psi}+\frac{\psi
z}{a+2},\quad Y=y,\quad Z=-\frac{\psi x}{a+2}+\frac{\omega z}{\psi},$ (2.8)
$\displaystyle{p_{X}}={}$
$\displaystyle\frac{\omega{p_{x}}}{\psi}+\frac{\psi{p_{z}}}{a+2},\quad{p_{Y}}={p_{y}},\quad{p_{Z}}=-\frac{\psi{p_{x}}}{a+2}+\frac{\omega
p_{z}}{\psi},$

where we have defined a real nonvanishing quantity (within our range for $a$)

$\psi=\sqrt{-a^{2}-a+2},$ (2.9)

yields a simpler form with a better physical interpretation. This
transformation is consistent for any $a$ in the range (2.7), including $a=0$,
and represents a rotation by the angle $\varphi$ satisfying
$\cos(\varphi)=\frac{\omega}{\psi}$.

In the rotated Cartesian coordinates, which we hereafter denote again by
$(x,y,z)$ for simplicity, the integral $X_{2}$ simplifies significantly if we
take its linear combination with the Hamiltonian and rescale it. We have

$\displaystyle X_{1}={}$
$\displaystyle-3a\ell_{z}^{A}p_{z}^{A}-2\omega\ell_{z}^{A}p_{x}^{A}+\frac{a-1}{a+2}c(p_{x}^{A})^{2}+c(p_{y}^{A})^{2}-\frac{2a+1}{a+2}c(p_{z}^{A})^{2}$
(2.10) $\displaystyle-\frac{2\omega
c}{a+2}p_{x}^{A}p_{z}^{A}+\frac{2\psi}{\omega}bp_{x}^{A}p_{y}^{A}-\frac{2b\psi}{a+2}p_{y}^{A}p_{z}^{A}+\sum_{j}s_{1}^{x_{j}}(\vec{x})p_{j}^{A}+m_{1}(\vec{x}),$
$\displaystyle X_{2}={}$
$\displaystyle(p_{z}^{A})^{2}+\sum_{j}s_{2}^{x_{j}}(\vec{x})p_{j}^{A}+m_{2}(\vec{x}),$
(2.11)

where the lower order terms were without any loss of generality redefined
appropriately after the rotation.

Given the form of the first term in $X_{1}$, it is convenient to treat the
cases with $a\neq 0$ and $a=0$ separately, as we shall do in the following two
sections.

##  3 Parameter $a$ nonvanishing

When $a\neq 0$, we simplify the integrals by an additional shift of
coordinates,

$\tilde{x}=x+\frac{2b\psi}{3a\left(a+2\right)},\quad\tilde{y}=y-\frac{2c\omega}{3a\left(a+2\right)},\quad\tilde{z}=z,$
(3.1)

which eliminates some of the mixed terms in $X_{1}$ while leaving $X_{2}$
intact (up to redefinition of lower order terms)

$\displaystyle X_{1}={}$
$\displaystyle-3a\ell_{\tilde{z}}^{A}p_{\tilde{z}}^{A}-2\omega
l_{\tilde{z}}^{A}p_{\tilde{x}}^{A}-\frac{\left(5a^{2}-a-4\right)c(p_{\tilde{x}}^{A})^{2}}{3a(a+2)}+c(p_{\tilde{y}}^{A})^{2}$
$\displaystyle-\frac{2a+1}{a+2}c(p_{\tilde{z}}^{A})^{2}-\frac{2b\omega\psi}{3a(a-1)}p_{\tilde{x}}^{A}p_{\tilde{y}}^{A}+s_{1}^{\tilde{x}}p_{\tilde{x}}+s_{1}^{\tilde{y}}p_{\tilde{y}}+s_{1}^{\tilde{z}}p_{\tilde{z}}+m_{1},$
(3.2)

where $s_{1}^{x_{i}}$ and $m_{1}$ are now functions of
$(\tilde{x},\tilde{y},\tilde{z})$.

We proceed by the standard direct method, i.e. we impose the Poisson
commutativity (1.1), separate the equations into different powers in momenta
$p_{i}$, e.g. $p_{i}p_{j}$, and impose that all the coefficients vanish. In
this way we obtain an overdetermined system of PDEs determining the form of
the integrals of motion $X_{1},X_{2}$ as well as the magnetic field $\vec{B}$
and electrostatic potential $W$, hence we call them determining equations.

We present the equations corresponding to the quadratic terms in momenta in
(1.1), first for the integral $X_{1}$

$\displaystyle\partial_{\tilde{x}}s_{1}^{\tilde{x}}=3a\tilde{y}B^{\tilde{y}}+2\omega\left(\tilde{x}+\frac{b\psi}{3a(a-1)}\right)B^{\tilde{z}},$
(3.3a)
$\displaystyle\partial_{\tilde{y}}s_{1}^{\tilde{y}}=3a\tilde{x}B^{\tilde{x}}-2\omega\left(\tilde{x}+\frac{b\psi}{3a(a-1)}\right)B^{\tilde{z}},$
(3.3b)
$\displaystyle\partial_{\tilde{z}}s_{1}^{\tilde{z}}=-3a(\tilde{x}B^{\tilde{x}}+\tilde{y}B^{\tilde{y}}),$
(3.3c)
$\displaystyle\partial_{\tilde{x}}s_{1}^{\tilde{y}}+\partial_{\tilde{y}}s_{1}^{\tilde{x}}=-3a(\tilde{y}B^{\tilde{x}}+\tilde{x}B^{\tilde{y}})+\left(4\omega\tilde{y}-\frac{2c\left(8a^{2}+5a-4\right)}{3a(a+2)}\right)B^{\tilde{z}},$
(3.3d)
$\displaystyle\partial_{\tilde{x}}s_{1}^{\tilde{z}}+\partial_{\tilde{z}}s_{1}^{\tilde{x}}=-2\omega\left(\tilde{x}+\frac{b\psi}{3a\left(a-1\right)}\right)B^{\tilde{x}}-2\left(2\omega\tilde{y}+\frac{c\left(a+2\right)}{3a}\right)B^{\tilde{y}}+3a\tilde{x}B^{\tilde{z}},$
(3.3e)
$\displaystyle\partial_{\tilde{y}}s_{1}^{\tilde{z}}+\partial_{\tilde{z}}s_{1}^{\tilde{y}}=\frac{6c(a+1)}{a+2}B^{\tilde{x}}+2\omega\left(\tilde{x}+\frac{b\psi}{3a(a-1)}\right)B^{\tilde{y}}+3a\tilde{y}B^{\tilde{z}},$
(3.3f)

and for the integral $X_{2}$

$\displaystyle\partial_{\tilde{x}}s_{2}^{\tilde{x}}=0,\quad\partial_{\tilde{y}}s_{2}^{\tilde{y}}=0,\quad\partial_{z}s_{2}^{\tilde{y}}=0,$
(3.4)
$\displaystyle\partial_{\tilde{y}}s_{2}^{\tilde{x}}+\partial_{\tilde{x}}s_{2}^{\tilde{y}}=0,$
$\displaystyle\partial_{\tilde{z}}s_{2}^{\tilde{y}}+\partial_{\tilde{y}}s_{2}^{\tilde{z}}=-2B^{\tilde{x}},$
$\displaystyle\partial_{\tilde{z}}s_{2}^{\tilde{x}}+\partial_{\tilde{x}}s_{2}^{\tilde{z}}=2B^{\tilde{y}}.$

The equations corresponding to the involutivity condition
$\\{X_{1},X_{2}\\}=0$ are rather long, we thus introduce them later when they
are needed and some of the functions above already get reduced by a partial
solution.

Regardless of the shifted coordinates (3.1), equations (3.4) have a standard
form and have therefore been solved e.g. in [8]. We add to these equations
also the condition coming from Maxwell equations,

$\nabla\cdot\vec{B}=\partial_{\tilde{x}}B^{\tilde{x}}+\partial_{\tilde{y}}B^{\tilde{y}}+\partial_{\tilde{z}}B^{\tilde{z}}=0,$
(3.5)

which we solve for $\partial_{z}B^{\tilde{z}}.$ We arrive at

$\begin{split}B^{\tilde{x}}={}&\tfrac{1}{2}(\partial_{\tilde{z}}S^{\tilde{x}}_{21}(\tilde{z})\tilde{x}-\partial_{\tilde{y}}S^{\tilde{z}}_{2}(\tilde{x},\tilde{y})-\partial_{\tilde{z}}S^{\tilde{y}}_{2}(\tilde{z})),\\\
B^{\tilde{y}}={}&\tfrac{1}{2}(\partial_{\tilde{z}}S^{\tilde{x}}_{21}(\tilde{z})\tilde{y}+\partial_{\tilde{x}}S^{\tilde{z}}_{2}(\tilde{x},\tilde{y})+\partial_{\tilde{z}}S^{\tilde{x}}_{22}(\tilde{z})),\\\
B^{\tilde{z}}={}&-S^{\tilde{x}}_{21}(\tilde{z})+B^{\tilde{z}}_{1}(\tilde{x},\tilde{y}),\\\
s_{2}^{\tilde{x}}={}&S^{\tilde{x}}_{21}(\tilde{z})\tilde{y}+S^{\tilde{x}}_{22}(\tilde{z}),\quad
s_{2}^{\tilde{y}}=-S^{\tilde{x}}_{21}(\tilde{z})\tilde{x}+S^{\tilde{y}}_{2}(\tilde{z}),\quad
s_{2}^{\tilde{z}}=S^{\tilde{z}}_{2}(\tilde{x},\tilde{y}).\end{split}$ (3.6)

We next look at the involutivity equation $\\{X_{1},X_{2}\\}=0$. The
coefficient of $p_{\tilde{x}}^{2}$ gives

$-3a^{2}(a-1)\tilde{y}^{2}\partial_{\tilde{z}}S^{\tilde{x}}_{21}(\tilde{z})-3a^{2}(a-1)\tilde{y}\partial_{\tilde{z}}S^{\tilde{x}}_{22}(\tilde{z})+2\omega\left(a(a-1)S^{\tilde{y}}_{2}(\tilde{z})+\frac{b\psi}{3}S^{\tilde{x}}_{21}(\tilde{z})\right)=0,$
(3.7)

which implies that all three functions
$S^{\tilde{x}}_{21}(\tilde{z}),S^{\tilde{x}}_{22}(\tilde{z})$ and
$S^{\tilde{y}}_{2}(\tilde{z})$ are constants. Combining this with the
coefficients of $p_{\tilde{y}}^{2}$ and $p_{\tilde{x}}p_{\tilde{y}}$, we
arrive at

$S^{\tilde{x}}_{21}(\tilde{z})=0,\quad S^{\tilde{y}}_{2}(\tilde{z})=0,\quad
S^{\tilde{x}}_{22}(\tilde{z})=0.$ (3.8)

These results mean that the magnetic field does not depend on $\tilde{z}$,
namely it has the form

$B^{\tilde{x}}=-\tfrac{1}{2}\partial_{\tilde{y}}S^{\tilde{z}}_{2}(\tilde{x},\tilde{y}),\quad
B^{\tilde{y}}=\tfrac{1}{2}\partial_{\tilde{x}}S^{\tilde{z}}_{2}(\tilde{x},\tilde{y}),\quad
B^{\tilde{z}}=B^{\tilde{z}}_{1}(\tilde{x},\tilde{y}).$ (3.9)

The remaining second order equations from $\\{X_{1},X_{2}\\}=0$ are equivalent
to the conclusion that $s_{1}^{x_{i}}$ are independent of $\tilde{z}$ as well.

This suggests that $\tilde{z}$ might be a cyclic coordinate if it also
disappears from the electrostatic potential
$W(\tilde{x},\tilde{y},\tilde{z})$. Indeed, the first order equations for
$X_{2}$,

$\displaystyle\partial_{\tilde{x}}m_{2}=\tfrac{1}{2}S^{\tilde{z}}_{2}(\tilde{x},\tilde{y})\partial_{\tilde{x}}S^{\tilde{z}}_{2}(\tilde{x},\tilde{y}),\quad\partial_{\tilde{y}}m_{2}=\tfrac{1}{2}S^{\tilde{z}}_{2}(\tilde{x},\tilde{y})\partial_{\tilde{y}}S^{\tilde{z}}_{2}(\tilde{x},\tilde{y}),\quad\partial_{\tilde{z}}m_{2}=2\partial_{\tilde{z}}W,$
(3.10)

imply

$m_{2}(\tilde{x},\tilde{y},\tilde{z})=\tfrac{1}{4}\left(S^{\tilde{z}}_{2}(\tilde{x},\tilde{y})\right)^{2}+W_{3}(\tilde{z}),\quad
W(\tilde{x},\tilde{y},\tilde{z})=W_{1}(\tilde{x},\tilde{y})+W_{3}(\tilde{z}).$
(3.11)

$W_{3}(\tilde{z})$ is however immediately eliminated (up to an irrelevant
additive constant) using the first order involutivity equations, because
coefficients of $p_{\tilde{x}}$ and $p_{\tilde{y}}$ there read

$a\tilde{x}\partial_{\tilde{z}}W_{3}(\tilde{z})=0,\quad-a\tilde{y}\partial_{\tilde{z}}W_{3}(\tilde{z})=0.$
(3.12)

With $W_{3}$ eliminated, all equations for $X_{2}$ are satisfied. We only have
to set $m_{1}=m_{1}(\tilde{x},\tilde{y})$ to satisfy the involutivity
condition $\\{X_{1},X_{2}\\}=0$ as well.

We therefore have a cyclic coordinate $\tilde{z}$, which implies the existence
of a first-order integral $\tilde{X}_{2}=p_{\tilde{z}}$ in a suitably chosen
gauge, while $X_{2}=\tilde{X}_{2}^{2}$. We stress that this result holds
whenever $a\neq 0$ regardless of the values of the other parameters $b$ and
$c$.

We have yet to solve the equations for $X_{1}$, starting with (3.3). The left-
hand sides of the first 3 equations are vanishing. The third one, (3.3c),
gives a simple compatibility

$\tilde{x}\partial_{\tilde{y}}S^{\tilde{z}}_{2}(\tilde{x},\tilde{y})-\tilde{y}\partial_{\tilde{x}}S^{\tilde{z}}_{2}(\tilde{x},\tilde{y})=0.$
(3.13)

This implies that we may introduce the shifted polar coordinates, which we
denote $(r,\phi)$,

$\tilde{x}=r\cos(\phi),\quad\tilde{y}=r\sin(\phi),$ (3.14)

in which (3.13) implies $S^{\tilde{z}}_{2}(r,\phi)=S^{\tilde{z}}_{2}(r).$

The combination
$\tilde{y}\cdot\eqref{conds2ord15}-\tilde{x}\cdot\eqref{conds2ord16}$
determines the dependence of $S_{1}^{\tilde{z}}(r,\phi)$ on $\phi$ explicitly,
introducing an arbitrary function $S^{\tilde{z}}_{11}(r)$. We now eliminate
derivatives of $S_{1}^{\tilde{x}}$ and $S^{\tilde{y}}_{1}$ from the equations
(3.3a), (3.3b), and (3.3d) and subsequently use (3.3e) to eliminate
$B^{\tilde{z}}_{1}$, which yields an equation with explicit trigonometric
functions of $\phi$ and unknown functions of $r$, $S^{\tilde{z}}_{11}(r)$ and
$S_{2}^{\tilde{z}}(r)$. Each of the coefficients of $\sin(\phi)$ and
$\cos(\phi)$ must vanish separately, we therefore obtain ODEs for the
functions of $r$, some of them multiplied by $b$ and/or $c$.

For example the coefficients of $\sin(\phi)$ and $\cos(\phi)^{2}\sin(\phi)$
are up to nonvanishing constant factors, respectively,

$b\left(r^{3}\partial_{r}^{4}S^{\tilde{z}}_{2}\\!\left(r\right)+r^{2}\partial_{r}^{3}S^{\tilde{z}}_{2}\\!\left(r\right)-2r\partial_{r}^{2}S^{\tilde{z}}_{2}\\!\left(r\right)+2\partial_{r}S^{\tilde{z}}_{2}\\!\left(r\right)\right)=0$
(3.15)

and

$b\left(r^{3}\partial_{r}^{4}S^{\tilde{z}}_{2}\\!\left(r\right)-3r\partial_{r}^{2}S^{\tilde{z}}_{2}\\!\left(r\right)+3\partial_{r}S^{\tilde{z}}_{2}\\!\left(r\right)\right)=0,$
(3.16)

and we get similar equations including the parameter $c$ together with those
for $S^{\tilde{z}}_{11}(r)$.

This leads to a split: On the one hand $b\neq 0$ or $c\neq 0$ for which
(3.15)–(3.16) and/or similar equations with the parameter $c$ remain ODEs and
we always arrive at the solution (3.17), and on the other hand $b=c=0$ for
which most of the equations for $S_{11}^{\tilde{z}}(r)$ and
$S^{\tilde{z}}_{2}(r)$ vanish identically, cf. (3.25)–(3.26) below.

###  3.1 Parameter $b$ or $c$ nonvanishing

For concreteness, let us consider $b\neq 0$; the case $b=0,$ $c\neq 0$ is very
similar and leads to the same final result.

We continue with solving equations (3.3) as outlined above.

The final form of $S_{11}^{\tilde{z}}$ and $S^{\tilde{z}}_{2}$ obtained from
(3.15)–(3.16) and other ODEs determining also $S^{\tilde{z}}_{11}$ is

$S^{\tilde{z}}_{11}(r)=\tfrac{1}{2}as^{\tilde{z}}_{11}r^{2}+s^{\tilde{z}}_{13},\quad
S^{\tilde{z}}_{2}(r)=s^{\tilde{z}}_{24}r^{2}+s^{\tilde{z}}_{21}.$ (3.17)

This result together with the solution of (3.3e) gives us the magnetic field

$B^{\tilde{x}}(\tilde{x},\tilde{y},\tilde{z})=-s^{\tilde{z}}_{24}\tilde{y},\quad
B^{\tilde{y}}(\tilde{x},\tilde{y},\tilde{z})=s^{\tilde{z}}_{24}\tilde{x},\quad
B^{\tilde{z}}(\tilde{x},\tilde{y},\tilde{z})=\frac{2\omega
s^{\tilde{z}}_{24}}{a}\tilde{y}+\frac{c\left(10a^{2}+13a+4\right)s^{\tilde{z}}_{24}}{9a^{2}\left(a+2\right)}+\frac{s^{\tilde{z}}_{11}}{3}.$
(3.18)

Using the remaining equations from (3.3), we find also $S^{\tilde{x}}_{1}$ and
$S^{\tilde{y}}_{1}$.

We now use the first-order terms in momenta from the involution
$\\{X_{1},X_{2}\\}=0$ again in the polar coordinates (3.14). This determines
the dependence of $W_{1}(r,\phi)$ on $\phi$, leaving an undetermined function
$W_{1}(r)$.

Next, we consider the zeroth order equation for $X_{1}$. The relevant term
reads

$bs^{\tilde{z}}_{24}\left(a+\tfrac{1}{2}\right)\psi
r^{2}\left(a^{4}(a+2)^{2}W_{1}^{\prime}(r)-s^{\tilde{z}}_{24}P(r)\right)=0$
(3.19)

with $P(r)$ a polynomial in $r$ stemming from the dependence of $W$ on $\phi$,
not explicit for brevity.

Given that $s^{\tilde{z}}_{24}=0$ implies constant magnetic field $\vec{B}$
and vanishing electrostatic potential $W=0$, we continue with
$s^{\tilde{z}}_{24}\neq 0$ and solve for $W_{1}(r)$. The last undetermined
function, $m_{1}$, is simply integrated from the first order terms in momenta
in $\\{X_{1},H\\}=0$.

We hereby summarize our solution to the problem with $a\neq 0,b\neq 0$ with a
suitable redefinition of constants for better readability. The magnetic field
$\vec{B}$ depends on 2 arbitrary constants, $b_{x}=s^{\tilde{z}}_{24}$ and
$b_{z}=\frac{c\left(10a^{2}+13a+4\right)}{9a^{2}\left(a+2\right)}s^{\tilde{z}}_{24}+\frac{1}{3}s^{\tilde{z}}_{11}$,

$B^{\tilde{x}}(\tilde{x},\tilde{y},\tilde{z})=-b_{x}\tilde{y},\quad
B^{\tilde{y}}(\tilde{x},\tilde{y},\tilde{z})=b_{x}\tilde{x},\quad
B^{\tilde{z}}(\tilde{x},\tilde{y},\tilde{z})=\frac{2\omega
b_{x}}{a}\tilde{y}+b_{z},$ (3.20)

and the electrostatic potential $W$ has 3 more constants, $w_{1},w_{2}$ and
$w_{3}$,

$\begin{split}W(\tilde{x},\tilde{y},\tilde{z})={}&\Bigg{[}-\frac{\left(a+2\right)^{2}b_{x}\tilde{x}^{4}}{72a^{2}}+\frac{2\left(2a+1\right)b\psi
b_{x}\tilde{x}^{3}}{27a^{3}}+\left(\frac{\left(5a^{2}-4a-4\right)b_{x}\tilde{y}^{2}}{12a^{2}}-\frac{b_{z}\omega\tilde{y}}{3a}\right.\\\
{}&-\left.\frac{b_{z}\left(8a^{2}+5a-4\right)c}{9a^{2}\left(a+2\right)}\right)\tilde{x}^{2}+\left(\frac{2\left(2a+1\right)b\psi
b_{x}\tilde{y}^{2}}{9a^{3}}-\frac{{2bb_{z}\omega\psi\tilde{y}}}{9a^{2}\left(a-1\right)}\right)\tilde{x}\\\
{}&+\frac{\left(31a^{2}-20a-20\right)b_{x}\tilde{y}^{4}}{72a^{2}}+\left(\frac{2\left(8a^{2}+5a-4\right)c\omega
b_{x}}{27\left(a+2\right)a^{3}}-\frac{b_{z}\omega}{3a}\right)\tilde{y}^{3}\\\
{}&+w_{3}(\tilde{x}^{2}+\tilde{y}^{2})+w_{1}\tilde{x}+w_{2}\tilde{y}\Bigg{]}b_{x}.\end{split}$
(3.21)

Given that the coordinate $\tilde{z}\equiv z$ is cyclic (in a suitable gauge),
we have a first order integral

$\tilde{X}_{2}=p_{z}^{A}+\frac{b_{x}}{2}(\tilde{x}^{2}+\tilde{y}^{2})$ (3.22)

instead of $X_{2}$ (which is a quadratic polynomial in $\tilde{X}_{2}$). The
integral $X_{1}$ (3.2) remains quadratic, and we only present the functions in
it for better readability:

$\displaystyle s_{1}^{\tilde{x}}(\tilde{x},\tilde{y},\tilde{z})={}$
$\displaystyle\left(\frac{\left(-5a^{2}+4a+4\right)b_{x}\tilde{y}}{2a}+b_{z}\omega\right)\tilde{x}^{2}-\left(\frac{4b_{x}b\psi\left(2a+1\right)\tilde{y}}{3a^{2}}-\frac{2\omega
b\psi b_{z}}{3a\left(a-1\right)}\right)\tilde{x}$
$\displaystyle-\frac{\left(31a^{2}-20a-20\right)b_{x}\tilde{y}^{3}}{6a}-\left(\frac{2\left(8a^{2}+5a-4\right)c\omega
b_{x}}{3\left(a+2\right)a^{2}}-3b_{z}\omega\right)\tilde{y}^{2}-6aw_{3}\tilde{y}-3aw_{2},$
$\displaystyle s_{1}^{\tilde{y}}(\tilde{x},\tilde{y},\tilde{z})={}$
$\displaystyle-\frac{\left(a+2\right)^{2}b_{x}\tilde{x}^{3}}{6a}+\frac{2b_{x}b\psi\left(2a+1\right)\tilde{x}^{2}}{3a^{2}}+\left(\frac{\left(5a^{2}-4a-4\right)b_{x}\tilde{y}^{2}}{2a}-{2b_{z}\omega\tilde{y}}\right.$
(3.23)
$\displaystyle-\left.\frac{2\left(8a^{2}+5a-4\right)b_{z}c}{3a\left(a+2\right)}+6aw_{3}\right)\tilde{x}+\frac{2b_{x}b\psi\left(2a+1\right)\tilde{y}^{2}}{3a^{2}}-\frac{2\omega
b\psi b_{z}\tilde{y}}{3a\left(a-1\right)}+3aw_{1},$ $\displaystyle
s_{1}^{\tilde{z}}(\tilde{x},\tilde{y},\tilde{z})={}$
$\displaystyle\left(2\omega
b_{x}\tilde{y}-\frac{c\left(a+2\right)b_{x}}{3a}+\frac{3ab_{z}}{2}\right)\tilde{x}^{2}+\frac{2\omega
b\psi b_{x}\tilde{y}\tilde{x}}{3a\left(a-1\right)}+2\omega b_{x}\tilde{y}^{3}$
$\displaystyle-\left(\frac{3\left(a+1\right)cb_{x}}{a+2}-\frac{3ab_{z}}{2}\right)\tilde{y}^{2}.$

The function $m_{1}$ can be obtained by integrating the first order equations,
which in the shifted coordinates (3.1) but without substitution of
(3.20)–(3.23) read

$\displaystyle\partial_{\tilde{x}}m_{1}={}$ $\displaystyle
s^{\tilde{z}}_{1}B^{\tilde{y}}-s^{\tilde{y}}_{1}B^{\tilde{z}}+2\left(2\omega\tilde{y}-\frac{c(5a^{2}-a-4)}{3a(a+2)}\right)\partial_{\tilde{x}}W-2\left(\omega\tilde{x}+\frac{b\psi}{3a(a-1)}\right)\partial_{\tilde{y}}W+3a\tilde{y}\partial_{\tilde{z}}W,$
$\displaystyle\partial_{\tilde{y}}m_{1}={}$ $\displaystyle
s^{\tilde{x}}_{1}B^{\tilde{z}}-s^{\tilde{z}}_{1}B^{\tilde{x}}-2\omega\left(\tilde{x}+\frac{b\psi}{3a(a-1)}\right)\partial_{\tilde{x}}W+2c\partial_{\tilde{y}}W-3a\tilde{x}\partial_{\tilde{z}}W,$
(3.24) $\displaystyle\partial_{\tilde{z}}m_{1}={}$ $\displaystyle
s^{\tilde{y}}_{1}B^{\tilde{x}}-s^{\tilde{x}}_{1}B^{\tilde{y}}+3a\tilde{y}\partial_{\tilde{x}}W-3a\tilde{x}\partial_{\tilde{y}}W-\frac{2c(2a+1)}{a+2}\partial_{\tilde{z}}W.$

The very long explicit solution for $m_{1}(\tilde{x},\tilde{y},\tilde{z})$ can
be found in Appendix A.

There is no additional first order integral commuting with $X_{1}$ or $X_{2}$.
A more thorough search for additional integrals, i.e. superintegrability, is
out of scope of this article.

###  3.2 Parameters $b$ and $c$ vanishing

We note that under these assumptions the shift (3.1) is just the identity
transformation, hence we drop the tildes in this subsection.

To recapitulate, we have already determined that the integral $X_{2}$ from
(2.11) reduces to a first order integral $\tilde{X}_{2}=p_{z}$ in a suitable
gauge.

We proceed by solving (3.3), the second order equations for $X_{1}$.
Proceeding as described under (3.14), we obtain conditions on $S^{z}_{11}(r)$
and $S^{z}_{2}(r)$, which now reduce to only two equations:

$r^{2}\partial_{r}^{3}S^{z}_{11}(r)+r\partial_{r}^{2}S^{z}_{11}(r)-\partial_{r}S^{z}_{11}(r)=0$
(3.25)

and

$r^{3}\partial_{r}^{4}S^{z}_{2}(r)+5r^{2}\partial_{r}^{3}S^{z}_{2}(r)+2r\partial_{r}^{2}S^{z}_{2}(r)-2\partial_{r}S^{z}_{2}(r)=0.$
(3.26)

The solution has more arbitrary constants than before in (3.17), namely

$S^{z}_{11}(r)=s^{z}_{11}+s^{z}_{12}\ln(r)+s^{z}_{13}r^{2},\quad
S^{z}_{2}(r)=s^{z}_{21}r^{2}+s^{z}_{22}\ln(r)+s^{z}_{23}+\frac{s^{z}_{24}}{r}.$
(3.27)

We solve the equations for $B^{z}_{1}$ and subsequently for $S^{x}_{1}$ and
$S^{y}_{1}$ in Cartesian coordinates analogously to the previous case, this
time the solution includes $\ln(x^{2}+y^{2})$ and $\arctan(\frac{y}{x}).$ Thus
the second order equations are solved and the preliminary form of the magnetic
field reads

$\begin{split}B^{x}(\vec{x})={}&-3a\left(a+2\right)\left(s^{z}_{21}+\frac{s^{z}_{22}}{2r^{2}}-\frac{s^{z}_{24}}{2r^{3}}\right)y,\\\
B^{y}(\vec{x})={}&3a\left(a+2\right)\left(s^{z}_{21}+\frac{s^{z}_{22}}{2r^{2}}-\frac{s^{z}_{24}}{2r^{3}}\right)x,\\\
B^{z}(\vec{x})={}&\left(a+2\right)\left(\omega\left(6s^{z}_{21}+\frac{s^{z}_{22}}{r^{2}}\right)y+3a\left(a+2\right)\left(2s^{z}_{13}+\frac{s^{z}_{12}}{r^{2}}\right)\right).\end{split}$
(3.28)

We now use the first-order terms in momenta from $\\{X_{1},X_{2}\\}=0$ again
in the polar coordinates (3.14). This determines the dependence of
$W_{1}(r,\phi)$ on $\phi$, leaving an undetermined function $W_{11}(r)$.

We proceed to the zeroth order equation from $\\{X_{1},H\\}=0$ and separate
again into coefficients of $\sin(\phi)^{i}\cos(\phi)^{j}$ with different
powers $i,j$ ($i=0,1$). One of the coefficients is

$\omega^{2}\left(2rs^{z}_{22}-3s^{z}_{24}\right)\left(2s^{z}_{21}r^{2}+s^{z}_{22}\right)^{2}=0,$
(3.29)

which (by collecting powers of $r$) implies $s^{z}_{22}=0$. Two other
coefficients yield

$s^{z}_{24}(s^{z}_{21})^{2}=0,\quad\omega(s^{z}_{12})^{2}s^{z}_{24}=0,$ (3.30)

implying the split: $s^{z}_{24}=0$ or $s^{z}_{24}\neq 0$ with
$s^{z}_{12}=s^{z}_{21}=0$.

When $s^{z}_{24}=0$, the compatibility condition for the mixed derivatives of
$m_{1}(x,y,z)=M_{1}(x,y)$ implies that also $s_{12}^{z}=0$ and we find the
system (3.20) with $b=c=0$. The other case, $s^{z}_{24}\neq 0$, gives a
cylindrically separable system: Case II in [17, eq. (7.12)] with constants
$A_{2}=A_{3}=0$, i.e.

$B^{x}(\vec{x})=-A_{1}\frac{y}{r^{3}},\quad
B^{y}(\vec{x})=A_{1}\frac{x}{r^{3}},\quad B^{z}(\vec{x})=0,\quad
W(\vec{x})=\frac{w_{1}}{r}-\frac{A_{1}^{2}}{2r^{2}}.$ (3.31)

The system (3.31) admits $p_{z}$ and $\ell_{z}$ as integrals and $X_{1}$
reduces to a parabolic integral $\ell_{z}p_{x}$, i.e. it is a superintegrable
system separable in cylindrical as well as parabolic cylindrical coordinates.

##  4 Parameter $a$ vanishing

We note that $a=0$ implies $\omega=1$. Applying the assumptions on the
parameters, the integral $X_{1}$ of (2.10) reads

$X_{1}=-2\ell_{z}^{A}p_{x}^{A}-\frac{c}{2}\left(p_{x}^{A}\right)^{2}+c\left(p_{y}^{A}\right)^{2}-cp_{x}^{A}p_{z}^{A}+2\beta
p_{y}^{A}(p_{x}^{A}-p_{z}^{A})-\frac{c}{2}\left(p_{z}^{A}\right)^{2}+s_{1}^{x}p_{x}+s_{1}^{y}p_{y}+s_{1}^{z}p_{z}+m_{1},$
(4.1)

where we have introduced $\beta=b/\sqrt{2}$, which simplifies many formulas
below. The other integral $X_{2}$ retains the form (2.11).

As was observed in [20], the choice $a=0,\omega=1,c=b=0$ leads to the standard
parabolic cylindrical case. We are mainly interested in the integrals of the
generalized form; we will therefore assume that at least one of the parameters
$\beta,c$ is nonvanishing.

The form of $X_{2}$ implies that we can use the solution to the second order
equations (3.4) and $\nabla\cdot\vec{B}=0$ (3.5) from the $a\neq 0$ case,
namely (3.6), with $\tilde{x}=x,\tilde{y}=y,\tilde{z}=z$.

Looking at three of second order involutivity equations $\\{X_{1},X_{2}\\}=0$,
namely

$\begin{split}c\partial_{z}S^{x}_{21}\\!\left(z\right)x-2\left(\beta\partial_{z}S^{x}_{21}\\!\left(z\right)+S^{x}_{21}\\!\left(z\right)\right)y-2\beta\partial_{z}S^{x}_{22}\\!\left(z\right)+2S^{x}_{22}\\!\left(z\right)-c\partial_{z}S^{y}_{2}\\!\left(z\right)+3cS^{x}_{21}\\!\left(z\right)=0,\\\
(\beta\partial_{z}S^{x}_{21}\\!\left(z\right)+S^{x}_{21}\\!\left(z\right))x-\beta\partial_{z}S^{y}_{2}\\!\left(z\right)-\beta
S^{x}_{21}\\!\left(z\right)=0,\\\
c\partial_{z}S^{x}_{21}\\!\left(z\right)y-2\beta
S^{x}_{21}\\!\left(z\right)+c\partial_{z}S^{x}_{22}\\!\left(z\right)+2S^{y}_{2}\\!\left(z\right)=0,\end{split}$
(4.2)

it is clear that we have to split, this time into 3 subcases: both $\beta\neq
0$ and $c\neq 0$, $\beta\neq 0$ with $c=0,$ and $c\neq 0$ with $\beta=0$. The
first and second cases yield a new system, so we will consider them in detail.
The last case, $\beta=0,c\neq 0$, only contains a reduced version of the
system (4.6)–(4.7) from $\beta\neq 0$ with $c\neq 0$ together with known
superintegrable systems (the constant magnetic field with $W=0$, the Cartesian
superintegrable system [21, Case A.2] in the form of eq. (50) therein with
$\beta=-\alpha$ and $\Omega_{2}=0$, and the helical undulator from [16]).

###  4.1 Both $\beta$ and $c$ nonvanishing

This assumption implies that all functions in (4.2) vanish,

$S^{x}_{21}(z)=0,\quad S^{x}_{22}(z)=0,\quad S^{y}_{2}(z)=0$ (4.3)

and the remaining second order equations in $\\{X_{1},X_{2}\\}=0$ imply that
also the functions $s_{1}^{i}$ do not depend on $z$,
$s_{1}^{x_{i}}=s_{1}^{x_{i}}(x,y)$.

As a consequence of the first order terms in the involutivity equations
$\\{X_{1},X_{2}\\}=0,$ the electrostatic potential separates
$W=W_{1}(x,y)+W_{3}(z)$, but the zeroth order term in $\\{X_{1},X_{2}\\}=0$
eliminates $W_{3}(z)$. Therefore, the coordinate $z$ is cyclic and the
integral $X_{2}$ reduces to a first order one,

$\tilde{X}_{2}=p_{z}^{A}+\tfrac{1}{4}\left(S^{z}_{2}(x,y)\right)^{2}.$ (4.4)

We continue with equations for the integral $X_{1}$, starting with the
quadratic ones, i.e. the analogues of (3.3) but without the shift (3.1). The
key equation at this moment is

$2\beta\partial_{y}S^{z}_{2}+c\partial_{x}S^{z}_{2}=0,$ (4.5)

which implies $S^{z}_{2}=S^{z}_{2}(y-2\frac{\beta x}{c})$. Next, we integrate
equations with $S^{x}_{1}$ and $S^{z}_{1}$ differentiated by $y$, then find
$B^{z}_{1}$, introducing arbitrary functions $S^{x}_{11}(x),S^{z}_{11}(x)$ and
$B^{z}_{11}(y-2\frac{\beta x}{c})$ in the integration.

The remaining leading order equations from $\\{H,X_{1}\\}=0$ can be solved for
derivatives of $S^{x}_{1}(x,y)$. Imposing their compatibility
$\partial_{xy}S^{x}_{1}=\partial_{yx}S^{x}_{1}$ and differentiating with
respect to $y$ to eliminate functions of $x$, a coordinate change enables
reduction of $S^{z}_{2}(y-2\frac{\beta x}{c})$ and $B^{z}_{11}(y-2\frac{\beta
x}{c})$ to polynomials in their variable. The solution of the remaining second
order equations is tedious but straightforward.

The first order involutivity conditions $\\{X_{1},X_{2}\\}=0$ introduce a
function of $(y-2\frac{\beta x}{c})$ into $W$, but the zeroth order equation
from $\\{X_{1},H\\}=0$ reduces it to a polynomial or eliminates the magnetic
field (we do not consider the latter further).

The remaining calculations are again straightforward and we get (after
elimination of constants corresponding to the first order integral
$\tilde{X}_{2}$) a magnetic field $\vec{B}$ depending on 2 constants
$b_{y},b_{z}$

$B^{x}(\vec{x})=b_{y}c,\quad B^{y}(\vec{x})=2b_{y}\beta,\quad
B^{z}(\vec{x})=6b_{y}y+b_{z},$ (4.6)

and the electrostatic potential $W$ depending on 3 additional parameters
$w_{1},w_{2},w_{3},$

$\begin{split}W(\vec{x})={}&\left(-\frac{x^{4}}{2}+2\beta
x^{3}-\frac{\left(6y^{2}+4\beta^{2}-c^{2}\right)x^{2}}{2}+6\left(y+\frac{c}{3}\right)\beta
xy-\frac{5y^{4}}{2}+3cy^{3}\right)b_{y}^{2}\\\
{}&-\left(\left(\frac{\left(2y+3c\right)x^{2}}{2}-2\beta
xy+y^{3}\right){b_{z}}+\frac{{w_{1}}}{2}\left(x^{2}+y^{2}\right)+\beta{w_{2}}x+{w_{3}}y\right){b_{y}}.\end{split}$
(4.7)

The integral $X_{2}$, retaining the form (2.11), reduces to the first order
one,

$\tilde{X}_{2}=p_{z}^{A}+b_{y}(2\beta x-cy),$ (4.8)

while $X_{1}$ from (4.1) remains quadratic. Its lower order terms are defined
by the following functions:

$\displaystyle s_{1}^{x}(\vec{x})={}$
$\displaystyle\left(6b_{y}y+b_{z}\right)x^{2}-2\beta\left(6b_{y}y+b_{y}c+b_{z}\right)x+10b_{y}y^{3}-3\left(3b_{y}c-b_{z}\right)y^{2}+{w_{1}}y+{w_{3}},$
$\displaystyle s_{1}^{y}(\vec{x})={}$
$\displaystyle-2{b_{y}}x^{3}+6{b_{y}}\beta
x^{2}-\left(6{b_{y}}y^{2}+2{b_{z}}y+\left(4\beta^{2}-c^{2}\right){b_{y}}+3c{b_{z}}+w_{1}\right)x$
(4.9)
$\displaystyle+6\beta{b_{y}}y^{2}+2\beta\left({b_{y}}c+{b_{z}}\right)y-\beta{w_{2}},$
$\displaystyle s_{1}^{z}(\vec{x})={}$ $\displaystyle-
b_{y}cx^{2}+2\beta\left(2b_{y}y+b_{y}c+b_{z}\right)x-3b_{y}cy^{2}-\left(4\beta^{2}b_{y}-3b_{y}c^{2}+cb_{z}\right)y+\tfrac{1}{2}c{w_{2}},$

$\displaystyle m_{1}(\vec{x})={}$
$\displaystyle\left(4b_{y}^{2}y+\tfrac{1}{2}cb_{y}^{2}+b_{z}b_{y}\right)x^{4}-4\left[4b_{y}^{2}\beta
y+b_{y}\beta\left(b_{y}c+b_{z}\right)\right]x^{3}$
$\displaystyle+\left[16b_{y}^{2}y^{3}-\left(6cb_{y}^{2}-8b_{z}b_{y}\right)y^{2}+\left(20b_{y}^{2}\beta^{2}-c^{2}b_{y}^{2}+4b_{z}b_{y}c+2b_{y}w_{1}+b_{z}^{2}\right)y\right.$
$\displaystyle+\left.\frac{\left(12\beta^{2}c-c^{3}\right)b_{y}^{2}}{2}+\frac{\left(12\beta^{2}b_{z}+2b_{z}c^{2}+cw_{1}+2w_{3}\right)b_{y}}{2}+\frac{3cb_{z}^{2}}{2}+\frac{b_{z}w_{1}}{2}\right]x^{2}$
$\displaystyle-\left[32b_{y}^{2}\beta
y^{3}-2\left(b_{y}\left(b_{y}c-8b_{z}\right)\beta
y+2\beta\left(c^{2}-2\beta^{2}\right)b_{y}^{2}-\left(\left(3cb_{z}+w_{1}\right)\beta-\beta
w_{2}\right)b_{y}-\beta b_{z}^{2}\right)y\right.$
$\displaystyle-\left.\left(2b_{y}c+b_{z}\right)\beta
w_{2}+2b_{y}w_{3}\beta\right]x+12b_{y}^{2}y^{5}+\left(7b_{z}b_{y}-\tfrac{37}{2}cb_{y}^{2}\right)y^{4}$
(4.10)
$\displaystyle+\left(4b_{y}^{2}\beta^{2}+7c^{2}b_{y}^{2}-5b_{z}b_{y}c+2b_{y}w_{1}+b_{z}^{2}\right)y^{3}$
$\displaystyle+\left(\frac{\left(8\beta^{2}c-3c^{3}\right)b_{y}^{2}}{2}+\frac{\left(4\beta^{2}b_{z}+b_{z}c^{2}-2cw_{1}+6w_{3}\right)b_{y}}{2}+\frac{b_{z}w_{1}}{2}\right)y^{2}$
$\displaystyle-\frac{\left(4w_{2}b_{y}\beta^{2}+\left(4b_{y}c-2b_{z}\right)w_{3}+b_{y}c^{2}w_{2}\right)y}{2}.$

The quadratic terms in $X_{1}$ from (4.1) can be simplified by a shift of the
coordinates to

$X_{1}=\ell_{z}^{A}p_{x}^{A}-2\beta
p_{y}^{A}p_{z}^{A}-2cp_{x}^{A}p_{z}^{A}+\ldots$ (4.11)

Subsequently reducing the system to the $xy$-plane by setting
$X_{2}=\text{const.},$ the last two terms of $X_{1}$ become first order ones,
which leads to a parabolic structure of the integral [33]. The system can be
therefore interpreted as an extension of the 2D parabolic case [29]. On the
other hand, the motion in the plane is not completely decoupled unlike e.g.
the systems in [8]. (See [5] for the related notion of partial integrability
applied to various systems including one with a magnetic field.)

###  4.2 Parameter $\beta$ nonvanishing, $c$ vanishing

With this assumption on the parameters, the solution to (4.2) is no longer
trivial,

$S^{x}_{21}(z)=s^{x}_{21}{\mathrm{e}}^{-\frac{z}{\beta}},\quad
S^{x}_{22}(z)=s^{x}_{22}{\mathrm{e}}^{\frac{z}{\beta}},\quad
S^{y}_{2}(z)=\beta s^{x}_{21}{\mathrm{e}}^{-\frac{z}{\beta}},$ (4.12)

and the functions $s_{1}^{x_{i}}$ get some terms depending on $z$.

We proceed with the quadratic equations for $X_{1}$. Eq. (4.5) with $c=0$
implies $S^{z}_{2}=S^{z}_{2}(x)$, so we avoid the technical difficulties
encountered in previous subsection. The other quadratic equations can be
solved step by step. To identify splitting of the computation, we need to fix
the constants in the solution for $S^{z}_{2}(x)$,

$S^{z}_{2}(x)=s^{z}_{20}+2s^{z}_{21}(\beta-x)+\frac{s^{z}_{22}}{16(\beta-x)^{4}}+\frac{s^{z}_{23}}{16(\beta-x)^{2}}.$
(4.13)

When $s^{x}_{21}=s^{x}_{22}=0$, we do not get anything new, only a shifted
version of a known superintegrable system [22, Cases I(b) and I(c)]
corresponding to $s^{z}_{21}=0$,

$\begin{split}B^{x}(\vec{x})={}&0,\quad
B^{y}(\vec{x})=\frac{b_{1}}{\left(\beta-x\right)^{3}},\quad
B^{z}(\vec{x})=0,\quad
W(\vec{x})=-\frac{b_{1}^{2}}{8\left(\beta-x\right)^{4}}+\frac{w_{1}}{8\left(\beta-x\right)^{2}},\end{split}$
(4.14)

and $s^{z}_{22}=s^{z}_{23}=0$ leads to the special case of (4.6)–(4.7) with
$c=0$.

When at least one of the constants $s^{x}_{21}$ and $s^{x}_{22}$ is
nonvanishing, which implies that the integral $X_{2}$ does not reduce to a
first order one, we also have to go through several subcases depending on
constants in $S^{z}_{2}(z).$ However, only the branch with $s^{x}_{21}\neq 0$
and $s^{x}_{22}\neq 0$ is relevant in the end, all other give reduced version
of this branch, i.e. its limits with either $s^{x}_{21}\to 0$ or
$s^{x}_{22}\to 0$. The final solution for the fields has 4 constants in the
magnetic field $\vec{B}$,
$b_{1}=-\frac{s^{x}_{21}}{2\beta},b_{2}=-\frac{s^{x}_{22}}{2\beta},b_{y}=-\frac{s^{z}_{21}}{2\beta}$,
and $b_{z}$ which is the constant term in $B^{z}_{1}(x,y)$ (the other
constants in (4.13) vanish),

$B^{x}(\vec{x})=b_{1}{\mathrm{e}}^{-\frac{z}{\beta}}\left(x-\beta\right),\quad
B^{y}(\vec{x})=b_{1}{\mathrm{e}}^{-\frac{z}{\beta}}y-b_{2}{\mathrm{e}}^{\frac{z}{\beta}}+2\beta
b_{y},\quad B^{z}(\vec{x})=2\beta
b_{1}{\mathrm{e}}^{-\frac{z}{\beta}}+6b_{y}y+b_{z},$ (4.15)

and the electrostatic potential $W$ has 2 more constants, $w_{1},w_{3}$ (the
choice of the notation will be clear soon),

$\displaystyle W(\vec{x})={}$
$\displaystyle-\tfrac{1}{2}\beta^{2}b_{1}^{2}\left(x^{2}-2\beta
x+y^{2}\right){\mathrm{e}}^{-\frac{2z}{\beta}}-\beta b_{1}\left[(x^{2}-2\beta
x+y^{2})(2b_{y}y+\tfrac{1}{2}b_{z})+\tfrac{1}{2}w_{3}\right]{\mathrm{e}}^{-\frac{z}{\beta}}$
$\displaystyle-\beta b_{2}\left[(x^{2}-2\beta
x+3y^{2})b_{y}+b_{z}y+\tfrac{1}{2}w_{1}\right]{\mathrm{e}}^{\frac{z}{\beta}}-\tfrac{1}{2}\beta^{2}b_{2}^{2}{\mathrm{e}}^{\frac{2z}{\beta}}-\beta^{2}b_{1}b_{2}y$
$\displaystyle+\left[-\tfrac{1}{2}x^{4}+2\beta
x^{3}-\left(3y^{2}+2\beta^{2}\right)x^{2}+6\beta
xy^{2}-\tfrac{5}{2}y^{4}\right]b_{y}^{2}$ (4.16)
$\displaystyle-\left[\left(x^{2}-2\beta
x+y^{2}\right)b_{z}y+\frac{w_{1}}{2}\left(x^{2}+y^{2}\right)-\beta{w_{1}}x+{w_{3}}y\right]{b_{y}}.$

If we eliminate the exponentials in (4.15) by setting $b_{1}=b_{2}=0$, the
electrostatic potential $W$ goes to (4.7) where we set $c=0$, $w_{2}=-w_{1}$
and $b_{1}=b_{y}$.

This time, both integrals $X_{1}$ of the form (4.1) and $X_{2}$ as in (2.11)
are quadratic, the functions determining their lower order terms read

$\displaystyle s_{1}^{x}(\vec{x})={}$ $\displaystyle 4\beta
b_{2}y{\mathrm{e}}^{\frac{z}{\beta}}+2\beta b_{1}\left(x^{2}-2\beta
x+2y^{2}\right){\mathrm{e}}^{-\frac{z}{\beta}}+2\left(3x^{2}-6\beta
x+5y^{2}\right)b_{y}y$ $\displaystyle+\left(x^{2}-2\beta
x+3y^{2}\right)b_{z}+w_{1}y+w_{3},$ (4.17) $\displaystyle
s_{1}^{y}(\vec{x})={}$ $\displaystyle-\left[2\beta
b_{2}{\mathrm{e}}^{\frac{z}{\beta}}+2\beta
b_{1}y{\mathrm{e}}^{-\frac{z}{\beta}}+2b_{y}(x^{2}-2\beta
x+3y^{2})+2b_{z}y+w_{1}\right]\left(x-\beta\right),$ $\displaystyle
s^{z}_{1}(\vec{x})={}$ $\displaystyle
2\beta\left[b_{1}\beta{\mathrm{e}}^{-\frac{z}{\beta}}+2b_{y}y+b_{z}\right]\left(x-\beta\right),$

$\displaystyle m_{1}(\vec{x})={}$ $\displaystyle
2\beta^{2}b_{2}^{2}y{\mathrm{e}}^{\frac{2z}{\beta}}+b_{2}\beta\left[2\left(3x^{2}-6\beta
x+5y^{2}+2\beta^{2}\right)b_{y}y+b_{z}\left(x^{2}-2\beta
x+3y^{2}+2\beta^{2}\right)+w_{1}y+w_{3}\right]{\mathrm{e}}^{\frac{z}{\beta}}$
$\displaystyle+\beta b_{1}\left[2b_{y}(x^{2}-2\beta x+y^{2})(x^{2}-2\beta
x+5y^{2}+\beta^{2})+(x^{2}-2\beta x+y^{2})3b_{z}y+\right.$
$\displaystyle+\left.w_{1}\left((x-\beta)^{2}+y^{2}\right)+w_{3}y\right]{\mathrm{e}}^{-\frac{z}{\beta}}+2y\beta^{2}b_{1}^{2}\left(x^{2}+y^{2}-2\beta
x\right){\mathrm{e}}^{-\frac{2z}{\beta}}$ (4.18)
$\displaystyle+\left(x^{2}-2\beta
x+y^{2}\right)\left[4b_{y}^{2}y\left((x-\beta)^{2}+3y^{2}\right)+\left(x^{2}-2\beta
x+7y^{2}+2\beta^{2}\right)b_{y}b_{z}\right.$
$\displaystyle+\left.b_{z}^{2}y+\tfrac{1}{2}b_{z}w_{1}\right]+2\beta^{2}b_{1}b_{2}\left(x^{2}-2\beta
x+2y^{2}\right)$
$\displaystyle+2b_{y}w_{1}\left((x-\beta)^{2}+y^{2}\right)y+[(x^{2}-2\beta
x+3y^{2})b_{y}+b_{z}y]w_{3},$

$\displaystyle s_{2}^{x}(\vec{x})={}$ $\displaystyle-2\beta
b_{1}{\mathrm{e}}^{-\frac{z}{\beta}}y-2\beta
b_{2}{\mathrm{e}}^{\frac{z}{\beta}},\quad s_{2}^{y}(\vec{x})=2\beta
b_{1}{\mathrm{e}}^{-\frac{z}{\beta}}\left(x-\beta\right),\quad
s_{2}^{z}(\vec{x})=4\beta b_{y}\left(x-\beta\right),$ (4.19)

$\begin{split}m_{2}(\vec{x})={}&-2\beta^{2}b_{2}^{2}{\mathrm{e}}^{\frac{2z}{\beta}}-2\beta
b_{2}\left[\left(x^{2}-2\beta
x+3y^{2}-2\beta^{2}\right)b_{y}+b_{z}y+\tfrac{1}{2}w_{1}\right]{\mathrm{e}}^{\frac{z}{\beta}}\\\
{}&-4\beta
b_{1}\left[\left((x-\beta)^{2}+y^{2}\right)b_{y}y+\tfrac{1}{4}\left(x^{2}-2\beta
x+y^{2}\right)b_{z}+\tfrac{1}{4}w_{3}\right]{\mathrm{e}}^{-\frac{z}{\beta}}\\\
{}&-\beta^{2}b_{1}^{2}\left(\beta^{2}-4\beta
x+2x^{2}+2y^{2}\right){\mathrm{e}}^{-\frac{2z}{\beta}}+{4\beta^{2}b_{y}^{2}\left(x^{2}-2\beta
x\right)}-4\beta^{2}b_{1}b_{2}y.\end{split}$ (4.20)

##  5 Conclusions

In this article on quadratically integrable systems with magnetic field
possessing integrals of motion not connected to separability, we have
classified such systems in the generalized parabolic cylindrical class [20,
class (j)].

Depending on the values of the parameters $a,b,c$, there are 3 families of
such systems not known before: When $a\neq 0$ and the integral of motion
$X_{1}$ (3.2) has 2 terms depending on the angular momentum, we have found the
system (3.20)–(3.21) with linear magnetic field $\vec{B}$ and quartic scalar
potential $W$. When $a=0$, we have 2 families: The first one allowing both
remaining parameters $b$ and $c$ nonvanishing, (4.6)–(4.7), also has a rather
simple linear $\vec{B}$ and quartic potential, while the second one with
$b\neq 0$ and $c=0$, (4.15)–(4.16), has the fields depending on exponential
functions of the coordinate $z$.

Let us consider the limit of vanishing magnetic field. All three systems
reduce to the free motion, thus they are obviously separable in every choice
of orthogonal coordinates in the Euclidean space $\mathbb{E}_{3}$. This
behavior is consistent with previous results on nonstandard integrable
systems, which in the limit $\vec{B}\to 0$ always separate in Cartesian or
cylindrical coordinates, i.e. in subgroup-type coordinate systems, whether
they extend this type of coordinates or not [16, 12, 20]. On the other hand,
the limit for standard systems typically does not eliminate the potential $W$
and leads to separation only in the corresponding coordinates, see e.g. the
spherical case [23] (which does not admit extensions [16]). We therefore
confirm the not well understood disparity between subgroup- and non-subgroup-
type coordinates, which shall be investigated further.

Given the computational complexity in classifying the parabolic case in 2D
[29], the 3D cylindrical parabolic systems, here corresponding to the limit
$a,b,c\to 0$, have not yet been classified. We therefore investigate this
limit of the systems found above expecting to find yet unknown integrable
systems of that type.

In the case $a\neq 0$, (3.20)–(3.21), we have to take the limit $b,c\to 0$
first, redefine $b_{x}=a\tilde{b}_{x},w_{i}=\frac{\tilde{w}_{i}}{a}$ and only
then $a\to 0$. This leads to the following fields

$\begin{split}B^{x}(\vec{x})={}&0,\quad B^{y}(\vec{x})=0,\quad
B^{z}(\vec{x})=2\tilde{b}_{x}y+b_{z},\\\
W(\vec{x})={}&\left[-\tfrac{1}{18}\tilde{b}_{x}\left(x^{4}+6x^{2}y^{2}+5y^{4}\right)-\tfrac{1}{3}b_{z}\left(x^{2}+y^{2}\right)y+\tilde{w}_{3}\left(x^{2}+y^{2}\right)+\tilde{w}_{1}x+\tilde{w}_{2}y\right]\tilde{b}_{x},\end{split}$
(5.1)

yielding a parabolic cylindrical system with magnetic fields that can be seen
as a trivial extension of a parabolic 2D system. (In this limit, the shift
(3.1) is vacuous.) As far as we know this is a new system.

The limit $b,c\to 0$ (thus $\beta\to 0$) for the system (4.6)–(4.7) does not
entail such technical difficulties. We get the system (5.1) but with
$\tilde{w}_{1}=0,\tilde{w}_{2}=\tfrac{1}{3}w_{3},\tilde{w}_{3}=\tfrac{1}{6}w_{1}$
and $\tilde{b}_{x}\to 3b_{y}$. On the other hand, the limit for system
(4.15)–(4.16) with $b_{1}\neq 0$ or $b_{2}\neq 0$ diverges due to the
exponentials. When the exponentials of (4.15)–(4.16) vanish, i.e.
$b_{1}=b_{2}=0$, the system differs from (4.6)–(4.7) only by the absence of
the terms involving $c$, hence the limit $b,c\to 0$ again yields the system
(5.1) with the redefinitions above.

## Acknowledgment

Computations in this paper were performed using computer algebra software
MapleTM by Maplesoft, a division of Waterloo Maple, Inc., Waterloo, Ontario,
and MathematicaTM by Wolfram Research, Inc., Champaign, IL.

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

##  Appendix A Appendix

In this brief appendix we present the explicit solution of equations (3.24).

$\displaystyle m_{1}(\vec{x})={}$
$\displaystyle\left(\frac{b_{x}^{2}\left(a+2\right)\omega\tilde{y}}{9a^{2}}-\frac{c\left(4a^{2}+a+4\right)b_{x}^{2}}{108a^{3}}+\frac{b_{z}\left(a+2\right)b_{x}}{12a}\right)\left(a+2\right)\tilde{x}^{4}$
$\displaystyle+\left(-\frac{\left(23a^{2}-16a-16\right)b\psi\omega
b_{x}^{2}\tilde{y}}{27a^{3}\left(a-1\right)}+\frac{4\psi\omega^{2}bc\left(5a+4\right)b_{x}^{2}}{81a^{4}\left(a+2\right)}-\frac{8b\psi\left(a+\frac{1}{2}\right)b_{z}b_{x}}{9a^{2}}\right)\tilde{x}^{3}$
$\displaystyle+\left[-\frac{2\left(7a^{2}-8a-8\right)\omega
b_{x}^{2}\tilde{y}^{3}}{9a^{2}}+\left(\frac{c\left(2a+1\right)\left(a-2\right)b_{x}^{2}}{3a^{2}}-\frac{b_{z}\left(23a^{2}-16a-16\right)b_{x}}{6a}\right)\tilde{y}^{2}\right.$
$\displaystyle+\left.\left(-4\omega
b_{x}w_{3}+\frac{8\left(a+2\right)\left(a+\frac{1}{2}\right)b^{2}\omega
b_{x}^{2}}{27a^{4}}+\frac{2c\left(13a^{2}+4a-8\right)\omega
b_{x}b_{z}}{9a^{2}\left(a+2\right)}+b_{z}^{2}\omega\right)\tilde{y}\right.$
$\displaystyle+\left.\left(-\frac{2\left(5a^{2}-a-4\right)cb_{x}}{3a\left(a+2\right)}-3ab_{z}\right)w_{3}-b_{x}\omega
w_{2}\right.$
$\displaystyle+\left.\frac{2\left(2a^{4}b^{2}+40a^{4}c^{2}+13a^{3}b^{2}+17a^{3}c^{2}+30a^{2}b^{2}-57a^{2}c^{2}+28ab^{2}-16ac^{2}+8b^{2}+16c^{2}\right)b_{x}b_{z}}{27a^{3}\left(a+2\right)^{2}}\right.$
$\displaystyle+\left.\frac{b_{z}^{2}\left(8a^{2}+5a-4\right)c}{3a\left(a+2\right)}\right]\tilde{x}^{2}$
$\displaystyle+\left[{-\frac{\psi\left(55a^{2}-32a-32\right)b\omega
b_{x}^{2}\tilde{y}^{3}}{27a^{3}\left(a-1\right)}}-\frac{4\psi\left(2a+1\right)bb_{x}\left(4ab_{z}-cb_{x}\right)\tilde{y}^{2}}{9a^{3}}\right.$
$\displaystyle+\left.\left(-\frac{4b\psi\omega
b_{x}w_{3}}{3a\left(a-1\right)}-2b_{x}\omega
w_{1}+\frac{20\left(a+\frac{4}{5}\right)b\psi c\omega
b_{x}b_{z}}{27a^{3}\left(a+2\right)}+\frac{2b\psi
b_{z}^{2}\omega}{3a\left(a-1\right)}\right)\tilde{y}\right.$
$\displaystyle-\left.\frac{2bb_{x}\psi\omega
w_{2}}{3a\left(a-1\right)}+\left(\frac{2cb_{x}\left(-5a^{2}+a+4\right)}{3\left(a+2\right)a}-3ab_{z}\right)w_{1}\right]\tilde{x}-\frac{\left(5a^{2}-4a-4\right)\omega
b_{x}^{2}\tilde{y}^{5}}{3a^{2}}$
$\displaystyle+\left(\frac{c\left(196a^{4}+39a^{3}-312a^{2}-52a+48\right)b_{x}^{2}}{36a^{3}\left(a+2\right)}-\frac{b_{z}\left(47a^{2}-28a-28\right)b_{x}}{12a}\right)\tilde{y}^{4}$
$\displaystyle+\left(-4b_{x}\omega
w_{3}+\frac{4\left(2a^{3}b^{2}+24a^{3}c^{2}+9a^{2}b^{2}+15a^{2}c^{2}+12ab^{2}-12ac^{2}+4b^{2}\right)\omega
b_{x}^{2}}{81\left(a+2\right)a^{4}}\right.$
$\displaystyle-\left.\frac{22\left(a^{2}+a-\frac{4}{11}\right)c\omega
b_{x}b_{z}}{9a^{2}\left(a+2\right)}+b_{z}^{2}\omega\right)\tilde{y}^{3}$
$\displaystyle+\left((2cb_{x}-3ab_{z})w_{3}-3b_{x}\omega
w_{2}+\frac{4\left(a+2\right)\left(a+\frac{1}{2}\right)b^{2}b_{x}b_{z}}{27a^{3}}\right)\tilde{y}^{2}$
$\displaystyle-\left(\left(3ab_{z}-2cb_{x}\right)w_{2}+\frac{2bb_{x}\psi\omega
w_{1}}{3a\left(a-1\right)}\right)\tilde{y},$

where we recall that $\omega=\sqrt{1+a-2a^{2}}$ and $\psi=\sqrt{-a^{2}-a+2}$.

## References

  * [1] N. Balal, V. L. Bratman, and E. Magory.  New varieties of helical undulators.  Nucl. Instrum. Methods Phys. Res., 971:163895, 2020.
  * [2] S. Benenti, C. Chanu, and G. Rastelli.  Variable separation for natural Hamiltonians with scalar and vector potentials on Riemannian manifolds.  J. Phys. A, 42(5):2065–2091, 2001.
  * [3] J. Bérubé and P. Winternitz.  Integrable and superintegrable quantum systems in a magnetic field.  J. Math. Phys., 45(5):1959–1973, 2004.
  * [4] R. Campoamor-Stursberg, J. F. Cariñena, and M. F. Rañada.  Higher-order superintegrability of a Holt related potential.  J. Phys. A: Math. Theor., 46(43):435202, 2013.
  * [5] A. M. Escobar-Ruiz and R. Azuaje.  On particular integrability in classical mechanics.  J. Phys. A: Math. Theor., 57(10):105202, 2024.
  * [6] A. M. Escobar-Ruiz, J. C. López Vieyra, and P. Winternitz.  Fourth order superintegrable systems separating in polar coordinates. I. Exotic potentials.  J. Phys. A: Math. Theor., 50(49):495206, 2017.
  * [7] A. M. Escobar-Ruiz, P. Winternitz, and İ. Yurduşen.  General Nth-order superintegrable systems separating in polar coordinates.  J. Phys. A: Math. Theor., 51(40):40LT01, 2018.
  * [8] F. Fournier, L. Šnobl, and P. Winternitz.  Cylindrical type integrable classical systems in a magnetic field.  J. Phys. A, 53(8):085203, 2020.
  * [9] J. Friš, V. Mandrosov, Ya. A. Smorodinsky, M. Uhlíř, and P. Winternitz.  On higher symmetries in quantum mechanics.  Phys. Lett., 16(3):354–356, 1965.
  * [10] T. Heinzl and A. Ilderton.  Superintegrable relativistic systems in spacetime-dependent background fields.  J. Phys. A, 50(34):345204, 14, 2017.
  * [11] F. Hoque and L. Šnobl.  Family of nonstandard integrable and superintegrable classical Hamiltonian systems in non-vanishing magnetic fields.  J. Phys. A, 56:165203, 2023.
  * [12] M. F. Hoque, O. Kubů, A. Marchesiello, and L. Šnobl.  New classes of quadratically integrable systems with velocity dependent potentials: non-subgroup type cases.  Eur. Phys. J. Plus, 138(9):845, 2023.
  * [13] M. F. Hoque, I. Marquette, and Y.-Z. Zhang.  On superintegrable monopole systems.  J. Phys. Conf. Ser., 965:012018, 2018.
  * [14] E. G. Kalnins, W. Miller, Jr., and P. Winternitz.  The group ${\rm O}(4)$, separation of variables and the hydrogen atom.  SIAM J. Appl. Math., 30:630–664, 1976.
  * [15] J. Kress, K. Schöbel, and A. Vollmer.  An algebraic geometric foundation for a classification of second-order superintegrable systems in arbitrary dimension.  J. Geom. Anal., 33(11):360, 2023.
  * [16] O. Kubů, A. Marchesiello, and L. Šnobl.  New classes of quadratically integrable systems in magnetic fields: The generalized cylindrical and spherical cases.  Ann. Phys., 451:169264, 2023.
  * [17] O. Kubů, A. Marchesiello, and L. Šnobl.  Superintegrability of separable systems with magnetic field: the cylindrical case.  J. Phys. A: Math. Theor., 54(42):425204, 2021.
  * [18] A. A. Makarov, Ya. A. Smorodinsky, Kh. Valiev, and P. Winternitz.  A systematic search for nonrelativistic systems with dynamical symmetries.  Nuovo Cimento A Series, 10:1061–1084, 1967.
  * [19] A. Marchesiello, S. Post, and L. Šnobl.  Third-order superintegrable systems with potentials satisfying only nonlinear equations.  J. Math. Phys., 56(10):102104, 2015.
  * [20] A. Marchesiello and L. Šnobl.  Pairs of commuting quadratic elements in the universal enveloping algebra of Euclidean algebra and integrals of motion.  J. Phys. A, 55(14):145203, 2022.
  * [21] A. Marchesiello and L. Šnobl.  Superintegrable 3D systems in a magnetic field corresponding to Cartesian separation of variables.  J. Phys. A, 50(24):245202, 24, 2017.
  * [22] A. Marchesiello and L. Šnobl.  Classical superintegrable systems in a magnetic field that separate in Cartesian coordinates.  SIGMA Symmetry Integrability Geom. Methods Appl., 16, 2020.  015, 35 pages.
  * [23] A. Marchesiello, L. Šnobl, and P. Winternitz.  Spherical type integrable classical systems in a magnetic field.  J. Phys. A, 51(13):135205, 24, 2018.
  * [24] E. McSween and P. Winternitz.  Integrable and superintegrable Hamiltonian systems in magnetic fields.  J. Math. Phys., 41(5):2957–2967, 2000.
  * [25] W. Miller, Jr., E. G. Kalnins, and J. M. Kress.  Separation of Variables and Superintegrability: The Symmetry of Solvable Systems.  Bristol: IOP Publ., 2018.
  * [26] W. Miller, Jr., J. Patera, and P. Winternitz.  Subgroups of Lie groups and separation of variables.  J. Math. Phys., 22(2):251–260, 1981.
  * [27] W. Miller, Jr., S. Post, and P. Winternitz.  Classical and quantum superintegrability with applications.  J. Phys. A, 46(42):423001, 97, 2013.
  * [28] S. Post and P. Winternitz.  A nonseparable quantum superintegrable system in 2D real Euclidean space.  J. Phys. A: Math. Theor., 44(16):162001, 2011.
  * [29] G. Pucacco and K. Rosquist.  Integrable Hamiltonian systems with vector potentials.  J. Math. Phys., 46(1):012701, 2005.
  * [30] D. Reyes, P. Tempesta, and G. Tondo.  Classical multiseparable Hamiltonian systems, superintegrability and Haantjes geometry.  Commun. Nonlinear Sci., 104:106021, 2022.
  * [31] V. N. Shapovalov, V. G. Bagrov, and A. G. Meshkov.  Separation of variables in the stationary Schrödinger equation.  Sov. Phys. J., 15(8):1115–1119, 1972.
  * [32] P. Tempesta and G. Tondo.  Haantjes algebras of classical integrable systems.  Ann. Mat. Pur. Appl., 201(1):57–90, 2021.
  * [33] P. Winternitz, Ya. A. Smorodinsky, M. Uhlíř, and I. Friš.  Symmetry groups in classical and quantum mechanics.  Sov. J. Nucl. Phys., 4:444–450, 1967.
  * [34] P.-M. Zhang, L.-P. Zou, P. A. Horvathy, and G. W. Gibbons.  Separability and dynamical symmetry of quantum dots.  Ann. Phys., 341:94–116, 2014.