# Quantum Integrable Systems arising from Separation of Variables on $S^{3}$

Sean R. Dawson, Holger R. Dullin111Emails<EMAIL_ADDRESS><EMAIL_ADDRESS>
School of Mathematics and Statistics,
The University of Sydney, Australia

###### Abstract

We study the family of quantum integrable systems that arise from separating
the Schrödinger equation in all $6$ separable orthogonal coordinates on the
$3-$sphere: ellipsoidal, prolate, oblate, Lamé, spherical and cylindrical. On
the one hand each separating coordinate system gives rise to a quantum
integrable system on $S^{2}\times S^{2}$, on the other hand it also leads to
families of harmonic polynomials in $\mathbb{R}^{4}$. We show that separation
in ellipsoidal coordinates yields a generalised Lamé equation - a Fuchsian ODE
with 5 regular singular points. We seek polynomial solutions so that the
eigenfunctions are analytic at all finite singularities. We classify
eigenfunctions by their discrete symmetry and compute the joint spectrum for
each symmetry class. The latter 5 separable coordinate systems are all
degenerations of the ellipsoidal coordinates. We perform similar analyses on
these systems and show how the ODEs degenerate in a fashion akin to their
respective coordinates. For the prolate system we show that there exists a
defect in the joint spectrum which prohibits a global assignment of quantum
numbers: the system has quantum monodromy. This is a companion paper to [33]
where the respective classical systems were studied.

##  1 Introduction

The study of separation of variables yields a vast array of classical and
quantum integrable systems that are ripe for exploration. Separation of
variables originated with Jacobi and Stäckel [42]. On the quantum side,
seminal work was done in the early $20$th century by Robertson [36] who
showed, for a given orthogonal coordinate system, that the Schrödinger
equation is separable if both the classical Hamilton-Jacobi equation is
separable and the Robertson condition is satisfied. Work done soon after by
Eisenhart [19] showed that all systems that arise from separation of variables
are Stäckel systems and so the corresponding classical integrals obtained will
be quadratic in the momenta. Work by Kalnins and Miller in the late $20$th
century yielded a complete classification of all orthogonal separable
coordinates on both $\mathbb{R}^{n}$ and $S^{n}$ [26]. This work was further
summarised in [23], where connections to superintegrability were also drawn.
Separation for systems with potential are also well studied. In [25]
invariants for a system that admit separation of variables were constructed.
Schöbel and Veselov gave a topology to the space of separable coordinates on
the sphere [40, 39] and identified it with Stasheff polytopes. Previously, our
motivation was to study the corresponding integrable systems, which inherit
the same topology. In our paper [33], we showed that the space of orthogonally
separable coordinates on $S^{3}$ induces a family of classical integrable
systems on $S^{2}\times S^{2}$ after reduction. Here we seek to broaden this
analysis to the space of quantum integrable systems, their corresponding
separated ODEs and special functions. In particular, we study 5 families of
Fuchsian equations: the generalised Lamé, Heun, Jacobi, Gegenbauer and
Legendre differential equations. Similar to our previous work, we aim to
understand how degenerations in the coordinate systems descend to the space of
integrable systems, and correspondingly ODEs and special functions.

The general idea to start with a super-integrable system and exploit its
multiseparability to define interesting Liouville integrable systems by
reduction has been exploited in the case of the Kepler problem [18] and the
harmonic oscillator [12]. This paper follows the same approach for the
geodesic flow on the three dimensional sphere $S^{3}$. In this case the space
of separable coordinate systems is a Stasheff polytope which is the pentagon
shown in Fig. 1. This is the parameter space of the family of quantum
integrable systems on $S^{2}\times S^{2}$. While the previous examples (Kepler
and Harmonic oscillator) also lead to reduced systems on compact symplectic
manifolds and hence quantised systems on finite dimensional Hilbert spaces,
starting with $S^{3}$ separation of variables in the Schrödinger equation
leads to Fuchsian equations, while the presence of the potential in the
earlier examples lead to confluent Fuchsian equations. Degeneration of
parameters in the general separable coordinate system on $S^{3}$ to the so-
called prolate case (the bottom edge of the pentagon in Fig. 1) then leads to
a semi-toric family whose quantum mechanics is described by polynomial
solutions of the Heun equation. The semi-toric systems on $S^{2}\times S^{2}$
studied in [38, 28, 4] do not appear to include the prolate system studied
here, because the non-trivial integral is a homogeneous polynomial. In fact,
all our Liouville integrable systems have integrals that are in general
homogeneous quadratics in the momenta, and may degenerate into squares of
linear functions of the momenta. Because of this all our integrable systems
are related to spherical harmonics on $S^{3}$. The ellipsoidal reduced system
is related to the Manakov top and has been studied classically and quantum
mechanically in [41].

Figure 1: Stasheff polytope of orthogonal separable coordinates on $S^{3}$
with corresponding ODEs arising from separation.

In our study of these quantum systems, we perform a variety of numerical
computations to find the joint spectrum. For all degenerate systems, a power
series ansatz yields a three term recurrence relation. Given the necessary
quantisation conditions, this recurrence terminates and yields a polynomial as
an eigenfunction. The joint spectrum is then computed from the associated tri-
diagonal matrix. The ellipsoidal coordinates on $S^{3}$ gives the more
difficult generalised Lamé equation which has $5$ regular singularities.
Consequently, the aforementioned methodology falls short and we use the method
of computing Heine-Stieltjes polynomials as described in [3, 2, 45]. It is
interesting to note that for all systems on $S^{3}$, the joint spectrum is
recovered from finite matrices, with corresponding polynomial eigenfunction.
This contrasts with our previous work [15, 14] where the spectrum is
determined by infinite matrices and the eigenfunctions are given by infinite
series.

Another motivation for this work is to suggest a natural family of integrable
systems that contains the semi-toric case, but includes deformation to more
typical quantum integrable systems. The classification of semi-toric systems
is now fully understood classically [35] and quantum mechanically [29]. This
line of research was initiated by Duistermaat [17] and the first example of
quantum monodromy is due to Cushman and Duistermaat [13]. Monodromy, both
classical and quantum, is caused by the existence of a non-degenerate focus-
focus critical value in the classical system [32]. Many important examples of
integrable systems in physics do have monodromy, see, e.g. [13, 38, 14, 18,
15]. In this paper, we obtain a semi-toric system through separation of the
geodesic flow on $S^{3}$ in prolate coordinates. This system thus arises as a
degeneration of a more general 2 parameter family corresponding to the inside
of the pentagon in Fig. 1.

The outline of this paper is as follows. In Section 2 we describe the
ellipsoidal-spherical coordinates, the most general orthogonal separable
coordinates on $S^{3}$. The ODEs that arise from separating the Schrödinger
equation in these coordinates are known as the generalised Lamé equation. The
joint spectrum is computed, along with the corresponding quantised actions and
their symmetry classes. The first degenerate system is analysed in Section 3.1
with the prolate coordinates. Separation yields a semi-toric system whose
monodromy matrix we compute. This contrasts with the oblate system in Section
3.2 which, while similar to prolate coordinates, does not possess monodromy
and instead has hyperbolic singularities. Section 3.3 studies the Lamé
coordinates (an extension of ellipsoidal coordinates on $S^{2}$ onto $S^{3}$)
while Sections 3.4 and 3.5 focus on the spherical and cylindrical systems
respectively. All these results are visualised by the Stasheff polytope in
Fig. 1. In the appendix we recall some results on $S^{2}$ for the only two
separable coordinates: ellipsoidal and spherical.

##  2 Ellipsoidal Coordinates and the Generalised Lamé Equation

Let $\bm{x}\coloneqq(x_{1},x_{2},x_{3},x_{4})$ be Cartesian coordinates on
$\mathbb{R}^{4}$ and consider the unit sphere $S^{3}\subset\mathbb{R}^{4}$
where
$S^{3}=\left\\{(x_{1},x_{2},x_{3},x_{4})\in\mathbb{R}^{4}\,|\,x_{1}^{2}+x_{2}^{2}+x_{3}^{2}+x_{4}^{2}=1\right\\}$.
The geodesic flow on $S^{3}$ is a constrained system on $T^{*}\mathbb{R}^{4}$
with corresponding momentum $\bm{y}=(y_{1},y_{2},y_{3},y_{4})$ such that
$\bm{x}\cdot\bm{x}=1$ and $\bm{x}\cdot\bm{y}=0$. Set

$H=\frac{1}{2}\sum_{i>j}\ell_{ij}^{2}$

to be the Hamiltonian of the geodesic flow where the angular momentum in the
$(i,j)$ plane is given by $\ell_{ij}\coloneqq x_{i}y_{j}-x_{j}y_{i}$ and the
symplectic structure is $\sum_{i}dx_{i}\wedge dy_{i}$.

Let a homogeneous polynomial $\Psi_{D}(x_{1},x_{2},\dots,x_{n})$ of degree $D$
be a solution to Laplace’s equation in $\mathbb{R}^{n}$,

$\Delta\Psi_{D}=\sum\partial_{i}^{2}\Psi_{D}=0\,.$

By definition $\Psi_{D}$ is called a _harmonic_ polynomial. Introduce the
radius $r=\sqrt{\sum x_{i}^{2}}$ and some local coordinates $s$ on the sphere
$S^{n-1}$. The Laplacian in such a coordinate system is

$\Delta_{\mathbb{R}^{n}}=\frac{\partial^{2}}{\partial
r^{2}}+\frac{n-1}{r}\frac{\partial}{\partial
r}+\frac{1}{r^{2}}\Delta_{S^{n-1}}\,.$

Now rewrite the harmonic polynomial as $\Psi_{D}(\bm{x})=r^{D}\psi_{D}(s)$ and
thus define the spherical harmonic $\psi_{D}:S^{n-1}\to\mathbb{R}$. Inserting
this into the above Laplacian gives

$-\Delta_{S^{n-1}}\psi_{D}=D(D-2+n)\psi_{D}\,,$

so we see that spherical harmonics are eigenfunctions of the Laplace-Beltrami
operator on the sphere. For $n=3$ on $S^{2}$ these are the usual spherical
harmonics typically denoted by $Y_{lm}$, for the case of general $n$ see,
e.g., [6]. The above Hamiltonian can be canonically quantised which leads to
$\hat{\ell}_{ij}=-i\hbar(x_{i}\partial_{j}-x_{j}\partial_{i})$. This again
gives the Laplace-Beltrami operator on $S^{n-1}$ and hence

$2\hat{H}\Psi_{D}=D(D-2+n)\Psi_{D}\,.$

In the following we will speak of both the spherical harmonic $\psi_{D}$ in
local coordinates $s$ on the sphere $S^{n-1}$, and the corresponding harmonic
polynomial $\Psi_{D}$ in global Cartesian coordinates $\bm{x}$ on
$\mathbb{R}^{n}$, as eigenfunctions of the Laplace-Beltrami operator. These
considerations apply to any dimension $n$, but in the following we will
consider $n=4$ only.

Ellipsoidal-spherical (also known as sphero-conal or spherical ellipsoidal)
coordinates $\bm{s}\coloneqq(s_{1},s_{2},s_{3})$ on $S^{3}$ are defined as the
roots of

$\displaystyle T(s)$ $\displaystyle=\sum_{i=1}^{4}\frac{x_{i}^{2}}{s-e_{i}}$

where $0\leq e_{j}\leq s_{j}\leq e_{j+1}$ for all $j=1,2,3$ and the $e_{i}$
are called the semi-major axes. This is by analogy with ellipsoidal
coordinates on the tri-axial ellipsoid which are different and defined through
$T(s)=1$. We will only use ellipsoidal-spherical coordinates in the following
and call them ellipsoidal for short since no confusion can arise. The
transformation to ellipsoidal coordinates is found by solving $T(s_{j})=0$
where $j=1,2,3$ together with $r^{2}=\sum_{i=1}^{4}x_{i}^{2}=1$ to give

$\displaystyle x_{i}^{2}$
$\displaystyle=\frac{\prod_{j=1}^{3}(s_{j}-e_{i})}{\prod_{k\neq
i}(e_{k}-e_{i})}.$ (1)

From [33] we have the following result.

######  Lemma 1.

The Hamiltonian $H$ for the geodesic flow on $S^{3}$ separates in ellipsoidal
coordinates (1) to give integrals

$\displaystyle\eta_{1}=$
$\displaystyle\sum_{i<j}\left(\ell_{ij}^{2}\sum_{k\neq i,j}e_{k}\right),$
$\displaystyle\eta_{2}=\sum_{i<j}\left(\ell_{ij}^{2}\prod_{k\neq
i,j}e_{k}\right).$ (2)

The triple $(H,\eta_{2},\eta_{2})$ is a Liouville integrable system on
$T^{*}S^{3}$.

Consider the following linear coordinate transformation

$\displaystyle\begin{aligned}
X_{1}&=\frac{1}{\sqrt{2}}(\ell_{12}+\ell_{34})&&&Y_{1}&=\frac{1}{\sqrt{2}}(\ell_{12}-\ell_{34})\\\
X_{2}&=\frac{1}{\sqrt{2}}(\ell_{13}-\ell_{24})&&&Y_{2}&=-\frac{1}{\sqrt{2}}(\ell_{13}+\ell_{24})\\\
X_{3}&=\frac{1}{\sqrt{2}}(\ell_{14}+\ell_{23})&&&Y_{3}&=\frac{1}{\sqrt{2}}(\ell_{14}-\ell_{23}).\end{aligned}$

Since $\sqrt{H}$ generates an $S^{1}$ flow with constant period, we are able
to reduce by the flow of $\sqrt{H}$ to $S^{2}\times S^{2}$ with local
coordinates $(\bm{X},\bm{Y})=(X_{1},X_{2},X_{3},Y_{1},Y_{2},Y_{3})$ and
Casimirs $\bm{X}\cdot\bm{X}=\bm{Y}\cdot\bm{Y}=1$. Again, from [33] we have the
following result.

######  Lemma 2.

The original integrable system $(H,\eta_{1},\eta_{2})$ on $T^{*}S^{3}$
descends to a two degree of freedom integrable system
$(\eta_{1}(\bm{X},\bm{Y}),\eta_{2}(\bm{X},\bm{Y}))$ on $S^{2}\times S^{2}$
with Poisson structure

$\displaystyle B_{\bm{X},\bm{Y}}$
$\displaystyle=\begin{pmatrix}\hat{\bm{X}}&\bm{0}\\\
\bm{0}&\hat{\bm{Y}}\end{pmatrix}$ (3)

where, for a vector $\bm{v}\in\mathbb{R}^{3}$, the corresponding antisymmetric
hat matrix $\hat{\bm{v}}$ is given by

$\displaystyle\hat{\bm{v}}\bm{u}$ $\displaystyle=\bm{v}\times\bm{u},$
$\displaystyle\forall\bm{u}\in\mathbb{R}^{3}.$ (4)

When written in terms of $\bm{X}$ and $\bm{Y}$ on $S^{2}\times S^{2}$ the
system can be quantised through Berezin-Toeplitz quantisation in the same way
it was done in [28, 27]. The quantisation is the same, but the operators are
different. In their case the global $S^{1}$ integral is linear and the
commuting second integral has linear and quadratic terms, while the integrals
$\eta_{1}$, $\eta_{2}$ are both homogeneous quadratic. In later chapters on
degenerate coordinate systems (e.g. the prolate case) linear integrals will
also appear, in particular the global $S^{1}$ actions are linear in
$\bm{X},\bm{Y}$. We will not use Berezin-Toeplitz quantisation, but instead
employ separation of variables, so that the connection to spherical harmonics
and special functions remains.

Quantising through separation of variables reduces the quantisation problem to
the problem of finding polynomial solutions of a single 2nd order Fuchsian
ODE. It would be interesting to study in detail the relationship between the
two approaches. The one thing that we do use from the Berezin-Toeplitz
approach is the fact that $\hbar$ is the inverse of an integer. This also
appears naturally through the separation of variables, but would be harder to
justify without reference to the Berezin-Toeplitz approach, because separation
of variables does not clearly execute the reduction to a system with compact
phase space.

To construct a quantised version of the Hamiltonian we derive the diagonal
metric tensor $g$ in ellipsoidal coordinates on the $3-$sphere given by (1).
This is an orthogonal coordinate systems, so the off-diagonal entries of the
metric vanish. The diagonal entries of the metric are

$g_{jj}=\frac{-(s_{j}-s_{k})(s_{j}-s_{m})}{4(s_{j}-e_{1})(s_{j}-e_{2})(s_{j}-e_{3})(s_{j}-e_{4})}$

where $k$ and $m$ are the two unique distinct indices different from $j$.
Denote the determinant of the metric by $g=\prod g_{jj}$. The inverse metric
with upper indices has entries $g^{jj}=1/g_{jj}$. Thus the Laplace-Beltrami
operator on $S^{3}$ in these coordinates is given by

$\Delta=-\frac{1}{\sqrt{g}}\sum\partial_{i}(\sqrt{g}\,g^{jj}\partial_{i})\,.$

When quantising with the semi-classical parameter $\hbar$ the corresponding
stationary Schrödinger equation is given by

$\displaystyle-\hbar^{2}\Delta\psi$ $\displaystyle=\tilde{E}\psi\,.$ (5)

As already remarked, it is possible to quantise in the original Euclidean
coordinates and replace $\ell_{ij}$ with $x_{i}\partial_{j}-x_{j}\partial i$,
see [21, 44]. We need separation of variables, however, and hence obtain
commuting operators via the Stäckel matrix of the system. The classical
commuting integrals (2) are obtained from separation of variables in the
Hamilton-Jacobi equation. The same separation works for the Schrödinger
equation because the sphere has constant curvature, and hence the Robertson
condition is satisfied, see, e.g. [24].

The reduction to a two degree of freedom integrable system is achieved on the
quantum level by fixing the eigenvalue of the Laplace-Beltrami operator, which
for $S^{3}$ is $E=D(D+2)$ for non-negative integer $D$.

######  Theorem 1.

There are three commuting 2nd order differential operators
$(\hat{H},\hat{\eta}_{1},\hat{\eta}_{2})$ on $S^{3}$ with rational
coefficients, whose classical limit are the integrals $(H,\eta_{1},\eta_{2})$
given in (2).

###### Proof.

A Stäckel matrix for ellipsoidal coordinates is given by

$\sigma_{\text{el}}=\frac{1}{4}\begin{pmatrix}-\frac{s_{1}^{2}}{A(s_{1})}&-\frac{s_{1}}{A(s_{1})}&-\frac{1}{A(s_{1})}\\\
-\frac{s_{2}^{2}}{A(s_{2})}&-\frac{s_{2}}{A(s_{2})}&-\frac{1}{A(s_{2})}\\\
-\frac{s_{3}^{2}}{A(s_{3})}&-\frac{s_{3}}{A(s_{3})}&-\frac{1}{A(s_{3})}\end{pmatrix}.$
(6)

