# Quantum three-rotor problem in the identity representation

Govind S. Krishnaswamia<EMAIL_ADDRESS>Himalaya Senapatia,b
<EMAIL_ADDRESS>aPhysics Department, Chennai Mathematical Institute, SIPCOT
IT Park, Siruseri 603103, India
bHSBC, RMZ Ecospace Business Park, Bellandur, Bengaluru 560103, India

(July 21, 2024)

###### Abstract

The quantum three-rotor problem concerns the dynamics of three equally massive
particles moving on a circle subject to pairwise attractive cosine potentials
and can model coupled Josephson junctions. Classically, it displays order-
chaos-order behavior with increasing energy. The quantum system admits a
dimensionless coupling with semiclassical behavior at strong coupling. We
study stationary states with periodic ‘relative’ wave functions. Perturbative
and harmonic approximations capture the spectrum at weak coupling and that of
low-lying states at strong coupling. More generally, the cumulative
distribution of energy levels obtained by numerical diagonalization is well-
described by a Weyl-like semiclassical estimate. However, the system has an
$S_{3}\times\mathbb{Z}_{2}$ symmetry that is obscured when working with
relative angles. By exploiting a basis for invariant states, we obtain the
purified spectrum in the identity representation. To uncover universal quantum
hallmarks of chaos, we partition the spectrum into energy windows where the
classical motion is regular, mixed or chaotic and unfold each separately. At
strong coupling, we find striking signatures of transitions between regularity
and chaos: spacing distributions morph from Poisson to Wigner-Dyson while the
number variance shifts from linear to logarithmic behavior at small lengths.
Some nonuniversal features are also examined. For instance, the number
variance saturates and oscillates at large lengths for strong coupling while
deviations from Poisson spacings at asymptotically low and high energies are
well-explained by purified quantum harmonic and free-rotor spectra at strong
and weak coupling. Interestingly, the degeneracy of free-rotor levels admits
an elegant formula that we deduce using properties of Eisenstein primes.

###### Contents

  1. I Introduction
  2. II Formulation
    1. II.1 Center of mass and relative variables
    2. II.2 Nondimensional coupling and Hamiltonian
    3. II.3 Separation of variables and boundary conditions
    4. II.4 Symmetries of the three-rotor system
  3. III Some approximate results
    1. III.1 Variational approach
    2. III.2 Perturbative approach
    3. III.3 Harmonic approximation around ground state
  4. IV Numerical diagonalization and purification of the spectrum
    1. IV.1 Total spectrum via Fourier basis
      1. IV.1.1 Density of states and Weyl’s law
    2. IV.2 Purification of the spectrum
      1. IV.2.1 Unitary irreducible representations of $S_{3}\times{\mathbb{Z}}_{2}$
      2. IV.2.2 Basis for the identity representation
      3. IV.2.3 Matrix elements of $\hat{H}$ in $s_{m,n}$ basis
      4. IV.2.4 Purified spectrum in the identity representation
  5. V Spacing distributions
    1. V.1 Semiclassical regime of large $\tilde{g}$
    2. V.2 Harmonic approximation: large $\tilde{g}$ and low $\hat{E}$
    3. V.3 Spacings at low and moderate $\tilde{g}$
  6. VI Number variance
  7. VII Discussion
  8. A Free rotor degeneracy and Eisenstein primes
  9. B Basis labels for identity representation

## I Introduction

The study of quantum manifestations of chaos and the transition from regular
to chaotic behavior in Hamiltonian systems has been a fertile area of research
with several applications such as to the electronic states of a donor impurity
in a semiconductor crystal [1] and to quantum dynamical localization and
transport [2]. The development of quantum chaology has benefited from the
study of several model systems such as billiard problems [3], the anisotropic
Kepler problem [1], free particle motion on constant negative curvature
Riemann surfaces [4], kicked rotors and tops [5] and the planar elastic
pendulum [6]. Some themes that have emerged from these studies in quantifying
quantum chaos are energy level spacing statistics (level repulsion and avoided
crossings) [7], trace formulae and sums over periodic orbits as a bridge
between classical and quantum mechanics [1], the nature of Wigner and Husimi
functions [8], nodal patterns of wave functions [9] and out-of-time-ordered
correlators [10, 11]. These investigations have revealed both universal
features (modeled by random matrix ensembles) depending only on the symmetry
class of the model as well as system-dependent features (such as oscillations
in the spectral rigidity) on various scales [12].

In this paper, we focus on the quantum three-rotor system that models a chain
of coupled Josephson Junctions [13, 14, 15, 16]. Interestingly, the classical
version has also been found to display rich features [17, 18]. In the
classical three-rotor problem, point particles of mass $m$ move on a circle of
radius $r$ subject to attractive cosine interparticle potentials. If the rotor
angles are denoted $\theta_{1,2,3}$, then their total energy is

$\displaystyle E_{\rm cl}$ $\displaystyle=$
$\displaystyle\frac{mr^{2}}{2}(\dot{\theta}_{1}^{2}+\dot{\theta}_{2}^{2}+\dot{\theta}_{3}^{2})+g[3-\cos(\theta_{1}-\theta_{2})$
(2)
$\displaystyle-\cos(\theta_{2}-\theta_{3})-\cos(\theta_{3}-\theta_{1})].\;\quad$

In terms of the center of mass angle
$\varphi_{0}=(\theta_{1}+\theta_{2}+\theta_{3})/3$ and relative angles
$\varphi_{1}=\theta_{1}-\theta_{2}$ and $\varphi_{2}=\theta_{2}-\theta_{3}$,
the energy of relative motion is

$\displaystyle E$ $\displaystyle=$
$\displaystyle\frac{mr^{2}}{3}\left(\dot{\varphi}_{1}^{2}+\dot{\varphi}_{2}^{2}+\dot{\varphi}_{1}\dot{\varphi}_{2}\right)$
(3) $\displaystyle+$ $\displaystyle
g[3-\cos\varphi_{1}-\cos\varphi_{2}-\cos(\varphi_{1}+\varphi_{2})].$ (4)

In the application to coupled Josephson junctions, rotor angles correspond to
superconducting phases of metallic segments [16]. Quite apart from this
physical realization, the three-rotor problem is interesting for its rich
dynamics. For instance, it displays order-chaos-order behavior with increasing
relative energy $E$ and a globally chaotic phase ($5.33g\lesssim E\lesssim
5.6g$) where the dynamics has been argued to be ergodic and mixing [18].
Intriguingly, the transition to widespread chaos around $E=4g$ seems to be
associated with (i) an accumulation of a doubly infinite sequence of fork-like
isochronous and period-doubling bifurcations [16] in the pendulum family of
periodic orbits and (ii) a loss of strict positivity in the curvature of the
Jacobi-Maupertuis metric on the configuration space [17]. Several other
remarkable phenomena concerning the classical three-rotor problem have been
reported in Refs. [17, 18, 16].

A natural question concerns the quantum dynamics of three rotors and in
particular, quantum manifestations of its order-chaos-order behavior. In this
paper, we initiate the study of the quantum three-rotor problem and address
this question by examining the spectrum of stationary states and investigate
both universal as well as system-specific properties. We begin in Section II
by formulating the quantum theory via canonical quantization. The parameters
$m,r$ and $g$ of the classical system do not admit any dimensionless
combination. Thus, the energy $E$ of a state in units of $g$ served as a
dimensionless classical control parameter, with order-chaos-order behavior
with increasing energy [17]. A distinctive feature of the quantum system is
the possibility of using Planck’s constant to define a dimensionless coupling
$\tilde{g}=gmr^{2}/\hbar^{2}$, which helps in organizing the quantum theory.
The weak coupling limit $\tilde{g}\to 0$ is an extreme quantum limit while
semiclassical behavior may be expected when $\tilde{g}$ is large. In Section
II.3, we show that the passage from rotor angles to center of mass and
relative angles also entails new features following from the requirement that
the wave function be singlevalued. The center of mass and relative wave
functions arising from separation of variables in the Schrödinger eigenvalue
problem must satisfy any one of three ‘periodicity up to a phase’ boundary
conditions, labelled by the cube roots of unity. In the sequel, we restrict to
the simplest type of boundary condition and the corresponding sector with
center of mass angular momentum quantum number $\ell$ divisible by $3$.
Although the use of relative angles reduces the number of degrees of freedom,
it obscures an $S_{3}\times\mathbb{Z}_{2}$ permutation-parity symmetry of the
model. In Section II.4, we identify the action of this group on relative
angles. It will be seen to play a crucial role in ‘purifying’ the energy
spectrum into irreducible representations.

Section III contains results that employ variational, perturbative and
harmonic approximations. Our variational ansatz accurately captures the energy
of the ground state over the entire range of values of $\tilde{g}$ from weak
to strong coupling. Perturbation theory reveals the pattern of energy level
degeneracies and their breaking at weak coupling. By exploiting number
theoretic properties of the ring of Eisenstein integers, we derive a closed
form expression for the unperturbed level degeneracies in terms of prime
factorization and explain why (with the exception of the ground state) they
are multiples of $6$ (details are provided in Appendix A). On the other hand,
our harmonic approximation around the ground state retains the
$S_{3}\times\mathbb{Z}_{2}$ symmetry and works well for a significant number
of low lying states at strong coupling.

In Section IV, we begin by describing a spectral method for numerical
diagonalization of the relative Hamiltonian in a Fourier basis. We find that,
aside from the very low-lying states, the cumulative distribution of energy
levels is closely captured by a Weyl-like semiclassical approximation that we
derive. To examine spectral statistics and its universality, we need to
‘purify’ the spectrum, i.e., decompose it into unitary irreducible
representations of $S_{3}\times\mathbb{Z}_{2}$. After enumerating these
representations in Section IV.2, we restrict attention to the identity
representation. A basis for the invariant states carrying this representation
is constructed and the matrix elements of the Hamiltonian are found in this
basis. Numerical diagonalization of a truncated Hamiltonian of size $\sim
3\times 10^{4}$ in this basis allows us to reliably determine the lowest $\sim
10^{4}$ levels and their spacings.

Having obtained the purified spectrum in the identity representation, we
examine the statistics of nearest neighbor spacings in Section V. To look for
quantum manifestations of the classical order-chaos-order behavior with
increasing energy, we first consider the semiclassical regime of large
$\tilde{g}\gtrsim 1000$. Here, we partition the spectrum into energy windows
where the classical dynamics is regular, mixed and chaotic, and unfold each
separately to have approximately unit mean spacing. This enables us to glimpse
striking universal signatures of the transitions between regularity and chaos
that were hidden in the unpurified or unpartitioned spectrum: spacing
distributions morph from Poisson to Wigner-Dyson (Gaussian orthogonal
ensemble, since the Hamiltonian is time-reversal invariant and real symmetric)
and back. However, even in the semiclassical regime, the low-lying spectrum
does not show Poisson statistics, although the classical dynamics is nearly
integrable. Interestingly, the behavior at low energies is accurately captured
by a purified version of our harmonic approximation (Section V.2). On the
other hand, in the weak coupling limit of small $\tilde{g}$, we see departures
from universal semiclassical expectations. In fact, in Section V.3, we see
signatures of the purified free-rotor spectrum superposed on the semiclassical
Poisson distributions of spacings at high energies.

Section VI is devoted to another spectral statistic: the number variance
$\Sigma(L)$. It is the variance in the number of unfolded levels in a spectral
window of length $L$. The variance is computed via an ensemble average over
the centers of such spectral windows. As with level spacings, we find that
evaluating the number variance in energy windows where the classical dynamics
is chaotic or regular facilitates interpretation and comparison with universal
expectations for small $L$ from the Gaussian orthogonal ensemble and a Poisson
process. We find that for relatively large $\tilde{g}$, the number variance in
classically regular (for $\tilde{g}=300$) and chaotic (for $\tilde{g}=6000$)
energy windows shows the expected linear and logarithmic behavior for small
$L$ and then displays saturation and nonuniversal oscillations when $L$ is
larger. We conclude in Section VII with a brief discussion of our results and
open questions and conjectures inspired by our work.

## II Formulation

The classical Hamiltonian for the equal-mass three rotor system is

$H_{\rm
cl}=\sum_{i=1}^{3}\left\\{\frac{\pi_{i}^{2}}{2mr^{2}}+g[1-\cos\left(\theta_{i}-\theta_{i+1}\right)]\right\\}$
(5)

with $\theta_{1}\equiv\theta_{4}$ and $g>0$ the strength of the attractive
coupling. Here, $\theta_{i}$ are the rotor angles and $\pi_{i}$ their
conjugate angular momenta: $\\{\theta_{i},\pi_{j}\\}=\delta_{ij}$ [17].
Quantizing canonically by taking $\pi_{i}=-i\hbar\partial_{\theta_{i}}$, we
get the quantum Hamiltonian

$H_{\rm
tot}=\sum_{i=1}^{3}-\frac{\hbar^{2}}{2mr^{2}}\frac{\partial^{2}}{\partial\theta_{i}^{2}}+g[1-\cos(\theta_{i}-\theta_{i+1})].$
(6)