where $A(z)=\prod_{k=1}^{4}(z-e_{k})$. Note that the metric and Stäckel matrix
(6) are related by

$\frac{1}{g_{jj}}=\frac{\det(\Omega_{j})}{\det(\sigma_{\text{el}})}$

where $\Omega_{j}$ is the minor formed by deleting the $j^{th}$ row and first
column of $\sigma_{\text{el}}$.

Denote the rows of the inverse of the Stäckel matrix by $r^{1},r^{2},r^{3}$.
By definition of the Stäckel matrix the first row contains the diagonal
entries of the metric $g^{jj}=r^{1j}$ with upper indices. Then the commuting
operators are

$\hat{\eta}_{k-1}=-\frac{1}{\sqrt{g}}\sum_{i}\partial_{i}(\sqrt{g}\,r^{ki}\partial_{i}),\quad
k=1,2,3.$

where $\hat{H}=\hat{\eta}_{0}$ corresponding to the first row of the inverse
of the Stäckel matrix is the Laplace-Beltrami operator. The corresponding
classical integrals are $\eta_{k-1}=\sum r^{ki}p_{i}^{2}$ as defined in (2) in
Cartesian coordinates. Here $p_{i}$ are the canonical curvilinear momenta
conjugate to the coordinates $s_{i}$. We denote the eigenvalues of
$(\hat{H},\hat{\eta}_{1},\hat{\eta}_{2})$ by $(E,\lambda_{1},\lambda_{2})$. ∎

The system of ellipsoidal coordinates on the sphere is used in the separation
of the quantised C. Neumann system, see e.g. [21, 22] and also [10, 8, 9]. Our
system before reduction is the Neumann system without potential, and hence
formulas in [21] can be specialised to our case. The reduction step (or fixing
the eigenvalue of the Schrödinger operator) does not make sense for the
Neumann system, because unlike the geodesic flow on $S^{3}$ it is not a
superintegrable system.

It is also possible to obtain the operators $\hat{\eta}_{i}$ by canonical
quantisation of the angular momentum operators expressed in ellipsoidal
coordinates using

$\ell_{ij}=2x_{i}(s)x_{j}(s)(e_{i}-e_{j})\sum_{k=1}^{3}c_{ij}(s_{k})p_{k},\quad
c_{ij}(s_{k})=\frac{\prod_{n\not=i,j}(s_{k}-e_{n})}{\prod_{m\not=k}(s_{k}-s_{m})}$

where $x_{i}(s)$ as given in (1). This equation is linear in momenta and can
be quantised by replacing $p_{k}\to-i\hbar\partial_{s_{k}}$. In this way the
classical commuting integrals expressed in terms of $\ell_{ij}$ can be
directly converted into the corresponding quantum operators. However, the
calculation using the Stäckel matrix is much simpler.

The inverse of the Stäckel matrix is used to define the commuting operators,
while the Stäckel matrix itself determines how the spectral parameters enter
the separated equation(s), which is the content of

######  Lemma 3.

In ellipsoidal coordinates (1) on $S^{3}$, the Schrödinger equation (5)
separates into

$\psi_{j}^{{}^{\prime\prime}}+\frac{1}{2}\left(\frac{1}{s_{j}-e_{1}}+\frac{1}{s_{j}-e_{2}}+\frac{1}{s_{j}-e_{3}}+\frac{1}{s_{j}-e_{4}}\right)\psi_{j}^{{}^{\prime}}+\frac{-Es_{j}^{2}+\lambda_{1}s_{j}-\lambda_{2}}{4(s_{j}-e_{1})(s_{j}-e_{2})(s_{j}-e_{3})(s_{j}-e_{4})}\psi_{j}=0$
(7)

for all $j=1,2,3$ with $E\coloneqq\frac{\tilde{E}}{\hbar^{2}}$, $\lambda_{1}$
and $\lambda_{2}$ are spectral parameters and $\psi=\psi_{1}\psi_{2}\psi_{3}$
where $\psi_{j}=\psi(s_{j})$.

###### Proof.

Let

$\displaystyle f_{j}\coloneqq$ $\displaystyle\frac{\partial}{\partial
s_{j}}\log(g^{jj}\sqrt{\det(g)}).$

Following [21] the equation (5) separates into

$\displaystyle\psi_{j}^{{}^{\prime\prime}}+f_{j}\psi_{j}^{{}^{\prime}}+\sum_{k=1}^{3}c_{k-1}(\sigma_{\text{el}})_{jk}\psi_{j}=0$

where $(c_{0},c_{1},c_{2})=(-E,\lambda_{1},-\lambda_{2})$ are spectral
parameters. Computing the relevant quantities gives (7). ∎

Equation (7) is known as the generalised Lamé equation. This is a second order
Fuchsian equation with $4$ finite regular singular points which can be writen
as

$\displaystyle\frac{d^{2}w}{dz^{2}}+\left(\sum_{j=1}^{4}\frac{\gamma_{j}}{z-e_{j}}\right)\frac{dw}{dz}+\left(\sum_{j=1}^{4}\frac{q_{j}}{z-e_{j}}\right)w$
$\displaystyle=0$ (8)

where $\sum_{j=1}^{4}q_{j}=0$ and the $q_{j}$ are known as the accessory
parameters. Often, we write the multiplicative term of (8) over a common
denominator, i.e.

$\sum_{j=1}^{4}\frac{q_{j}}{z-e_{j}}=\frac{c_{0}z^{2}+c_{1}z+c_{2}}{4\Pi_{j=1}^{4}(z-e_{j})}.$
(9)

The exponents at the pole $e_{j}$ are $(0,1-\gamma_{j})$ while those at
infinity are $(\alpha,\beta)$ defined by

$\displaystyle\alpha+\beta+1$ $\displaystyle=\sum_{j=1}^{4}\gamma_{j},$
$\displaystyle\alpha\beta$ $\displaystyle=\sum_{j=1}^{4}e_{j}q_{j}.$ (10)

From (10) it is clear that there are two free (so-called) accessory
parameters, while the more well known Lamé equation only has one, see [16]
31.14 for more detail. These two free parameters turn out to be the spectral
parameters, i.e. the eigenvalues of $\hat{\eta}_{1}$ and $\hat{\eta}_{2}$. For
some other examples of occurrences of the generalised Lamé equation, see [34,
31, 11] and the references therein.

In this paper, we seek polynomial solutions of (5) giving harmonic polynomials
on $S^{3}$ in Cartesian coordinates. Ellipsoidal coordinates are singular
along co-dimension one submanifolds of $S^{3}$, and these coordinate
singularities translate to the singularities of the Fuchsian ODE (8).
Solutions to the separated equations that are analytic at all finite
singularities correspond to polynomial solutions in the local coordinates
$s_{i}$, possibly after factoring out single powers of $x_{k}(s)$ to achieve
all possible discrete symmetries in the original coordinates. Thus polynomial
solutions of the Fuchsian ODEs in local coordinates lead to polynomial
solutions in the original cartesian coordinates. The quantisation conditions
at the level of the Fuchsian equation _are_ the conditions that enforce a
polynomial separated eigenfunction $\psi(s)$. In general these are conditions
for each factor $\psi_{j}(s_{j})$ of the separated eigenfunction
$\psi(s)=\prod\psi_{j}(s_{j})$, but it turns out that in the present highly
symmetric case the three equations are the same equation, just evaluated on
different intervals $e_{j}\leq s_{j}\leq e_{j+1}$, and hence a single
polynomial gives three $\psi_{j}(s_{j})$ and hence the complete eigenfunction
$\psi(s)$. In cases of more degenerate coordinate systems discussed in the
next section there will actually be different equations for the different
coordinates.

It is known from 31.15 of [16] and [45] that for every vector
$\bm{n}=(n_{1},n_{2},n_{3})$ of non negative integers, there are uniquely
determined real values of $(c_{0},c_{1},c_{2})$ such that (8) has polynomial
solutions with $n_{j}$ zeros in the interval $(e_{j},e_{j+1})$ for all
$j=1,2,3$. Such solutions $S_{d}(z)$ of degree $d$ where $d=n_{1}+n_{2}+n_{3}$
are known as Heine-Stieltjes polynomials.

From [45] we have the following Lemma for the spectral parameter $c_{0}$.

######  Lemma 4.

Let $S_{d}(z)$ be a Heine-Stieltjes polynomial solution to (8) of degree $d$
over the interval $[e_{1},e_{4}]$. Then

$\displaystyle\frac{1}{4}c_{0}$
$\displaystyle=-d(d-1+\sum_{i=1}^{4}\gamma_{i})=\alpha\beta$ (11)

and $(\alpha,\beta)=(-d,(-1+d+\sum\gamma_{i}))$. Further, let $K(d)$ be the
space of all such polynomials that satisfy (8). Then, we have

$\displaystyle\text{dim}(K(d))$ $\displaystyle\coloneqq\begin{pmatrix}d+2\\\
2\end{pmatrix}.$ (12)

Note that $\alpha$ is a negative integer in order to obtain polynomial
solutions. We also have from [45] the following Lemma relating polynomial
solutions of the Schrödinger equation in Cartesian coordinates to products of
Heine-Stieltjes polynomials.

######  Proposition 5.

The product $S_{d}(s_{1})S_{d}(s_{2})S_{d}(s_{3})$, where $S_{d}(s)$ is a
Heine-Stieltjes polynomial solution to (7) of degree $d$, when expressed in
the original Cartesian coordinates using (1), is a homogeneous polynomial of
degree $\tilde{D}\coloneqq 2d$ in the variables $x_{i}$ given by

$\displaystyle\Phi_{\tilde{D}}(\bm{x})\coloneqq\left(\prod_{i=1}^{4}\prod_{k=1}^{d}(z_{k}-e_{i})\right)\prod_{k=1}^{d}\sum_{i=1}^{4}\frac{x_{i}^{2}}{z_{k}-e_{i}}\,,$
(13)

where $z_{k}$ are the roots of $S_{d}(z)$.

In Cartesian coordinates there are 16 discrete symmetry classes of
eigenfunctions corresponding to parity about each of the $x_{i}$ axes. Let
$\bm{\mu}=(\mu_{1},\mu_{2},\mu_{3},\mu_{4})$ with $\mu_{i}\in\\{0,1\\}$ where
$\mu_{i}=1$ denotes a wavefunction odd about the respective axis and even
otherwise. By (13) the wave function $\Phi_{\tilde{D}}(\bm{x})$ is quadratic
in all $x_{i}$, and hence cannot give any odd symmetry. In order to find these
other symmetries, consider the change of dependent variable in (7) given by

$\phi_{j}=\prod_{i}(s_{j}-e_{i})^{\mu_{i}/2}\psi_{j}.$ (14)

This corresponds to multiplication by
$x_{1}^{\mu_{1}}x_{2}^{\mu_{2}}x_{3}^{\mu_{3}}x_{4}^{\mu_{4}}$ in Cartesian
coordinates. This yields the transformed generalised Lamé equation

$\phi_{j}^{{}^{\prime\prime}}+\left(\frac{\tilde{\gamma}_{1}}{z-e_{1}}+\frac{\tilde{\gamma}_{2}}{z-e_{2}}+\frac{\tilde{\gamma}_{3}}{z-e_{3}}+\frac{\tilde{\gamma}_{4}}{z-e_{4}}\right)\phi_{j}^{{}^{\prime}}+\frac{(u_{0}-E)z^{2}+(u_{1}+\lambda_{1})z+(u_{2}-\lambda_{2})}{4(z-e_{1})(z-e_{2})(z-e_{3})(z-e_{4})}\phi_{j}=0$
(15)

where

$\displaystyle\tilde{\gamma}_{m}$
$\displaystyle=\begin{cases}\frac{3}{2}&\text{if $\mu_{m}=1$}\\\
\frac{1}{2}&\text{if $\mu_{m}=0$}\end{cases}$ (16)

and the $u_{j}$ are shown in Table 1 for a given symmetry class $\bm{\mu}$.
Note that in the table we adopted the notation where $(\mu_{i})$ represents
symmetry classes where only one $\mu_{i}$ is $1$ (i.e. $(1,0,0,0),(0,1,0,0)$
etc). Similarly, $(\mu_{i},\mu_{j})$ represents all classes which are odd
about two axes and so forth for $(\mu_{i},\mu_{j},\mu_{k})$.

Denote by $S_{d}^{\bm{\mu}}(z)$ solutions to (15) for the symmetry class
$\bm{\mu}$. Repeating the same reasoning used to obtain (13), we let the
product
$S_{d}^{\bm{\mu}}(s_{1})S_{d}^{\bm{\mu}}(s_{2})S_{d}^{\bm{\mu}}(s_{3})$, when
converted to Cartesian coordinates, be given by
$\Phi_{\tilde{D}}^{\bm{\mu}}(\bm{x})$. Finally, we set

$\Psi_{D}^{\bm{\mu}}(\bm{x})\coloneqq
x_{1}^{\mu_{1}}x_{2}^{\mu_{2}}x_{3}^{\mu_{3}}x_{4}^{\mu_{4}}\Phi_{\tilde{D}}^{\bm{\mu}}(\bm{x})$
(17)

where $D=\tilde{D}+\sum_{i=1}^{4}\mu_{i}$. Note that
$\Phi_{\tilde{D}}(\bm{x})$ in (13) is the special case of (17) with
$\bm{\mu}=(0,0,0,0)$.

Symmetry | $u_{2}$ | $u_{1}$ | $u_{0}$
---|---|---|---
$(0,0,0,0)$ | $0$ | $0$ | $0$
$(\mu_{i})$ | $e_{j}e_{k}+e_{j}e_{m}+e_{k}e_{m}$ | $-2(e_{j}+e_{k}+e_{m})$ | $3$
$(\mu_{i},\mu_{j})$ | $e_{i}e_{k}+e_{i}e_{m}+e_{j}e_{k}+e_{j}e_{m}+4e_{k}e_{m}$ | $-2(e_{i}+e_{j})-6(e_{k}+e_{m})$ | $8$
$(\mu_{i},\mu_{j},\mu_{k})$ | $e_{i}e_{j}+e_{i}e_{k}+e_{j}e_{k}+4(e_{i}e_{m}+e_{j}e_{m}+e_{k}e_{m})$ | $-6(e_{i}+e_{j}+e_{k}+2e_{m})$ | $15$
$(1,1,1,1)$ | $4(e_{1}e_{2}+e_{1}e_{3}+e_{2}e_{3}+e_{1}e_{4}+e_{2}e_{4}+e_{3}e_{4})$ | $-12(e_{1}+e_{2}+e_{3}+e_{4})$ | $24$
Table 1: Coefficients $u_{j}$ in (15) for each symmetry class.

######  Theorem 2.

The function $\Psi_{D}^{\bm{\mu}}(\bm{x})$ given by (17) is an eigenfunction
of the Schrödinger operator (5). The energy eigenvalue $E$ is given by

$\displaystyle E$ $\displaystyle=D(D+2)$ (18)

where $D=2d+\sum_{i=1}^{4}\mu_{i}$ is the degree of the harmonic polynomial
$\Psi_{D}$.

###### Proof.

We have

$\Delta(\Psi_{D}^{\bm{\mu}}(\bm{x}))=\Delta(\Phi_{\tilde{D}}^{\bm{\mu}}(\bm{x}))+\sum_{j=1}^{4}\frac{2\mu_{j}}{x_{j}}\frac{\Phi_{\tilde{D}}^{\bm{\mu}}(\bm{x})}{\partial
x_{j}}.$ (19)

From [45], we know that since $\Phi_{\tilde{D}}^{\bm{\mu}}(\bm{x})$ satisfy
(15) for the appropriate parameter choices in Table 1, they also satisfy the
right hand side of (19). Hence, $\Psi_{D}^{\bm{\mu}}(\bm{x})$ are solutions of
(5). To compute the energy for a given symmetry class $\bm{\mu}$, we note that
the shift $u_{2}$ in Table 1 is simply $U(U+2)$ where
$U=\sum_{i=1}^{4}\mu_{i}$. Combining this with (11) and using $E=-c_{0}$ gives

$E=U(U+2)+4d(d+1+U)$ (20)

which simplifies to $E=D(D+2)$ for each class. ∎

From (17) we observe the following

######  Proposition 6.

For polynomial solutions to (5) of fixed degree $D$ corresponding to energy
$E=D(D+2)$, only 8 of the 16 total discrete symmetry classes can be present.
If $E$ is even, then the $(0,0,0,0)$, $(1,1,1,1)$ and all symmetry classes
even about 2 of the $x_{i}$ axes are present. Similarly, if $E$ is odd, then
the remaining 8 symmetry classes odd about and odd number of $x_{i}$ are
present.

We refer to solutions which are even about an even (odd) number of axes (with
even (odd) energies) as “even” (“odd”) solutions. In particular, the even
symmetry classes are

$\mathcal{S}_{E}=\\{(0,0,0,0),(1,1,0,0),(1,0,1,0),(1,0,0,1),(0,1,1,0),(0,1,0,1),(0,0,1,1),(1,1,1,1)\\}$
(21)

and the odd symmetry classes are

$\mathcal{S}_{O}=\\{(0,0,0,1),(0,0,1,0),(0,1,0,0),(1,0,0,0),(1,1,1,0),(1,1,0,1),(1,0,1,1),(0,1,1,1)\\}.$
(22)

Using the above results, we also have the following Lemma, which shows that
all eigenfunctions are obtained in this way.

######  Lemma 7.

The total number of eigenstates $N$ for a polynomial solution to (5) of degree
$D$ is given by

$N=(D+1)^{2}$ (23)

###### Proof.

Without loss of generality, we prove this for $D$ even (i.e. “even”
solutions). From (12), we know that the total number of eigenstates for the
$(0,0,0,0)$ symmetry class is $\begin{pmatrix}d+2\\\
2\end{pmatrix}=\begin{pmatrix}\frac{D+4}{2}\\\ 2\end{pmatrix}$. For a fixed
$D$, we similarly have the total number of states for an element of the
$(\mu_{i},\mu_{j})$ class as $\begin{pmatrix}\frac{D+2}{2}\\\ 2\end{pmatrix}$
and for the $(1,1,1,1)$ class we obtain $\begin{pmatrix}\frac{D}{2}\\\
2\end{pmatrix}$. The result follows since

$N=\begin{pmatrix}\frac{D+4}{2}\\\
2\end{pmatrix}+6\begin{pmatrix}\frac{D+2}{2}\\\
2\end{pmatrix}+\begin{pmatrix}\frac{D}{2}\\\ 2\end{pmatrix}=(D+1)^{2}.$

In Table 2 we show the number of eigenstates per symmetry class and energy
parity. ∎

Symmetry | Number of States | Energy Parity
---|---|---
$(0,0,0,0)$ | $\begin{pmatrix}(D+4)/{2}\\\ 2\end{pmatrix}$ | Even
$(\mu_{i})$ | $\begin{pmatrix}(D+3)/{2}\\\ 2\end{pmatrix}$ | Odd
$(\mu_{i},\mu_{j})$ | $\begin{pmatrix}(D+2)/{2}\\\ 2\end{pmatrix}$ | Even
$(\mu_{i},\mu_{j},\mu_{k})$ | $\begin{pmatrix}(D+1)/{2}\\\ 2\end{pmatrix}$ | Odd
$(1,1,1,1)$ | $\begin{pmatrix}D/{2}\\\ 2\end{pmatrix}$ | Even
Table 2: Number of states for each symmetry class and energy parity.

Recall that the basis of our study is separation of variables on the
$3-$sphere. It is no surprise then that the total energy (18) and number of
eigenstates (23) precisely coincide with standard results for spherical
harmonics on the sphere, see, e.g. [6]. To compute the joint spectrum of
$(\lambda_{1},\lambda_{2})$ in (7) for all 16 discrete symmetry classes, we
employ the following Lemma from [2].

######  Lemma 8.

Let $S_{d}$ be a Heine-Stieltjes polynomial of degree $d$ and denote its
(real) roots by $z_{1},\dots,z_{d}$. If, in (8), the $\gamma_{j}>0$ and
$e_{j}\in\mathbb{R}$ with $e_{j}<e_{j+1}$ for all $j$, then every root is a
solution to the system of equations

$\displaystyle\sum_{j=1}^{4}\frac{\gamma_{j}/2}{z_{k}-e_{j}}+\sum_{j=1,j\neq
k}^{d}\frac{1}{z_{k}-z_{j}}$ $\displaystyle=0,\quad k=1,\dots,d\,.$ (24)

For each solution $z_{1},\dots,z_{d}$ the accessory parameters $q_{j}$ are
given by

$\displaystyle q_{j}$
$\displaystyle=\gamma_{j}\sum_{k=1}^{d}\frac{1}{z_{k}-e_{j}},\quad
j=1,2,3,4\,,$

and in terms of these the spectral parameters $(\lambda_{1},\lambda_{2})$ are

$\displaystyle\lambda_{1}$ $\displaystyle=-4\sum e_{i}e_{j}(q_{k}+q_{m})$
(25a) $\displaystyle\lambda_{2}$ $\displaystyle=4\sum e_{i}e_{j}e_{k}q_{m}$
(25b)

where the sums are over all $i,j,k,m$ that are pairwise distinct.

This well known Lemma is surprising in the sense that the system of equations
(24) determines the eigenfunction $S_{d}$ and joint spectrum
$(\lambda_{1},\lambda_{2})$ simultaneously. Comparing (25a), (25b) and the
classical equation (2), the degrees in $e_{j}$ are different. The reason for
this is that $q_{j}$ for fixed energy is of degree $-1$ in $e_{j}$, as can be
seen from (11). Note that the solution to the nonlinear equations (24) is such
that automatically $\sum q_{i}=0$ and $\sum
q_{i}e_{i}=\alpha\beta=-d(d-1+\sum\gamma_{i})$, combining (11) and (10).

We can use the value of $\hbar$ to scale the energy eigenvalue to 1 in the
semiclassical limit $\hbar\to 0$, and similar scalings for the other
eigenvalues. After the scaling, changing $\hbar$ changes the number of states,
but not the size of the image of the momentum map. Recall the definition of
$E=\frac{\tilde{E}}{\hbar^{2}}$ and the discrete values of $E=D(D+2)$ given in
(18). To scale such that $\tilde{E}$ is exactly 1 would require
$\hbar=1/\sqrt{D(D+2)}$, which is, however, not an integer. The reduced
symplectic manifold is $S^{2}\times S^{2}$ and hence compact, and therefore we
require $\hbar$ to be the inverse of an integer. Thus

$\displaystyle\hbar$ $\displaystyle=\frac{1}{\sqrt{D(D+2)+1}}=\frac{1}{D+1}$
(26)

so that $1/\hbar^{2}$ is exactly the number of states in the reduced compact
system. With this definition of $\hbar$ we find

$\tilde{E}=E\hbar^{2}=\frac{D(D+2)}{(D+1)^{2}}=1-\hbar^{2}$

where the last equality is exact.

Using (24) with the 16 forms of the generalised Lamé equation, we produce
examples of the joint spectrum shown in Fig. 2 with $D=18$ (a) and $D=19$ (b).
We also show that the boundary of the joint spectrum is given by the classical
momentum map (solid black lines). For their derivation, see [33].

In our numerics, we seed roots $z_{k}$ randomly in the appropriate intervals
$(e_{j},e_{j+1})$ as outlined in [45]. Standard root finding techniques like
Newton’s method are then used to find accurate solutions to (24). Larger
values of $D$ result in substantially longer computation time.

Since we have the freedom to choose the quadratic integrals via an affine
transformation, we can give a clearer representation of the spectrum, as shown
in Fig. 3 a) and b) with the hyperbolic-hyperbolic point centred at the
origin. This clearly shows three regions with a $\mathbb{Z}^{2}$ lattice.

Note that lattice points in these regions are actually two points,
representing period doubling. Classically, this corresponds to two tori in the
pre-image of the classical momentum map. Moving from these chambers into the
fourth (bounded by the curve), we see a halving of visible dots as we see
period doubling occurring again. This phenomena has been described in more
detail in [41].

Figure 2: a) Joint spectrum $(\lambda_{1},\lambda_{2})$ where
$(e_{1},e_{2},e_{3},e_{4})=(1,2,5,8)$ with $D=18$ and b) $D=19$. A direct
correspondence between the coloured dots and symmetry classes is shown in
Table 3.

Figure 3: Spectra corresponding to Fig. 2 a) and b) after performing an affine transformation to map the hyperbolic-hyperbolic point to the origin.  $\bm{\mu}$ | Colour a) | $\bm{\mu}$ | Colour b)
---|---|---|---
$(0,0,0,0)$ | Blue | $(1,0,0,0)$ | Purple
$(1,1,0,0)$ | Red | $(0,1,0,0)$ | Yellow
$(1,0,1,0)$ | Magenta | $(0,0,1,0)$ | Grey
$(1,0,0,1)$ | Green | $(0,0,0,1)$ | Pink
$(0,1,1,0)$ | Brown | $(1,1,1,0)$ | Light Blue
$(0,1,0,1)$ | Black | $(1,1,0,1)$ | Dark Green
$(0,0,1,1)$ | Cyan | $(0,1,1,1)$ | Light Red
$(1,1,1,1)$ | Orange | $(1,0,1,1)$ | Dark Brown
Table 3: Symmetry class corresponding to each coloured eigenstate shown in
Figure 2.

Equipped with the joint spectrum $(\lambda_{1},\lambda_{2})$, we can
investigate the corresponding action variables. Using the Stäckel matrix (6)
and the results of [21] and [33], we have the following Lemma.

######  Lemma 9.

The three continuous, classical actions $(J_{1},J_{2},J_{3})$ are given by

$\displaystyle J_{1}=\frac{1}{\pi}\int_{e_{1}}^{\min(R_{1},e_{2})}p(z)dz$
$\displaystyle
J_{2}=\frac{1}{\pi}\int_{\max(R_{1},e_{2})}^{\min(R_{2},e_{3})}p(z)dz$
$\displaystyle J_{3}=\frac{1}{\pi}\int_{\max(R_{2},e_{3})}^{e_{4}}p(z)dz$ (27)

where the momentum $p(z)$ is given by

$p^{2}=\frac{-\tilde{E}z^{2}+\lambda_{1}z-\lambda_{2}}{4(z-e_{1})(z-e_{2})(z-e_{3})(z-e_{4})},$

and $R_{1},R_{2}$ denote the roots of $p^{2}$ with $0\leq R_{1}\leq R_{2}$.

It was shown in [33] that the 3 actions are linearly dependent, continuous and
satisfy the following Lemma.

######  Lemma 10.

The actions are constrained by the following relation

$\displaystyle J_{1}+J_{2}+J_{3}$ $\displaystyle=\sqrt{\tilde{E}}.$ (28)

The fact that $H$ can be written as a function of the sum of the actions in
(28) follows from the fact that $H$ is superintegrable. In Fig. 4 a) and b) we
show the actions of the even and odd spectra in Fig. 2 a) and b) respectively,
computed using (27) with $\tilde{E}=1$. Similar plots for each symmetry class
can be found in Appendix A.3.

The classical action variables even though continuous are complicated
functions, and defining corresponding quantum mechanical operators may be
difficult. Thus presenting the joint spectrum in the space of action variables
is not done by quantising these functions. Instead we take the eigenvalues of
the operators $(\hat{H},\hat{\eta}_{1},\hat{\eta}_{2})$ and map them to
$(J_{1},J_{2},J_{3})$ with the previous Lemma. In this way a striking
representation of the joint spectrum inside an equilateral triangle is
obtained. Is shown in [33] this triangle is rigid even when passing to
degenerate coordinates systems. Note that this triangle is not Delzant (our
system is not toric to begin with), but it seems to serve a similar role for
our class of quantum integrable systems. Near each corner of the triangle the
lattice is approximately $\mathbb{Z}^{2}$ with basis vectors given by the
sides of the triangle. The action triangle is a natural global representation
because there are three elliptic-elliptic equilibrium points in the reduced
system, so that it represents three different $\mathbb{Z}^{2}$ lattices in one
moment map.

Among the degenerate systems discussed in the next section the cylindrical
coordinate system does lead to a toric system on $S^{2}\times S^{2}$ with a
square as Delzant polytope. The triangle appears as the quotient of the square
by discrete symmetry group $\mathbb{Z}^{2}\times\mathbb{Z}^{2}$ acting by
reflection across the diagonals of the square, see below.

Figure 4: Joint spectrum in action variables in the ellipsoidal case
corresponding to the joint spectra shown in Fig 2 a) and b) respectively.

##  3 Degenerate Systems

We now study the degenerate systems arising from the generalised ellipsoidal
coordinates $(1234)$. We originally focus on the prolate, oblate and Lamé
coordinates; these correspond to the edges of the Stasheff polytope shown in
Fig. 1. Further degenerations of these yield the spherical and cylindrical
coordinates, represented by the corners of the polytope.

We will show that the separated ODEs for these coordinates can be obtained by
smoothly degenerating those of the generalised ellipsoidal coordinates (i.e.
the generalised Lamé equation). In doing this, we prove the following theorem
which is the classical analogue of Theorem 3 of [33]:

######  Theorem 3.

For each pair of parentheses in the labelled Stasheff polytope (see Fig. 1)
that enclose two adjacent numbers, the corresponding quantum integrable system
has an $SO(2)$ symmetry and a corresponding ODE that can be transformed to a
trivial ODE with trigonometric solutions.

For each pair of parentheses that enclose three adjacent numbers, the
corresponding quantum integrable system has a global $SO(3)$ symmetry and a
corresponding ODE that can be transformed to the Gegenbauer equation
(hypergeometric type). This applies also when there is a pair of parentheses
inside.

If there are two pairs of parentheses (corresponding to the edges of the
Stasheff polytope) then one ODE is of Heun type (Fuchsian with 4 regular
singular points).

If there are three pairs of parentheses (corresponding to the corners of the
Stasheff polytope) then one ODE is of hypergeometric type (Fuchsian with 3
regular singular points). For spherical coordinates this can be transformed to
the associated Legendre equation and for cylindrical coordinates this is the
Jacobi equation.

Note that a single ODE of Heun type gives solutions for _two_ factors of the
eigenfunction $\psi$. In the general ellipsoidal case treated in the previous
section the single Fuchsian ODE with 4 finite regular singular points gives
_three_ such factors. Any ODE of hypergeometric type (including the trivial
trigonometric one) gives a solution for a single factor of the eigenfunction.
Hence another way of stating which ODEs occur for the various cases is this:
On the inside/edges/corners of the Stasheff polytope the highest Fuchsian ODE
has 5/4/3 regular singular points, accounting for 3/2/1 of the factors of the
eigenfunctions. The remaining 0/1/2 ODEs are of hypergeometric type.

###  3.1 Prolate Coordinates and the Heun Equation

Prolate coordinates on $S^{3}$ are a degeneration of ellipsoidal coordinates
arising from limiting the middle two semi major axes to each other, i.e.
$e_{3}\to e_{2}^{+}$. We normalise the $e_{i}$ according to
$(e_{1},e_{2}=e_{3},e_{4})=(0,1,a)$.

From [26], an explicit representation of prolate coordinates is

$\displaystyle x_{1}^{2}$ $\displaystyle=\frac{s_{1}s_{3}}{a}$ $\displaystyle
x_{2}^{2}=-\frac{\left(s_{1}-1\right)s_{2}\left(s_{3}-1\right)}{a-1}$ (29)
$\displaystyle x_{3}^{2}$
$\displaystyle=\frac{\left(s_{1}-1\right)\left(s_{2}-1\right)\left(s_{3}-1\right)}{a-1}$
$\displaystyle
x_{4}^{2}=\frac{\left(a-s_{1}\right)\left(a-s_{3}\right)}{(a-1)a}$

where $0\leq s_{1},s_{2}\leq 1\leq s_{3}\leq a$ and a possible Stäckel matrix
is given by

$\sigma_{\text{pro}}=\frac{1}{4}\left(\begin{array}[]{ccc}-\frac{1}{\left(s_{1}-1\right)\left(s_{1}-a\right)}&-\frac{1}{\left(s_{1}-1\right)s_{1}\left(s_{1}-a\right)}&-\frac{1-a}{\left(s_{1}-1\right){}^{2}s_{1}\left(s_{1}-a\right)}\\\
0&0&-\frac{1}{\left(s_{2}-1\right)s_{2}}\\\
-\frac{1}{\left(s_{3}-1\right)\left(s_{3}-a\right)}&-\frac{1}{\left(s_{3}-1\right)s_{3}\left(s_{3}-a\right)}&-\frac{1-a}{\left(s_{3}-1\right){}^{2}s_{3}\left(s_{3}-a\right)}\end{array}\right).$
(30)

From [33] we have the following result.

######  Lemma 11.

Separating the Hamilton-Jacobi equation in prolate coordinates and
subsequently reducing by $\sqrt{H}$ gives a semi-toric integrable system on
$S^{2}\times S^{2}$ with integrals

$\displaystyle G=a\ell_{12}^{2}+a\ell_{13}^{2}+\ell_{14}^{2}$ $\displaystyle
M=\ell_{23}^{2}$

and Poisson structure $B_{\bm{X},\bm{Y}}$ as given in (3).

To obtain the separated equations, we use the Stäckel matrix (30) and the
metric tensor for prolate coordinates (29), and follow the same methodology as
described in Lemma 3. For all future systems studied in this paper, this same
technique is used to obtain the separated ODEs. This gives the following.

######  Lemma 12.

Separating the Schrödinger equation (5) in prolate coordinates (29) gives the
following separated equations

$\displaystyle\psi_{j}^{{}^{\prime\prime}}+\frac{1}{2}\left(\frac{1}{s_{j}}+\frac{2}{s_{j}-1}+\frac{1}{s_{j}-a}\right)\psi_{j}^{{}^{\prime}}+\frac{-Es_{j}^{2}+\left(\lambda+E+(a-1)m^{2}\right)s_{j}-\lambda}{4s_{j}(s_{j}-1)^{2}(s_{j}-a)}\psi_{j}$
$\displaystyle=0$ for $j=1,3$ (31a)
$\displaystyle\psi_{j}^{{}^{\prime\prime}}+\frac{1}{2}\left(\frac{1}{s_{j}}+\frac{1}{s_{j}-1}\right)\psi_{j}^{{}^{\prime}}+\frac{m^{2}}{4s_{j}(1-s_{j})}\psi_{j}$
$\displaystyle=0$ for $j=2$ (31b)

where $(m,\lambda)$ are spectral parameters.

We observe that the separated equation in (31b) is of hypergeometric type, but
via the coordinate transformation $s_{2}=\cos^{2}\phi$ we recover the trivial
ODE

$\displaystyle\psi^{{}^{\prime\prime}}_{2}(\phi)+m^{2}$ $\displaystyle=0.$
(32)

To ensure a smooth globally defined solution on $S^{3}$, we enforce periodic
boundary conditions. This yields a discrete spectral parameter
$m\in\mathbb{Z}$ representing the quantised angular momentum in the
$(x_{2},x_{3})$ plane with solutions $e^{im\phi}$ to (32).

In [33] we showed that the classical integrals for the degenerate systems are
obtained by a smooth degeneration of the coordinates. The following Lemma
shows how the ODEs (31a) and (31b) naturally arise from the generalised Lamé
equation (7) in a similar fashion.

######  Lemma 13.

The separated equations for $\psi_{j}$ in (31a) and (31b) can be obtained by
smoothly degenerating the generalised Lamé equation (7).

###### Proof.

As noted in [33], the canonical transformation from ellipsoidal to prolate
coordinates is given by

$\displaystyle e_{3}$ $\displaystyle=e_{2}+\epsilon,$ $\displaystyle s_{2}$
$\displaystyle=e_{2}+\epsilon\tilde{s}_{2},$ $\displaystyle p_{2}$
$\displaystyle=\frac{\tilde{p}_{2}}{\epsilon},$ (33)

in the limit $\epsilon\to 0$ where $\tilde{s}_{2}\in[0,1]$ and $\tilde{p}_{2}$
is its conjugate momenta. We normalise the semi-major axes with
$(e_{1},e_{2}=e_{3},e_{4})=(0,1,a)$ and set

$\displaystyle\lambda=\lambda_{2},$ $\displaystyle
m^{2}=\frac{1}{a-1}(\lambda_{1}-\lambda_{2}-E).$

For the equations in $s_{1}$ and $s_{3}$, the first derivative terms trivially
degenerate while for the multiplicative term, substituting (33) into (7) and
taking the limit $\epsilon\to 0$ gives (31a).

Repeating the same process for the equations in $s_{2}$ gives the desired
result. ∎

Note that the Stäckel matrix used to obtain the separation constants is not
unique and so neither are the arising ODEs. Another choice of Stäckel matrix
would result in a different (yet equivalent) set of ODEs since the separation
constants would be different. The first derivative term is specific to the
coordinates and so is independent of the choice of Stäckel matrix.