### II.1 Center of mass and relative variables

As in the classical theory[17], it is convenient to introduce center of mass
(CM) and relative angles

$\begin{pmatrix}\varphi_{0}\\\ \varphi_{1}\\\
\varphi_{2}\end{pmatrix}=\begin{pmatrix}1/3&1/3&1/3\\\ 1&-1&0\\\
0&1&-1\end{pmatrix}\begin{pmatrix}\theta_{1}\\\ \theta_{2}\\\
\theta_{3}\end{pmatrix},$ (7)

in terms of which $H_{\text{tot}}=H_{\rm CM}+H$ where

$\displaystyle H_{\rm CM}$ $\displaystyle=$
$\displaystyle-\frac{\hbar^{2}}{6mr^{2}}\partial^{2}_{\varphi_{0}}\quad\text{and}\quad
H=\frac{\hbar^{2}}{mr^{2}}T+gV(\varphi_{1},\varphi_{2})$ (8)
$\displaystyle\text{with}\;\;T$ $\displaystyle=$
$\displaystyle-\left(\partial^{2}_{\varphi_{1}}-\partial_{\varphi_{1}}\partial_{\varphi_{2}}+\partial^{2}_{\varphi_{2}}\right)$
(9) $\displaystyle\text{and}\;\;V$ $\displaystyle=$
$\displaystyle-\left[\cos\varphi_{1}+\cos\varphi_{2}+\cos(\varphi_{1}+\varphi_{2})-3\right].\quad$
(10)

Here, the conjugate momenta are represented by
$p_{j}=-i\hbar\partial_{\varphi_{j}}$ for $j=0,1$ and $2$.

### II.2 Nondimensional coupling and Hamiltonian

The classical three-rotor model has three parameters $m,g$ and $r$. Since
their dimensions are independent, the classical theory does not admit any
dimensionless parameter. However, physical quantities may be
nondimensionalized by measuring them in units of suitable combinations of
$m,r$ and $g$. For instance, energies of states may be measured in units of
$g$ and time in units of $\sqrt{mr^{2}/g}$. On the other hand, the quantum
theory admits one dimensionless parameter due to the introduction of Planck’s
constant. We may take this dimensionless parameter to be the coupling constant
$\tilde{g}=mr^{2}g/\hbar^{2}$. Thus, upon nondimensionalization, measured
values of physical quantities (eigenvalues of observables) become
dimensionless functions of $\tilde{g}$. Two limiting values of $\tilde{g}$ are
noteworthy. The weak coupling limit or extreme quantum limit is the one where
$\tilde{g}\to 0$. The strong coupling limit $\tilde{g}\to\infty$ is a
semiclassical limit.

Bearing this in mind, it is convenient to define a dimensionless Hamiltonian.
In fact, there are two obvious possibilities

$\tilde{H}=\frac{H}{g}=\frac{1}{\tilde{g}}T+V\quad\text{and}\quad\hat{H}=\frac{H}{\hbar^{2}/mr^{2}}=T+\tilde{g}V,\quad$
(11)

which are related by $\hat{H}=\tilde{g}\tilde{H}=(\tilde{g}/g)H$. We will
denote the eigenvalues of $H$, $\hat{H}$ and $\tilde{H}$ by $E$, $\hat{E}$ and
$\tilde{E}$. Similarly, $\hat{E}_{\rm CM}$ denotes eigenvalues of
$\hat{H}_{\rm CM}$.

### II.3 Separation of variables and boundary conditions

The boundary conditions (BC) on the wave function $\Psi$ in the CM and
relative variables follow from the requirement that it be single-valued when
viewed as a function of the $2\pi$-periodic rotor angles $\theta_{1,2,3}$. For
instance, from (7), $\theta_{1}\mapsto\theta_{1}+2\pi$ implies that
$\varphi_{0}\mapsto\varphi_{0}+2\pi/3$ and
$\varphi_{1}\mapsto\varphi_{1}+2\pi$. Thus, we have the following ‘boundary’
conditions on $\Psi(\varphi_{0,1,2})$:

$\displaystyle\Psi(\varphi_{0,1,2})$ $\displaystyle=$
$\displaystyle\Psi(\varphi_{0}\pm 2\pi/3,\varphi_{1}\pm 2\pi,\varphi_{2})$
(12) $\displaystyle=$ $\displaystyle\Psi(\varphi_{0}\pm
2\pi/3,\varphi_{1},\varphi_{2}\mp 2\pi)$ (13) $\displaystyle=$
$\displaystyle\Psi(\varphi_{0}\pm 2\pi,\varphi_{1},\varphi_{2}).$ (14)

Since $H_{\rm tot}$ is a sum of CM and relative energies, we may factorize
$\Psi(\varphi_{0},\varphi_{1},\varphi_{2})=\Upsilon(\varphi_{0})\psi(\varphi_{1},\varphi_{2})$
in the Schrödinger eigenvalue problem $H_{\rm tot}\Psi=E_{\rm tot}\Psi$. We
now derive BCs on $\psi(\varphi_{1},\varphi_{2})$ and $\Upsilon(\varphi_{0})$
by imposing (14). To begin with, (14) implies

$\displaystyle\Upsilon(\varphi_{0})$ $\displaystyle=$
$\displaystyle\Upsilon(\varphi_{0}\pm 2\pi)\quad\text{and}$ (15)
$\displaystyle\Upsilon(\varphi_{0})\psi(\varphi_{1},\varphi_{2})$
$\displaystyle=$ $\displaystyle\Upsilon(\varphi_{0}\pm
2\pi/3)\psi(\varphi_{1}\pm 2\pi,\varphi_{2})$ (16) $\displaystyle=$
$\displaystyle\Upsilon(\varphi_{0}\pm 2\pi/3)\psi(\varphi_{1},\varphi_{2}\mp
2\pi).$ (17)

Now, without loss of generality, we may write

$\psi(\varphi_{1}+2\pi,\varphi_{2})=\varpi_{1}(\varphi_{1},\varphi_{2})\psi(\varphi_{1},\varphi_{2}),$
(18)

for some function $\varpi_{1}(\varphi_{1},\varphi_{2})$. This by (17) implies

$\Upsilon(\varphi_{0}+2\pi/3)=\varpi_{1}(\varphi_{1},\varphi_{2})^{-1}\Upsilon(\varphi_{0}).$
(19)

Thus $\varpi_{1}$ must be a constant which again by (17) must be a cube root
of unity:

$\Upsilon(\varphi_{0})=\Upsilon(\varphi_{0}+2\pi)\quad\Rightarrow\quad\varpi_{1}^{3}=1.$
(20)

Similarly, we may write

$\psi(\varphi_{1},\varphi_{2}-2\pi)=\varpi_{2}(\varphi_{1},\varphi_{2})\psi(\varphi_{1},\varphi_{2})$
(21)

which by (17) implies

$\Upsilon(\varphi_{0}+2\pi/3)=\varpi_{2}^{-1}\Upsilon(\varphi_{0}).$ (22)

From (19) and (22), $\varpi_{1}=\varpi_{2}\equiv\varpi$. In summary, the BCs
on the CM and relative wave functions are

$\displaystyle\Upsilon(\varphi_{0}+2\pi/3)$ $\displaystyle=$
$\displaystyle\varpi^{-1}\Upsilon(\varphi_{0})\quad\text{and}$ (23)
$\displaystyle\psi(\varphi_{1}+2\pi,\varphi_{2})$ $\displaystyle=$
$\displaystyle\psi(\varphi_{1},\varphi_{2}-2\pi)=\varpi\psi(\varphi_{1},\varphi_{2})$
(24)

where $\varpi$ must be a cube root of unity. Hence, we arrive at the CM and
relative Schrödinger eigenvalue problems which must be solved subject to (24):

$-\frac{1}{6}\partial^{2}_{\varphi_{0}}\Upsilon=\hat{E}_{\rm
CM}\Upsilon\quad\text{and}\quad\left[T+\tilde{g}V\right]\psi=\hat{E}\psi$ (25)

with $\hat{E}_{\rm tot}=\hat{E}_{\rm CM}+\hat{E}$. Although the BCs for
$\varpi\neq 1$ render $\Upsilon$ and $\psi$ discontinuous (say at 0, $\pm
2\pi/3$ on the $\varphi_{0}$-circle and along the two cycles $\varphi_{1,2}=0$
on the $\varphi_{1}$-$\varphi_{2}$ torus), the total wave function $\Psi$ is
nevertheless continuous.

The center of mass eigenvalue problem is particularly simple with the
solutions

$\Upsilon(\varphi_{0})=e^{i\ell\varphi_{0}}\quad\text{and}\quad\hat{E}_{\rm
CM}=\ell^{2}/6$ (26)

where $\ell$ is an integer with $\ell\equiv 0,1,2$ $(\text{mod~{}}3)$
according as $\varpi=1$, $e^{4\pi i/3}$ and $e^{2\pi i/3}$. All CM energy
levels other than $E_{\rm CM}=0$ are doubly degenerate (clockwise and
counterclockwise rotation of the CM corresponding to $\ell$ and $-\ell$) and
is a consequence of the parity symmetry $\theta_{i}\mapsto-\theta_{i}$. For
simplicity, we will henceforth restrict attention to the case $\varpi=1$ and
take $\ell=0$ so that the CM is at rest. In this case, the BC on
$\psi(\varphi_{1},\varphi_{2})$ reduces to $2\pi$ periodicity in each relative
angle.

### II.4 Symmetries of the three-rotor system

Simultaneous rotation of rotor angles $\theta_{j}\to\theta_{j}+\theta$ by any
fixed $\theta$ is a symmetry of $H_{\rm tot}=H_{\rm CM}+H$. It takes
$\varphi_{0}\to\varphi_{0}+\theta$ but acts trivially on the relative angles
$\varphi_{1,2}$ so that it is a symmetry of $H_{\rm CM}$ but not $H$. The
system also admits discrete symmetries. $H_{\rm tot}$ is symmetric under
$G=S_{3}\times\mathbb{Z}_{2}$ where
$S_{3}=\\{e,\tau_{12},\tau_{23},\tau_{31},\sigma_{123},\sigma_{132}\\}$ is the
group of permutations of the rotor angles and $\mathbb{Z}_{2}=\\{e,\pi\\}$ is
the parity group with $\pi$ taking $\theta_{j}\to-\theta_{j}$ (see Table 1).
These continue to be symmetries of the relative Hamiltonian $H$ with the
action on relative angles given in Table 1. The CM angle $\varphi_{0}$ is
invariant under rotor permutations but reverses sign under parity. As a
consequence, $H_{\rm CM}$ is also invariant under $G$. These discrete
symmetries will play a crucial role in ‘purifying’ the spectrum of the
Hamiltonian into irreducible representations of $G$.

g | acts on $(\theta_{1},\theta_{2},\theta_{3})$ | acts on $(\varphi_{1},\varphi_{2})$ | acts on $(m,n)$
---|---|---|---
$e$ | $(\theta_{1},\theta_{2},\theta_{3})$ | $(\varphi_{1},\varphi_{2})$ | $(m,n)$
$\pi$ | $(-\theta_{1},-\theta_{2},-\theta_{3})$ | $(-\varphi_{1},-\varphi_{2})$ | $(-m,-n)$
$\tau_{12}$ | $(\theta_{2},\theta_{1},\theta_{3})$ | $(-\varphi_{1},\varphi_{1}+\varphi_{2})$ | $(n-m,n)$
$\tau_{23}$ | $(\theta_{1},\theta_{3},\theta_{2})$ | $(\varphi_{1}+\varphi_{2},-\varphi_{2})$ | $(m,m-n)$
$\tau_{31}$ | $(\theta_{3},\theta_{2},\theta_{1})$ | $(-\varphi_{2},-\varphi_{1})$ | $(-n,-m)$
$\sigma_{123}$ | $(\theta_{2},\theta_{3},\theta_{1})$ | $(\varphi_{2},-\varphi_{1}-\varphi_{2})$ | $(-n,m-n)$
$\sigma_{132}$ | $(\theta_{3},\theta_{1},\theta_{2})$ | $(-\varphi_{1}-\varphi_{2},\varphi_{1})$ | $(n-m,-m)$
Table 1: Action of the elements of the parity ${\mathbb{Z}}_{2}=\\{e,\pi\\}$
and permutation $S_{3}=\\{e,\tau,\sigma\\}$ groups on rotor angles
$\theta_{1,2,3}$, relative angles $\varphi_{1,2}$ and the Fourier basis
element $(m,n)\equiv e^{i(m\varphi_{1}+n\varphi_{2})}/2\pi$. The action of
ordered pairs of elements of ${\mathbb{Z}}_{2}$ and $S_{3}$ is obtained by
composition. Under this action, ${\mathbb{Z}}_{2}\times S_{3}$ is a symmetry
of $H_{\rm tot}$ as well as of $H$.

## III Some approximate results

Here we discuss variational, perturbative and harmonic approximations that
allow us to understand the spectrum of the relative Hamiltonian $\hat{H}$ in
limiting cases. They also help to validate our numerical approach to
diagonalizing the Hamiltonian in Section IV.