Equation (31a) is an example of a Heun equation; a Fuchsian $4$ ODE with
regular singularities at $0,1,a$ and $\infty$:

$\displaystyle
W^{{}^{\prime\prime}}+\left(\frac{\gamma}{z}+\frac{\delta}{z-1}+\frac{\epsilon}{z-a}\right)W^{{}^{\prime}}+\frac{\alpha\beta
z-q}{z(z-1)(z-a)}W$ $\displaystyle=0$ (34)

where $\alpha+\beta+1=\gamma+\delta+\epsilon$. The corresponding Riemann
symbol is given by

$\displaystyle\mathcal{S}_{\text{Heun}}$ $\displaystyle=\left(\begin{aligned}
0&&1&&a&&\infty\\\ 0&&0&&0&&\alpha&;z\\\
1-\gamma&&1-\delta&&1-\epsilon&&\beta\end{aligned}\right).$ (35)

######  Lemma 14.

By the change of dependent variable $\psi_{j}=(z-1)^{|m|/2}W_{j}$, (31a) can
be written in the form of (34) where
$\gamma=\epsilon=\frac{1}{2},\delta=1+|m|,q=\frac{1}{4}\left(a|m|-\lambda\right)$
and

$\displaystyle\alpha,\beta$
$\displaystyle=\frac{1}{2}\left(1\pm\sqrt{1+E}+\left|m\right|\right).$ (36)

Specifically, we have

$W_{j}^{{}^{\prime\prime}}+\left(\frac{1}{2z}+\frac{1+|m|}{2(z-1)}+\frac{1}{2(z-a)}\right)W_{j}^{{}^{\prime}}+\frac{(-E+|m|(|m|+2))z-(a|m|-\lambda)}{4z(z-1)(z-a)}W_{j}$
(37)

The solution $H\ell(a,q;\alpha,\beta,\gamma,\delta;z)$ that corresponds to
exponent $0$ about $z=0$ such that $H\ell(a,q;\alpha,\beta,\gamma,\delta;0)=1$
is called a Heun function. We use the notation from [16]. Assuming
$\gamma\notin-\mathbb{N}$, the Heun function can be expanded as an infinite
series

$\displaystyle H\ell(a,q;\alpha,\beta,\gamma,\delta;z)$
$\displaystyle=\sum_{i=0}^{\infty}c_{i}z^{i}$ (38)

where $|z|<1$ and $c_{0}=1$. Substituting (38) into (34) gives the following
Lemma taken from [37].

######  Lemma 15.

The coefficients $c_{i}$ in (38) satisfy the following three term recurrence
relation

$\displaystyle A_{i}c_{i+1}-(B_{i}+q)c_{i}+C_{i}c_{i-1}$ $\displaystyle=0$
(39)

where, for $i\geq 1$:

$\displaystyle A_{i}$ $\displaystyle=a(i+1)(i+\gamma)$ $\displaystyle B_{i}$
$\displaystyle=i\left[(i-1+\gamma)(a+1)+a\delta+\epsilon\right]$ (40)
$\displaystyle C_{i}$ $\displaystyle=(i-1+\alpha)(i-1+\beta)$

subject to the conditions $c_{0}=1$ and

$a\gamma c_{1}-qc_{0}=0.$ (41)

From (40) it is clear that $\alpha=-d$ where $d\in\mathbb{N}$ forces
$C_{d+1}=0$ and the three term recurrence (39) truncates (38) to a polynomial
of degree $d$. These solutions, analytic at all three finite singularities
$z=0,1,a$ are known as Heun polynomials and will be denoted by
$Hp(a,q;\alpha,\beta,\gamma,\delta;z)$. In the general ellipsoidal case Heine-
Stieltjes polynomials were computed by solving a non-linear system of
equations for the roots of the eigenfunctions. In the Heun case the above
recursion gives a different and simpler method for computing the
eigenfunctions and eigenvalues. In the remainder of this section we show how
to obtain the corresponding harmonic polynomials in $\mathbb{R}^{4}$ and how
to compute the various discrete symmetry classes of eigenfunctions.

Using a technique similar to the ellipsoidal case and outlined by [45], we
have the following Lemma.

######  Lemma 16.

Let $Hp_{d}(z)$ be a Heun polynomial solution to (37) of degree $d$ and denote
by $z_{1},\dots,z_{d}$ its roots. The product $Hp_{d}(s_{1})Hp_{d}(s_{3})$,
expressed in the original Cartesian coordinates, is given by

$\displaystyle Hp_{d}(s_{1})Hp_{d}(s_{3})$
$\displaystyle=\prod_{k=1}^{d}z_{k}(z_{k}-a)\left[\frac{x_{1}^{2}}{z_{k}}-r^{2}+x_{4}^{2}\frac{1-a}{z_{k}-a}\right]$
(42)

where $r^{2}=x_{1}^{2}+x_{2}^{2}+x_{3}^{2}+x_{4}^{2}$.

###### Proof.

Firstly, we claim that

$\displaystyle(\theta-s_{1})(\theta-s_{3})$
$\displaystyle=\theta(\theta-a)\left(r^{2}-\frac{x_{1}^{2}}{\theta}+\frac{(a-1)}{\theta-a}x_{4}^{2}\right)$
(43)

for all $\theta\neq 0,a$. Both the left and right hand sides of (43) are monic
quadratic functions of $\theta$. Further, both expressions have roots at
$\theta=s_{1},s_{3}$. This is easily verified for the right hand side by using
the coordinate transform (29) and the condition $r^{2}=1$.

Let $Hp_{d}(z)$ be a Heine-Stieltjes polynomial of degree $d$ with roots
$z_{1},z_{2},\dots,z_{d}$ in factored form, i.e.

$Hp_{d}(z)=(z_{1}-z)(z_{2}-z)\dots(z_{d}-z).$ (44)

Combining (43) with (44) gives the desired result. ∎

Note that we have solved the canonical Heun equation (34). To write the
combined solution to (31a) and (31b), we include the normalisation factor by
recalling $\psi_{1}=(1-s_{1})^{|m|/2}W_{1}$. For the $s_{3}$ equation, since
$s_{3}\geq 1$ we have $\psi_{3}=(s_{3}-1)^{|m|/2}W_{3}$. Hence, the combined
solution to (5) in prolate coordinates is the product

$\psi_{1}\psi_{2}\psi_{3}=(1-s_{1})^{|m|/2}(s_{3}-1)^{|m|/2}Hp_{d}(s_{1})Hp_{d}(s_{2})e^{im\phi}.$
(45)

This yields the following Lemma.

######  Lemma 17.

Expressed in Cartesian coordinates, the product $\psi_{1}\psi_{2}\psi_{3}$ in
(45) can be written as the following $\tilde{D}$ degree homogeneous polynomial
where $\tilde{D}=2d+|m|$:

$\Phi_{\tilde{D}}(\bm{x})=(a-1)^{|m|/2}\left(x_{2}+i\operatorname{sign}(m)x_{3}\right)^{|m|}\prod_{k=1}^{d}z_{k}(z_{k}-a)\left[\frac{x_{1}^{2}}{z_{k}}-r^{2}+x_{4}^{2}\frac{1-a}{z_{k}-a}\right]$
(46)

and $m\in[-\tilde{D},\tilde{D}]$ is an integer.

###### Proof.

Using the definition of prolate coordinates in (29) we obtain

$(1-s_{1})^{m/2}(s_{3}-1)^{m/2}=\left[(a-1)(x_{2}^{2}+x_{3}^{2})\right]^{m/2}.$

We also know that
$e^{im\phi}=\left(\frac{x_{2}+ix_{3}}{\sqrt{x_{2}^{2}+x_{3}^{2}}}\right)^{m}$.
Combining this with (42) and (45) gives the product in (46). For homogeneity,
given that $r^{2}=\sum_{i=1}^{4}x_{i}^{2}$, it is clear that all terms in (46)
are of degree $\tilde{D}=2d+|m|$. ∎

Like in the ellipsoidal case, there are $16$ discrete symmetry classes of
$\Phi_{\tilde{D}}$ corresponding to parity about each of the $x_{i}$ axes.
However, we note that a simple parity flip of either $x_{2}$ or $x_{3}$ for
the the complex phase $(x_{2}+ix_{3})^{m}$ would not necessarily yield a
materially different wave function depending on the parity of $m$.

To address this, we note that for odd $m$ we have
$\Re(x_{2}+ix_{3})^{m}=x_{2}P(x_{2}^{2},x_{3}^{2})$ and
$\Im(x_{2}+ix_{3})^{m}=x_{3}Q(x_{2}^{2},x_{3}^{2})$ where $P$ and $Q$ are some
polynomials of degree $(m-1)/2$. These therefore give solutions which are odd
about the $x_{2}$ $(x_{3})$ axes and even about $x_{3}$ $(x_{2})$
respectively. Similarly, we note that for even $m$ we have
$\Re(x_{2}+ix_{3})^{m}=\tilde{P}(x_{2}^{2},x_{3}^{2})$ and
$\Im(x_{2}+ix_{3})^{m}=x_{2}x_{3}\tilde{Q}(x_{2}^{2},x_{3}^{2})$ where
$\tilde{P}$ and $\tilde{Q}$ are polynomials of degree $m/2$ and $(m-2)/2$
respectively. The former give solutions which are even about both $x_{2}$ and
$x_{3}$ and the latter those which are odd about both of these axes.

In the following discussion, we only focus on the discrete symmetries about
the $x_{1}$ and $x_{4}$ axes. This is because symmetries about the $x_{2}$ and
$x_{3}$ axes are encoded in taking the real and imaginary components of the
wave function, in conjunction with a parity choice of $m$, as just described.
This is similar in other degenerate systems (oblate, spherical, etc).

Let $\bm{\mu}=(\mu_{1},\mu_{4})$ where $\mu_{i}\in\\{0,1\\}$ and $\mu_{i}$
being $1$ ($0$) denotes a solution odd (even) about the $x_{i}$ axis. For
these symmetries, we consider the change of dependent variable in (37) given
by $\phi_{1}=s_{1}^{\mu_{1}/2}(a-s_{1})^{\mu_{1}/2}W_{1}$ and
$\phi_{3}=s_{3}^{\mu_{4}/2}(a-s_{3})^{\mu_{4}/2}W_{3}$ which corresponds to
multiplication of the total wavefunction by $x_{1}^{\mu_{1}}x_{4}^{\mu_{4}}$.
Doing so gives a resulting Heun equation with parameters shown in Table 4.

Denote by $Hp_{d}^{\bm{\mu}}(z)$ solutions to (34) corresponding to the
$\bm{\mu}$ symmetry class. Let the product
$Hp_{d}^{\bm{\mu}}(s_{1})Hp_{d}^{\bm{\mu}}(s_{3})e^{im\phi}$, when converted
back to Cartesian coordinates, be given by
$\Phi_{\tilde{D}}^{\bm{\mu}}(\bm{x})$. As with the ellipsoidal case, we set

$\Psi_{D}^{\bm{\mu}}(\bm{x})\coloneqq
x_{1}^{\mu_{1}}x_{4}^{\mu_{4}}\Phi_{\tilde{D}}^{\bm{\mu}}(\bm{x})$ (47)

where $D\coloneqq\tilde{D}+\mu_{1}+\mu_{4}$ and $m\in[-D,D]$ is an integer for
all symmetry classes. Note that (46) is a special case of
$\Psi_{D}^{\bm{\mu}}(\bm{x})$ with $\bm{\mu}=(0,0)$.

######  Lemma 18.

The $\Psi_{D}^{\bm{\mu}}(\bm{x})$ are $D$ degree, homogeneous harmonic
polynomials and the energy eigenvalue is given by $E=D(D+2)$.

###### Proof.

The proof that $\Psi_{D}^{\bm{\mu}}(\bm{x})$ satisfy (5) is identical to the
ellipsoidal case. For the energy eigenvalue, we note that a necessary
condition for polynomial solutions is for the three-term recurrence relation
(39) to vanish. This forces $C_{i}=0$ in (40), leading to the condition

$\alpha=-d+\frac{1}{2}(\mu_{1}+\mu_{4}).$ (48)

Substituting (48) into the expression for $\alpha$ in (36) gives the desired
result. ∎

From (48) we observe that for the $(0,0)$ and $(1,1)$ symmetry classes,
$\alpha$ is a negative integer while for the $(0,1)$ and $(1,0)$ classes,
$\alpha$ is a negative half integer. Further, we note that fixing a value of
$E$ (in turn fixing $D$) enforces a parity relationship between $E$ and $m$
for a given symmetry class. For even (odd) $E$, $m$ must be even (odd) for the
$(0,0)$ and $(1,1)$ classes while $m$ must be odd (even) for the others.

$(\mu_{1},\mu_{4})$ | $\tilde{\alpha}$ | $\tilde{\beta}$ | $\tilde{\gamma}$ | $\tilde{\delta}$ | $\tilde{\epsilon}$ | $\tilde{q}$
---|---|---|---|---|---|---
$(0,0)$ | $\alpha$ | $\beta$ | $\gamma$ | $\delta$ | $\epsilon$ | $q$
$(1,0)$ | $\alpha+1-\gamma$ | $\beta+1-\gamma$ | $2-\gamma$ | $\delta$ | $\epsilon$ | $(1-\gamma)(a\delta+\epsilon)+q$
$(0,1)$ | $\alpha+1-\epsilon$ | $\beta+1-\epsilon$ | $\gamma$ | $\delta$ | $2-\epsilon$ | $\gamma(1-\epsilon)+q$
$(1,1)$ | $2+\alpha-\gamma-\epsilon$ | $2+\beta-\gamma-\epsilon$ | $2-\gamma$ | $\delta$ | $2-\epsilon$ | $q+2+a\delta-\epsilon-\gamma(1+a\delta)$
Table 4: Parameters of (34) for the various symmetry classes.

As in the ellipsoidal case, we can compute the total number of states for a
fixed energy $E$ as well as that for each symmetry class.

######  Lemma 19.

Fix an energy $E=D(D+2)$. Then the total number of eigenstates is given by

$\displaystyle N$ $\displaystyle=(D+1)^{2}.$ (49)

The total number of eigenstates for each symmetry class is shown in Table 5.

###### Proof.

Assume $D$ is even and consider the $(\mu_{1},\mu_{4})=(0,0)$ symmetry class;
the other cases are proven in a similar manner. We observe that $|m|$ ranges
over the $D+1$ values $-D,-D+2,\dots,0,\dots,D-2,D$ since $|m|$ is required to
be an even integer based on our condition for $\alpha$. For a fixed energy and
the relationship between $d$ and $|m|$ this enforces, we have a total of
$\frac{D-|m|}{2}+1$ states for a given $m$.

Summing the total number of states over all possible values of $m$ gives

$\displaystyle N_{(0,0)}^{even}=$
$\displaystyle\left(\frac{D}{2}+1\right)+2\times\left(\frac{D}{2}+\left(\frac{D}{2}-1\right)+\dots+1\right)=\left(\frac{D}{2}+1\right)^{2}.$

Repeating a similar procedure for all other symmetry classes gives the desired
result. Regardless of the parity of $D$, summing over each column of Table 5
gives (49). ∎

$(\mu_{1},\mu_{4})$ |  $E$ even |  $E$ odd
---|---|---
$(0,0)$ | $\left(\frac{D}{2}+1\right)^{2}$ | $\left(\frac{D+1}{2}\right)\left(\frac{D+3}{2}\right)$
$(1,0)$ | $\left(\frac{D}{2}+1\right)\frac{D}{2}$ | $\left(\frac{D+1}{2}\right)^{2}$
$(0,1)$ | $\left(\frac{D}{2}+1\right)\frac{D}{2}$ | $\left(\frac{D+1}{2}\right)^{2}$
$(1,1)$ | $\left(\frac{D}{2}\right)^{2}$ | $\left(\frac{D+1}{2}\right)\left(\frac{D-1}{2}\right)$
Table 5: Number of states for each symmetry class and each parity of $E$.

As with the ellipsoidal case, we have $\hbar$ as the inverse of an integer,
namely $\hbar=\frac{1}{D+1}$. The joint spectrum $(m\hbar,\lambda\hbar^{2})$
is obtained for fixed $\tilde{E}$ and total degree $d$ by numerically
computing the eigenvalues of the matrix (40) to obtain $d$ eigenvalues $q$ and
hence $\lambda$ for given $m$. We use Mathematica’s Eigenvalue[] command and
the numerics is fast and accurate for moderately large $D=2d+|m|\approx 20$.
Since the matrices (40) are not self adjoint, substantially larger $d$ results
in non-trivial numerical error.

An example of the joint spectrum is shown in Fig. 5 a) below. We have chosen
$D=20$ giving $\hbar=\frac{1}{21}$ and a total of $(20+1)^{2}=441$ total
eigenstates. The blue, orange, red and cyan points correspond to the
$(0,0),(1,0),(0,1),(1,1)$ symmetry classes respectively. Since we have chose
an even value of $E$ (and thus $D$) we observe states of the $(0,0)$ and
$(1,1)$ symmetry classes only along even $m$. Similarly, along odd $m$ only
the $(1,0)$ and $(0,1)$ states are present. Had we chosen an odd $E$ then
orange/red dots would be seen along even $m$ and vice versa. We also show the
classical boundary of the momentum map in black. As mentioned, this is a semi-
toric system and an isolated rank-2 critical point of focus-focus variety is
located at $(m\hbar,\lambda\hbar^{2})=(0,1)$. This is shown as the magenta
dot.

From the same reasoning used to obtain (27), we obtain the following formulae
for the actions

$\displaystyle J_{1,pro}=\frac{1}{\pi}\int_{0}^{\min(r_{1},1)}p_{1}ds_{1}$
$\displaystyle J_{2,pro}$ $\displaystyle=\frac{1}{\pi}\int_{0}^{1}p_{2}ds_{2}$
$\displaystyle J_{3,pro}=\frac{1}{\pi}\int_{\max(1,r_{2})}^{a}p_{3}ds_{3}$
(50)

where, using (30), for $j=1,3$

$p_{j}^{2}=\frac{\tilde{E}s_{j}^{2}+\left(\lambda+\tilde{E}+(a-1)m^{2}\right)s_{j}-\lambda}{4s_{j}(s_{j}-1)^{2}(s_{j}-a)}$
(51)

and $(r_{1},r_{2})$ denote the roots of $p_{j}^{2}$ with $0\leq r_{1}\leq
1\leq r_{2}\leq a$. Note that $J_{2}$ simplifies to
$m=\left|\ell_{23}\right|$. Like for the ellipsoidal system, the prolate
actions also satisfy (28). The action map for the prolate system is shown in
Fig. 5 b) where the magenta dot corresponding to the focus-focus point is
located at
$\frac{2}{\pi}(\sin^{-1}(\frac{1}{a}),0,\frac{\pi}{2}-\sin^{-1}(\frac{1}{a}))$.
Eigenstates with $m=0$ have $J_{1}=0$ and are hence on the boundary of the
triangle. The reason for this is that the discrete symmetry reduced action
variables in the quantum setting would require to impose Neumann boundary
conditions.

Figure 5: a) Prolate spectrum with $D=20$, $a=2.4$ and a total of $21^{2}$
eigenstates. Different symmetry classes are represented with different
colours. Note the focus-focus critical value shown in magenta. b) Action map
corresponding to the joint spectrum.

Since the classical integrable system is semi-toric and has a focus-focus
point, we expect the corresponding joint spectrum to exhibit quantum
monodromy. This is shown in Fig. 5 a). A unit cell is parallel transported
around the focus-focus critical value. Denote the basis vectors of this unit
cell as $v_{1}$ (vertical) and $v_{2}$ (horizontal). As the cell completes a
full loop, we observe a basis transformation where $v_{1}$ remains unchanged
but $v_{2}$ is updated to $v_{2}+2v_{1}$. In a similar vein to our previous
work [14], we have the following basis transformation

$\begin{pmatrix}v_{1}^{{}^{\prime}}\\\
v_{2}^{{}^{\prime}}\end{pmatrix}=\begin{pmatrix}1&0\\\
\omega&1\end{pmatrix}\begin{pmatrix}v_{1}\\\ v_{2}\end{pmatrix}$ (52)

where $\omega=2$. As was observed in [14], when considering the symmetry
classes one at a time, $\omega=1$ whereas combining them two at a time gives
$\omega=2$. Finally, in Fig. 6 we show the projection of the joint spectrum in
the action variables onto either the $(m,J_{1})$ or $(m,J_{3})$ axes (and
plotting the signed $m$ instead of $|m|$). The resulting polygons are two
different representations of the semi-toric polygon invariant with the cut
above and below the focus-focus point.

Figure 6: a) The semi-toric polygon in the prolate case $D=20$, $a=2.4$;
projection of the actions onto the $(m,J_{1})$ axes, b) Projection onto the
$(m,J_{3})$ axes.

###  3.2 Oblate Coordinates

Degenerating either the two largest $(e_{3}=e_{4})$ or smallest
$(e_{1}=e_{2})$ semi major axes gives rise to the oblate coordinates on
$S^{3}$. They are equivalent under sending $e_{i}\to-e_{i}$ and then
translating and scaling. In this paper we only address the former case and
normalise the $e_{i}$ with $(e_{1},e_{2},e_{3}=e_{4})=(0,1,a)$.

From [26], an explicit representation of oblate coordinates is

$\displaystyle x_{1}^{2}$ $\displaystyle=\frac{s_{1}s_{2}}{a}$ $\displaystyle
x_{2}^{2}=\frac{-\left(s_{1}-1\right)\left(s_{2}-1\right)}{a-1}$ (53)
$\displaystyle x_{3}^{2}$
$\displaystyle=\frac{\left(s_{1}-a\right)\left(s_{2}-a\right)s_{3}}{a\left(a-1\right)}$
$\displaystyle
x_{4}^{2}=\frac{\left(s_{1}-a\right)\left(s_{2}-a\right)\left(1-s_{3}\right)}{a\left(a-1\right)}$

where $0\leq s_{1},s_{2}\leq 1\leq s_{3}\leq a$. A possible Stäckel matrix is
given by

$\Phi_{\text{obl}}=\frac{1}{4}\left(\begin{array}[]{ccc}-\frac{1}{\left(s_{1}-1\right)\left(s_{1}-a\right)}&-\frac{1}{\left(s_{1}-1\right)s_{1}\left(s_{1}-a\right)}&\frac{(1-a)a}{\left(s_{1}-1\right)s_{1}\left(s_{1}-a\right){}^{2}}\\\
-\frac{1}{\left(s_{2}-1\right)\left(s_{2}-a\right)}&-\frac{1}{\left(s_{2}-1\right)s_{2}\left(s_{2}-a\right)}&\frac{(1-a)a}{\left(s_{2}-1\right)s_{2}\left(s_{2}-a\right){}^{2}}\\\
0&0&-\frac{1}{\left(s_{3}-1\right)s_{3}}\end{array}\right).$ (54)

From [33] and similar to the prolate system, we have the following.

######  Lemma 20.

Separating the Hamilton Jacobi equation in oblate coordinates using the
Stäckel matrix (54) and subsequently reducing by $H$ gives a two degree of
freedom integrable system on $S^{2}\times S^{2}$ with integrals

$\displaystyle G_{\text{obl}}=a\ell_{12}^{2}+\ell_{13}^{2}+\ell_{14}^{2}$
$\displaystyle M=\ell_{34}^{2}$

and Poisson structure $B_{\bm{X},\bm{Y}}$ as given in (3).

Using (54) and a computation similar to the prolate case we have the following
Lemma.

######  Lemma 21.

Separating the Schrödinger equation (5) in oblate coordinates (53) gives the
following ODEs

$\displaystyle\psi_{j}^{{}^{\prime\prime}}+\frac{1}{2}\left(\frac{1}{s_{j}}+\frac{1}{s_{j}-1}+\frac{2}{s_{j}-a}\right)\psi_{j}^{{}^{\prime}}+\frac{-Es_{j}^{2}+\left(\lambda+aE-(a-1)m^{2}\right)s_{j}-a\lambda}{4s_{j}(s_{j}-1)(s_{j}-a)^{2}}\psi_{j}$
$\displaystyle=0$ for $j=1,2$ (55a)
$\displaystyle\psi_{j}^{{}^{\prime\prime}}+\frac{1}{2}\left(\frac{1}{s_{j}}+\frac{1}{s_{j}-1}\right)\psi_{j}^{{}^{\prime}}+\frac{m^{2}}{4s_{j}(1-s_{j})}\psi_{j}$
$\displaystyle=0$ for $j=3$ (55b)

where $(m,\lambda)$ are spectral parameters.

Again, the equation for $\psi_{3}$ in (55b) can be converted to the trivial
equation (32) using the same transformation as for prolate. This yields a
discrete spectral parameter $m\in\mathbb{Z}$ which represents quantised
angular momentum in the $(x_{3},x_{4})$ plane with solutions $e^{im\phi}$. As
with the prolate system, the separated ODEs can also be obtained by
degenerating the generalised Lamé equation.

######  Lemma 22.

The separated equations for the oblate system (55a) and (55b) smoothly
degenerate from those of the ellipsoidal system.

###### Proof.

The proof is identical to that of prolate but uses the transformation

$\displaystyle(e_{4},s_{3},p_{3})$
$\displaystyle=(e_{3}+\epsilon,e_{3}+\epsilon\tilde{s}_{3},\frac{\tilde{p}_{3}}{\epsilon})$

and taking the limit as $\epsilon\to 0$. For more detail, see [33]. ∎

Equation (55a) is an example of a Heun equation.

######  Lemma 23.

By the change of dependent variable $\psi_{j}=(z-a)^{|m|/2}W_{j}$, (55a) can
be written in the form of (34) where
$\gamma=\delta=\frac{1}{2},\epsilon=1+|m|,q=\frac{1}{4}\left(|m|-\lambda\right)$
and $\alpha,\beta$ as in (36). In particular, we have

$W_{j}^{{}^{\prime\prime}}+\left(\frac{1}{2z}+\frac{1}{2(z-1)}+\frac{1+|m|}{z-a}\right)W_{j}^{{}^{\prime}}+\frac{(-E+|m|(|m|+2))z-(|m|-\lambda)}{4z(z-1)(z-a)}W_{j}.$
(56)

Let $Hp_{d}(z)$ be a Heun polynomial solution to (56) of degree $d$ and let
$z_{1},\dots,z_{d}$ be its roots. Many of the results for prolate apply for
the oblate system which we state in the following Lemma.

######  Lemma 24.

The combined solution to (55a) and (55b) is the product

$\psi_{1}\psi_{2}\psi_{3}=(a-s_{1})^{m/2}(a-s_{2})^{m/2}Hp_{d}(s_{1})Hp_{d}(s_{2})e^{im\phi}.$
(57)

Expressed in the original Cartesian coordinates using (53), (57) is a
homogeneous polynomial of total degree $\tilde{D}=2d+|m|$ given by

$\Phi_{\tilde{D}}=\left(a(a-1)\right)^{m/2}\left(x_{3}+i\text{sign}(m)x_{4}\right)^{|m|}\prod_{k=1}^{d}z_{k}(z_{k}-1)\left(r^{2}-\frac{x_{1}^{2}}{z_{k}}a-\frac{a-1}{z_{k}-1}x_{2}^{2}\right).$
(58)

###### Proof.

The proof follows the same logic as used for the prolate system. To obtain
(58), we note the use of the following identity:

$(\theta-s_{1})(\theta-
s_{2})=\theta(\theta-1)\left(r^{2}-\frac{x_{1}^{2}}{\theta}a-\frac{a-1}{\theta-1}x_{2}^{2}\right),$

which holds since both sides are monic quadratic polynomials in $\theta$ that
vanish at $\theta=s_{1},s_{2}$. Using the definition of oblate coordinates in
(53) and observing that
$e^{im\phi}=\left(\frac{x_{3}+ix_{4}}{\sqrt{x_{3}^{2}+x_{4}^{2}}}\right)^{m}$,
the result follows. ∎

As with the prolate system, we only consider the 4 discrete symmetry classes
about the $x_{1}$ and $x_{2}$ planes. Let $\bm{\mu}=(\mu_{1},\mu_{2})$ where
$\mu_{i}\in\\{0,1\\}$ and $\mu_{1}$ ($\mu_{2}$) being $1$ ($0$) denotes a
solution odd (even) about the $x_{1}$ ($x_{2}$) axis. To investigate these
symmetries, consider the change of dependent variable in (56) given by
$\phi_{1}=s_{1}^{\mu_{1}/2}(a-s_{1})^{\mu_{1}/2}W_{1}$ and
$\phi_{2}=s_{2}^{\mu_{2}/2}(a-s_{2})^{\mu_{2}/2}W_{2}$ which corresponds to
multiplication by $x_{1}^{\mu_{1}}x_{2}^{\mu_{2}}$. This change of variables
leads to a Heun equation (34) with parameters given in Table 6. As in the
prolate example, denote by $Hp_{d}^{\bm{\mu}}(z)$ solutions to (34) for the
symmetry class $\bm{\mu}$ and let the product
$Hp_{d}^{\bm{\mu}}(s_{1})Hp_{d}^{\bm{\mu}}(s_{2})e^{im\phi}$, when converted
back into Cartesian coordinates, be given by
$\Phi_{\tilde{D}}^{\bm{\mu}}(\bm{x})$. Set

$\Psi_{D}^{\bm{\mu}}(\bm{x})\coloneqq
x_{1}^{\mu_{1}}x_{2}^{\mu_{2}}\Phi_{\tilde{D}}^{\bm{\mu}}(\bm{x})$ (59)

where $D=\tilde{D}+\sum_{i=1}^{2}\mu_{i}$. Again, (58) is a special case with
$\bm{\mu}=(0,0)$. Following the same argument used in prolate, we have the
following Lemma.

######  Lemma 25.

The $\Psi_{D}^{\bm{\mu}}(\bm{x})$ are $D$ degree homogeneous, harmonic
polynomials and the energy eigenvalue is given by $E=D(D+2)$.

$(\mu_{1},\mu_{2})$ | $\tilde{\alpha}$ | $\tilde{\beta}$ | $\tilde{\gamma}$ | $\tilde{\delta}$ | $\tilde{\epsilon}$ | $\tilde{q}$
---|---|---|---|---|---|---
$(0,0)$ | $\alpha$ | $\beta$ | $\gamma$ | $\delta$ | $\epsilon$ | $q$
$(1,0)$ | $\alpha+1-\gamma$ | $\beta+1-\gamma$ | $2-\gamma$ | $\delta$ | $\epsilon$ | $q+(1-\gamma)(a\delta+\epsilon)$
$(0,1)$ | $\alpha+1-\delta$ | $\beta+1-\delta$ | $\gamma$ | $2-\delta$ | $\epsilon$ | $q+(1-\gamma)(a\delta+\epsilon)$
$(1,1)$ | $\alpha+2-\gamma-\delta$ | $\beta+2-\gamma-\delta$ | $2-\gamma$ | $2-\delta$ | $\epsilon$ | $q+2a(1-\gamma)+2(1-\gamma)\epsilon$
Table 6: Parameters of (34) for the various symmetry classes of the oblate
system.

For a given energy the total number of eigenstates is computed using the same
methodology as for the prolate case.

######  Lemma 26.

For a fixed value of $D$, the total number of eigenstates is given by
$N=(D+1)^{2}$. The number of states for each symmetry class given in Table 6
is the same as that shown in Table 5.

Using the three term recurrence (40) and the same methodology for prolate, we
find the joint spectrum. In Fig. 7 a) we show the joint spectrum
$(m\hbar,\lambda\hbar^{2})$ for $D=20$ and $a=2.4$. As in the prolate system,
we have $\hbar=\frac{1}{D+1}$ and a total of $441$ total eigenstates. The
blue, red, orange and cyan points correspond to the $(0,0),(1,0),(0,1)$ and
$(1,1)$ symmetry classes respectively.

We also show the classical outline of the momentum map in black. Unlike the
prolate system, this is not a semi-toric system. One interesting feature of
the joint spectrum is the existence of period doubling. A quantum consequence
of this is the near-degeneracy in the upper chamber.

Figure 7: a) Joint spectrum of the oblate system $(m\hbar,\lambda\hbar^{2})$
for $D=20$ and $a=2.4$ showing a total of 441 states. b) Corresponding action
map. Note that light blue (cyan) and orange states are hidden by the nearby
degenerate red and blue states in the upper chamber.

To compute the actions, we use the same methodology as for the prolate system.
Doing so gives the following three continuous actions:

$\displaystyle J_{1}$
$\displaystyle=\frac{2}{\pi}\int_{0}^{\min(r_{1},1)}p_{1}ds$ $\displaystyle
J_{2}$
$\displaystyle=\frac{2}{\pi}\int_{\max(r_{1},1)}^{\min(r_{2},a)}p_{2}ds$
$\displaystyle J_{3}$ $\displaystyle=|\ell_{34}|$

where

$p_{i}^{2}=\frac{-\tilde{E}s_{i}^{2}+(\tilde{E}a+g-(a-1)m^{2})s_{i}-ag}{4s_{i}(s_{i}-1)(s_{i}-a)^{2}}$

for $i=1,2$, $r_{2}\geq 1$ and $0\leq r_{1}\leq r_{2}\leq a.$ The action map
corresponding to the joint spectrum in Fig. 7 a) is shown in b). Again, we
note that $J_{1}+J_{2}+J_{3}=1$.

###  3.3 Lamé Coordinates

The Lamé coordinates are unique compared to the other singly degenerate
coordinate systems (oblate, prolate) since they are an extension of
ellipsoidal coordinates from $S^{2}$ onto $S^{3}$. Consequently, while prolate
and oblate coordinates have symmetry group $SO(2)$, the Lamé coordinates have
the larger symmetry group $SO(3)$.

As with oblate, there are two possible (yet equivalent) ways to define Lamé
coordinates. Either the largest three $e_{i}$ become equal, or the smallest
three. We choose the former and define Lamé coordinates by

$\displaystyle x_{1}^{2}$ $\displaystyle=s_{1}$ $\displaystyle x_{2}^{2}$
$\displaystyle=\frac{(1-s_{1})(s_{2}-f_{1})(s_{3}-f_{1})}{(f_{2}-f_{1})(f_{3}-f_{1})}$
$\displaystyle x_{3}^{2}$
$\displaystyle=\frac{(s_{1}-1)(f_{2}-s_{2})(f_{2}-s_{3})}{(f_{1}-f_{2})(f_{2}-f_{3})}$
$\displaystyle x_{4}^{2}$
$\displaystyle=\frac{(s_{1}-1)(f_{3}-s_{2})(f_{3}-s_{3})}{(f_{2}-f_{3})(f_{3}-f_{1})}$
(60)

where $0\leq s_{1}\leq 1$ and $0\leq f_{1}\leq s_{2}\leq f_{2}\leq s_{3}\leq
f_{3}$ with real parameters $f_{i}$.

A possible Stäckel matrix is given by

$\sigma_{\text{Lamé}}=\frac{1}{4}\begin{pmatrix}\frac{1}{s_{1}(1-s_{1})}&\frac{-1}{s_{1}(1-s_{1})^{2}}&0\\\
0&\frac{1}{(f_{3}-s_{2})(s_{2}-f_{2})}&\frac{1}{(f_{3}-s_{2})(s_{2}-f_{2})(s_{2}-f_{1})}\\\
0&\frac{1}{(f_{3}-s_{3})(s_{3}-f_{2})}&\frac{1}{(f_{3}-s_{3})(s_{3}-f_{2})(s_{3}-f_{1})}\end{pmatrix}$
(61)

with corresponding integrals obtained by separation as

$\displaystyle F$ $\displaystyle=\ell_{12}^{2}+\ell_{13}^{2}+\ell_{14}^{2},$
$\displaystyle G$
$\displaystyle=f_{1}\ell_{34}^{2}+f_{2}\ell_{24}^{2}+f_{3}\ell_{23}^{2}.$

As in the other coordinates, reduction by $H$ gives a two degree of freedom
integrable system $(F,G)$ on $S^{2}\times S^{2}$ with Poisson structure
$B_{\bm{X},\bm{Y}}$. The integral $F$ is a result of the $SO(3)$ symmetry.
Using (61) and the same methodology used in the previous systems, we obtain
the following separated ODEs

$\displaystyle\psi_{1}^{{}^{\prime\prime}}+\frac{1}{2}\left(\frac{1}{s_{1}}+\frac{3}{s_{1}-1}\right)\psi_{1}^{{}^{\prime}}+\frac{f-Es_{1}}{4s_{1}(1-s_{1})^{2}}\psi_{1}$
$\displaystyle=0$ (62a)
$\displaystyle\psi_{k}^{{}^{\prime\prime}}+\frac{1}{2}\left(\frac{1}{z-f_{1}}+\frac{1}{z-f_{2}}+\frac{1}{z-f_{3}}\right)\psi_{k}^{{}^{\prime}}+\frac{(f-E)z+g}{4(z-f_{1})(z-f_{2})(z-f_{3})}\psi_{k}$
$\displaystyle=0$ $\displaystyle k=2,3$ (62b)