### III.1 Variational approach

We will obtain a variational upper bound on the ground state energy

$\hat{E}_{0}=\min_{\psi}\langle\psi|\hat{H}|\psi\rangle/\langle\psi|\psi\rangle$
(27)

by minimizing over a class of wave functions obeying appropriate boundary
conditions. For $\varpi=1$, from (24), the wave functions and their first
partial derivatives in $\varphi_{1,2}$ must be $2\pi$-periodic. A constant
trial function $\psi$ leads to the crude upper bound $\hat{E}_{0}\leq
3\tilde{g}$. A better ansatz with dimensionless variational parameter $a$ is
given by $\psi(\varphi_{1},\varphi_{2})=e^{-aV}$, which automatically
satisfies the BCs as $V(\varphi_{1},\varphi_{2})$ (10) is periodic. The
expectation values and optimal values of $a$ are found numerically for
$10^{-3}\leq\tilde{g}\leq 100$. The resulting variational estimates for
$\hat{E}_{0}(\tilde{g})$ and optimal $a(\tilde{g})$ are plotted in Fig. 1.
They compare favorably with the results of numerical diagonalization in
Section IV. For instance, for $\tilde{g}=1$, the variational estimate
$\hat{E}_{0}^{\rm var}\approx 1.62$ for the optimal value $a\approx.68$ is
quite close to the numerically obtained $\hat{E}_{0}\approx 1.593$. By fitting
the variational estimates over the above range of $\tilde{g}$, we find the
following asymptotic behaviors:

$\displaystyle\hat{E}^{\rm var}_{0}$ $\displaystyle\approx$
$\displaystyle\begin{cases}3\tilde{g}-1.5{\tilde{g}}^{2}&\text{as}\quad\tilde{g}\to
0\\\
\sqrt{3\tilde{g}}&\text{as}\quad\tilde{g}\to\infty\end{cases}\quad\text{and}$
(28) $\displaystyle a$ $\displaystyle\approx$
$\displaystyle\begin{cases}\tilde{g}&\text{as}\quad\tilde{g}\to 0\\\
\sqrt{\tilde{g}/e}&\text{as}\quad\tilde{g}\to\infty.\end{cases}$ (29)

In the weak coupling limit ($\tilde{g}\to 0$), this variational estimate
matches the ground state energy obtained from second order perturbation theory
in Section III.2. Moreover, in this limit, $a\to 0$ and the optimal trial wave
function approaches the constant ground state wave function of noninteracting
rotors. On the other hand, in the semi-classical ($\tilde{g}\to\infty$) limit,
$a\to\infty$, and $\psi$ is concentrated at the minimum of $V$ with
$\tilde{E}^{\rm var}_{0}\sim\sqrt{3/\tilde{g}}\to 0$, in agreement with the
vanishing energy (in units of $g$) of the classical ground state [17].

(a)

(b)

Figure 1: (a) Variational ground state energies of $\hat{H}$ (blue dots) as a
function of $\tilde{g}$ based on the ansatz $\psi\propto e^{-aV}$ are very
close to the numerical values (red squares). Optimal values of the variational
parameter $a$ are also shown (black stars). (b) The variational ground state
normalized probability density on the $-\pi\leq\varphi_{1},\varphi_{2}\leq\pi$
torus for $\tilde{g}\approx 0.3$ where $a=0.3$. The state is localized around
the minimum $(0,0)$ of the potential $V$ (10).

### III.2 Perturbative approach

$x$$y$$-\omega^{2}\;(1,1)$$-\omega\;(0,-1)$$1\;(1,0)$$\quad\;-1\;(-1,0)$$\omega\;(0,1)$$\omega^{2}\;(-1,-1)$$(1,2)$$(1,-1)$$(2,1)$
Figure 2: Eigenstates of kinetic energy $T$, represented as $(m,n)$ pairs on
the $\mathbb{Z}+\omega\mathbb{Z}$ lattice on the complex plane for
$\omega=e^{2\pi i/3}$. Rotations by elements of the cyclic group
$C_{6}=\langle e^{i\pi/3}\rangle$ are symmetries of the lattice that preserve
the distance $|m+n\omega|$ from the origin (and hence preserve the eigenvalues
of $T$). So the cardinality of each orbit must be divisible by $6$ implying
that the degeneracies are multiples of 6. Dashed (red) and solid (blue) arrows
point to two such orbits corresponding to energies $T=1$ and $T=3$.

Here we use first order perturbation theory for the Hamiltonian
$\hat{H}=T+\tilde{g}V$ (11) to estimate the energies of several low-lying
states. The eigenvalues of $T$ are of the form $m^{2}+n^{2}-mn$ where $m$ and
$n$ are integers corresponding to the Fourier eigenmodes
$e_{m,n}=e^{i(m\varphi_{1}+n\varphi_{2})}/2\pi$. The first few eigenenergies
$T$ along with their degeneracies $d$ are $T^{d}=0^{1}$, $1^{6}$, $3^{6}$,
$4^{6}$, $7^{12}$, $9^{6}$, $12^{6}$ and $13^{12}$. Aside from the ground
state, the degeneracies must be multiples of six. To see this, notice that
$m^{2}+n^{2}-mn=|m+n\omega|^{2}$ where $\omega=e^{2\pi i/3}$ (or $e^{4\pi
i/3}$) is a nontrivial cube root of unity. Hence, we may view the energy
eigenvalue $m^{2}+n^{2}-mn$ for $m,n\in\mathbb{Z}$ as the square of the
distance from the origin to the point $m+n\omega$ lying on the triangular
integral lattice spanned by $(1,\omega)$ in the complex plane, as shown in
Fig. 2. The vertices of this lattice labelled $(m,n)$ represent the linearly
independent eigenstates $e_{m,n}$ of the kinetic energy $T$. The cyclic group
$C_{6}$ of rotations generated by $e^{i\pi/3}$ (sixth roots of unity) is a
symmetry of the lattice that preserves the distance from the origin. Thus, the
number of lattice points at a fixed distance from the origin must be a
multiple of six. In fact, using further properties of the Unique Factorization
Domain $\mathbb{Z}+\mathbb{Z}\omega$, one may obtain an explicit formula for
the degeneracies. Suppose the eigenenergy $T$ has the prime factorization
$3^{\gamma}p_{1}^{\alpha_{1}}\cdots p_{r}^{\alpha_{r}}q_{1}^{\beta_{1}}\cdots
q_{s}^{\beta_{s}}$, where $p_{i}\equiv 1\;\text{mod~{}}3$ and $q_{i}\equiv
2\;\text{mod~{}}3$ with $\beta_{i}$ even. Then the degeneracy of the
eigenvalue $T$ is (see Appendix A)

$\text{degen}(T)=6(\alpha_{1}+1)(\alpha_{2}+1)\cdots(\alpha_{r}+1).$ (30)

If any of the $\beta_{i}$ is odd, then $T$ cannot arise as an eigenvalue and
the ‘degeneracy’ vanishes. For example,
$T\in\\{1,3,4=2^{2},9=3^{2},12=2^{2}3\\}$ have no prime factors congruent to 1
modulo 3 and the $\beta$’s are even so that all these degeneracies are equal
to 6. $T=7$ and $T=13$ are both primes $\equiv 1$ $(\text{mod~{}}3)$ so that
$r=1$, $\alpha_{1}=1$ and $\gamma=s=0$ leading to a degeneracy of 12 for both
levels. These multiplicities agree with the ones given at the beginning of
this section.

The ground state $e_{0,0}$ of the noninteracting system $T$ is nondegenerate
with $\hat{E}_{0}^{(0)}=0$ and $\psi_{0}^{(0)}=1/2\pi$. This implies that

$\hat{E}_{0}^{(1)}=\langle 1/2\pi|V|1/2\pi\rangle=3.$ (31)

Thus, to the first order,
$\hat{E}_{0}\approx\hat{E}_{0}^{(0)}+\tilde{g}\hat{E}_{0}^{(1)}=3\tilde{g}$.
At second order, $\hat{E}_{0}^{(2)}=-1.5$, so that $\hat{E}_{0}\approx
3\tilde{g}-1.5\tilde{g}^{2}$. For small $\tilde{g}$, this agrees well with the
numerically obtained ground state energy as well as our variational
approximation (see (29) and Fig. 1).

The first excited state of $T$ is 6-fold degenerate with energy
$\hat{E}_{1}^{(0)}=1$ and the eigenspace is spanned by $e^{\pm
i\varphi_{1}}/2\pi$, $e^{\pm i\varphi_{2}}/2\pi$ and $e^{\pm
i(\varphi_{1}+\varphi_{2})}/2\pi$. We find that the matrix representation of
$V$ has eigenvalues $2$, $2.5$, $2.5$, $3.5$, $3.5$ and $4$. Thus, to first
order in degenerate perturbation theory, we have

$\displaystyle\hat{E}_{1}=1+2\tilde{g},\quad\hat{E}_{2,3}=1+2.5\tilde{g}$ (32)
$\displaystyle\hat{E}_{4,5}=1+3.5\tilde{g}\quad\text{and}\quad\hat{E}_{6}=1+4\tilde{g}.$
(33)

We find that this agrees well with the low-lying spectrum derived numerically
for small $\tilde{g}\lesssim 0.15$. For the second and third excited states of
$T$, first order perturbation theory does not break the 6-fold degeneracies,
and we get $\hat{E}_{7-12}=3+3\tilde{g}$ and $\hat{E}_{13-18}=4+3\tilde{g}$.
However, the 12-fold degenerate $4^{\rm th}$ excited state of $T$ splits into
two 6-fold degenerate levels $\hat{E}_{19-24}=7+2.5\tilde{g}$ and
$\hat{E}_{25-30}=7+3.5\tilde{g}$.

### III.3 Harmonic approximation around ground state

In Jacobi-like coordinates $\varphi_{\pm}=\varphi_{1}\pm\varphi_{2}$, the
kinetic energy operator (10) becomes diagonal:

$\displaystyle\hat{H}$ $\displaystyle=$
$\displaystyle-\left(\partial^{2}_{\varphi_{+}}+3\partial^{2}_{\varphi_{-}}\right)-\tilde{g}\bigg{[}\cos\frac{\varphi_{+}+\varphi_{-}}{2}$
(35)
$\displaystyle\quad+\cos\frac{\varphi_{+}-\varphi_{-}}{2}+\cos\varphi_{+}-3\bigg{]}.$

In the harmonic approximation around the classical ground state,
$\varphi_{\pm}$ are small and the Hamiltonian

$\hat{H}^{\rm
harm}=-\left(\partial^{2}_{\varphi_{+}}+3\partial^{2}_{\varphi_{-}}\right)+\frac{\tilde{g}}{4}(3\varphi_{+}^{2}+\varphi_{-}^{2})$
(36)

splits into a pair of uncoupled linear oscillators so that

$\hat{E}^{\rm harm}_{m,n}=\sqrt{3\tilde{g}}(m+n+1)\quad\text{for}\quad
m,n=0,1,2,\cdots.$ (37)

Pleasantly, the $S_{3}\times\mathbb{Z}_{2}$ symmetry of the three-rotor system
survives the harmonic approximation. The action of group elements on
$\varphi_{\pm}$ deduced from Table 1 are given in Table 2. It turns out that
$\hat{H}^{\rm harm}$ (36) is invariant under this group action as we now
argue. The kinetic energy is unaffected by the approximation. To see the
invariance of the potential, we note that up to an additive constant, it is
proportional to
$\varphi_{1}^{2}+\varphi_{2}^{2}+(\varphi_{1}+\varphi_{2})^{2}$ and observe
that elements of $S_{3}$ permute the three terms while parity just reverses
the signs of $\varphi_{1,2}$. The implementation of this symmetry in the
harmonic approximation will be relevant to the purification of the spectrum to
be discussed in Section V.2.

For large $\tilde{g}\approx 300$, we find that energies of low-lying states
are well captured by the harmonic approximation as seen in Fig. 3. Although
degeneracies of the 2d SHO are somewhat broken, the multiplet structure
($n^{\rm th}$ energy level being $n$-fold degenerate) is approximately
retained for the lowest few levels. This is to be expected because for large
$\tilde{g}$ (stiff spring), the low-lying states are strongly localized around
the classical ground state where the harmonic part of the potential dominates.
As $\tilde{g}$ is increased, a larger number of low-lying levels (roughly when
$m+n+1\ll\sqrt{\tilde{g}/3}$) follow the harmonic approximation. This
observation suggests that we should omit these low-lying states when comparing
spectral statistics with Poisson or Wigner-Dyson spacing distributions as the
spacing of an oscillator is quite different from the latter ones, see Fig. 5.