The separated equations can also be obtained by smoothly degenerating the
generalised Lamé equation.

######  Lemma 27.

Equations (62a) and (62b) arise from smoothly degenerating the separated
generalised Lamé equations in (7).

###### Proof.

From [33], it is known that the transformation

$\displaystyle\left(e_{1},e_{2},e_{3},e_{4}\right)$
$\displaystyle=\left(-\frac{1}{\epsilon},f_{1},f_{2},f_{3}\right)$ (63)
$\displaystyle\left(s_{1},p_{1}\right)$
$\displaystyle=\left(f_{2}-\frac{s_{1}}{\epsilon},\epsilon p_{1}\right)$ (64)

defines a canonical transformation that smoothly degenerates the classical
ellipsoidal system into that for Lamé. Further, applying (63) and taking the
limit as $\epsilon\to 0$ gives

$(\eta_{1},\eta_{2})=\left(-\frac{2H-F}{\epsilon},-\frac{G}{\epsilon}\right).$
(65)

Substituting (63), (64) and (65) into (7) and and considering dominating terms
of $\epsilon$ gives the degenerated ODEs (63) and (64). ∎

Equation (62a) has regular singularities at $s_{1}=0,1,\infty$ meaning that it
is of the hypergeometric type. We observe that the form as written in (62a) is
not canonical. The roots of the indicial equation at $s_{1}=1$ are

$r_{\pm}=\frac{1}{4}\left(-1\pm\sqrt{1+4E-4f}\right),$ (66)

where $r_{-}<0<r_{+}$. To change (62a) into canonical form, we make the change
of dependent variable $\psi_{1}\to y(1-s_{1})^{r_{+}}$ followed by the change
of independent variable $s_{1}=x_{1}^{2}$, yielding the well known Gegenbauer
equation

$(1-x_{1}^{2})y^{{}^{\prime\prime}}-(2u+1)x_{1}y^{{}^{\prime}}+n(n+2u)y=0$
(67)

where $u=\frac{1}{2}(1+\sqrt{1+4E-4f})$ and

$n=-\frac{1}{2}\left(\sqrt{4E-4f+1}-\sqrt{4E+4}+1\right).$ (68)

From [1, 30], it is known that polynomial solutions to (67) occur when $n$ is
an integer. These are known as the Gegenbauer polynomials $C_{n}^{u}(x_{1})$
where $n$ denotes the degree of the polynomial.

The equation in the $s_{2},s_{3}$ variables is known as the Lamé equation, a
special case of the Heun equation. As with the previous degenerate cases, we
simplify by mapping $(f_{1},f_{2},f_{3})$ to $(0,1,a)$ where $a>1$. The
parameters $(f_{1},f_{2},f_{3})$ in (62b) can be mapped to $(0,1,a)$ by a
Möbius transformation, then the Riemann symbol is given by

$\mathcal{S}_{2,3}=\left(\begin{aligned} 0&&1&&a&&\infty\\\
0&&0&&0&&\alpha&;z\\\ 1/2&&1/2&&1/2&&\beta\end{aligned}\right)$

where $\alpha,\beta=\frac{1}{4}\left(1\pm\sqrt{1+4E-4f}\right)$ and $q=g/4$.
The Lamé equation is typically written as follows [7]:

$W^{{}^{\prime\prime}}+\frac{1}{2}\left(\frac{1}{\xi}+\frac{1}{\xi-1}+\frac{1}{\xi-k^{-2}}\right)W^{{}^{\prime}}+\frac{bk^{-2}-\nu(\nu+1)\xi}{4\xi(\xi-1)(\xi-k^{-2})}W=0$
(69)

where $k=\frac{1}{\sqrt{a}}$, $b=g/a$ and
$\nu=\frac{1}{2}\left(\sqrt{4E-4f+1}-1\right)=2\beta-1$ with $\nu\geq-1/2$ and
$0<k<1$. Polynomial solutions of the Lamé equation are known as Lamé
polynomials [5] and we denote them as $Lp_{d}(z)$ where the subscript $d$
signifies the solution’s degree. Note that $\nu$ being a non negative integer
is a necessary condition for polynomial solutions.

######  Lemma 28.

Solutions to (62b) are polynomial iff $E=\tilde{D}(\tilde{D}+2)$ where
$\tilde{D}=2d+n$ is an integer.

###### Proof.

From [45] we know that polynomial solutions are only obtainable for (62b) if

$\frac{f-E}{4}=-d(d+\frac{1}{2}).$ (70)

Combining (70) and (68) gives the result. ∎

Using (70), we know that $r_{+}$ in (66) simplifies to $d$,
$\nu=2d,\alpha=-d,\beta=2d+\frac{1}{2}$ and $u=2d+1$. Denote the full solution
$\psi=\psi_{1}\psi_{2}\psi_{3}$ to (62a) and (62b) as

$\psi=(1-s_{1})^{d}C_{n}^{2d+1}(x_{1})Lp_{d}(s_{2})Lp_{d}(s_{3}).$ (71)

We have the following Lemma.

######  Lemma 29.

Expressed in the original Cartesian coordinates, $\psi$ is a homogeneous
polynomial of degree $\tilde{D}$. Specifically, we can rewrite (71) as

$\Phi_{\tilde{D}}(\bm{x})=r^{n}C_{n}^{2d+1}\left(\frac{x_{1}}{r}\right)\prod_{i=1}^{3}\prod_{k=1}^{d}(z_{k}-f_{i})\prod_{k=1}^{d}\sum_{i=1}^{3}\frac{x_{i+1}^{2}}{z_{k}-f_{i}}$
(72)

where the $z_{k}$ are the roots of the Lamé polynomials $Lp_{d}$,
$(f_{1},f_{2},f_{3})=(0,1,a)$ and $r^{2}=\sum_{i}x_{i}^{2}$.

###### Proof.

From the definition of the Lamé coordinates (60), it is easily seen that

$(z_{k}-s_{2})(z_{k}-s_{3})=\frac{(z_{k}-f_{1})(z_{k}-f_{2})(z_{k}-f_{3})}{1-x_{1}^{2}}\left(\frac{x_{2}^{2}}{z_{k}-f_{1}}+\frac{x_{3}^{2}}{z_{k}-f_{2}}+\frac{x_{4}^{2}}{z_{k}-f_{3}}\right).$

Taking the product over all $d$ roots and using $r^{2}=1$ along with
$s_{1}=x_{1}^{2}$ gives (72). To prove homogeneity, we note that

$C^{2d+1}_{n}(x_{1})=r^{n}C^{2d+1}_{n}(\frac{x_{1}}{r})$ (73)

and so the argument of the Gegenbauer polynomial will be of degree 0. We
observe that $r^{n}$ will be of degree $n$ in the $x_{i}$ and
$C^{2d+1}_{n}(\frac{x_{1}}{r})$ will only consist of terms whose powers have
the same parity as $n$. Thus, we have shown that (72) is homogeneous and
polynomial. ∎

As with the previous systems, let $\bm{\mu}=(\mu_{2},\mu_{3},\mu_{4})$ where
$\mu_{i}\in\\{0,1\\}$ denote a symmetry class where $\mu_{i}=0$ represent the
wavefunction odd about the $x_{i}$ axis and even otherwise. Consider the
change of dependent variable

$\displaystyle\phi_{j}$
$\displaystyle=\prod_{i=1}^{3}(s_{j}-f_{i})^{\mu_{i+1}/2}\psi_{j}$ (74)

where $j=2,3$ in (62b). This is equivalent to multiplying the original
solution by $x_{2}^{\mu_{2}}x_{3}^{\mu_{3}}x_{4}^{\mu_{4}}$. Note that
performing this change of variables in the Lamé equation (62b) results in a
Heun equation of the form (34) with parameters for each symmetry class
$\bm{\mu}$ as shown in Table 7.

For the Gegenbauer equation (67), we note that the independent variable is
$x_{1}$ and so parity about this axis is determined simply by the parity of
the Gegenbauer polynomial $C_{n}^{u}(x_{1})$, which is determined by $n$. We
let $\mu_{1}=1$ if $n$ is odd and $\mu_{1}=0$ otherwise.

For a given symmetry class $\bm{\mu}$, denote by $Lp_{d}^{\bm{\mu}}(z)$
solutions to (34) with parameters given in Table 7. We consider the product
$\Phi_{\tilde{D}}^{\bm{\mu}}(\bm{x})=C_{n}^{2d+1+U}(x_{1})Lp_{d}^{\bm{\mu}}(s_{2})Lp_{d}^{\bm{\mu}}(s_{3})$
where $U=\sum_{i=1}^{3}\mu_{i+1}$. Finally, we set

$\Psi_{D}^{\bm{\mu}}(\bm{x})\coloneqq
x_{2}^{\mu_{2}}x_{3}^{\mu_{3}}x_{4}^{\mu_{4}}\Phi_{\tilde{D}}^{\bm{\mu}}(\bm{x})$
(75)

where $D=n+\tilde{D}+U$. As before, (72) is a special case of (75) with
$\bm{\mu}=(0,0,0)$.

######  Lemma 30.

The $\Psi_{D}^{\bm{\mu}}(\bm{x})$ are homogeneous polynomial eigenfunctions of
the original Laplacian (5) and the energy is given by $E=D(D+2)$.

###### Proof.

Proving $\Psi_{D}^{\bm{\mu}}(\bm{x})$ is an eigenfunction follows the same
argument as the previous cases. To compute the energy for a given symmetry
class $\bm{\mu}$, we combine the parameters of Table 7 with the results of
[45] to get a more general form of (70), namely

$\frac{f-E+U(U+1)}{4}=-d(d+1/2+U).$ (76)

Combining (76) and (68) gives the result for $E$. Similarly, we have that
$u=2d+1+U$. ∎

$\bm{\mu}$ |  $\tilde{\alpha}$ |  $\tilde{\beta}$ |  $\tilde{\gamma}$ |  $\tilde{\delta}$ |  $\tilde{\epsilon}$ |  $\tilde{q}$
---|---|---|---|---|---|---
$(0,0,0)$ |  $\alpha$ |  $\beta$ |  $\gamma$ |  $\delta$ |  $\epsilon$ |  $q$
$(1,0,0)$ |  $\alpha+1-\gamma$ |  $\beta+1-\gamma$ |  $2-\gamma$ |  $\delta$ |  $\epsilon$ |  $q+(1-\gamma)(a\delta+\epsilon)$
$(0,1,0)$ |  $\alpha+1-\delta$ |  $\beta+1-\delta$ |  $\gamma$ |  $2-\delta$ |  $\epsilon$ |  $q-a\gamma(\delta-1)$
$(0,0,1)$ |  $\alpha+1-\epsilon$ |  $\beta+1-\epsilon$ |  $\gamma$ |  $\delta$ |  $2-\epsilon$ |  $q+\gamma(1-\epsilon)$
$(1,1,0)$ |  $\alpha-\gamma-\delta+2$ |  $\beta-\gamma-\delta+2$ |  $2-\gamma$ |  $2-\delta$ |  $\epsilon$ |  $q-a(\gamma+\delta-2)-\gamma\epsilon+\epsilon$
$(1,0,1)$ |  $\alpha-\gamma-\epsilon+2$ |  $\beta-\gamma-\epsilon+2$ |  $2-\gamma$ |  $\delta$ |  $2-\epsilon$ |  $q-\gamma(a\delta+1)+a\delta-\epsilon+2$
$(0,1,1)$ |  $\alpha-\delta-\epsilon+2$ |  $\beta-\delta-\epsilon+2$ |  $\gamma$ |  $2-\delta$ |  $2-\epsilon$ |  $q+\gamma(-a\delta+a-\epsilon+1)$
$(1,1,1)$ |  $\alpha-\gamma-\delta-\epsilon+3$ |  $\beta-\gamma-\delta-\epsilon+3$ |  $2-\gamma$ |  $2-\delta$ |  $2-\epsilon$ |  $q-a(\gamma+\delta-2)-\gamma-\epsilon+2$
Table 7: Symmetry class $\bm{\mu}$ and corresponding Heun (Lamé) equation
parameters for which the resulting polynomial is an eigenfunction.

From the above, it is clear that for a fixed energy $E$, only 8 of the 16
discrete symmetry classes for the Lamé system can be present. If $E$ is even,
then the symmetry classes $\mathcal{S}_{E}$ given in (21) are present.
Conversely, if $E$ is odd, we have the remaining 8 symmetry classes
$\mathcal{S}_{O}$ from (22). Note that while we have omitted $\mu_{1}$ from
the definition of $\bm{\mu}$ in our discussion of the Lamé systems so far,
this symmetry is accounted for in the parity of $n$ which is set by the value
of $D$. When addressing all symmetry classes below, we will write
$\tilde{\bm{\mu}}=(\mu_{1};\mu_{2},\mu_{3},\mu_{4})$ to highlight this fact.

We have the following observation.

######  Lemma 31.

Chebyshev polynomials of the second kind are the eigenfunction corresponding
to the joint spectrum point $(f,g)=(1-\hbar^{2},0)$.

###### Proof.

For $\bm{\mu}=(0,0,0)$, we have from (76) the relation $f=E-4d(d+1/2)$. When
$d=0$ one obtains $f=E$ which gives $f=1-\hbar^{2}$ (after normalisation by
$\hbar$) and $g=0$. This gives Chebyshev polynomials of the second kind
$U_{n}(x_{1})$ since $U_{n}(x_{1})=C_{n}^{1}(x_{1}).$ Note that in the limit
$\hbar\to 0$ this joint spectrum point becomes the degenerate point $(1,0)$
for the classical system. ∎

The Lamé equation (69) can be transformed by the change of independent
variable $\xi=\text{sn}^{2}(z,k)$ into the following form:

$W^{{}^{\prime\prime}}+(b-\nu(\nu+1)k^{2}\text{sn}^{2}(z,k))W=0.$ (77)

This form of the Lamé equation has regular singular points at
$2pK+(2q+1)iK^{{}^{\prime}}$ where $K,K^{{}^{\prime}}$ are the complete
elliptic integrals of the first kind with moduli $k$ and $\sqrt{1-k^{2}}$
respectively and $p,q\in\mathbb{Z}$.

From [5], it is known that there are $8$ polynomial solutions to (77) if $\nu$
is a non-negative integer. These are shown in Table 8 where $m=0,1,\dots,\nu$
and $Ec^{m}_{d},Es^{m}_{d}$ are known as the Lamé functions [20]. Since
$\xi=s_{i}=\text{sn}^{2}(z,k)$, we are able to connect these eigenfunctions to
$Lp_{d}^{\bm{\mu}}(z)$ as given in Table 7.

The total number of states for each symmetry class is shown in Table 9 below.
Here we use the notation $(\mu_{1};\mu)$ to represent a state with parity
given by $\mu_{1}$ about the $x_{1}$ axis and odd about one remaining axis
given by $\mu$. E.g. $(1;\mu)$ represents the classes
$(1;0,0,1),(1;0,1,0),(1;1,0,0)$. Similarly, we have $(\mu_{1};\mu,\nu)$ and so
forth. The number of states are computed using similar methods to those shown
in the prolate case.

Summing the relevant entries for a given parity of the energy shows that the
total number of states for a given energy $E=D(D+2)$ is given by $(D+1)^{2}$.
To compute the joint spectrum for a given $E$, we use the parameters shown in
Table 7 and the recurrence relationship given in Lemma 15. An example for even
and odd energies is shown in Fig. 8. As with the other cases such as
ellipsoidal and prolate, we have $\hbar=\frac{1}{D+1}$. As with the prolate
and oblate systems, since the matrices involved are not self adjoint, larger
values of $D$ lead to numerical instability.

Lamé Polynomial | Corresponding Symmetry $\bm{\mu}$
---|---
$Ec_{2d}^{2m}(z,k^{2})$ | $(0,0,0)$
$Ec_{2d+1}^{2m+1}(z,k^{2})$ | $(1,0,0)$
$Es_{2d+1}^{2m+1}(z,k^{2})$ | $(0,1,0)$
$Ec_{2d+1}^{2m}(z,k^{2})$ | $(0,0,1)$
$Es_{2d+1}^{2m+2}(z,k^{2})$ | $(1,1,0)$
$Ec_{2d+2}^{2m+1}(z,k^{2})$ | $(1,0,1)$
$Es_{2d+2}^{2m+1}(z,k^{2})$ | $(0,1,1)$
$Es_{2d+3}^{2m+2}(z,k^{2})$ | $(1,1,1)$
Table 8: Lamé polynomials solutions to (77) and their corresponding symmetries given in Table 7. | Number of States | $E$
---|---|---
$(1;\mu)$ | $D(D+2)/8$ | Even
$(1;1,1,1)$ | $D(D-2)/8$ | Even
$(0;0,0,0)$ | $(D+2)(D+4)/8$ | Even
$(0;\mu,\nu)$ | $D(D+2)/8$ | Even
$(1;0,0,0)$ | $(D+1)(D+3)/8$ | Odd
$(1;\mu,\nu)$ | $(D+1)(D-1)/8$ | Odd
$(0;\mu)$ | $(D+1)(D+3)/8$ | Odd
$(0;1,1,1)$ | $(D+1)(D-1)/8$ | Odd
Table 9: Number of states per symmetry class and corresponding energy parity.

Figure 8: Joint spectrum $(f,g)$ in the Lamé case with
$(f_{1},f_{2},f_{3})=(0,1,2.4)$ for a) even energy with $D=20$ and b) odd
energy with $D=21$ thereby representing all $16$ symmetry classes.
Correspondence between the coloured dots and symmetry class are the same for
the ellipsoidal system shown in Table 3.

The actions of the Lamé system are given by the following

$\displaystyle J_{1}=\frac{2}{\pi}\int_{0}^{f}p_{1}ds,$ $\displaystyle
J_{2}=\frac{2}{\pi}\int_{f_{1}}^{\min(r_{2},f_{2})}p_{2}ds,$ $\displaystyle
J_{3}=\frac{2}{\pi}\int_{\max(f_{2},r_{2})}^{f_{3}}p_{3}ds$

where $r_{2}=\frac{g}{1-f}$ and

$\displaystyle p_{1}^{2}=\frac{f-\tilde{E}s_{1}}{4(s_{1}-1)^{2}s_{1}}$
$\displaystyle
p_{i}^{2}=\frac{(f-\tilde{E})s_{i}+g}{4(s_{i}-f_{3})(s_{i}-f_{2})(s_{i}-f_{1})}$

for $i=2,3.$ The first action simplifies to

$J_{1}=1-\sqrt{1-f}.$ (78)

The image of spectra shown in Fig. 8, under the action map is shown in Fig. 9.
Unlike the prolate and oblate cases, there are no states on the boundary.
There is a single state located near the top corner $(f,g)=(1-\hbar^{2},0)$.

Figure 9: Joint spectrum in the Lamé case in action variables corresponding to
the joint spectra shown in Fig. 8 a) and b).

###  3.4 Spherical Coordinates

Spherical coordinates (also known as poly-spherical coordinates) are a further
degeneration of prolate, oblate or Lamé coordinates and possess both an
$SO(2)$ and $SO(3)$ symmetry. Unlike prolate coordinates, spherical
coordinates - despite having a global $S^{1}$ action - are not semi-toric
since they have a degenerate point. This point corresponds to the critical
values at which the global action is not differentiable, and the pre-image of
the critical value is a sphere. We will see that there is a simple
eigenfunction built from Chebyshev polynomials that is related to this sphere.

Spherical coordinates can be obtained from the Lamé case by setting either
$f_{1}=f_{2}$ or $f_{2}=f_{3}$. These result in different coordinate systems,
respectively called the $12$ and $23$ spherical systems and are related to
each other via a permutation of coordinates. Similarly, setting $a=1$ in both
the prolate and oblate cases gives the $23$ and $12$ spherical systems
respectively. See [33] for more details.

In this paper we focus on the $23$ spherical system, defined by

$\displaystyle x_{1}^{2}$ $\displaystyle=s_{1}$ $\displaystyle x_{2}^{2}$
$\displaystyle=(1-s_{1})s_{2}$ $\displaystyle x_{3}^{2}$
$\displaystyle=(1-s_{1})(1-s_{2})s_{3}$ $\displaystyle x_{4}^{2}$
$\displaystyle=(1-s_{1})(1-s_{2})(1-s_{3}),$ (79)

where $0\leq s_{k}\leq 1$ for all $k=1,2,3$. Unlike the previously described
systems, the spherical coordinates are relatively simple and separation can be
done "by hand", without the use of a Stäckel matrix. The classical integrals
obtained by separation are $(F,\ell_{34}^{2})$ where
$F=\ell_{12}^{2}+\ell_{13}^{2}+\ell_{14}^{2}$. Note this is the same $F$ as in
the Lamé system. The separated ODEs are

$\displaystyle\psi_{1}^{{}^{\prime\prime}}+\frac{1}{2}\left(\frac{3}{s_{1}-1}+\frac{1}{s_{1}}\right)\psi_{1}^{{}^{\prime}}+\frac{\lambda-
Es_{1}}{4s_{1}(1-s_{1})^{2}}\psi_{1}$ $\displaystyle=0$ (80a)
$\displaystyle\psi_{2}^{{}^{\prime\prime}}+\left(\frac{1}{2s_{2}}+\frac{1}{s_{2}-1}\right)\psi_{2}^{{}^{\prime}}+\frac{(\lambda-E)(s_{2}-1)-m}{4s_{2}(1-s_{2})^{2}}\psi_{2}$
$\displaystyle=0$ (80b)
$\displaystyle\psi_{3}^{{}^{\prime\prime}}+\frac{1}{2}\left(\frac{1}{s_{3}}+\frac{1}{s_{3}-1}\right)\psi_{3}^{{}^{\prime}}+\frac{m}{4s_{3}(1-s_{3})}\psi_{3}$
$\displaystyle=0.$ (80c)

The equation for $\psi_{3}$ (80c) can be converted to the trivial equation
(32) with the coordinate transformation $s_{3}=\cos^{2}(\phi)$. Using this and
periodic boundary conditions, we find that $m\in\mathbb{Z}$ and
$\psi_{3}=e^{im\phi}$ where $\phi=\arctan(\frac{x_{3}}{x_{4}})$.

As for the previous degenerate systems, we have the following Lemma.

######  Lemma 32.

The system of separated equations in (80a), (80b), (80c) can be obtained from
either those of the Lamé system (62a), (62b) or those of the oblate system
(55a), (55b).

###### Proof.

For the Lamé system (62a), (62b), the transformation
$(s_{3},p_{3})\to(f_{2}+\tilde{s_{3}}\epsilon,\tilde{p_{3}}/\epsilon)$ is
canonical. This, coupled with $f_{3}\to f_{2}+\epsilon$ and a corresponding
transformation of the integrals followed by the normalization
$(f_{1},f_{2})=(0,1)$ yields the result. For the oblate direction, a similar
method with
$(a,s_{3},p_{3})=(1+\epsilon,1+\epsilon\tilde{s_{3}},\tilde{p_{3}}/\epsilon)$
in (55a) and (55b) gives the result. ∎

Next, we turn our attention to the non trivial equations. For $\psi_{1}$ in
(80a), we observe that it is of hypergeometric type. It has 3 regular
singularities - two finite poles at $s_{1}=0,1$ and the third at infinity. To
assist in finding polynomial solutions, we transform the equation for
$\psi_{1}$ using the change of dependent variable
$\psi_{1}=(1-s_{1})^{\frac{1}{4}(-1+X)}W_{1}$ where $X=\sqrt{1+4E-4\lambda}$,
followed by the change of independent variable $s_{1}=x_{1}^{2}$. This yields
the Gegenbauer equation

$(1-x_{1}^{2})W_{1}^{{}^{\prime\prime}}-2x_{1}(2u+1)W_{1}^{{}^{\prime}}+n(n+2u)W_{1}=0$
(81)

where $u=\frac{1}{2}(1+X)$ and

$n=-\frac{1}{2}(1-2\sqrt{1+E}+X).$ (82)

Gegenbauer polynomial solutions to (81), denoted by $C_{n}^{u}(x_{1})$, are
obtained when $n_{1}$ is a non negative integer.

Repeating a similar process for $\psi_{2}$ in (80b), we again note the
equation is hypergeometric with regular singularities at $0,1$ and infinity.
The change of independent variable
$s_{2}=x^{2}\coloneqq\frac{x_{2}^{2}}{\tilde{r}^{2}}$ where
$\tilde{r}\coloneqq x_{2}^{2}+x_{3}^{2}+x_{4}^{2}$ gives the Associated
Legendre equation written as follows

$(1-x^{2})\psi_{2}^{{}^{\prime\prime}}-2x\psi_{2}^{{}^{\prime}}+\left(\ell(\ell+1)-\frac{m^{2}}{1-x^{2}}\right)\psi_{2}=0$
(83)

where $\ell(\ell+1)=E-\lambda$, i.e.

$\ell=\frac{1}{2}(-1+\sqrt{1+4E-4\lambda}).$ (84)

In (83), $\ell$ and $m$ are referred to as the degree and order respectively
of the Associated Legendre equation. Non trivial non singular solutions
$P_{\ell}^{m}(x)$ are yielded when both $\ell$ and $m$ are integers. Further,
said solutions are polynomial only if $m$ is an even integer. For polynomial
solutions, both $n$ and $\ell$ from (82) and (84) respectively must be
integers. Solving these simultaneously for $E$ gives

$E=D(D+2)$ (85)

where $D=\ell+n$. Using (85) with (84) gives analytic values of the eigenvalue
$\lambda$ as

$\lambda=D(D+2)-\ell(\ell+1)$ (86)

Using (86), we also have

$X=2\ell+1$ (87)

along with $u=\ell+1$. Combining the above, we have the following Lemma.

######  Lemma 33.

For a fixed value of $D=n+\ell$, the combined wavefunction solution to (5),
expressed in Cartesian coordinates, is the homogeneous polynomial given by the
product

$\Psi_{D}(\bm{x})=r^{n}\tilde{r}^{\ell}\left(\frac{x_{3}+ix_{4}}{\sqrt{x_{3}^{2}+x_{4}^{2}}}\right)^{m}C_{n}^{\ell+1}\left(\frac{x_{1}}{r}\right)P_{\ell}^{m}\left(\frac{x_{2}}{\tilde{r}}\right)$
(88)

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

###### Proof.

First, we note as before
$e^{im\phi}=(\frac{x_{3}+ix_{4}}{\sqrt{x_{3}^{2}+x_{4}^{2}}})^{m}$. The
associated Legendre polynomials $P_{\ell}^{m}(x)$ can be written in the form
(see, e.g., [16, 14.7])

$P^{m}_{\ell}(x)=\frac{(-1)^{m}}{2^{\ell}\ell!}(1-x^{2})^{m/2}\frac{d^{\ell+m}}{dx^{\ell+m}}(x^{2}-1)^{\ell}\,.$
(89)

Thus the denominator $\sqrt{x_{3}^{2}+x_{4}^{2}}^{m}$ in the Lemma cancels
with $\tilde{r}^{m}(1-x_{2}^{2}/\tilde{r}^{2})^{m/2}$. The remaining terms
produce a homogeneous degree $\ell-m$ polynomial in $x_{2}$ and
$\tilde{r}^{2}-x_{2}^{2}=x_{3}^{2}+x_{4}^{2}$. We note that
$r^{n}C_{n}^{\ell+1}(\frac{x_{1}}{r})$ is homogeneous of degree $n$ and so it
follows that (88) is homogeneous of degree $D=n+\ell$. ∎

Let $\bm{\mu}=(\mu_{1},\mu_{2})$ where $\mu_{i}\in\\{0,1\\}$ and $\mu_{1}$
$(\mu_{2})$ being $1$ $(0)$ denotes a solution odd (even) about the $x_{1}$
$(x_{2})$ axis. As with the Lamé system, we note that both the Gegenbauer and
Associated Legendre equations have independent variables $x_{1}$ and $x_{2}$
respectively and thus their parity about these axes will be determined by that
of $n$ and $\ell-m$ respectively. Consequently, a fixed value of $D$ enforces
a relationship between $n,\ell$ and $m$.

We note that $\ell=0$ forces $D=n$ and $m=0$. Consequently, we recover
Chebyshev polynomials of the second kind as $C_{n}^{1}(x_{1})=U_{n}(x_{1})$.
The corresponding state on the joint spectrum is
$(m,\lambda)=(0,1-\hbar^{2})$, which becomes the classical degenerate point in
the limit of $\hbar\to 0$. Depending on the parity of $E$, this corresponds to
either the $(0,0)$ symmetry (if $E$ is even) or $(1,0)$ if $E$ is odd.

Fixing an energy $E$ allows us to compute the joint spectrum for the spherical
system. An example of the joint spectrum $(m,\lambda)$ is shown in Fig. 10 a).
The degenerate point is shown in magenta at $(0,1)$.

From [33], it is known that the actions of the spherical system are given as
follows:

$\displaystyle J_{1}=1-\sqrt{1-f}$ $\displaystyle J_{2}=\sqrt{1-f}-m$
$\displaystyle J_{3}=m$ (90)

where $f$ is a value of the integral $F$. In Fig. 10 b) we show the
corresponding action map to the joint spectrum shown in a).

Figure 10: a) The joint spectrum in the spherical case with $D=30$. Blue, red,
orange and cyan represent the $(0,0),(1,0),(0,1)$ and $(1,1)$ symmetries
respectively. There are $31^{2}=961$ states in total. b) The corresponding
spectrum in action variables (90).

###  3.5 Cylindrical Coordinates

Cylindrical (or Hopf, or doubly cylindrical) coordinates are a further
degeneration of oblate coordinates, found by setting $e_{1}=e_{2}$ in addition
to $e_{3}=e_{4}$. The relationship between Cartesian and cylindrical
coordinates is given by

$\displaystyle x_{1}^{2}$ $\displaystyle=s_{1}s_{2}$ $\displaystyle x_{2}^{2}$
$\displaystyle=s_{2}(1-s_{1})$ (91) $\displaystyle x_{3}^{2}$
$\displaystyle=s_{3}(1-s_{2})$ $\displaystyle x_{4}^{2}$
$\displaystyle=(1-s_{2})(1-s_{3}).$

The separated equations are given by

$\displaystyle\psi_{1}^{{}^{\prime\prime}}+\frac{1}{2}\left(\frac{1}{s_{1}}+\frac{1}{s_{1}-1}\right)\psi_{1}^{{}^{\prime}}+\frac{m_{1}^{2}}{4s_{1}(1-s_{1})}\psi_{1}$
$\displaystyle=0$ (92a)
$\displaystyle\psi_{2}^{{}^{\prime\prime}}+\left(\frac{1}{s_{2}-1}+\frac{1}{s_{2}}\right)\psi_{2}^{{}^{\prime}}+\frac{m_{1}(s_{2}-1)-(m_{2}+E(s_{2}-1))s_{2}}{4s_{2}^{2}(1-s_{2})^{2}}\psi_{2}$
$\displaystyle=0$ (92b)
$\displaystyle\psi_{3}^{{}^{\prime\prime}}+\frac{1}{2}\left(\frac{1}{s_{3}}+\frac{1}{s_{3}-1}\right)\psi_{3}^{{}^{\prime}}+\frac{m_{2}^{2}}{4s_{3}(1-s_{3})}\psi_{3}$
$\displaystyle=0$ (92c)

where $m_{1},m_{2}$ are spectral parameters. The equations in $\psi_{1}$ and
$\psi_{3}$ in (92a) and (92c) are reducible to the trivial equation (32) with
the trigonometric transformation used for the previous systems. This, along
with periodic boundary conditions, give $m_{1},m_{2}\in\mathbb{Z}$ with
solutions $\psi_{1}=e^{im_{1}\phi_{1}}$ and $\psi_{3}=e^{im_{2}\phi_{3}}$
where $\phi_{1}=\arctan(\frac{x_{2}}{x_{1}})$ and
$\phi_{3}=\arctan(\frac{x_{4}}{x_{3}})$.

In (92b), the equation in $\psi_{2}$ is of hypergeometric type. The change of
dependent variable
$\psi_{2}=s_{2}^{\frac{\left|m_{1}\right|}{2}}(1-s_{2})^{\frac{\left|m_{2}\right|}{2}}y$
and independent $s_{2}=\frac{1-x}{2}$ variables yields

$y^{{}^{\prime\prime}}+\left(\frac{1+\left|m_{1}\right|}{x-1}+\frac{1+\left|m_{2}\right|}{x+1}\right)y^{{}^{\prime}}+\frac{-E+2\left|m_{1}\right|\left(\left|m_{2}\right|+1\right)+m_{1}^{2}+m_{2}^{2}+2\left|m_{2}\right|}{4(x^{2}-1)}y=0$

whose Riemann symbol is given by

$\mathcal{S}_{\text{Cyl}}=\left(\begin{aligned} -1&&1&&\infty\\\
0&&0&&\frac{1}{2}\left(-\sqrt{E+1}+\left|m_{1}\right|+\left|m_{2}\right|+1\right)&;x\\\
-\left|m_{2}\right|&&-\left|m_{1}\right|&&\frac{1}{2}\left(\sqrt{E+1}+\left|m_{1}\right|+\left|m_{2}\right|+1\right)\end{aligned}\right).$

Polynomial solutions to $\psi_{2}$ are obtained when an index at infinity
vanishes. This leads to the condition

$\frac{1}{2}\left(-\sqrt{E+1}+\left|m_{1}\right|+\left|m_{2}\right|+1\right)=-d$
(93)

which gives the quantised energy

$E=D(D+2)$

where $D=2d+\left|m_{1}\right|+\left|m_{2}\right|$. These solutions, denoted
by $P_{d}^{(|m_{1}|,|m_{2}|)}(x)$, are the well known Jacobi polynomials.

The combined wavefunction solution to (5) is therefore

$\psi=P_{d}^{(|m_{1}|,|m_{2}|)}(x)e^{im_{1}\phi_{1}}e^{im_{2}\phi_{3}}.$ (94)

######  Lemma 34.

For a given $D=2d+|m_{1}|+|m_{2}|$, the total wave function (94) in the
original Cartesian coordinates is a harmonic homogeneous polynomial of degree
$D$ given by

$\Psi_{D}=r^{2d}(x_{1}+i\operatorname{sign}(m_{1})x_{2})^{|m_{1}|}(x_{3}+i\operatorname{sign}(m_{2})x_{4})^{|m_{2}|}P_{d}^{(\left|m_{1}\right|,\left|m_{2}\right|)}\left(\frac{x_{3}^{2}+x_{4}^{2}-x_{1}^{2}-x_{2}^{2}}{r^{2}}\right).$
(95)

###### Proof.

Using the definition of cylindrical coordinates in (91), we note that
$s_{2}=x_{1}^{2}+x_{2}^{2}$. Since $x=1-2s_{2}$ and $r^{2}=\sum_{i}x_{i}^{2}$,
the argument of the Jacobi polynomial follows. Similarly, rewriting the
exponential terms in (95) using this definition and cancelling terms gives the
final result. To prove homogeneity, we note that the argument of the Jacobi
polynomial is of degree $0$ and so the final degree of $\Psi_{D}$ will be
$|m_{1}|+|m_{2}|+2d=D$. ∎

Note that for $m_{1}=m_{2}=0$ we recover the Legendre polynomials
$P_{d}^{0}(x)=P_{d}^{(0,0)}(x)$. For $m_{1}=m_{2}$ we recover Gegenbauer
polynomials $C_{d}^{l}(x)/P_{d}^{(l+1/2,l+1/2)}=const$, i.e. they are the same
up to normalisation, see, e.g. [43].

Since there are two continuous $S^{1}$ symmetries, the discrete symmetries of
(95) are represented by the parity of $m_{1}$ and $m_{2}$. Note that
$D-2d=|m_{1}|+|m_{2}|$ has the same parity as $D$ and $E$. In particular this
implies that $m_{1}=m_{2}=0$ is only possible for even $D$ and $E$.

An example of the joint spectrum is given in Fig. 11 a) with $D=20$. The total
number of eigenstates is $(D+1)^{2}=441$.

Two of the spectral parameters are trivial angular momenta:
$(J_{1},J_{3})=(\left|m_{1}\right|,\left|m_{2}\right|)$. As
$J_{1}+J_{2}+J_{3}=1$, the remaining action $J_{2}$ is easily found. In Fig.
11 b) an example of the action map for the corresponding spectrum in a) is
shown.

Figure 11: a) Joint spectrum with $D=20$. b) Corresponding action map.

##  Appendix A Appendix

###  A.1 Ellipsoidal Coordinates on $S^{2}$

In these remaining sections we compare our results on $S^{3}$ to the more
familiar $S^{2}$. Ellipsoidal coordinates on $S^{2}$ are defined as follows:

$\displaystyle x_{1}^{2}$
$\displaystyle=\frac{(s_{1}-e_{1})(s_{2}-e_{1})}{(e_{2}-e_{1})(e_{3}-e_{1})}$
$\displaystyle
x_{2}^{2}=\frac{(s_{1}-e_{2})(s_{2}-e_{2})}{(e_{1}-e_{2})(e_{3}-e_{2})}$
$\displaystyle
x_{3}^{2}=\frac{(s_{1}-e_{3})(s_{2}-e_{3})}{(e_{1}-e_{3})(e_{2}-e_{3})}.$ (96)

where $0\leq e_{1}\leq s_{1}\leq e_{2}\leq s_{2}\leq e_{3}$. Separating (5) in
these coordinates gives the following separated equations

$\psi_{i}^{{}^{\prime\prime}}+\frac{1}{2}\left(\frac{1}{s_{i}-e_{1}}+\frac{1}{s_{i}-e_{2}}+\frac{1}{s_{i}-e_{3}}\right)\psi_{i}^{{}^{\prime}}+\frac{-Es_{i}+\lambda}{4(s_{i}-e_{1})(s_{i}-e_{2})(s_{i}-e_{3})}\psi_{i}=0$
(97)

for $i=1,2$. Equations (97) for $i=1,2$ are the Lamé equations studied in
section 3.3 and so we denote their polynomial solutions by $Lp_{d}(s_{i})$.
Classically, the integrals obtained by separation are

$(I_{1},I_{2})=(\ell_{12}^{2}+\ell_{13}^{2}+\ell_{23}^{2},e_{3}\ell_{12}^{2}+e_{2}\ell_{13}^{2}+e_{1}\ell_{23}^{2}).$

A combined solution to (5) in ellipsoidal coordinates is given by

$\psi=Lp_{d}(s_{1})Lp_{d}(s_{2})$ (98)

where there are $d=n_{1}+n_{2}$ total roots of $\psi$ with $n_{1}$ occurring
in the interval $[e_{1},e_{2}]$ and $n_{2}$ in $[e_{2},e_{3}]$.

######  Lemma 35.

Written in Cartesian coordinates, $\psi$ in (98) is a homogeneous
$\tilde{\ell}\coloneqq 2d$ degree polynomial

$\Phi_{\tilde{l}}=\left(\prod_{i=1}^{3}\prod_{k=1}^{d}(z_{k}-e_{i})\right)\prod_{k=1}^{d}\sum_{i=1}^{3}\frac{x_{i}^{2}}{z_{k}-e_{i}}$
(99)

where the $z_{k}$ are solutions of

$\sum_{i=1}^{3}\frac{1/4}{z_{k}-e_{i}}+\sum_{j=1,j\neq
k}^{d}\frac{1}{z_{k}-z_{j}}=0.$

Inspecting (99), it is clear there are 8 symmetry classes in total (even/odd
parity about each of the $x_{i}$ axes). Let
$\bm{\mu}\coloneqq(\mu_{1},\mu_{2},\mu_{3})$ where $\mu_{i}\in\\{0,1\\}$
represent a symmetry class of the solution (99) where $\mu_{i}=1$ signifies a
solution odd about the $x_{i}$ axis and even otherwise. As with the $S^{3}$
ellipsoidal case, we consider these symmetries by performing the change of
variables

$\phi_{j}=\prod_{i=1}^{3}(s_{j}-e_{i})^{\mu_{i}/2}\psi_{j}$ (100)

in (97). Let solutions to the transformed (97) for a given symmetry class
$\bm{\mu}$ be given by $\Phi_{\tilde{l}}^{\bm{\mu}}(\bm{x})$ and set

$\Psi_{\ell}^{\bm{\mu}}(\bm{x})\coloneqq
x_{1}^{\mu_{1}}x_{2}^{\mu_{2}}x_{3}^{\mu_{3}}\Phi_{\tilde{l}}^{\bm{\mu}}(\bm{x})$
(101)

We then have the following Lemma.

######  Lemma 36.

The $\Psi_{\ell}^{\bm{\mu}}(\bm{x})$ are eigenfunctions of the Schrödinger
equation (5) and the energy is given by $E=\ell(\ell+1)$ where $\ell\coloneqq
2d+\sum_{i=1}^{3}\mu_{i}$.

Similar to the $S^{3}$ ellipsoidal and Lamé system, for a fixed value of
$\ell$, only four of these symmetries appear at any time. When $\ell$ is even,
these are the $(0,0,0),(1,1,0),(1,0,1),(0,1,1)$ classes (blue, purple, orange,
gray respectively in Fig. 12 a)) while for odd $\ell$ we have the
$(1,0,0),(0,1,0),(0,0,1),(1,1,1)$ classes (red, green, cyan, brown). Note that
$E$ is always even.

We demonstrate an example of this in Fig. 12 a) and b) with $\ell=20$ and
$\ell=19$ respectively and $(e_{1},e_{2},e_{3})=(0,1,2.4)$. The joint spectrum
is obtained using the same three term recurrence relations from Lemma 15 used
to find the prolate, oblate and Lamé spectra.

Figure 12: Example spectra with a) $\ell=20$ and b) $\ell=19$ respectively and
$(e_{1},e_{2},e_{3})=(0,1,2.4)$.

For ellipsoidal coordinates on the $2-$sphere, consider a fixed value of the
energy on $(E_{2}),$ the corresponding eigenvalues $\lambda$, as well as an
energy on the $3-$sphere $(E_{3})$ and eigenvalues for the Lamé system
$(f,g)$. Comparing (97) and (62a) for $k=1,2$, it is clear that

$\displaystyle E_{2}=f-E_{3}$ $\displaystyle\lambda$ $\displaystyle=g.$

This can be geometrically interpreted as follows: each joint spectrum for the
ellipsoidal $S^{2}$ system (such as those shown in Fig. 12) corresponds to a
horizontal slice of the Lamé joint spectrum (see Fig. 8), where the slice is
specified for a given combination of $(f,E_{3})$.

###  A.2 Spherical Coordinates on $S^{2}$

Spherical coordinates on $S^{2}$ are defined as follows:

$\displaystyle x_{1}^{2}=s_{1}$ $\displaystyle x_{2}^{2}=(1-s_{1})s_{2}$
$\displaystyle x_{3}^{2}=(1-s_{1})(1-s_{2})$ (102)

where $0\leq s_{1}\leq s_{2}\leq 1$. Separating the Hamilton-Jacobi equation
gives classical integrals
$(I_{1},I_{2})=(\ell_{12}^{2}+\ell_{13}^{2}+\ell_{23}^{2},\ell_{12}^{2})$
while doing so for the Schrödinger equation (5) yields the following ODEs

$\displaystyle\psi_{1}^{{}^{\prime\prime}}+\frac{1}{2}(\frac{1}{s_{1}}+\frac{2}{s_{1}-1})\psi_{1}^{{}^{\prime}}+\frac{-Es_{1}+(E-m^{2})}{4s_{1}(s_{1}-1)^{2}}\psi_{1}$
$\displaystyle=0$ (103a)
$\displaystyle\psi_{2}^{{}^{\prime\prime}}+\frac{1}{2}(\frac{1}{s_{2}}+\frac{1}{s_{2}-1})\psi_{2}^{{}^{\prime}}+\frac{m}{4s_{2}(s_{2}-1)}\psi_{2}$
$\displaystyle=0$ (103b)

where $(E,m)$ are spectral parameters. The change of independent variable
$s_{2}^{2}=\cos^{2}\phi$ transforms (103b) to the trivial equation (32) with
solution $e^{im\phi}$ and the associated boundary conditions force
$m\in\mathbb{Z}$.

As with the previous degenerate systems, we have the following Lemma.

######  Lemma 37.

The spherical ODEs (103a), (103b) can be obtained by degenerating those
arising from ellispoidal coordinates (99)

For the equation in $\psi_{1}$ (103a), the change of independent variable
$s_{1}=x_{1}^{2}$ transforms (103a) into the form of the associated Legendre
equation given in (83) where $E=\ell(\ell+1)$ is always even and $\ell$ is an
integer. Polynomial solutions are this given by the Associated Legendre
polynomials $P^{m}_{\ell}(x_{1})$. We have the following Lemma.

######  Lemma 38.

For a given energy $E=\ell(\ell+1)$ the wave function of the Schrödinger
equation (5) in spherical coordinates on $S^{2}$ is one of the following
degree $\ell$ harmonic homogeneous polynomials

$\Psi_{\ell}=r^{\ell}\left(\frac{x_{2}+ix_{3}}{\sqrt{x_{2}^{2}+x_{3}^{2}}}\right)^{m}P_{\ell}^{m}\left(\frac{x_{1}}{r}\right),$
(104)

where $m=-\ell,\dots,\ell$.

###### Proof.

The form of (104) is obtained by recognising
$e^{im\phi}=\left(\frac{x_{2}+ix_{3}}{\sqrt{x_{2}^{2}+x_{3}^{2}}}\right)^{m}$.
Harmonicity is clear and homogeneity follows by noting that $P_{\ell}^{m}$ is
degree $\ell$ in its argument but the argument $x_{1}/r$ is degree 0 in terms
of the Cartesian coordinates. Moreover, $x_{2}^{2}+x_{3}^{2}=r^{2}-x_{1}^{2}$
and so the denominator $(x_{2}^{2}+x_{3}^{2})^{m/2}$ cancels with a factor in
$r^{m}P_{\ell}^{m}(x_{1}/r)$ and the remaining polynomial in $x_{1}/r$ is of
degree $\ell-m$. ∎

For discrete symmetries reflecting parity about the $x_{1}$ axis, we let
$\mu\in\\{0,1\\}$ denote the solution which is even $(0)$ or odd $(1)$ about
the $x_{1}$ axis. The symmetry is even when $\ell-m$ is even, and odd
otherwise. We note that, like for the spherical $S^{3}$ system, for a fixed
$\ell$ and fixed symmetry class only every 2nd $m$ is allowed so that the
parity of $\ell-m$ remains fixed.

###  A.3 Ellipsoidal $S^{3}$ Quantised Actions

Here, we show the quantised actions for each symmetry class $\bm{\mu}$
originally presented in Fig. 4 a) and b). The number of eigenstates per
symmetry class is given in Table 2. For a legend to convert between colour and
symmetry class, see Table 3.

Figure 13: Actions for each symmetry class shown in Fig. 2 a).

Figure 14: Actions for each symmetry class shown in Fig. 2 b).

## References

  * [1] M. Abramowitz and I. A. Stegun, editors.  Handbook of mathematical functions with formulas, graphs, and mathematical tables.  Dover Publications Inc., New York, 1992.  Reprint of the 1972 edition.
  * [2] A. M. Al-Rashed and N. Zaheer.  Zeros of Stieltjes and Van Vleck polynomials and applications.  Journal of Mathematical Analysis and Applications, 110(2):327–339, sep 1985.
  * [3] M. Alam.  Zeros of Stieltjes and Van Vleck polynomials.  Transactions of the American Mathematical Society, 252(0):197–204, 1979.
  * [4] J. Alonso, H. R. Dullin, and S. Hohloch.  Symplectic classification of coupled angular momenta.  Nonlinearity, 33:417–468, 2019.
  * [5] F. Arscott, I. Sneddon, M. Stark, and S. Ulam.  Periodic Differential Equations: An Introduction to Mathieu, Lamé, and Allied Functions.  International series of monographs in pure and applied mathematics. Elsevier Science, 2014.
  * [6] S. Axler, P. Bourdon, and W. Ramey.  Graduate Texts in Mathematics: Harmonic Function Theory.  Springer New York, 2001.
  * [7] H. Bateman and A. Erdélyi.  Higher transcendental functions Vol. III.  California Institute of Technology. Bateman Manuscript project. McGraw-Hill, New York, NY, 1955.
  * [8] M. Bellon and M. Talon.  Spectrum of the quantum Neumann model.  Physics Letters A, 337(4):360–368, 2005.
  * [9] M. Bellon and M. Talon.  The quantum Neumann model: Refined semiclassical results.  Physics Letters A, 356(2):110–114, 2006.
  * [10] M. P. Bellon and M. Talon.  The Quantum Neumann model: Asymptotic analysis.  Phys. Lett. A, 351:283–289, 2006.
  * [11] Z. Chen, T.-J. Kuo, and C.-S. Lin.  The geometry of generalized Lamé equation, i.  Journal de Mathématiques Pures et Appliquées, 127:89–120, Jul 2019.
  * [12] I. Chiscop, H. R. Dullin, K. Efstathiou, and H. Waalkens.  A Lagrangian fibration of the Isotropic 3-Dimensional Harmonic Oscillator with monodromy.  Journal of Mathematical Physics, 60(3):032103, 2019.
  * [13] R. H. Cushman and J. J. Duistermaat.  The quantum mechanical spherical pendulum.  Bull. Amer. Math. Soc., 19:475–479, 1988.
  * [14] S. R. Dawson, H. R. Dullin, and D. M. H. Nguyen.  Monodromy in Prolate Spheroidal harmonics.  Studies in Applied Mathematics, 146(4):953–982, 2021.
  * [15] S. R. Dawson, H. R. Dullin, and D. M. H. Nguyen.  The harmonic Lagrange top and the confluent Heun equation.  Regular and Chaotic Dynamics, 27(4):443–459, 2022.
  * [16] NIST Digital Library of Mathematical Functions.  http://dlmf.nist.gov/, Release 1.0.24 of 2019-09-15.  F. W. J. Olver, A. B. Olde Daalhuis, D. W. Lozier, B. I. Schneider, R. F. Boisvert, C. W. Clark, B. R. Miller, B. V. Saunders, H. S. Cohl, and M. A. McClain, eds.
  * [17] J. J. Duistermaat.  On global action-angle coordinates.  Comm. Pure Appl. Math., 33:687–706, 1980.
  * [18] H. R. Dullin and H. Waalkens.  Defect in the joint spectrum of hydrogen due to monodromy.  Phys. Rev. Lett., 120:020507, 2018.
  * [19] L. P. Eisenhart.  Separable Systems in Euclidean 3-Space.  Annals of Mathematics, 35:284, 1934.
  * [20] A. Erdélyi, W. Magnus, F. Oberhettinger, and F. G. Tricomi.  Higher Transcendental Functions. Vol. III.  McGraw-Hill Book Company, Inc., New York-Toronto-London, 1955.
  * [21] D. Gurarie.  Quantized Neumann problem, separable potentials on ${S}^{n}$ and the Lame equation.  Journal of Mathematical Physics, 36(10):5355–5391, 1995.
  * [22] D. Gurarie.  Symmetries and Laplacians: Introduction to Harmonic Analysis, Group Representations and Applications.  Springer, 2008.
  * [23] E. G. Kalnins, J. M. Kress, and W. Miller.  Second-order superintegrable systems in conformally flat spaces. IV. The classical 3D Stäckel transform and 3D classification theory.  Journal of mathematical physics, 47(4):043514, 2006.
  * [24] E. G. Kalnins, J. M. Kress, and W. Miller.  Separation of Variables and Superintegrability.  Bristol: IOP Publishing, 2018.
  * [25] E. G. Kalnins, J. M. Kress, W. Miller, and G. S. Pogosyan.  Complete sets of invariants for dynamical systems that admit a separation of variables.  Journal of Mathematical Physics, 43(7):3592–3609, 2002.
  * [26] E. G. Kalnins and W. Miller.  Separation of variables on the $n$-dimensional Riemannian manifolds. i. the $n$-sphere $S^{n}$ and Euclidean $n$-space $\mathbb{R}^{n}$.  Journal of Mathematical Physics, 27(7):1721–1736, 1986.
  * [27] Y. Le Floch.  Berezin-Toeplitz quantization and semitoric systems.  Habilitation à diriger des recherches, Université de Strasbourg, 2023.
  * [28] Y. Le Floch and A. Pelayo.  Symplectic geometry and spectral properties of classical and quantum coupled angular momenta.  Journal of Nonlinear Science, 29(2):655–708, 2018.
  * [29] Y. Le Floch and S. Vũ Ngoc.  The inverse spectral problem for quantum semitoric systems.  arXiv preprint, 2021.
  * [30] P. Morse and H. Feshbach.  Methods of Theoretical Physics.  McGraw-Hill, New York, 1953.
  * [31] G. Mussardo, V. Riva, G. Sotkov, and G. Delfino.  Kink scaling functions in 2d non-integrable quantum field theories.  Nuclear Physics B, 736(3):259–287, feb 2006.
  * [32] S. V. Ngọc.  Quantum monodromy in integrable systems.  Comm. Math. Phys., 203(2):465–479, 1999.
  * [33] D. M. H. Nguyen, S. R. Dawson, and H. R. Dullin.  Integrable systems arising from separation of variables on ${S}^{3}$, 2023.
  * [34] M. Pawellek.  Quasi-doubly periodic solutions to a generalized Lamé equation.  Journal of Physics A: Mathematical and Theoretical, 40(27):7673–7686, jun 2007.
  * [35] A. Pelayo and S. V. Ngọc.  Semitoric integrable systems on symplectic 4-manifolds.  Inventiones mathematicae, 177(3):571–597, 2009.
  * [36] H. P. Robertson.  Bemerkung über separierbare Systeme in der Wellenmechanik.  Mathematische Annalen, 98:749, 1927.
  * [37] A. Ronveaux and F. Arscott.  Heun’s differential equations.  Oxford University Press, 1995.
  * [38] D. A. Sadovskií and B. I. Zĥilinskií.  Monodromy, diabolic points, and angular momentum coupling.  Phys. Lett. A, 256(4):235–244, 1999.
  * [39] K. Schöbel.  Are orthogonal separable coordinates really classified?  Symmetry, Integrability and Geometry: Methods and Applications (SIGMA), 12:No 041, 16, 2016.
  * [40] K. Schöbel and A. P. Veselov.  Separation coordinates, moduli spaces and Stasheff polytopes.  Communications in Mathematical Physics, 337(3):1255–1274, 2015.
  * [41] E. Sinitsyn and B. Zhilinskii.  Qualitative analysis of the Classical and Quantum Manakov Top.  SIGMA. Symmetry, Integrability and Geometry: Methods and Applications, 3:Paper 046, 23 p., 2007.
  * [42] P. Stäckel.  Über die Bewegung eines Punktes in einer n-fachen Mannigfaltigkeit.  Math. Ann., 42:537–563, 1893.
  * [43] G. Szegő.  Orthogonal Polynomials.  American Math. Soc: Colloquium publ. American Mathematical Society, 1975.
  * [44] J. A. Toth.  On a class of spherical harmonics associated with rigid body motion.  Mathematical Research Letters, 1(2):203–210, 1994.
  * [45] H. Volkmer.  Expansions in Products of Heine-Stieltjes Polynomials.  Constructive Approximation, 15(4):467–480, oct 1999.