g | action on $(\varphi_{+},\varphi_{-})$
---|---
$e$ | $(\varphi_{+},\varphi_{-})$
$\pi$ | $(-\varphi_{+},-\varphi_{-})$
$\tau_{12}$ | $((\varphi_{+}-\varphi_{-})/2,-(3\varphi_{+}+\varphi_{-})/2$
$\tau_{23}$ | ($(\varphi_{+}+\varphi_{-})/2,(3\varphi_{+}-\varphi_{-})/2$)
$\tau_{31}$ | $(-\varphi_{+},\varphi_{-})$
$\sigma_{123}$ | $(-(\varphi_{+}+\varphi_{-})/2,(3\varphi_{+}-\varphi_{-})/2)$
$\sigma_{132}$ | $(-(\varphi_{+}-\varphi_{-})/2,-(3\varphi_{+}+\varphi_{-})/2)$
Table 2: Action of the elements of the parity ${\mathbb{Z}}_{2}=\\{e,\pi\\}$
and permutation $S_{3}=\\{e,\tau,\sigma\\}$ groups on Jacobi angles
$\varphi_{\pm}$. Figure 3: Lowest 100 numerically obtained eigenvalues of
$\hat{H}$ (red, lower) follow a pattern similar to that of the 2d harmonic
oscillator (blue, higher) provided $\tilde{g}$ is large.

## IV Numerical diagonalization and purification of the spectrum

In this section, we first describe a scheme to determine the spectrum of the
relative Hamiltonian by numerical diagonalization. The cumulative distribution
function of this ‘total spectrum’ is then compared with a Weyl-type
semiclassical asymptotic formula. We then address the question of purifying
the spectrum i.e., decomposing the total spectrum into various irreducible
representations of the discrete symmetry group $G=S_{3}\times\mathbb{Z}_{2}$.
Finally, an explicit basis for $G$-invariant states is constructed to compute
the purified spectrum in the identity representation.

### IV.1 Total spectrum via Fourier basis

As mentioned in II.3, for zero center of mass angular momentum ($\ell=0$), the
wave function $\psi$ must be $2\pi$-periodic in both $\varphi_{1}$ and
$\varphi_{2}$. Such wave functions may be expanded in an orthonormal Fourier
basis

$\langle\varphi_{1},\varphi_{2}|m,n\rangle=e_{m,n}=e^{i(m\varphi_{1}+n\varphi_{2})}/2\pi\equiv(m,n).$
(38)

Noting that the potential $V$ in (10) can be rewritten as

$\displaystyle V$ $\displaystyle=$
$\displaystyle-\frac{1}{2}\left(e^{i\varphi_{1}}+e^{-i\varphi_{1}}+e^{i\varphi_{2}}+e^{-i\varphi_{2}}\right)$
(40)
$\displaystyle-\frac{1}{2}\left(e^{i(\varphi_{1}+\varphi_{2})}+e^{-i(\varphi_{1}+\varphi_{2})}\right)+3,$

the matrix elements of $\hat{H}$ are

$\displaystyle\hat{H}_{a,b;m,n}$ $\displaystyle=$ $\displaystyle\langle
a,b|\hat{H}|m,n\rangle$ (41) $\displaystyle=$
$\displaystyle\delta^{a}_{m}\delta^{b}_{n}\left(3\tilde{g}+[m^{2}+n^{2}-mn]\right)$
(44)
$\displaystyle-\frac{\tilde{g}}{2}\left[\delta^{b}_{n}(\delta^{a}_{m+1}+\delta^{a}_{m-1})+\delta^{a}_{m}(\delta^{b}_{n+1}+\delta^{b}_{n-1})\right]$
$\displaystyle-\frac{\tilde{g}}{2}\left[\delta^{a}_{m+1}\delta^{b}_{n+1}+\delta^{a}_{m-1}\delta^{b}_{n-1}\right].$

To find the eigenvalues of $\hat{H}$, we proceed as follows. First, we
truncate the Fourier basis $e_{m,n}$ to $0\leq|m|,|n|\leq m_{\rm max}$
corresponding to the $(2m_{\rm max}+1)\times(2m_{\rm max}+1)$ integer lattice
on the plane centered at the origin. To get an ordered basis, we traverse this
lattice columnwise upwards starting from the south-west corner. This leads us
to an ordered orthonormal Fourier basis $e_{i}$ labeled by the single index
$i(m,n)=(2m_{\rm max}+1)m+n+2m_{\rm max}^{2}+2m_{\rm max}+1$ with $1\leq
i\leq(2m_{\rm max}+1)^{2}$. The matrix elements $\hat{H}_{ij}=\langle
e_{i}|\hat{H}|e_{j}\rangle$ are computed in this basis resulting in a
truncated Hamiltonian matrix of dimension $N=(2m_{\rm max}+1)^{2}$ that we
diagonalize numerically. However, the larger eigenvalues suffer from
truncation errors. By requiring that the eigenvalues do not change appreciably
when $N$ is increased, we isolate the reliably computed portion of the
spectrum.

#### IV.1.1 Density of states and Weyl’s law

Figure 4: Numerically computed number of energy levels $N(\hat{E},\tilde{g})$
below energy $\hat{E}$ for $\tilde{g}=300$ (blue dots) agrees well with
semiclassical estimate of Eq. (51) (black dashes for $\hat{E}>4.5\tilde{g}$
and black circles for $\hat{E}<4.5\tilde{g}$). Asymptotically,
$N\sim(2\pi/\sqrt{3})(\hat{E}-3\tilde{g})$.

A version of Weyl’s semiclassical law [19] for systems with a potential says
that the asymptotic number $N(E,\hbar)$ of energy eigenvalues below $E$ is
given by

$N(E,\hbar)\sim\frac{1}{(2\pi\hbar)^{2}}\int_{H<E}d\varphi_{1}dp_{1}d\varphi_{2}dp_{2}.$
(45)

The RHS is simply the energetically-allowed phase space volume in units of
$h^{2}=(2\pi\hbar)^{2}$. Since $N(E,\hbar)$ is dimensionless, it can depend on
$E$ and $\hbar$ only via the dimensionless energy
$\hat{E}=E/(\hbar^{2}mr^{2})$ and coupling $\tilde{g}=gmr^{2}/\hbar^{2}$.
Introducing the dimensionless angular momenta
$\tilde{p}_{i}=p_{i}/\sqrt{gmr^{2}}$ allows us to write

$N(\hat{E},\tilde{g})\sim\frac{\tilde{g}}{(2\pi)^{2}}\int_{\hat{H}\leq\hat{E}}d\varphi_{1}d\tilde{p}_{1}d\varphi_{2}d\tilde{p}_{2}.$
(46)

Here, from (10),

$\displaystyle\hat{H}$ $\displaystyle=$ $\displaystyle
T+\tilde{g}V=\tilde{g}[\tilde{p}_{1}^{2}-\tilde{p}_{1}\tilde{p}_{2}+\tilde{p}_{2}^{2}$
(48)
$\displaystyle+3-\cos\varphi_{1}-\cos\varphi_{2}-\cos(\varphi_{1}+\varphi_{2})].$

Formula (46) is expected to hold in a semiclassical limit where $\hat{E}\gg 1$
and/or $\tilde{g}\gg 1$. To evaluate (46), the angular momentum integral is
performed first for fixed $(\varphi_{1},\varphi_{2})$ followed by an integral
over the energetically allowed Hill region of the configuration torus. For
each point $(\varphi_{1},\varphi_{2})$, the energetically allowed region
$R_{\varphi_{1},\varphi_{2}}$ in angular momentum space is given by
$T\leq\hat{E}-\tilde{g}V$. Defining
$K(\varphi_{1},\varphi_{2})=\hat{E}/\tilde{g}-V(\varphi_{1},\varphi_{2})$, the
boundary of this region is the ellipse

$\tilde{p}_{1}^{2}-\tilde{p}_{1}\tilde{p}_{2}+\tilde{p}_{2}^{2}=K(\varphi_{1},\varphi_{2}).$
(49)

It may be brought to the standard form
$\tilde{p}_{+}^{2}/(2K)+\tilde{p}_{-}^{2}/(2K/3)=1$ upon defining
$\tilde{p}_{\pm}=(\tilde{p}_{1}\pm\tilde{p}_{2})/\sqrt{2}$, whence it encloses
an area $2\pi K/\sqrt{3}$. Consequently,

$\displaystyle N$ $\displaystyle\sim$
$\displaystyle\frac{\tilde{g}}{(2\pi)^{2}}\iint_{K(\varphi_{1},\varphi_{2})\geq
0}d\varphi_{1}d\varphi_{2}\int_{R_{\varphi_{1},\varphi_{2}}}d\tilde{p}_{1}d\tilde{p}_{2}$
(50) $\displaystyle=$ $\displaystyle\frac{\tilde{g}}{(2\pi)^{2}}\iint_{K\geq
0}d\varphi_{1}d\varphi_{2}\>(2\pi/\sqrt{3})K(\varphi_{1},\varphi_{2}).\qquad$
(51)

The angular integrals simplify if the rotors can explore the whole of the
configuration torus. This happens if $\hat{E}>4.5\tilde{g}$ because $T\geq 0$
and $0\leq V\leq 4.5$ leading to $N\sim(2\pi/\sqrt{3})(\hat{E}-3\tilde{g})$.
When $\hat{E}<4.5\tilde{g}$, the integral (51) extends only over the
energetically allowed region $V<\hat{E}/\tilde{g}$ of the
$\varphi_{1}-\varphi_{2}$ torus and is evaluated numerically. Fig. 4 shows
that these semiclassical estimates are in good agreement with the numerically
obtained cumulative distribution of levels for $\tilde{g}=300$.

### IV.2 Purification of the spectrum

The Hamiltonian possesses the discrete symmetry
$G=S_{3}\times{\mathbb{Z}}_{2}$ (see II.4). Thus the Hilbert space $\cal H$
spanned by the Fourier modes $e_{m,n}$ may be written as a direct sum of
subspaces that carry irreducible representations of $G$. The Hamiltonian is
block diagonal with respect to this decomposition. To examine spectral
statistics and its universality, we must ‘purify’ the spectrum by isolating
the energy levels corresponding to each irreducible block since the universal
Poisson and Wigner-Dyson spacing distributions expected in integrable and
chaotic systems do not incorporate any such discrete symmetries [7, 20].

#### IV.2.1 Unitary irreducible representations of
$S_{3}\times{\mathbb{Z}}_{2}$

The group $S_{3}$ (see Table 1) has two 1d irreducible representations (i)
identity or trivial: $\rho_{i}(g)=1$ $\forall g\in S_{3}$, (ii) sign:
$\rho_{s}(g)=-1$ for transpositions and $\rho_{s}(g)=1$ for the remaining
permutations and (iii) a 2d unitary representation coming from the dihedral
symmetries of an equilateral triangle centered at the origin with vertices
labeled counterclockwise and with vertex 2 on the upper vertical axis:

$\displaystyle\rho_{d}(e)=\left(\begin{smallmatrix}1&0\\\
0&1\end{smallmatrix}\right),$
$\displaystyle\quad\rho_{d}(\tau_{12})=\frac{1}{2}\left(\begin{smallmatrix}1&\sqrt{3}\\\
\sqrt{3}&-1\end{smallmatrix}\right),$ (52)
$\displaystyle\rho_{d}(\tau_{23})=\frac{1}{2}\left(\begin{smallmatrix}1&-\sqrt{3}\\\
-\sqrt{3}&-1\end{smallmatrix}\right),$
$\displaystyle\quad\rho_{d}(\tau_{31})=\left(\begin{smallmatrix}-1&0\\\
0&1\end{smallmatrix}\right),$ (53)
$\displaystyle\rho_{d}(\sigma_{123})=\frac{1}{2}\left(\begin{smallmatrix}-1&-\sqrt{3}\\\
\sqrt{3}&-1\end{smallmatrix}\right),$
$\displaystyle\quad\rho_{d}(\sigma_{132})=\frac{1}{2}\left(\begin{smallmatrix}-1&\sqrt{3}\\\
-\sqrt{3}&-1\end{smallmatrix}\right).\quad$ (54)

On the other hand, ${\mathbb{Z}}_{2}$ has only two irreducible
representations, both of dimension one: the identity representation
$\varrho_{i}(g)=1$ for all $g\in\mathbb{Z}_{2}$ and the sign representation:
$\varrho_{s}(e)=1$ and $\varrho_{s}(\pi)=-1$. Since the tensor product of
irreducible representations is an irreducible representation of the Cartesian
product (see Section 2.6 of [21]), we get six unitary irreducible
representations of $S_{3}\times\mathbb{Z}_{2}$:

$\displaystyle\text{(i)
}\;\rho_{ii}(g_{1},g_{2})=\rho_{i}(g_{1})\otimes\varrho_{i}(g_{2}),$ (55)
$\displaystyle\text{(ii)
}\;\rho_{is}(g_{1},g_{2})=\rho_{i}(g_{1})\otimes\varrho_{s}(g_{2}),$ (56)
$\displaystyle\text{(iii)
}\;\rho_{si}(g_{1},g_{2})=\rho_{s}(g_{1})\otimes\varrho_{i}(g_{2}),$ (57)
$\displaystyle\text{(iv)
}\;\rho_{ss}(g_{1},g_{2})=\rho_{s}(g_{1})\otimes\varrho_{s}(g_{2}),$ (58)
$\displaystyle\text{(v)
}\;\rho_{di}(g_{1},g_{2})=\rho_{d}(g_{1})\otimes\varrho_{i}(g_{2})\quad\text{and}$
(59) $\displaystyle\text{(vi)
}\;\rho_{ds}(g_{1},g_{2})=\rho_{d}(g_{1})\otimes\varrho_{s}(g_{2})$ (60)

where $g_{1}\in S_{3}$ and $g_{2}\in{\mathbb{Z}}_{2}$. These representations
have dimensions $1,1,1,1,2,2$. The sum of the squares of the dimensions is 12,
which is the same as the cardinality of $S_{3}\times\mathbb{Z}_{2}$. Thus, by
the dimension formula [21], these exhaust all the irreducible representations
of $S_{3}\times\mathbb{Z}_{2}$.

#### IV.2.2 Basis for the identity representation

In what follows, for simplicity, we will restrict attention to the identity
representation $\rho_{ii}$. To find the purified spectrum in this case, we
need to diagonalize $H$ in the invariant subspace of $\cal H$ corresponding to
$\rho_{ii}$. To this end, we seek a basis $s_{m,n}(\varphi_{1},\varphi_{2})$
for the identity representation. A simple way of obtaining such a basis is by
‘averaging’. Given a Fourier basis element $e_{m,n}(\varphi_{1},\varphi_{2})$,
we find its orbit under the action of $S_{3}\times{\mathbb{Z}}_{2}$ (see Table
1) and define $s_{m,n}$ to be the (normalized) sum:

$\displaystyle s_{m,n}$ $\displaystyle=$ $\displaystyle
c_{m,n}^{-1}(e_{m,n}+e_{n-m,n}+e_{m,m-n}+e_{-n,-m}$ (61) $\displaystyle+$
$\displaystyle e_{-n,m-n}+e_{n-m,-m}+e_{-m,-n}+e_{m-n,-n}$ (62)
$\displaystyle+$ $\displaystyle e_{-m,n-m}+e_{n,m}+e_{n,n-m}+e_{m-n,m})$ (63)

where $c_{m,n}$ are normalization constants defined in (65). To enumerate the
independent basis elements we need to restrict to a single $(m,n)$ pair from
each orbit. In Appendix B, we show that this is achieved by restricting to

$m\geq 0\quad\text{and}\quad 0\leq n\leq\lfloor m/2\rfloor,$ (64)

where $\lfloor\cdot\rfloor$ denotes the greatest integer part. We order them
in increasing order of $m$ and in decreasing order of $n$ for fixed $m$:
$s_{0,0}$, $s_{1,0}$, $s_{2,1}$, $s_{2,0}$, $s_{3,1},\cdots$ as indicated in
Fig 9. It is noteworthy that all the Fourier basis states $e_{p,q}$ that
appear in the formula (63) for a particular $s_{m,n}$ are eigenvectors of
kinetic energy $T=\hat{H}(\tilde{g}=0)$ with the same eigenvalue. Thus each
$s_{m,n}$ is an eigenstate of the kinetic energy:
$Ts_{m,n}=(m^{2}+n^{2}-mn)s_{m,n}$. Note that the above order does not
correspond to increasing kinetic energy although $T$ does increase with
decreasing $n$ for fixed $m$. Furthermore, each level of $T$ could contribute
more than one symmetrized basis state $s_{m,n}$. For instance, $s_{7,0}$ and
$s_{8,3}$ have the same kinetic energy $49$.

By virtue of the orthogonality of $e_{m,n}$, the $s_{m,n}$ are also
orthogonal. If the normalization factors are chosen as

$c_{m,n}=\begin{cases}12&\text{if}\quad m=n=0,\\\ 2\sqrt{6}&\text{if}\quad
m\neq n=0,\\\ 2\sqrt{6}&\text{if}\quad n=m/2\neq 0\quad\text{and}\\\
2\sqrt{3}&\text{otherwise,}\end{cases}$ (65)

then $s_{m,n}$ subject to (64) furnish an orthonormal basis for the identity
representation:

$\langle s_{a,b}|s_{m,n}\rangle=\delta^{a}_{m}\delta^{b}_{n}.$ (66)

#### IV.2.3 Matrix elements of $\hat{H}$ in $s_{m,n}$ basis

The matrix elements of $\hat{H}=T+\tilde{g}V$ in the $s_{m,n}$ basis for the
identity representation are

$\displaystyle T_{a,b;m,n}$ $\displaystyle=$
$\displaystyle(m^{2}+n^{2}-mn)\delta^{a}_{m}\delta^{b}_{n}\quad\text{and}$
(67) $\displaystyle V_{a,b;m,n}$ $\displaystyle=$ $\displaystyle
3\delta^{a}_{m}\delta^{b}_{n}-\frac{c^{-1}_{m,n}}{2}\big{(}c_{m+1,n}\delta^{a}_{m+1}\delta^{b}_{n}+c_{m,n-1}\delta^{a}_{m}\delta^{b}_{n-1}$
(70)
$\displaystyle+c_{m-1,n-1}\delta^{a}_{m-1}\delta^{b}_{n-1}+c_{m-1,n}\delta^{a}_{m-1}\delta^{b}_{n}$
$\displaystyle+c_{m,n+1}\delta^{a}_{m}\delta^{b}_{n+1}+c_{m+1,n+1}\delta^{a}_{m+1}\delta^{b}_{n+1}\big{)}.$

In this formula, if $(i,j)$ in $c_{i,j}\delta^{a}_{i}\delta^{b}_{j}$ does not
satisfy (64), then it is to be replaced by the appropriate orbit
representative that does satisfy (64). In other words, we replace

$\displaystyle(0,\pm 1)\to(1,0),\quad(-1,j)\to(1,0),\quad(i,-1)\to(i+1,1)$
(71) $\displaystyle\text{and}\;\;(i,\lfloor i/2\rfloor+1)\to(i,i-\lfloor
i/2\rfloor-1)\;\;\text{when}\;\;i>0.\qquad$ (72)

Taking these replacements into account, we get

$\displaystyle
V_{a,b;m,n}=3\delta^{a}_{m}\delta^{b}_{n}-\begin{cases}3(c_{1,0}/c_{0,0})\delta^{a}_{1}\delta^{b}_{0}&\text{if}\;\;m=n=0,\\\
\frac{c^{-1}_{1,0}}{2}(c_{2,0}\delta^{a}_{2}\delta^{b}_{0}+2c_{2,1}\delta^{a}_{2}\delta^{b}_{1}+2c_{1,0}\delta^{a}_{1}\delta^{b}_{0}+c_{0,0}\delta^{a}_{0}\delta^{b}_{0})&\text{if}\;\;m=1,n=0,\\\
\frac{c^{-1}_{m,0}}{2}(c_{m+1,0}\delta^{a}_{m+1}\delta^{b}_{0}+2c_{m+1,1}\delta^{a}_{m+1}\delta^{b}_{1}+c_{m-1,0}\delta^{a}_{m-1}\delta^{b}_{0}+2c_{m,1}\delta^{a}_{m}\delta^{b}_{1})&\text{if}\;\;m>1,n=0\\\
\frac{c^{-1}_{m,n}}{2}(c_{m+1,n}\delta^{a}_{m+1}\delta^{b}_{n}+c_{m,n-1}\delta^{a}_{m}\delta^{b}_{n-1}+c_{m-1,n-1}\delta^{a}_{m-1}\delta^{b}_{n-1}\cr+c_{m-1,m-n-1}\delta^{a}_{m-1}\delta^{b}_{m-n-1}+c_{m,m-n-1}\delta^{a}_{m}\delta^{b}_{m-n-1}+c_{m+1,m-n}\delta^{a}_{m+1}\delta^{b}_{m-n})&\text{if}\;\;\lfloor
m/2\rfloor=n\neq 0\;\;\text{and}\\\
\frac{c^{-1}_{m,n}}{2}(c_{m+1,n}\delta^{a}_{m+1}\delta^{b}_{n}+c_{m,n-1}\delta^{a}_{m}\delta^{b}_{n-1}+c_{m-1,n-1}\delta^{a}_{m-1}\delta^{b}_{n-1}\cr+c_{m-1,n}\delta^{a}_{m-1}\delta^{b}_{n}+c_{m,n+1}\delta^{a}_{m}\delta^{b}_{n+1}+c_{m+1,n+1}\delta^{a}_{m+1}\delta^{b}_{n+1})&\text{if}\;\;m\geq
4,0<n<\lfloor m/2\rfloor.\end{cases}$ (73)

#### IV.2.4 Purified spectrum in the identity representation

For numerical calculations, we limit $0\leq m\leq m_{\rm max}$. We define a
single index $i(m,n)=\sum_{j=0}^{m}(\lfloor j/2\rfloor+1)-n$ to label the
basis elements in the identity representation. The resulting truncated
Hamiltonian matrix $\hat{H}$ in this basis is of dimension $(\lfloor m_{\rm
max}/2\rfloor+1)(\lfloor m_{\rm max}/2\rfloor+2)$ when $m_{\rm max}$ is odd.
To avoid truncation errors and ensure convergence, we stipulate that the
fractional difference in energy eigenvalues be less than $10^{-10}$ when
$m_{max}$ is increased from $339$ to $349$ corresponding to matrix size
increasing from $28,000$ to $30,000$. This allows us to reliably determine the
lowest $\approx 12,000$ eigenvalues.

## V Spacing distributions

Nearest neighbor energy level spacing distributions are a standard way of
quantifying quantum manifestations of regular or chaotic dynamics [5, 3]. For
integrable and chaotic systems, the spacing distributions of purified levels
are expected to display universal Poisson and Wigner-Dyson statistics [7, 20].
To extract this universal behavior, we ‘unfold’ the purified spectrum to have,
on average, unit level spacing [5, 22]. Denoting energy levels by $E_{i}$ for
$i=1,2,3,\ldots,N$, the unfolded levels are defined as $\xi_{i}=\xi(E_{i})$
where $\xi(E)$ is the polynomial (say of degree 7) that best fits the
cumulative distribution function $n(E)=\sum_{i=1}^{N}\theta(E-E_{i})$ where
$\theta$ is the right-continuous Heaviside step function.

However, the classical three-rotor problem is neither fully chaotic nor
integrable but shows order-chaos-order behavior with increasing energy. To
interpret the spacing distributions, we partition the spectrum into energy
windows where the classical dynamics is regular, mixed or chaotic and do the
unfolding separately in each window. In the semiclassical regime of large
$\tilde{g}\gtrsim 1000$, this partitioning allows us to see striking quantum
signatures of the regular to chaotic and chaotic to regular transitions.

Interestingly, even in the semiclassical regime, we find that the low-lying
energy levels do not follow Poisson statistics although the classical dynamics
is nearly integrable. However, we find that the harmonic approximation of
Sections III.3 and V.2 is able to capture the spacing distribution in this
regime.

Finally, when $\tilde{g}$ is relatively small (e.g. $\tilde{g}=2,10$)
semiclassical expectations are not a useful guide and the order-chaos
transition is barely visible due to a paucity of levels with
$\hat{E}<6\tilde{g}$. However, as shown in Section V.3, we find that the
noninteracting rotor system provides a good approximation, especially at
relatively high energies.

### V.1 Semiclassical regime of large $\tilde{g}$

In the classical three rotor problem [17], we have regular behavior for
$E\lesssim 3.5g$, a transition to chaos for $3.5g\lesssim E\lesssim 5g$
followed by a globally chaotic band $5.33g\lesssim E\lesssim 5.6g$. There is
then a gradual return to regularity as $E\to\infty$. Thus, in the semi-
classical limit of large $\tilde{g}$, we expect to see quantum signatures of
this behavior in spacing distributions in the corresponding energy windows.
Working in the identity representation, we partition the energy spectrum and
investigate the nearest-neighbor spacings of unfolded levels in various
$\hat{E}$ windows. One expects to find an exponential spacing distribution
$P_{\rm e}(s)=e^{-s}$ when the classical behavior is close to integrable and
the Wigner surmise $P_{\rm w}(s)\approx(\pi s/2)e^{-\pi s^{2}/4}$ for
classically chaotic behavior [20, 5]. The latter corresponds to the Gaussian
Orthogonal Ensemble of random matrices that is appropriate to a real symmetric
and time-reversal invariant Hamiltonian [23]. When the classical dynamics is
mixed (islands of regular behavior interspersed with chaotic regions), we
expect to find a spacing distribution that may be crudely modeled by a Brody
distribution [24]

$\displaystyle P(s;\nu)$ $\displaystyle=$
$\displaystyle(1+\nu)f(\nu)s^{\nu}e^{-f(\nu)s^{1+\nu}}$ (74)
$\displaystyle\text{where}\;\;f(\nu)$ $\displaystyle=$
$\displaystyle\Gamma\left((2+\nu)/(1+\nu)\right)^{1+\nu}$ (75)

that interpolates between exponential ($\nu=0$) and Wigner ($\nu=1$).

(a) exact oscillator spacing

(b) oscillator-like

(c) regular: exponential

(d) mixed

(e) chaotic: Wigner

(f) mixed

(g) regular: exponential

(h) regular: exponential

Figure 5: Nearest neighbor spacing histograms (binwidth 0.1 for (a-b) and 0.3
for (c-g)) in various energy windows for $\tilde{g}=6000$ [in (g) and (h) we
take $\tilde{g}=1000$ and $\tilde{g}=300$ due to inadequate high energy data
for $\tilde{g}=6000$]. At low energies the spacing distribution in (b) is a
perturbation to that of a ‘purified’ 2d harmonic oscillator (see Section V.2)
shown in (a) which displays a large number of nearly degenerate levels. In
(c-h) the spacing distributions are compared with exponential (black dotted),
Wigner surmise (red dashed) and Brody (blue dot-dashed) distributions showing
an order-chaos-order transition. The fitted Brody parameter $\nu$ is
indicated. Semiclassical expectations are seen to hold even when $\tilde{g}$
is reduced to $300$.

As shown in Fig. 5, for $\tilde{g}=6000$, we observe (b) 2d harmonic
oscillator-like spacing for $\hat{E}<0.5\tilde{g}$ anticipated in Section
III.3 and further discussed in Section V.2, (c) exponential spacing when
$0.5\tilde{g}\leq\hat{E}\leq 2\tilde{g}$, (d) a Brody-like distribution with
$\nu\approx 0.3$ for $3.5\tilde{g}\leq\hat{E}\leq 4.5\tilde{g}$, (e) Wigner-
Dyson spacing for $5\tilde{g}\leq\hat{E}\leq 6\tilde{g}$ and (f) Brody with
$\nu=0.5$ for $6\tilde{g}\leq\hat{E}\leq 7\tilde{g}$. In (g) and (h) we see a
return to regularity at higher energies ($20\tilde{g}\leq\hat{E}\leq
29\tilde{g}$ when $\tilde{g}=1000$ and $120\tilde{g}<\hat{E}<130\tilde{g}$
when $\tilde{g}=300$). Thus, at least for large $\tilde{g}$, we see a quantum
analogue of the classical order-chaos-order behavior. Notably, while the
classical degree of chaos (as quantified by the fraction of chaos estimated in
Fig. 12 of [17]) showed a rather sharp transition around $E\approx 3.8g$, the
quantum spacing distribution shows a more gradual transition for large
$\tilde{g}$.

### V.2 Harmonic approximation: large $\tilde{g}$ and low $\hat{E}$

As mentioned in Section III.3, for large $\tilde{g}$, the states with low
$\hat{E}$ may be captured by the harmonic approximation (36). Prior to
purification, the harmonic spectrum (37) is

$\displaystyle\frac{\hat{E}^{\rm harm}_{m,n}}{\sqrt{3\tilde{g}}}$
$\displaystyle=$ $\displaystyle m+n+1\quad\text{for}\quad m,n=0,1,2,\ldots$
(76) $\displaystyle=$ $\displaystyle 1,2,2,3,3,3,4,4,4,4,5,5,5,5,5,\cdots.$
(77)

Evidently, the degeneracy is $m+n+1$. By comparing with the numerically
obtained purified spectrum, we empirically deduce that the states in the
identity representation have the energies

$\frac{\hat{E}^{\rm harm}_{\rm
identity}}{\sqrt{3\tilde{g}}}=1^{1},3^{1},5^{1},7^{2},9^{2},11^{2},13^{3},15^{3},17^{3},19^{4},21^{4},\ldots$
(78)

with superscripts denoting degeneracies. In other words, from the full
spectrum (77), we select only levels with ‘odd’ energies with the following
frequency: $\hat{E}/\sqrt{3\tilde{g}}=2i+1$ is chosen
$\lceil\frac{2i+1}{6}\rceil$ times. Here, $\lceil x\rceil$ is the smallest
integer greater than or equal to $x$. The resulting spacing distribution
(after unfolding) shown in Fig. 5(a) accurately captures the spacings obtained
from numerical diagonalization shown in Fig. 5(b).

It would be useful to find a basis for the identity representation of
$S_{3}\times\mathbb{Z}_{2}$ in the harmonic approximation to check our
proposed rule for purification. Prior to purification, the energy
eigenfunctions [with energy $\hat{E}_{m,n}=\sqrt{3\tilde{g}}(m+n+1)$] in this
approximation are

$\displaystyle\psi_{m,n}(\varphi_{+},\varphi_{-})$ $\displaystyle=$
$\displaystyle\psi^{+}_{m}(\varphi_{+})\psi^{-}_{n}(\varphi_{-})\quad\text{where}$
(79) $\displaystyle\psi^{+}_{m}(\varphi_{+})$ $\displaystyle\propto$
$\displaystyle
H_{m}((3\tilde{g}/4)^{1/4}\varphi_{+})\;e^{-\sqrt{3\tilde{g}/4}\varphi_{+}^{2}}\quad\text{and}$
(80) $\displaystyle\psi^{-}_{n}(\varphi_{-})$ $\displaystyle\propto$
$\displaystyle
H_{n}((\tilde{g}/12)^{1/4}\varphi_{-})\;e^{-\sqrt{\tilde{g}/12}\varphi_{-}^{2}}.$
(81)

Here, $\varphi_{\pm}$ are viewed as taking values on the real line although
they arose as angles. This is a reasonable approximation since the Gaussian
factors ensure that $\psi^{+}_{m}$ and $\psi^{-}_{n}$ decay rapidly. $H_{n}$
is the $n^{\rm th}$ Hermite polynomial. The exponent in the Gaussian factor in
$\psi_{m,n}$ is invariant under the symmetry group
$G=S_{3}\times\mathbb{Z}_{2}$ since it is proportional to the potential
$V(\varphi_{+},\varphi_{-})$ (36). Hence, the desired basis consists of those
linear combinations of $H_{m}H_{n}$ that are invariant under $G$. One may find
them by group averaging: add to $H_{m}H_{n}$ all its translates under $G$.
Since $G$ is a symmetry of the Hamiltonian, the orbit of $H_{m}H_{n}$ consists
of states with the same energy. From (78), we expect that

  1. (i)

for even $\hat{E}/\sqrt{3\tilde{g}}$, this average vanishes so that there is
no invariant combination, while

  2. (ii)

for odd $\hat{E}/\sqrt{3\tilde{g}}$ there should be
$\left\lceil\hat{E}/6\sqrt{3\tilde{g}}\right\rceil$ independent invariant
energy eigenstates.

We now furnish a proof of (i). The action of $S_{3}\times\mathbb{Z}_{2}$ on
$\varphi_{\pm}$ can be inferred from Table 2. Given the action of
$\tau_{31}\in S_{3}$ and $\pi\in\mathbb{Z}_{2}$, the translates of
$(\varphi_{+},\varphi_{-})$ may be organized as six pairs. The two translates
in each pair differ only in the sign of the image of $\varphi_{+}$:

$\displaystyle(\pm\varphi_{+},\varphi_{-}),\quad(\pm\varphi_{+},-\varphi_{-}),$
(82)
$\displaystyle(\pm(\varphi_{+}-\varphi_{-})/2,(3\varphi_{+}+\varphi_{-})/2),$
(83)
$\displaystyle(\pm(\varphi_{+}-\varphi_{-})/2,-(3\varphi_{+}+\varphi_{-})/2),$
(84)
$\displaystyle(\pm(\varphi_{+}+\varphi_{-})/2,(3\varphi_{+}-\varphi_{-})/2),$
(85)
$\displaystyle(\pm(\varphi_{+}+\varphi_{-})/2,-(3\varphi_{+}-\varphi_{-})/2).$
(86)

Evidently, they may also be arranged so that each pair differ only in the sign
of the image of $\varphi_{-}$. Now we consider the group average of

$\psi_{m,n}(\varphi_{+},\varphi_{-})=\psi^{+}_{m}(\varphi_{+})\psi^{-}_{n}(\varphi_{-}).$
(87)

Suppose $m$ is odd, so that $\psi^{+}_{m}$ is an odd function of
$\varphi_{+}$. Then the two terms from each pair in the above $G$-orbit
cancel, leading to a vanishing group average. Analogously, the group average
vanishes if $n$ is odd. Thus, there is no nontrivial $G$-invariant state in
the $\hat{E}^{\rm harm}_{m,n}$ eigenspace if either $m$ or $n$ is odd. From
(77), for $\hat{E}_{m,n}/\sqrt{3\tilde{g}}$ to be even, either $m$ or $n$ must
be odd. We thus conclude that in the identity representation, there cannot be
any energy levels with even $\hat{E}/\sqrt{3\tilde{g}}$ within the harmonic
approximation.

Although we have not yet found a general argument for (ii), we have verified
our conjecture on degeneracies for low-lying levels with odd
$\hat{E}/\sqrt{3\tilde{g}}$ in the identity representation. The first three
levels $\hat{E}/\sqrt{3\tilde{g}}=1,3,5$ are nondegenerate as we now show. The
ground state has a constant $G$-invariant eigenfunction $\propto
H_{0}H_{0}=12$. The group averages of $H_{2}H_{0}$ and $H_{0}H_{2}$ are both
equal to $4(3\sqrt{3}\varphi_{+}^{2}+\sqrt{3}\varphi_{-}^{2}-6)$ while the
group averages of $H_{4}H_{0},H_{2}H_{2}$ and $H_{0}H_{4}$ produce only one
independent $G$-invariant state given by 6, 2 and 6 times

$(24-8\sqrt{3}\varphi_{-}^{2}+\varphi_{-}^{4}-24\sqrt{3}\varphi_{+}^{2}+6\varphi_{-}^{2}\varphi_{+}^{2}+9\varphi_{+}^{4}).$
(88)

For the fourth level with $\hat{E}=7\sqrt{3\tilde{g}}$, we find that the group
averages of $H_{0}H_{6}\pm H_{6}H_{0}$ are equal to those of
$H_{4}H_{2}-H_{2}H_{4}$ and $5(H_{4}H_{2}+H_{2}H_{4})$, leading to two
linearly independent states, as expected from (78). It would of course be
satisfying to extend these special cases to a proof of the proposed degeneracy
formula $\left\lceil\hat{E}/6\sqrt{3\tilde{g}}\right\rceil$ for odd
$\hat{E}/\sqrt{3\tilde{g}}$.

### V.3 Spacings at low and moderate $\tilde{g}$

(a) free rotors: $\tilde{g}=0$

(b) low $\tilde{g}$, high $\hat{E}$

(c) moderate $\tilde{g}$, high $\hat{E}$

(d) moderate $\tilde{g}$ and $\hat{E}$

(e) $\tilde{g}=100$, moderate $\hat{E}$

(f) $\tilde{g}=100$, low $\hat{E}$

Figure 6: Exact free rotor ($\tilde{g}=0$) spacing in the identity
representation shown in (a) captures that of the highly excited levels for low
$\tilde{g}=2$ shown in (b). For $\tilde{g}=10$, the highly excited levels in
(c) show some departures from free rotor spacings while the moderately excited
levels in (d) show features of both exponential and free rotor spacings. When
$\tilde{g}=100$, spacings of moderately excited states in (e) is nearly
exponential while those of intermediate levels in the classically chaotic
regime is Wigner-like as shown in (f).

When $\tilde{g}$ is smaller, the main change is that there are very few energy
levels in the range $0\leq\hat{E}\leq 6\tilde{g}$ making it difficult to use
spacing distributions to examine the quantum analog of the classical
transition from regularity to chaos. For example, there are only 14 levels for
$\tilde{g}=10$ and 105 levels for $\tilde{g}=100$ in this energy window. On
the other hand, when $\hat{E}>25\tilde{g}$, where the classical behavior is
relatively regular, we may expect spacing distributions to display Poisson
statistics. However, there can be departures from this expectation for at
least two reasons: (i) the semiclassical approximation may break down for
small $\tilde{g}$ and (ii) when $\hat{E}\gg\tilde{g}$, the dynamics should be
a small perturbation to free rotor motion, for which the spacing distribution
is different from exponential. In fact, the free rotor ($\tilde{g}=0$) energy
eigenvalues in the identity representation are given by
$\hat{E}_{m,n}=m^{2}+n^{2}-mn$ for $m,n$ satisfying (64) with eigenfunctions
$s_{m,n}(\varphi_{1},\varphi_{2})$ leading to the spacing distribution shown
in Fig. 6(a). Pleasantly, for $\tilde{g}=2$, the spacings of highly excited
unfolded levels shown in Fig. 6(b) closely mimics the free rotor spacings. For
$\tilde{g}=10$, departures from free rotor at high energies are more
pronounced as seen in Fig. 6(c). By contrast, at lower energies in the
classically regular regime, we begin to observe signatures of the exponential
spacing distribution in addition to that of free rotors. This is visible in
Fig. 6(d) for $\tilde{g}=10$. However, when $\tilde{g}$ is increased to 100,
the free rotor effects seem to become insignificant for these energies with
spacing distributions approaching exponential as seen in Fig. 6(e). Moreover,
at $\tilde{g}=100$, we also begin to see Wigner-like spacing distributions in
the classically chaotic regime (see Fig. 6(f)).

## VI Number variance

Let $d(\xi)=\sum_{i=1}^{N}\delta(\xi-\xi_{i})$ denote the spectral density
function for the unfolded energy levels $\xi_{1},\ldots,\xi_{N}$. Then the
number of levels in a spectral window $[\xi-L/2,\xi+L/2]$ is given by

$n(\xi,L)=\int_{\xi-L/2}^{\xi+L/2}d(\xi^{\prime})\>d\xi^{\prime}.$ (89)

The number average $\langle n(\xi,L)\rangle$ is defined as the mean value of
$n(\xi,L)$ as $\xi$ ranges over a desired portion of the spectrum. In
practice, this ensemble average is performed by summing over an equally spaced
grid of values of $\xi$ (with a spacing of $0.1$ for $L>2$ and $0.05$ for
$L\leq 2$). Since the unfolded spectrum is constructed to have approximately
unit spacing on average, $\langle n(\xi,L)\rangle\approx L$. Fluctuations in
$n(\xi,L)$ are captured by the number variance [23, 3] defined as

$\Sigma(L)=\langle[n(\xi,L)-\langle n(\xi,L)\rangle]^{2}\rangle.$ (90)

(a)

(b)

Figure 7: Number variance $\Sigma(L)$ in classically regular regime of large
$\hat{E}\gg\tilde{g}$. (a) $\Sigma(L)$ shows linear growth for small $L$ with
a universal slope of unity. (b) For larger $L$, the number variance saturates
and oscillates. The universal linear behavior persists up to higher values of
$L$ at higher energies, where the dynamics approaches integrability.

(a)

(b)

Figure 8: Number variance for 1836 levels in the classically chaotic energy
band $(5\tilde{g}<\hat{E}<6\tilde{g})$ for $\tilde{g}=6000$ displays GOE
asymptotics for moderately small $L$ ($0.1\leq L\leq 4$) and then displays
oscillations for larger $L$.

To interpret the number variance, we evaluate it for portions of the three-
rotor spectrum in two representative regimes: (a) $1836$ levels with
$5\tilde{g}<\hat{E}<6\tilde{g}$ for $\tilde{g}=6000$ and (b) 921, 914 and 913
levels when $\hat{E}\in[30\tilde{g},40\tilde{g}$], $[70\tilde{g},80\tilde{g}]$
and $[120\tilde{g},130\tilde{g}]$ for $\tilde{g}=300$ where level spacing
distributions in Fig. 5 of Section V.1 indicate that the quantum dynamics are
predominantly chaotic and regular respectively. It is conjectured that in
these regimes, the number variance should show universal behavior for small
$L$ that is captured by the Gaussian orthogonal ensemble (GOE) of random
matrices and a Poisson process respectively. For the GOE, Mehta (see Appendix
A38 of [23]) derives the asymptotic behavior

$\Sigma(L)\sim\frac{2}{\pi^{2}}(\log(2\pi L)+\gamma+1-\pi^{2}/8)+{\cal
O}(L^{-1}),$ (91)

while for a Poisson process, $\Sigma(L)=L$. Note that, though Mehta’s
logarithmic GOE formula is derived for asymptotically large $L$, while
comparing with our number variance, it may be expected to hold in an
intermediate regime: the ${\cal O}(L^{-1})$ remainder term dominates when
$L\to 0$ while nonuniversal effects can become important for large $L$. We
compare the number variance in the above regular and chaotic regimes against
these expectations in Figs. 7 and 8. Although universal behavior is manifested
for small/moderately small $L$, we find that in both the regimes, the number
variance saturates and displays oscillations whose character depends on the
details of the three rotor system. In fact, it should be possible to capture
these oscillations in a semiclassical approximation by examining the short
periodic orbits of the three-rotor system [12].

## VII Discussion

Here, we discuss our results while highlighting some novel features of the
quantum three-rotor problem and mention some avenues for future research.

To begin with, the formulation of the quantum three-rotor problem involved
some subtleties: singlevaluedness of the wave function allowed for nontrivial
boundary conditions on the relative wave function (labelled by the cube roots
of unity) with consequences for the center of mass angular momentum quantum
number as well. The presence of a dimensionless coupling constant $\tilde{g}$
distinguished the quantum theory from its classical counterpart. This was
exploited in organizing the quantum theory into weak and strong coupling
regimes with semiclassical behavior observed at strong coupling allowing us to
interpret subsequent results on spectral statistics. However, even for large
$\tilde{g}$, making contact with universal semiclassical expectations in
spectral statistics was complicated by the need to partition the spectrum into
energy windows where the classical dynamics is regular, mixed or chaotic. By
contrast, in models such as the planar elastic pendulum [25, 26], which also
displays order-chaos-order behavior, there is a control parameter (other than
energy) whose variation shifts the entire spectrum from Poisson to Wigner-
Dyson.

Another interesting feature of the three-rotor model is the presence of an
$S_{3}\times\mathbb{Z}_{2}$ discrete symmetry of the relative Hamiltonian.
Although this symmetry is manifest in the rotor angle description, it is not
obvious while working with the relative or Jacobi angles. A key step was to
identify the implementation of this symmetry on relative angles in order to be
able to decompose the energy spectrum into irreducible representations where
universal behavior may be expected. In fact, incorrect purification via a more
obvious $\mathbb{Z}_{2}\times\mathbb{Z}_{2}$ symmetry would lead to spurious
Shnirelman-like peaks [27, 28] in spacing distributions.

Group averaging allowed us to construct a convenient basis for the subspace of
the Hilbert space that carries the identity representation of
$S_{3}\times\mathbb{Z}_{2}$, allowing us to purify the spectrum. However, we
suspect that there is a more elegant way of choosing orbit representatives to
label the basis elements. The analogous problem of finding suitable bases for
the remaining five irreducible representations and qualitative features of the
corresponding energy spectra has not been addressed in this paper. Pleasantly,
for the identity representation, we have a closed-form sparse matrix
representation of the Hamiltonian. Moreover, we demonstrate that numerically
determined eigenvalues are insensitive to an increase in the size of the
truncated Hamiltonian ensuring that spacing distributions are computed to
desired accuracy.

Remarkably, the above $S_{3}\times\mathbb{Z}_{2}$ symmetry continues to apply
even at the level of the harmonic approximation, which allowed us to use the
latter to predict the spacing distributions at low energies where there are
marked departures from the semiclassical Poisson distribution. However, while
we have conjectured which harmonic levels belong to the identity
representation (along with their degeneracies), we have not yet found general
expressions for suitable basis elements for the invariant subspace in the
harmonic approximation (unlike in the full problem where we have explicit
formulae for an orthonormal basis).

On the other hand, we showed that the free rotor approximation captures
deviations from Poisson statistics at asymptotically high energies and weak
coupling. Interestingly, finding the degeneracies of the free rotor spectrum
led us to a number theoretic problem involving the Eisenstein integers, whose
solution could be expressed in terms of the prime factorization of the energy
in suitable units.

There are several other interesting directions for further research such as to
understand (i) the nature of level crossings and avoided crossings, (ii)
higher order spacings and spacing ratios, especially as a tool for
purification of the spectrum [29] in other representations, (iii) the behavior
of wave functions, Husimi functions and their nodal curves in the regular and
chaotic regimes, (iv) the saturation and oscillation of the number variance
$\Sigma(L)$ for large $L$ using our knowledge of classical periodic orbits,
and (v) the connection of out-of-time-ordered correlators (OTOC) to classical
Lyapunov exponents.

###### Acknowledgements.

We thank M V Berry, A Lakshminarayan and M S Santhanam for helpful comments
and references and G D Patil for guiding us through the number theory behind
the degeneracy formula for free rotor energy levels given in Section III.2 and
Appendix A. Finally, although Ankit Yadav did not wish to be included as an
author, we are grateful to him for his valuable contributions at an
intermediate stage of this work. This work was supported in part by the
Infosys Foundation and grants (MTR/2018/000734, CRG/2018/002040) from the
Science and Engineering Research Board, Govt. of India.

##  Appendix A Free rotor degeneracy and Eisenstein primes

Recall from Section III.2 that the spectrum of the free rotor Hamiltonian
$\hat{H}=T$ is given by $m^{2}+n^{2}-mn$ where $m$ and $n$ are integers. Here,
we outline a proof of formula (30) for the degeneracy of free rotor energy
levels: if the prime factorization of the energy eigenvalue is of the form
$T=3^{\gamma}p_{1}^{\alpha_{1}}\cdots
p_{r}^{\alpha_{r}}q_{1}^{\beta_{1}}\cdots q_{s}^{\beta_{s}}$, where
$p_{i}\equiv 1$ $(\text{mod~{}}3)$ and $q_{i}\equiv 2$ $(\text{mod~{}}3)$,
then the number of integer solutions to $m^{2}+n^{2}-mn=T$ is

$\text{degen}(T)=6(\alpha_{1}+1)(\alpha_{2}+1)\cdots(\alpha_{r}+1)$ (92)

if all $\beta_{i}$ are even and $\text{degen}(T)=0$ otherwise.

To begin with, the factorization

$\displaystyle m^{2}+n^{2}-mn$ $\displaystyle=$
$\displaystyle|m+n\omega|^{2}=(m+n\omega)(m+n\omega^{2})$ (93)
$\displaystyle=$ $\displaystyle(m+n\omega)((m-n)-m\omega),$ (94)

where $\omega$ is a primitive cube root of unity suggests that the integral
lattice spanned by $1$ and $\omega$ is relevant to this problem. Thus, let
$\omega$ denote a complex number satisfying $\omega^{2}+\omega+1=0$ and let
$L=\mathbb{Z}[\omega]=\mathbb{Z}+\mathbb{Z}\omega=\\{\>m+n\omega|m,n\in\mathbb{Z}\\}$
be the triangular lattice in the complex plane displayed in Fig. 2. The set
$L$ is closed under addition and multiplication and forms a ring called the
Eisenstein integers, the ring of integers in $\mathbb{Q}(\omega)$, which is an
imaginary quadratic (same as $\mathbb{Q}+\mathbb{Q}(\sqrt{-3})$) and
cyclotomic number field [30]. The labels in Fig. 2 correspond to the choice
$\omega=e^{2\pi i/3}$, although one gets the same triangular lattice upon
choosing $\omega=e^{4\pi i/3}$.

Let us start with some terminology and basic facts about $L$. (i) Given
$r,s\in L$, we say that $s$ divides $r$ (denoted $s|r$) if there is some $q\in
L$ such that $r=sq$. (ii) The elements of the set

$\displaystyle U$ $\displaystyle=$
$\displaystyle\\{1,-\omega^{2}=\omega+1,\omega,-1,\omega^{2}=-\omega-1,-\omega\\}$
(95) $\displaystyle=$ $\displaystyle\\{(-\omega)^{k},k=0,1,\ldots,5\\}$ (96)

lying at the vertices of the hexagon marked in Fig. 2 of the lattice
constitute the group of units of $L$. They behave like $\pm 1$ in the ring of
integers $\mathbb{Z}$. More precisely, if $xy=1$ for $x,y\in L$ then $x,y\in
U$. Furthermore, units are the only elements of $L$ of unit magnitude: $|x|=1$
implies $x\in U$. (iii) An element $q$ that is not a unit ($q\in L\setminus
U$) is called irreducible if in any factorization $q=ab$, at least one of the
factors is a unit ($a\in U$ or $b\in U$). On the other hand, $q\in L\setminus
U$ is called a prime if whenever $q|rs$ either $q|r$ or $q|s$. While primes
are always irreducible, one may show that all irreducible elements of $L$ must
be primes as $L$ is a Unique Factorization Domain (see below).

Using these definitions, we state the following division algorithm. Given a
dividend $a\in L$ and a nonzero divisor $b\in L$, there are uniquely
determined quotient and remainder elements $q,r\in L$ with $a=bq+r$ where
$0\leq|r|\leq|b/2|$ or $r=b/2$. In other words, $L$ is a norm Euclidean ring.
Consequently, $L$ is a Principal Ideal Domain and hence a Unique Factorization
Domain (Corollary 4.4 of [30]). Moreover, from Proposition 4.7 of [30], the
primes in $L$ (called Eisenstein primes) must be of one of the following sorts
up to units: (i) $1+2\omega$ or (ii) $m+n\omega$ or $m+n\omega^{2}$ where
$m^{2}+n^{2}-mn$ is an ordinary prime $\equiv 1$ $(\text{mod~{}}3)$ in
$\mathbb{Z}$ or (iii) an ordinary prime $q\equiv 2$ $(\text{mod~{}}3)$ in
$\mathbb{Z}$.

We now return to the integers and explain the degeneracy formula (92). For
given kinetic energy eigenvalue $T$, suppose $(m,n)$ is one integer pair such
that $T=m^{2}+n^{2}-mn=|m+n\omega|^{2}$. Then since the six units $u\in U$
have unit magnitude, $T=|u(m+n\omega)|^{2}$ where
$u(m+n\omega)=m^{\prime}+n^{\prime}\omega\in L$. Thus we get five additional
distinct integer pairs with the same eigenvalue $T$. This explains the
prefactor $6$ in (92).

It remains to find the contribution to the degeneracy that does not arise from
the freedom to multiply by units. To do this, we will exploit the prime
factorization (in $\mathbb{Z}$) $T=3^{\gamma}p_{1}^{\alpha_{1}}\cdots
p_{r}^{\alpha_{r}}q_{1}^{\beta_{1}}\cdots q_{s}^{\beta_{s}}$, where
$p_{i}\equiv 1$ $(\text{mod~{}}3)$ and $q_{i}\equiv 2$ $(\text{mod~{}}3)$. The
general idea is that any representation of $T$ as
$|m+n\omega|^{2}=(m+n\omega)(m+n\omega^{2})$ must arise from a product of such
representations of the individual prime factors. The paired appearance of
$m+n\omega$ and $m+n\omega^{2}$ is because $\omega^{2}$ is the complex
conjugate of $\omega$. This is analogous to the nonreal roots of a polynomial
with real coefficients appearing in complex conjugate pairs. We will now deal
sequentially with the three sorts of primes in $\mathbb{Z}$: $3$, $p\equiv 1$
$(\text{mod~{}}3)$ and $q\equiv 2$ $(\text{mod~{}}3)$.

(a) We begin with the ordinary prime $3$ and note that
$3=|1+2\omega|^{2}=(1+2\omega)(1+2\omega^{2})=-(1+2\omega)^{2}$. Thus although
$3$ is a prime in $\mathbb{Z}$, it is not a prime in $L$. In fact, $3$
ramifies in $L$ with ramification index two since it is a unit times the
square of the prime $1+2\omega$ in $L$. Moreover, up to multiplication by
units, there is only one way of expressing $3$ as $|m+n\omega|^{2}$. Thus,
even if the prime factorization of $T$ includes $3^{\gamma}$ for some
$\gamma\geq 1$, there is only one way to write $3^{\gamma}$ as
$(m+n\omega)\times(m+n\omega^{2})$:
$(1+2\omega)^{\gamma}\times(1+2\omega^{2})^{\gamma}$. Hence, it will only
contribute a factor of one to $\text{degen}(T)$.

(b) Suppose $p\equiv 1$ $(\text{mod~{}}3)$ is a prime in $\mathbb{Z}$, then
there exist integers $m,n$ such that $p=m^{2}+n^{2}-mn$. For example, $p=7,13$
and $19$ are primes $\equiv 1$ $(\text{mod~{}}3)$ with
$7=|1+3\omega|^{2}=|3+\omega|^{2}$, $13=|1+4\omega|^{2}=|4+\omega|^{2}$ and
$19=|2+5\omega|^{2}=|5+2\omega|^{2}$. An inductive argument can be used to
show that the same applies to any other such prime. Moreover, in each case,
there are precisely two representations
$p=|m+n\omega|^{2}=|\tilde{m}+\tilde{n}\omega|^{2}$ that do not differ by a
unit but are related by $\tilde{m}+\tilde{n}\omega=(-\omega)^{k}(n+m\omega)$
for some $0\leq k\leq 5$. This is a reflection of the Galois symmetry
$\omega\to\omega^{2}$ that matches the complex conjugation $i\to-i$. More
generally, if $T=p^{\alpha}$ where $\alpha\geq 1$ and $p\equiv 1$
$(\text{mod~{}}3)$, then (aside from the freedom to multiply by units) there
are $\alpha+1$ distinct ways of writing $T=|m+n\omega|^{2}$. To see this, we
denote $A=m+n\omega$ and $B=n+m\omega$. Then
$p^{2}=|A^{2}|^{2}=|AB|^{2}=|B^{2}|^{2}$ and in general,
$p^{\alpha}=|A^{\alpha}|^{2}=|A^{\alpha-1}B|^{2}=\cdots=|AB^{\alpha-1}|^{2}=|B^{\alpha}|^{2}$.

(c) Finally, suppose $q$ is a prime number in $\mathbb{Z}$ and $q\equiv 2$
$(\text{mod~{}}3)$. Then $q$ is also an Eisenstein prime. Thus, if
$T=q^{\beta}$, then: (i) if $\beta$ is odd, there is no way of writing $T$ as
$|m+n\omega|^{2}$ but (ii) if $\beta$ is even, there is only way of writing
$T$ as $|m+n\omega|^{2}$ i.e., by taking $m+n\omega=q^{\beta/2}$ up to a unit.
This generalizes to $T=q_{1}^{\beta_{1}}\cdots q_{s}^{\beta_{s}}$.

Combining these, we arrive at the degeneracy formula (92).

##  Appendix B Basis labels for identity representation

As mentioned in Table 1, the symmetry group $G$ acts on the Fourier basis
states $e_{m,n}$ which may be viewed as points on the lattice of ordered pairs
of integers $(m,n)$. The orbit of $(m,n)$ is

$\displaystyle{\cal O}_{(m,n)}$ $\displaystyle=$
$\displaystyle\\{(m,n),(n-m,n),(m,m-n),(-n,-m),$ (99)
$\displaystyle(-n,m-n),(n-m,-m),(-m,-n),(m-n,-n),$
$\displaystyle(-m,n-m),(n,m),(n,n-m),(m-n,m)\\}.$

It is typically of length twelve except in the extreme cases [boundary of the
wedge defined in (100) and displayed in Fig. 9] when (a) $m=n=0$ and it
shrinks to a single point and when (b) ($m>0$, $n=0$ or $\lfloor m/2\rfloor$)
and it has length six. Now, we show that each orbit can be labelled uniquely
by a representative $(m,n)$ lying in the ‘wedge’ defined by

$m\geq 0\quad\text{and}\quad 0\leq n\leq\lfloor m/2\rfloor,$ (100)

where $\lfloor\cdot\rfloor$ is the greatest integer part. The proof proceeds
in two steps: (i) We use the group action to argue that every orbit enters the
wedge and (ii) show that an orbit cannot have more than one representative in
the wedge.

$m$$n$(3,1)(4,2)(1,2)(1,0)(5,2)= basis for identity representation= ordering
of basis elements$\text{{=} orbit of }(1,0):\mathcal{O}_{(1,0)}$$\text{{=}
orbit of }(2,1):\mathcal{O}_{(2,1)}$ Figure 9: Basis labels for identity
representation of $S_{3}\times\mathbb{Z}_{2}$ are orbit representatives
$(m,n)$ (denoted by black dots) lying in a wedge in the
$\mathbb{Z}\oplus\mathbb{Z}$ lattice. Each orbit intersects the wedge exactly
once.

Proof of (i). Of the points $(m,n)$ and $(-m,-n)$ lying in ${\cal O}_{(m,n)}$
[see (99)] one must lie in the right half plane. Without loss of generality we
suppose that $(m,n)$ lies in the right half plane, i.e., $m\geq 0$. Next, one
of the two orbit elements $(m,n)$ or $(m-n,-n)$ must lie in the upper half
plane. Taking $(m,n)$ to lie in the upper half plane, we get $n\geq 0$
allowing us to restrict $(m,n)$ to the first quadrant. In a similar manner,
either $(m,n)$ or $(n,m)$ must lie in the south-east half plane allowing us to
restrict to the $\pi/4$ wedge: $m\geq n\geq 0$. Finally, either $(m,m-n)$ or
$(m,n)$ must lie in the $\pi/8$ wedge of Eqn. (100). Thus we have shown that
every orbit has a representative $(m,n)$ satisfying (100). The Proof of (ii)
will show that there is only one such representative.

Proof of (ii). We will now argue that if $(m,n)$ satisfies (100) then all
other points on ${\cal O}_{(m,n)}$ either violate (100) or coincide with
$(m,n)$. Let us illustrate this with the point $(n-m,n)\in{\cal O}_{(m,n)}$.
If this point lies in the wedge, then (100) implies $n-m\geq 0$. But since
$(m,n)$ already lies in the wedge we must also have $n\leq\lfloor m/2\rfloor$.
Both conditions are satisfied only when $m=n=0$, in which case both points
coincide: $(m,n)=(n-m,n)=(0,0)$. A similar argument applies to all other
points in ${\cal O}_{m,n}$. This establishes that every orbit has a unique
representative in the wedge defined by (100). Thus, we may use $(m,n)$
satisfying (100) to label a basis for the identity representation of $G$.

## References

  * [1] M C Gutzwiller, Chaos in Classical and Quantum mechanics, Springer-Verlag, New York (1990).
  * [2] M S Santhanam, S Paul and J Bharathi Kannan, Quantum kicked rotor and its variants: Chaos, localization and beyond, Physics Reports, 956, 1 (2022).
  * [3] H J Stöckmann, Quantum chaos, an introduction, Cambridge University Press, Cambridge (1999).
  * [4] N L Balazs and A Voros, Chaos on the pseudosphere, Physics Reports, 143 (3), 109 (1986).
  * [5] F Haake, Quantum signatures of chaos, Springer-Verlag, Berlin (2010).
  * [6] J H Lowenstein, Essentials of Hamiltonian Dynamics, Cambridge University Press, Cambridge (2012).
  * [7] O Bohigas, M-J Giannoni and C Schmit, Characterization of Chaotic Quantum Spectra and Universality of Level Fluctuation Laws, Phys. Rev. Lett., 52, 1 (1984).
  * [8] M Tabor, Chaos and Integrability in Nonlinear Dynamics: An Introduction, Wiley-Interscience (1989).
  * [9] S R Jain and R Samajdar, Nodal portraits of quantum billiards: Domains, lines, and statistics, Rev. Mod. Phys., 89, 045005 (2017).
  * [10] A Larkin and Y N Ovchinnikov, Quasiclassical method in the theory of superconductivity, Sov. Phys. JETP, 28, 1200 (1969).
  * [11] S H Shenker and D Stanford, Black holes and the butterfly effect, J. High Energy Phys., 2014 (3), 67 (2014).
  * [12] M V Berry, Semi-classical theory of spectral rigidity, Proc. R. Soc. Lond. A, 400, 229 (1985).
  * [13] S L Sondhi, S M Girvin, J P Carini and D Shahar, Continuous quantum phase transitions, Rev. Mod. Phys., 69 (1), 315 (1997).
  * [14] A Shnirman, G Schoen and Z Hermon, Quantum Manipulations of Small Josephson Junctions, Phys. Rev. Lett., 79, 2371 (1997).
  * [15] T P Orlando, J E Mooij, L Tian, C H van der Wal, L S Levitov, S Lloyd and J J Mazo, Superconducting Persistent-Current Qubit, Phys. Rev. B, 60, 15398 (1999).
  * [16] G S Krishnaswami and A Yadav, Bifurcation cascade, self-similarity and duality in the 3-rotor problem, Chaos, 33, 083101 (2023).
  * [17] G S Krishnaswami and H Senapati, Classical three rotor problem: periodic solutions, stability and chaos, Chaos, 29 (12), 123121 (2019).
  * [18] G S Krishnaswami and H Senapati, Ergodicity, mixing and recurrence in the three rotor problem, Chaos, 30 (4), 043112 (2020).
  * [19] H Weyl, Das asymptotische Verteilungsgesetz der eigenwerte linearer partieller Differentialgleichungen (mit einer Anwendung auf die Theorie der Hohlraumstrahlung, Math Ann., 71, 441, (1912).
  * [20] M V Berry and M Tabor, Level clustering in the regular spectrum, Proc. R. Soc. Lond. A, 356, 375 (1977).
  * [21] S Sternberg, Group theory and physics, Cambridge University Press, Cambridge (1995).
  * [22] G S Krishnaswami and T R Vishnu, Screwon spectral statistics and dispersion relation in the quantum Rajeev-Ranken model, Physica D, 463, 134170 (2024).
  * [23] M L Mehta, Random Matrices, $3^{\rm rd}$ Ed., Elsevier, Amsterdam (2004).
  * [24] T A Brody, A statistical measure for the repulsion of energy levels, Lett. Nuovo Cimento, 7, 482 (1973).
  * [25] Anurag, B Mondal, J K Bhattacharjee and S Chakraborty, Understanding the order-chaos-order transition in the planar elastic pendulum, Physica D, 402, 132256 (2020).
  * [26] Anurag, B Mondal, T Shah and S Chakraborty, Chaos and order in librating quantum planar elastic pendulum, Nonlinear Dyn., 103, 2841 (2021).
  * [27] A I Shnirelman, The asymptotic multiplicity of the spectrum of the Laplace operator (in Russian), Uspekhi Mat. Nauk., 30:4(184), 265 (1975).
  * [28] B V Chirikov and D L Shepelyansky, Shnirelman peak in level spacing statistics, Phys. Rev. Lett., 74, 518 (1995).
  * [29] S H Tekur, U T Bhosale and M S Santhanam, Higher order spacing ratios in random matrix theory and complex quantum systems, Phys. Rev. B, 98, 104305 (2018).
  * [30] D A Cox, Primes of the Form $x^{2}+ny^{2}$: Fermat, Class Field Theory and Complex Multiplication, 2${}^{\text{nd}}$ Ed., John Wiley & Sons, Hoboken (2013).