Integrable families of hard-core particles with unequal masses in a one-dimensional harmonic trap

TL;DR

This paper identifies specific arrangements of unequal masses in a one-dimensional harmonic trap that lead to integrable and superintegrable dynamics, classified by Coxeter groups and solvable via Bethe ansatz-like methods.

Contribution

It introduces new integrable mass families in a harmonic trap with hard-core particles, classified by Coxeter groups, and suggests superintegrability and maximal superintegrability.

Findings

Existence of integrable mass families for any number of particles.

Classification of integrable families by Coxeter reflection groups.

Evidence and conjecture of superintegrability and maximal superintegrability.

Abstract

We show that the dynamics of particles in a one-dimensional harmonic trap with hard-core interactions can be solvable for certain arrangements of unequal masses. For any number of particles, there exist two families of unequal mass particles that have integrable dynamics, and there are additional exceptional cases for three, four and five particles. The integrable mass families are classified by Coxeter reflection groups and the corresponding solutions are Bethe ansatz-like superpositions of hyperspherical harmonics in the relative hyperangular coordinates that are then restricted to sectors of fixed particle order. We also provide evidence for superintegrability of these Coxeter mass families and conjecture maximal superintegrability.

Figures

Tables

Table 1. Table 1: This table provides information about the connected, non-branching, finite Coxeter reflection groups. For N 𝑁 N particles, each rank m = N − 1 𝑚 𝑁 1 m=N-1 Coxeter group defines a one-parameter family of masses for which the system is exactly solvable. For each group G m subscript 𝐺 𝑚 G_{m} , the following data are provided  Goodman ( 2004 ) : the Coxeter bracket [ q 1 , … , q m ] subscript 𝑞 1 … subscript 𝑞 𝑚 [q_{1},\ldots,q_{m}] from which one determines the angles of the integrable sector; the number of reflections λ 0 subscript 𝜆 0 \lambda_{0} in the group which determines the relative angular momentum of the ground state solution; and the order G 𝐺 G of the group which gives the number of integrable sectors required to tile the sphere. Note that there are two series of groups A m subscript 𝐴 𝑚 A_{m} and C m subscript 𝐶 𝑚 C_{m} that provide integrable mass families for any number of particles.

$N$	$A$-series	$C$-series	$H$-type	others
$3$	$A_2\equiv I_2\left(3\right)$	$C_2\equiv I_2\left(4\right)$	$H_2\equiv I_2\left(5\right)$	$I_2\left(q\right)$
	$\left[3\right]$	$\left[4\right]$	$\left[5\right]$	$\left[q\right]$
	${\lambda }_0=3$	${\lambda }_0=4$	${\lambda }_0=5$	${\lambda }_0=q$
	$G=6$	$G=8$	$G=10$	$G=2q$
$4$	$A_3$	$C_3$	$H_3$
	$[3,3]$	$[4,3]$	$[5,3]$
	${\lambda }_0=6$	${\lambda }_0=9$	${\lambda }_0=15$
	$G=24$	$G=48$	$G=120$
$5$	$A_4$	$C_4$	$H_4$	$F_4$
	$[3,3,3]$	$[4,3,3]$	$[5,3,3]$	$[3,4,3]$
	${\lambda }_0=10$	${\lambda }_0=16$	${\lambda }_0=60$	${\lambda }_0=24$
	$G=120$	$G=384$	$G=14400$	$G=1152$
$\ge 6$	$A_{N-1}$	$C_{N-1}$
	$\left[3^{N-1}\right]$	$[4,3^{N-2}]$
	${\lambda }_0=\frac{N\left(N-1\right)}{2}$	${\lambda }_0={\left(N-1\right)}^2$
	$G=\left(N-1\right)!$	$G=2^{N}N!$

Table 2. Table 2: This table provides information about the conjugacy classes of H 3 subscript 𝐻 3 H_{3} . The first column is the Schönflies notation for the elements in the class K i subscript 𝐾 𝑖 K_{i} . The second column gives the angle of rotation ϕ i subscript italic-ϕ 𝑖 \phi_{i} of the element realized in O ​ ( 3 ) O 3 \mathrm{O}(3) . The third column is whether it is generated by an even or odd number of reflections π i = ± 1 subscript 𝜋 𝑖 plus-or-minus 1 \pi_{i}=\pm 1 . All odd elements that are rotations can be considered as rotoreflections, i.e. a rotation followed by a reflection in the plane perpendicular to the rotation axis. The fourth and fifth columns are the order of the elements in the class K i subscript 𝐾 𝑖 K_{i} and the number of elements k i subscript 𝑘 𝑖 k_{i} in the class.

Elements	Angle	Even/odd	Order	Number
$E$	$0$	$+$	1	1
$C_5,C_5^4$	$2\pi /5$	$+$	5	12
$C_5^2,C_5^3$	$4\pi /5$	$+$	5	12
$C_3,C_3^2$	$2\pi /3$	$+$	3	20
$C_2$	$\pi$	$+$	2	15
$I$	$\pi$	$-$	2	1
$S_{10},S_{10}^9$	$\pi /5$	$-$	10	12
$S_{10}^3,S_{10}^7$	$3\pi /5$	$-$	10	12
$S_6,S_6^5$	$\pi /3$	$-$	6	20
$\sigma$	$0$	$-$	2	15

Equations

$$H=\sum _{i=1}^N\left(\frac{-{\hslash }^2}{2m_i}\frac{{\partial }^2}{\partial x_i^2}+\frac{1}{2}{m}_i{\omega }^2{x}_i^2\right)+g\sum _{i<j}\delta (x_i-x_j).$$

$$z_N=\sqrt{\frac{\omega }{\hslash M}}\sum _{i=1}^N{m}_i{x}_i,\text{mbox}and{\rho }^2=\sum _{j=1}^{N-1}{z}_j^2,$$

$${\omega }_{ijk}=\arctan \left(\sqrt{\frac{m_j(m_i+m_j+m_k)}{m_i{m}_k}}\right).$$

$$(R_{12}{R}_{34}{)}^2=(R_{34}{R}_{23}{)}^3=(R_{23}{R}_{12}{)}^q=1,$$

$${{\rm Y}}_{{\lambda }_0}(\hat{\mathbf{z}})=N_{{\lambda }_0}\prod _{R\in K}(\hat{\mathbf{\gamma }}(R)\cdot \hat{\mathbf{z}}),$$

$$A_3:$$

$$C_3:$$

$$H_3:$$

$$N(\tilde{E})=\frac{A}{4\pi }\tilde{E}-\frac{\ell }{4\pi }\sqrt{\tilde{E}},$$

$$\cos {\phi }_{ijk}=\frac{\cos {\omega }_{ijk}+\cos {\omega }_{jkl}\cos {\omega }_{ij,kl}}{\sin {\omega }_{jkl}\sin {\omega }_{ij,kl}}$$

$$q_m(z_1,z_2,z_3)=\sum _{\{\mathbf{\sigma }\}}(\mathbf{\sigma }\cdot (z_1,z_2,z_3){)}^m,$$

$$I_6\equiv q_6(a_1^{†},a_2^{†},a_3^{†})q_6(a_1,a_2,a_3),$$

$$H$$

$$\sum _{i<j}{\tilde{g}}_{ij}\delta \left({\hat{\mathbf{\beta }}}_{ij}\cdot \mathbf{y}\right).$$

$${\omega }_{ijk}=\arctan \left(\sqrt{\frac{m_j(m_i+m_j+m_k)}{m_i{m}_k}}\right).$$

$$H_{\mathrm{rel}}=\frac{\hslash \omega }{2}\sum _{i=1}^{N-1}\left(-\frac{{\partial }^2}{\partial z_i^2}+z_i^2\right)+\sum _{i<j}{\tilde{g}}_{ij}\delta \left({\hat{\mathbf{\gamma }}}_{ij}\cdot \mathbf{z}\right).$$

$${\rho }^2=\sum _{i=1}^{N-1}{z}_i^2=\sum _{i=1}^N{y}_i^2-Z^2,$$

$$H_{\mathrm{rel}}$$

$${\Phi }_{n,\nu ,\lambda ,\mu }(Z,\rho ,\Omega )={\zeta }_n(Z)R_{\nu \lambda }(\rho ){{\rm Y}}_{\lambda ,\mu }(\Omega ),$$

$$R_{\nu \lambda }(\rho )=A_{\nu ,\lambda }{\rho }^{\lambda }{L}_{\nu }^{\lambda +\frac{N-3}{2}}({\rho }^2)e^{-{\rho }^2/2},$$

$$L^2=\sum _{⟨i,j⟩}{L}_{ij}^2.$$

$$d(\lambda )=\frac{(m+2\lambda -2)(m+\lambda -3)!}{\lambda !(m-2)!}.$$

$$a_{\lambda }=\frac{1}{G}\sum _{K_i}{k}_i{\chi }^A(K_i){\chi }^{\lambda }(K_i).$$

$$P^A=\frac{1}{G}\sum _{g\in G_m}{\chi }^A(g)D^{\lambda }(g),$$

$${\chi }^{\lambda }(K_i)=\sum _{\mu =-\lambda }^{\lambda }\cos (m{\varphi }_i)({\pi }_i{)}^{\lambda -\mu },$$

$$\displaystyle Y_{\lambda\mu}(\theta,\phi)=\left\{\begin{array}[]{ll}\sqrt{2}N_{\lambda\mu}P^{\mu}_{\lambda}(\cos\theta)\cos(\mu\phi)&\mu>0\\ N_{\lambda 0}P^{0}_{\lambda}(\cos\theta)&\mu=0\\ \sqrt{2}N_{\lambda\mu}P^{|\mu|}_{\lambda}(\cos\theta)\sin(|\mu|\phi)&\mu<0\end{array}\right.$$

$$N_{\lambda \mu }=\sqrt{\frac{2\lambda +1}{4\pi }\frac{(\lambda -|\mu |)!}{(\lambda +|\mu |)!}}.$$

$$Y_{\lambda \mu }(\hat{\mathbf{z}})\equiv Y_{\lambda \mu }({\cos }^{-1}(z_3/\rho ),{\tan }^{-1}(z_2/z_1)).$$

$$\cos (\mu \varphi )$$

$$\sin (|\mu |\varphi )$$

$$P_{\lambda }^{\mu }(x)$$

$${\rho }^{\lambda }{P}^A{Y}_{\lambda \mu }(\hat{\mathbf{z}})=\frac{{\rho }^{\lambda }}{G}\sum _{g\in G_m}{\chi }^A(g)Y_{\lambda \mu }(O(g)\hat{\mathbf{z}}).$$

Full text

The complexity of interacting quantum systems can be partially tamed by extrapolating from solvable models, especially in one dimension Albeverio et al. (1988); Sutherland (2004); Gaudin (2014). Prominent examples include the Lieb-Liniger model with zero-range contact interactions in free space Lieb and Liniger (1963), the Tonks-Girardeau gas with hard-core contact interactions Girardeau (1960), the Calogero-Moser (CM) model with inverse square interactions either free or in a harmonic trap Calogero (1971, 2008), and the extended family of Calogero-Sutherland-Moser (CSM) models Olshanetsky and Perelomov (1983). Such models provide insights about the dynamics and thermodynamics of few-body and many-body physics, and they are proving grounds for inquiry into the nature of integrability, solvability, and chaos.

Interest in one-dimensional models has surged because of experiments with ultracold atoms trapped in tight wave guides with interactions controlled by Feshbach and confinement-induced resonances Cazalilla et al. (2011); Guan et al. (2013). These systems are well-described by a one-dimensional model with contact interactions Olshanii (1998). Dynamical effects predicted by this model like delayed thermalization due to integrability at the hard-core limit have been observed Kinoshita et al. (2006). Controllable dynamics and extended coherence times, possibly combined with internal degrees of freedom like spin or hyperfine structure, make such atomic systems suitable for exploring fundamental few-body and many-body quantum physics Serwane et al. (2011); Wenz et al. (2013) as well as for applications in quantum technologies DiVincenzo et al. (2000); Bose (2003); Volosniev et al. (2015).

However, the famous models mentioned above primarily consider equal-mass particles 111An exception is the free-space CSM model with mass-scaled interaction strengths in Ref. Sen (1996).. This article analyzes one-dimensional particles with different masses (but the same frequencies) in a harmonic trap with hard-core contact interactions. Our analysis shows that for particles with certain masses in a certain order, the mass-imbalanced hard-core system is integrable. Conversely, we provide numerical evidence that for other masses, or even the same masses but in a different order, the dynamics are quantum ergodic. Both of these limits possess potentially observable signatures in the energy level statistics Bohigas (1991); Whan (1997), particle correlations Brouzos and Foerster (2014), and thermalization dynamics Rigol et al. (2008); De Nardis et al. (2014). Measurements probing the relationship between ergodicity and entanglement in a closed quantum system have recently been performed in superconducting qubits Neill et al. (2016), showing the power for controllable quantum systems to test non-equilibrium thermodynamics.

The possibility for experimental implementation of mass-imbalanced atomic systems has driven multiple recent analyses. For general masses with contact interactions in an equal-frequency harmonic trap there are no exact solutions, so most previous approaches have relied on approximation schemes and numerical methods Mehta (2014); Dehkharghani et al. (2015); Pȩcak et al. (2016); Mehta and Morehead (2015); Pȩcak and Sowiński (2016); Dehkharghani et al. (2016); Pȩcak et al. (2017a, b). In this article, we show that for hard-core interactions there exist families of unequal masses for which there are exact solutions for the ground state and all excited states. This extends results first derived for hard-core interactions in free space Olshanii and Jackson (2015); Olshanii et al. (2016a, b). Exact solutions for hard-core interactions form the basis for approximation schemes for strongly-interacting systems, providing a valuable benchmark for testing numerical methods Volosniev et al. (2014); Deuretzbacher et al. (2014); Harshman (2016). This possibility is of special importance for mass imbalanced systems where ab initio calculations for strong interaction are especially challenging Roscher et al. (2014).

We derive the criteria for which sets of imbalanced masses are solvable and integrable using a geometrical approach. In the hard-core limit, configuration space is sectioned into $N!$ disconnected sectors, one for each order of particles. After separating out the center-of-mass and relative hyperradial degrees of freedom, the dynamics in each ordering sector of relative angular configuration space can be mapped onto hard-wall quantum billiards on a $(N\!-\!2)$-dimensional sphere. The domain is a simplex whose shape depends on the mass ratios. For three particles, the six ordering sectors are just arcs of a circles and every set of masses is therefore integrable by separation of variables. For four particles, the 24 ordering sectors are spherical triangles. Quantum billiards in a general spherical triangle cannot be solved by separability, nor can higher dimensional generalizations to $(N\!-\!2)$-simplexes on $(N\!-\!2)$-spheres. However, when a particular ordering sector tiles the sphere under reflections, then the problem is exactly solvable using something like the method of images. The possible spherical tilings are classified using Coxeter groups. Described in more detail below, Coxeter groups are point symmetry groups generated by reflections. They were originally developed for the purposes of analyzing symmetric polytopes Coxeter (1973). In hard-core contact interaction models, the $(N\!-\!1)$-dimensional hyperplanes where two particles coincide define planes of reflection. If the particle masses are correct, then these coincidence reflection planes generate a Coxeter group.

This logic can be reversed: we show that for every finite, connected, non-branching Coxeter group with rank $r$, there is one-parameter family of masses for which the dynamics of $r+1$ ordered particles is integrable. These models are non-trivial when $N>3$; we focus on the case $N=4$ where there are three families of solvable masses and the symmetries in relative configuration space are the same as the Platonic solids. We give special attention to the exceptional Coxeter group $H_{3}$ of icosahedral symmetries, and therefore this work is closely related to CSM models based on exceptional reflection groups Boreskov et al. (2005); García and Turbiner (2010, 2011). The role of Coxeter groups in providing integrability criteria for this model is perhaps not surprising because they have previously played an important role in the theory of classical and quantum dynamical systems. For example, extensions of the CSM have a closely-related classification scheme Olshanetsky and Perelomov (1983), and so do Gaudin models Gaudin (2014).

Besides applications to mass-imbalanced ultracold atomic gases, a motivating interest in this model is that it sits at the intersection of related notions of integrability and solvability. For sectors with Coxeter tiling symmetry, the energies can be calculated algebraically and all excited states can be expressed as orthogonal polynomials times the ground state; this property is called exact solvability Tempesta et al. (2001); Post et al. (2012). These solutions are constructed by Bethe-ansatz like superpositions, not of plane waves, but of spherical (or hyperspherical) harmonics. The conditions on the masses can be seen as the requirement that the scattering is non-diffractive, i.e. integrable in the Bethe-ansatz sense McGuire (1964); Sutherland (1980, 2004); Lamacraft (2013). On the other hand, we provide evidence that these models are also integrable in the classical, Liouvillian sense. Classically, Liouvillian integrability means there are $N$ functionally-independent invariant operators in involution that act continuously on $2N$ dimensional phase space. The quantum version of these operators commute with each other and therefore all states can be characterized by the spectrum of this set of operators (for a discussion of ambiguities defining integrability in quantum systems, see Caux and Mossel (2011)). As an example, we construct these operators for the four-particle case of $H_{3}$. Further, we conjecture that these solvable models are superintegrable, meaning they have more integrals of motion than degrees of freedom, and even maximally superintegrable with $2N-1$ integrals of motion (for a more complete discussion, see Evans (1990)). Superintegrable systems can have rich mathematical structure, like multi-separability and exact solvability. Further, perturbations from superintegrable models sometimes remain integrable. For example, the unitary (hard-core) limit of equal-mass particles in a harmonic trap with contact interactions is a maximally superintegrable system isomorphic to one limit of the CM model Wojciechowski (1983); Calogero (2008); Hakobyan et al. (2014). In the near-unitary limit, defined as a first order perturbation from the hard-core limit, maximal superintegrability is broken. However, the system still retains enough symmetry in the near-unitary limit that it can be mapped onto a spin chain model Volosniev et al. (2014); Deuretzbacher et al. (2014). At the other extreme from maximal superintegrability, the Coxeter group criteria could be used to identify mass-imbalanced systems where quantum ergodic dynamics should be expected. In the experimental outlook, we discuss the connections to quantum billiards and the consequences of integrability for thermalization in ultracold atomic gases.

Model and Symmetry

We consider the $N$-particle Hamiltonian with contact interactions. All particles are harmonically trapped with the same frequency. Written in terms of the particle coordinates ${\bf x}=(x_{1},x_{2},\ldots,x_{N})$, the Hamiltonian has the form

[TABLE]

In the limit $g\to\infty$ of hard-core contact interactions, the order of the particles is a dynamical invariant Harshman (2016). Configuration space is divided into $N!$ ordering sectors by $N(N-1)/2$ impenetrable coincidence hyperplanes, one for each pair of particles. The order $x_{p_{1}}\!\leq\!x_{p_{2}}\!\leq\!\cdots\!\leq\!x_{p_{N}}$ is labeled by a permutation $p=\{p_{1},p_{2},\ldots,p_{N}\}$, or more briefly $p=p_{1}p_{2}\cdots p_{N}$. Denote by $X_{ij}$ the coincidence hyperplane defined by $x_{i}-x_{j}=0$.

In Appendix A, we show that solving for the eigenstates of (1) in a particular ordering sector $p$ is equivalent to solving for the motion of a free quantum particle confined to an $(N\!-\!2)$-sphere and trapped inside an angular sector $\Omega_{p}$ bounded by $(N\!-\!1)$ hard walls. We establish this equivalence by making a mass-dependent transformation $T$ of the position coordinates ${\bf z}=T{\bf x}$, and then separating out the scaled center-of-mass $z_{N}$ and the scaled relative hyperradius $\rho$

[TABLE]

where $M$ is the total mass. The remaining relative coordinates are the $(N\!-\!2)$ hyperangles $\{\phi,\theta_{1},\ldots,\theta_{N\!-\!3}\}$ that cover the sphere $\mathcal{S}^{N\!-\!2}$. The transformation to these coordinates gives the harmonic potential a spherically symmetric form, but the coincidence hyperplanes, now transformed to $Z_{ij}=T(X_{ij})$, break that symmetry. For later convenience, denote by $\hat{\bm{\gamma}}_{ij}$ the unit normal vector to the $Z_{ij}$ coincidence hyperplane.

The specific angular ordering sector $\Omega_{p}$ is bounded by the intersection of the sphere with the $(N\!-\!1)$ coincidence hyperplanes $Z_{p_{1}p_{2}}$, $Z_{p_{2}p_{3}}$, …, $Z_{p_{N\!-\!1}p_{N}}$.Each sector $\Omega_{p}$ has $(N\!-\!2)$ angles $\omega_{ijk}$ of the form

[TABLE]

corresponding to the intersections of coincidence planes $Z_{ij}$ and $Z_{jk}$ that share a particle, and $(N\!-\!3)(N\!-\!2)/2$ angles of $\omega_{ij,kl}=\pi/2$ for the intersections of coincidence planes $Z_{ij}$ and $Z_{kl}$ that do not share a particle. For four particles each ordering sector $\Omega_{p}=\Omega_{ijkl}$ is a spherical triangle bounded by three great circles; see Fig. 1.

Solving the (hyper)spherical Helmholtz equation on an angular sector $\Omega_{p}$ with Dirichlet boundary conditions is an example of quantum billiards. The problem of quantum and classical billiards in planar triangles is well-studied Casati and Ford (1976); Richens and Berry (1981); Gutkin (1996); Glashow and Mittag (1997); Artuso et al. (1997); Casati and Prosen (1999); McGuire and Dirk (2001); Araújo Lima et al. (2013); Wang et al. (2014), and the integrability and solvability of the dynamics depends critically on the domain shape of the billiards. For example, the only three triangular billiards in a plane that have classically-integrable dynamics are the three triangles with distinguishable sides that tile the plane under reflections, without gaps or overlaps (see footnote 3 of Olshanii and Jackson (2015)). This serves as our guide for the following result for spherical quantum billiards. The dynamics in an angular sector $\Omega_{p}$ is integrable and exactly solvable when:

• •

The sector $\Omega_{p}$ tiles the $(N\!-\!2)$-sphere under reflections across its boundaries. The tiling covers the sphere with no gaps or overlaps and distinguishable sides. In other words, the $(N\!-\!2)(N\!-\!1)/2$ angles of a sector $\omega_{ijk}$ and $\omega_{ij,kl}$ define a spherical kaleidoscope.

• •

The $(N\!-\!1)$ reflections across the bounding hyperplanes $Z_{p_{1}p_{2}}$, $Z_{p_{2}p_{3}}$, …, $Z_{p_{N-1}p_{N}}$ generate a finite Coxeter reflection group. The $(N\!-\!1)$ reflection normals $\hat{\bm{\gamma}}_{p_{1}p_{2}}$, $\hat{\bm{\gamma}}_{p_{2}p_{3}}$, …, $\hat{\bm{\gamma}}_{p_{N-1}p_{N}}$ are the simple roots of the Coxeter group.

All finite reflection groups (in all dimensions) were classified by Coxeter Coxeter (1973); Humphreys (1992). Abstractly, a Coxeter group of rank $m$ is a finite group generated by $m$ reflections, where a reflection is a group element that squares to the identity. Every point symmetry group in $m$ dimensions is either a Coxeter group or a subgroup of a Coxeter group of rank $m$. For example, the three-dimensional point groups familiar from chemical and solid state physics are all subgroups of the the Coxeter groups $A_{3}$ (tetrahedral symmetry), $C_{3}$ (cubic symmetry), and $H_{3}$ (icosahedral symmetry, or they are subgroups of products of lower rank Coxeter groups.

The structure of the reflection group can be encoded by the Coxeter diagram, which can be branching or non-branching and connected or not connected. There is a family of $N$ masses that determine a ‘good’ sector for every non-branching and connected Coxeter reflection group with rank $N\!-\!1$. These groups are listed in Table 1. Only the non-branching Coxeter groups are relevant because in one-dimension each pair can have at most two adjacent pairs. Geometrically, this enforces that the coincidence planes of non-adjacent pairs like $Z_{ij}$ and $Z_{kl}$ are orthogonal $\omega_{ij,kl}=\pi/2$. We focus our attention on the connected Coxeter groups because they are relevant when all masses are finite. Disconnected graphs realize limiting cases of extreme mass imbalances. For example, for four particles, if the first or fourth particle is much more massive, i.e. like the ‘Born-Oppenheimer’ case of Mehta (2014), then the reducible Coxeter groups like $I_{2}(q)\times{A}_{1}$ could be employed.

For each group in Table 1, there exists a one-parameter family of mass sequences for which there are integrable sectors. The bracket notation for the Coxeter group $[q_{1},q_{2},\ldots,q_{N-2}]$ determines the sector angle by eq. 3, where $\omega_{i(i+1)(i+2)}=\pi/q_{i}$. In Figure 2, we show the integrable mass spectra for the three four-particle families $A_{3}$, $C_{3}$, and $H_{3}$. This can be reversed: given a set of $N$ masses in a particular order, one could check how close the sector angles derived from the masses come to the angles $\pi/q_{1},\ldots,\pi/q_{N-2}$ that define a rank $N-1$ Coxeter group.

For $N$ particles that define a ‘good’ sector, one that tiles the $(N\!-\!2)$-sphere, the generators of the Coxeter group can be chosen as the $m=N-1$ reflections $R_{ij}$ across the boundary hyperplanes $Z_{ij}$ of the sector $p={12\ldots N}$. For $N=4$, the Coxeter groups are rank $m=3$ and are generated by $R_{12}$, $R_{23}$, and $R_{34}$. All three generators square to the identity and the relations

[TABLE]

hold for $q=3$, $q=4$, and $q=5$ for $A_{3}$, $C_{3}$ and $H_{3}$, respectively. Generally, within each Coxeter group $G_{m}$, there is a conjugacy class $K\subset G_{m}$ of all reflections $R\in G_{m}$. We denote the number of reflection planes (and the order of $K$) by ${\lambda_{0}}$ and denote the normals to these planes by $\hat{\bm{\gamma}}(R)$, but remember that only $N\!-\!1$ of these planes and normals are “real”, i.e. they correspond to actual coincidence planes and normals. See Figure 3.

Note that for the Coxeter groups in the $A$-series, $C$-series and the exceptional group $F_{4}$, the sector angles are all $\pi/2$, $\pi/3$ or $\pi/4$. Inspecting eq. 3, we see that for these cases the masses are all rational fractions of each other. For example, for the group $C_{3}$ there are an countably-infinite number of rational mass sequences that give integrable sectors. The four with the lowest rational denominators are given by $(3m,m,2m,6m)$, $(10m,2m,3m,5m)$, $(12m,3m,5m,10m)$, and $(56m,7m,9m,12m)$. As discussed below, this allows the possibility of building integrable systems out of clusters of particles with the same mass.

Exact solvability and Bethe-ansatz integrability

For each Coxeter group, there is therefore a one-parameter family of masses such that the complete spectrum of energy eigenstates can be exactly solved in the ordering sector $\Omega_{1\cdots N}$ and its inverted sector $\Omega_{N\cdots 1}$. The ground state in each of these ‘Coxeter sectors’ is non-degenerate and its hyperangular wave function can be expressed as

[TABLE]

where $\hat{\bf z}=(z_{1},\ldots,z_{N-1})/\rho$ is a unit vector expressed in hyperspherical coordinates and $N_{\lambda_{0}}$ is a normalizing factor. For all $R\in K$, the function (5) is reflection antisymmetric $\Upsilon_{\lambda_{0}}(R\hat{\bf z})=-\Upsilon_{\lambda_{0}}(\hat{\bf z})$ and therefore vanishes on all reflections planes, including the coincidence planes. Note that $\rho^{\lambda_{0}}\Upsilon_{\lambda_{0}}(\hat{\bf z})$ is the lowest degree anti-invariant polynomial of the corresponding group Humphreys (1992). The function (5) is defined on the entire sphere, but its restriction to the ordering sectors $\Omega_{1\cdots N}$ or $\Omega_{N\cdots 1}$ provides that sector’s ground state with energy $\hbar\omega({\lambda_{0}}+N/2)$. Exploiting separability, a tower of states are laddered from the ground state manifold with energies $\hbar\omega(n+2\nu+{\lambda_{0}}+N/2)$, where $n$ is the center-of-mass excitation and $\nu$ is the relative hyperradial excitation [see Appendix A].

For $N$ equal masses, the Coxeter group is $A_{N-1}$ and the ground state corresponds to the lowest-energy fermionic state in a harmonic trap restricted to a sector (a la Girardeau) as expected Yukalov and Girardeau (2005); Harshman (2014). This equal-mass solution can also be seen as the limiting case of the ground state of the Calogero-Moser model with inverse-square interactions in a harmonic trap with zero coupling constant Brink et al. (1992); Vacek et al. (1994). However, unlike the equal-mass solutions, the non-equal mass solutions cannot be considered as restrictions of fermionic solutions to a single sector.

In addition to the ground state (5), the excited state relative hyperangular wave functions in a Coxeter sector is also constructed using a Bethe ansatz-like superposition of hyperspherical harmonics. Hyperspherical harmonics are homogenous polynomials in $\hat{z}_{i}$ that are eigenstates of the relative angular momentum $L^{2}_{rel}$ with eigenvalue $\lambda(\lambda+N-3)$ Avery (1989); Yáñez et al. (1994). The method takes advantage of the fact that reflections $R_{ij}$ commute with $L^{2}_{rel}$, or in other words the Coxeter group of rank $(N\!-\!1)$ is a subgroup of the orthogonal transformations $\mathrm{O}(N\!-\!1)$ [see Appendix B]. Like (5), excited solutions are first constructed over the whole sphere, and then restricted to the Coxeter sectors. By construction, the excited states are antisymmetric with respect to reflections in the Coxeter group. Not all values $\lambda>\lambda_{0}$ for the relative angular momentum support such solutions. For the three groups $A_{3}$, $C_{3}$ and $H_{3}$, the allowed spectra of $\lambda$ are

[TABLE]

In each case, the first number in the sum is $\lambda_{0}$, the relative angular momentum of the ground state and the number of reflections in the Coxeter group. Then the non-negative integers $n_{1}$ and $n_{2}$ label the excited states. Degeneracies in the hyperangular degrees of freedom arise when multiple pairs of integers provide the same $\lambda$, and the pattern of degeneracies matches the prediction of Weyl’s Law for a spherical triangle [see below]. The series of positive integers $3n_{1}+4n_{2}$, $4n_{1}+6n_{2}$, and $6n_{1}+10n_{2}$ in (6) correspond precisely to the orders of homogeneous polynomials that have definite relative angular momentum and are symmetric under the action of reflections in the groups $A_{3}$, $C_{3}$, and $H_{3}$, respectively Humphreys (1992). Incorporating the center-of-mass and hyperradial degrees of freedom, all energy eigenstates are uniquely identified by four quantum numbers: $\{n,\nu,n_{1},n_{2}\}$.

For $N\geq 4$ and general masses, or for Coxeter masses but in an arbitrary order, we believe the dynamics within sectors are not integrable in any sense. In the case of $H_{3}$, we provide evidence by numerically solving the spherical Laplacian for Coxeter masses in all sectors. We use the following method Dehkharghani et al. (2016): Each spherical triangle is flattened into an isosceles right triangle. The flattening coordinate transformation distorts the spherical Laplacian into a new operator whose spectrum must be solved inside the triangle with hard-wall boundary conditions. The spectrum is found by diagonalizing this transformed Laplacian in a basis of exact solutions for the right triangle. Details about this procedure and its convergence are described in Appendix C.

The level spacing statistics (after the standard unfolding Bohigas (1991)) for a set of four particles with $H_{3}$ mass ratios are depicted in Fig. 4. The first sector depicted is the integrable sector, whose numerical solution agrees with the prediction of (6c), the other five sectors are for the same masses arranged in other orders. For integrable sectors, the unfolded energy level statistics are expected to follow a Poissonian distribution. In our case, within the Coxeter sector the extra degeneracies of the system due to superintegrability distort the distribution so that it is peaked even more strongly at zero-energy gap Whan (1997). What we demonstrate in Fig. 4 is that even for Coxeter mass families, the incorrectly-ordered sectors show numerical evidence for quantum ergodicity in the form of Wigner-Dyson distributions for their eigenvalues.

Our numerical results for the spacing statistics open questions about the transition from integrability to ergodicity. We have performed numerical simulations on a variety of integrable mass families. While the characteristics of the Coxeter sectors remains stable, the other sectors sometimes look closer or further from Wigner-Dyson distributions, as is already visible in Fig. 4. We have investigated several possibilities for these intermediate distributions, such as integrable subclusters, but we have not arrived at any conclusive results. We have also investigated small random deviations from integrable mass sectors of the order of 5%. For this scale of deviation, the formerly integrable sectors still look far from Wigner-Dyson, but closer to Poissonian than the energy level statistics for exact Coxeter masses. Understanding the ragged edge between integrability and ergodicity using this model seems to be a productive avenue for future investigation.

To demonstrate that we find all the spectrum from this procedure, we have compared our results for $N(\tilde{E})$, the total number of energy eigenvalues below scaled energy $\tilde{E}$, to the prediction of Weyl’s Law Ivrii (2016) for a sphere:

[TABLE]

where $\tilde{E}=(2mR^{2}/\hbar^{2})E$ and $R$ is an arbitrary radius. The second term is the correction due to the Dirichlet boundary conditions proportional to the boundary length $\ell$. The area of spherical triangle $\Omega_{ijkl}$ is $A/R^{2}=\omega_{ijk}+\omega_{jkl}+\omega_{ij,kl}-\pi$ (Girard’s Theorem Brooks and Strantzen (2005)). The perimeter is $\ell/R=\varphi_{ijk}+\varphi_{jkl}+\varphi_{ij,kl}$, where the vertex angles ($\varphi_{ijk}$, $\varphi_{jkl}$, $\varphi_{ij,kl}$) satisfy Weisstein (2017)

[TABLE]

and cyclic permutations of ($\varphi_{ijk}$, $\varphi_{jkl}$, $\varphi_{ij,kl}$) and ($\omega_{ijk}$, $\omega_{jkl}$, $\omega_{ij,kl}$). Fig. 4 compares numerical solutions to this prediction.

Liouville integrability and superintegrability

The separability of the model (1) provides four functionally-independent integrals of the motion for $N\geq 3$ particles with any masses: the center-of-mass Hamiltonian $H_{\mathrm{COM}}$, the relative Hamiltonian $H_{\mathrm{rel}}$, the total angular momentum squared $L^{2}$, and the relative angular momentum squared $L_{\mathrm{rel}}^{2}$. Three of these integrals $\{H_{\mathrm{COM}},H_{\mathrm{rel}},L_{\mathrm{rel}}^{2}\}$ are in involution, but $L^{2}$ does not commute with them. This set of four integrals of motion is sufficient to prove integrability and superintegrability (but not maximal superintegrability) for $N=3$ with any masses. For $N=4$ this is not enough to prove integrability, which requires four integrals in involution, nor is it enough for superintegrability, requiring at least five total conserved quantities, and certainly not enough for maximal superintegrability, requiring seven conserved quantities.

Nonetheless, we conjecture that our model is maximally superintegrable for $N\geq 3$ when in a sector with Coxeter mass ratios. Our evidence is: (1) using the method of images Hobson (1975); Richens and Berry (1981), the classical problem can be shown to support closed orbits indicating maximal superintegrability; (2) the correspondence of our model as a limiting case of certain CSM models which are maximally superintegrable Wojciechowski (1983); Hakobyan et al. (2014). The general construction of a set of observables that provide maximal superintegrability for all reflection groups in any number of spatial dimensions requires the methods of algebraic geometry. This ongoing project will be the subject of a future publication.

Let us outline a potential strategy for searching for the missing integrals of motion for four particles using $H_{3}$ as an example. The key role is played by the invariant polynomials of the group $q_{m}(z_{1},z_{2},z_{3})$ with order $m$. These are the lowest-order, homogeneous, functionally-independent polynomials that remain unchanged under any of the group transformations. They are known and tabulated for all the reflection groups Mehta (1988). For $\mbox{H}_{3}$, there are three invariant polynomials $q_{m}$ with order $m=2$, $6$ and $10$ and (up to a normalization) they are constructed as

[TABLE]

where $\{\bm{\sigma}\}$ are the set of vectors describing the six five-fold rotation axes of $H_{3}$. From the three polynomials $q_{2}$, $q_{6}$, and $q_{10}$, we define the three operators $J_{m}\equiv q_{m}(L_{12},\,L_{23},\,L_{31})$. Here, $L_{ij}$ are components of the vector of the relative angular momentum in the $ij$-plane. Note that $J_{2}$ is proportional to $L_{\mathrm{rel}}^{2}$ and so it does not give an additional integral of motion.

However, the operator $J_{6}$ completes the commuting set $\{H_{\mathrm{COM}},H_{\mathrm{rel}},L_{\mathrm{rel}}^{2}\}$ to a Liouvillian set. Since $J_{6}$ commutes with mirror reflections of the $H_{3}$ group, by Schur’s Lemma it must act as a multiple of the identity on the antisymmetric states. It commutes with the previous three members of the Liouvillian set, and all four can be readily shown to be functionally-independent in the classical sense. The five-member set $\{H_{\mathrm{COM}},L^{2},H_{\mathrm{rel}},L_{\mathrm{rel}}^{2},J_{6}\}$ now establishes superintegrability for the $H_{3}$ mass family.

This set can be further extended to a maximally superintegrable set using the operator $J_{10}$ and another invariant operator $I_{6}$ defined by

[TABLE]

where $a_{j}=(-i\partial_{z_{j}}-iz_{j})/\sqrt{2}$ is an annihilation operator for the $j$-th component of the relative motion. The operator $I_{6}$ naturally commutes with the total Hamiltonian $H$. The resulting seven-member set is (classically) functionally independent and establishes maximal superintegrability for the $H_{3}$ model. The scheme can readily be generalized to the other two three-dimensional reflection groups, $A_{3}$ and $C_{3}$. However, no ready generalization to higher dimensions exists for the Liouvillian sets, because, a priori, the operators $J_{m}$ do not commute between themselves. Finding Liouvillian sets for higher-dimensional groups is a subject of future work. Identifying and classifying the maximal superintegrablty sets, and ideally connecting them to the known integrals for the Calogero-Moser model Wojciechowski (1983); Hakobyan et al. (2014); Saghatelian (2012), is another ongoing project.

Experimental outlook

At the moment, three possible experimental applications of the models considered in this article can be foreseen. The first possibility is the straightforward idea of finding a collection of atoms that naturally have the right mass ratios and seeking the signatures of integrability in the spectral, coherence, and thermalization properties of the system. Even if the particles are only close to Coxeter family, our numerical results for the energy spectrum suggest that traces of integrability should still be present. More generally, the Coxeter criteria can be used to measure how far from integrability particular arrangements of imbalanced masses are expected to be, or whether there are integrable subclusters possible within a multi-species ultracold atomic gas.

In the second scheme, if real masses with the correct ratio are not available, the atomic mass is controlled using optical lattices. Given sufficient laser power, the effective mass Kittel (2004) can be tuned from its “bare” value to almost zero Gadway et al. (2011). In particular the effective mass can be made $3$ times greater than its bare counterpart in a lattice of a depth $V_{0}=7E_{\text{R}}$, and $23$ times greater for $V_{0}=16E_{\text{R}}$ respectively. Here $E_{\text{R}}=\hbar^{2}k^{2}/(2m)$ is the so-called recoil energy, $k$ is the wave vector of light that creates the lattice, $m$ is the bare atomic mass. In both cases, harmonic confinement represents the most natural experimental environment, unlike the box and ring geometries traditionally studied using Bethe ansatz methods.

In the third scheme, described in more detail in Olshanii et al. (2016b), the role of massive particles is played by bosonic solitons in an atom waveguides Khaykovich et al. (2002); Strecker et al. (2002); Cornish et al. (2006). The solitons are made of atoms in two alternating internal states where the intraspecies interaction is attractive and the interspecies interaction is repulsive. The goal would be to engineer clusters of atoms whose combined masses satisfy the Coxeter criteria. For example, the clustering pattern $(3m,m,2m,6m)$ has $C_{3}$ symmetry. A mixture of ${}^{7}\mbox{Li}$ atoms with $m_{F}=-1$ and $m_{F}=0$ in a magnetic field of $855\,\mbox{G}$ constitutes an example 222R.G. Hulet, private communication..

Any implementation of the models considered in our article may constitute an efficient experimental realization of spherical triangular (or higher-dimensional, simplex-shaped) quantum billiards 333We note that in principle the relative Hamiltonian for any of the $N\!=\!4$ families also could be realized with a non-interacting Bose-Einstein condensate in a spherically-symmetric harmonic trap sliced by six sheets of intense laser light, or even just three sheets to create a sector.. The ergodicity of classical flat triangular billiards is conjectured to strongly depend on the rationality of the billiard angles Casati and Prosen (1999); Wang et al. (2014). Numerically, such questions about ergodicity are difficult, requiring long propagation for averages to converge to their infinite time limits. A study of the eigenstate-to-eigenstate variance of the expectation values of observables Deutsch (1991); Srednicki (1994); Rigol et al. (2008), which is a faithful quantum analogue of classical deviations from ergodicity (c.f. Zhang et al. (2015) for a comparison), may provide an efficient alternative to classical long-time averages. Experimentally, one may conjecture an appearance of a memory of initial conditions, if the billiard is not ergodic Yurovsky and Olshanii (2011). The mass mixtures considered in our paper could constitute a way to study multi-dimensional classical and quantum hard-wall billiards with continuously tunable geometry, a powerful extension of the existing experimental techniques Friedman et al. (2001).

Acknowledgments: This work was supported by the Aarhus University Research Foundation, the US National Science Foundation Grant No. PHY-1607221, the Office of Naval Research Grant N00014-12-1-0400, The Danish Council for Independent Research (Sapere Aude), and the Humboldt Foundation. Thanks to FF Bellotti, MES Andersen, L Rammelmüller, and F Werner for useful conversations.

Appendix A: Coordinate Transformations and the Map to Quantum Billiards

In the following section we establish the map from the model Hamiltonian (1) to a free particle in a bounded region on the $(N\!-\!2)$-sphere. Much of this is well-known Werner and Castin (2006); Dehkharghani et al. (2016), but we reproduce it here for the readers’ convenience and to establish notation.

The equipotentials of (1) are $N$-dimensional ellipsoids segmented into $N!$ sectors by $(N-1)$-dimensional hyperplanes $X_{ij}$ defined by the particle coincidences $x_{i}-x_{j}=0$; see Fig. 5. In the limit $g\to\infty$ these planes are impenetrable. The angle between coincidence hyperplanes $X_{ij}$ and $X_{jk}$ that share a particle is $\pi/3$ (only possible for $N\geq 3$); the angle between hyperplanes $X_{ij}$ and $X_{kl}$ that do no share a particle is $\pi/2$ (only possible for $N\geq 4$).

As a first step, we scale the position coordinates into unitless position variables $y_{i}=\sqrt{m_{i}\omega/\hbar}\,x_{i}$. Then the Hamiltonian (1) becomes

[TABLE]

where $\mu_{ij}=m_{i}m_{j}/(m_{i}+m_{j})$ and $\tilde{g}_{ij}=g_{ij}\sqrt{\mu_{ij}\omega/\hbar}$. This scaling transformation ${\bf y}=S{\bf x}$ has brought the harmonic potential into a form with $N$-spherical symmetry but at the cost of de-symmetrizing the coincidence planes. To describe the geometry, define the transformed coincidence planes $Y_{ij}\equiv SX_{ij}$ with normals $\hat{\bm{\beta}}_{ij}$. The contact interaction then has the form

[TABLE]

The angle $\omega_{ijk}$ between coincidence planes $Y_{ij}$ and $Y_{jk}$ with normals $\hat{\bm{\beta}}_{ij}$ and $\hat{\bm{\beta}}_{jk}$ is now

[TABLE]

Whereas in the equal mass case, $\omega_{ijk}$ is always $\pi/3$, for three arbitrary masses it can range from [math] ($m_{j}$ much lighter than other two masses) to $\pi/2$ ($m_{j}$ much heavier). The angle $\omega_{ij,kl}$ between coincidence planes $Y_{ij}$ and $Y_{kl}$ that do not share a particle remains $\pi/2$.

In the limit $g\to\infty$ the Hamiltonian (Appendix A: Coordinate Transformations and the Map to Quantum Billiards) separates in hyperspherical coordinates with radius $R^{2}=\sum_{i}y_{i}^{2}$. The interaction term [12] is proportional to $1/R$, so it is not separable for finite values of $\tilde{g}_{ij}$, but as $\tilde{g}_{ij}\to\infty$ there is no distinction between $1/R$ or $1/R^{2}$ times the sum of delta functions. So hyperspherical symmetry of (Appendix A: Coordinate Transformations and the Map to Quantum Billiards) emerges and there is $\mathrm{SO}(2,1)$ dynamical symmetry in the total hyperradial coordinate $R$ Werner and Castin (2006). To develop physical intuition, it is sometimes useful to imagine a single, classical particle bouncing around in this $N$-dimensional landscape. In this mass-rationalized geometry, the classical particle trajectory changes its direction of angular momentum when it bounces off of a coincidence hyperplane $Y_{ij}$, but it does not change its magnitude of “angular momentum” in configuration space. However, total angular momentum does not respect the center-of-mass separability and does not commute with the relative Hamiltonian or relative angular momentum [see below], so we do not exploit it here.

Next we rotate the coordinate system ${\bf z}=J{\bf y}$ so that the component $z_{N}\equiv Z$ is the scaled center-of-mass $Z=\sum_{i}y_{i}\sqrt{m_{i}/M}$ and $M$ is the total mass. The orthogonal transformation $J$ with this property is not unique and its selection determines a particular choice for Jacobi relative coordinates $z_{1}$ through $z_{N-1}$. The transformation $J$ also rotates the coincidence planes $Z_{ij}\equiv JY_{ij}$ and their normals $\hat{\bm{\gamma}}_{ij}\equiv J\hat{\bm{\beta}}_{ij}$, but leaves the angles between planes like $\omega_{ijk}$ and $\omega_{ij,kl}$ invariant. Since all the normal vectors $\hat{\bm{\gamma}}_{ij}$ have zero $Z$-components, the Hamiltonian in ${\bf z}$-coordinates separates into $H=H_{\mathrm{COM}}+H_{\mathrm{rel}}$, where $H_{\mathrm{COM}}$ is the Hamiltonian for a one-dimensional harmonic oscillator in the center-of-mass $Z$ coordinate and the relative Hamiltonian is

[TABLE]

Finally, we go to hyperspherical coordinates in the relative space, where the relative hyperradius $\rho$ is

[TABLE]

and there are $(N\!-\!2)$ angles charting the sphere $\mathcal{S}^{N\!-\!2}$, conventionally chosen as $\Omega=\{\phi,\theta_{1},\ldots,\theta_{N\!-\!3}\}$ with $\phi\in[0,2\pi)$ and $\theta_{i}\in[0,\pi]$. The relative Hamiltonian now becomes

[TABLE]

where $\Delta_{\Omega}$ is the angular part of the Laplacian in relative configuration space. As before, the relative Hamiltonian (16) is $\rho$-$\Omega$ separable in the limit $\tilde{g}_{ij}\to\infty$. A general energy eigenstate can be separated into a product of center-of-mass $\zeta(Z)$, relative hyperradial $R(\rho)$ and relative hyperangular $\Upsilon(\Omega)$ functions

[TABLE]

where $n$ is the center-of-mass quantum number, the function $\zeta_{n}(Z)$ is the one-dimensional harmonic oscillator wave function, and $\nu$ is the relative hyperradial quantum number. At this point, $\lambda$ is just derived from the judiciously-parametrized relative hyperangular separation constant $\lambda(\lambda+N-3)$ and $\mu$ is just an additional label to distinguish any possible degenerate states for a given $\lambda$. If there was no angular potential $\tilde{g}_{ij}=0$, then $\lambda$ would be a non-negative integer, the eigenfunctions on $\mathcal{S}^{N\!-\!2}$ would be the hyperspherical harmonics, and $\mu$ would be a collective index to label degeneracies Avery (1989); Yáñez et al. (1994). However since there is an angular potential in (16), the hyperangular solutions are unknown and we must explicitly solve for $\lambda$, including any possible degeneracies. Whatever value $\lambda$ takes (including non-integer values), the relative hyperradial function $R_{\nu\lambda}(\rho)$ is the standard solution for the radial factor of an $(N\!-\!1)$-dimensional isotropic harmonic oscillator in hyperspherical coordinates Yáñez et al. (1994) with energy $\hbar\omega(2\nu+\lambda+(N\!-\!1)/2)$:

[TABLE]

where $A_{\nu,\lambda}=\sqrt{2\,\nu!/\Gamma(\nu\!+\!\lambda\!+\!(N\!-\!1)/2)}$.

We have achieved our desired result: this series of coordinate transformations has reduced solving the $N$-particle Hamiltonian (1) with equal frequencies and infinite-strength contact interactions into solving hard-wall quantum billiards in $(N\!-\!2)$-simplexes on $(N\!-\!2)$-spheres.

Appendix B: Construction of Exact Solutions

In this section, we construct the wave functions within Coxeter sectors, i.e. sectors of the $\mathcal{S}^{N\!-\!2}$ hypersphere defined by the relative hyperangular coordinates that have the right shape to tile the sphere under reflection. For convenience, choose the Coxeter sector to be the ordering sector $\Omega_{12\cdots N}$ so that it is bounded by the coincidence hyperplanes $Z_{12}$ through $Z_{(N\!-\!1)N}$. These $(N\!-\!1)$ hyperplanes define the reflections $R_{12}$ though $R_{(N\!-\!1)N}$ that generate the Coxeter group.

The Coxeter group $G_{m}$ is generated by $m$ reflections in $m$ dimensions. As such, it can considered as a subgroup of $\mathrm{O}(m)$, orthogonal transformations in $m$ dimensions and the symmetry of the sphere $\mathcal{S}^{m\!-\!1}$. To summarize the method, there is a solution to the Hamiltonian in the Coxeter sectors $\Omega_{12\cdots N}$ and $\Omega_{N\cdots 21}$ whenever an irreducible representation (irrep) of $G_{N\!-\!1}$ that is antisymmetric under all reflections appears in the decomposition of an irrep of $\mathrm{O}(N\!-\!1)$. The irreps of $\mathrm{O}(m)$ generally are reducible with respect to the subgroup $G_{m}$. The method of characters can answer the question as to whether an irrep of a subgroup appears in the decomposition of the irrep of the group. When it does exist, the corresponding states can be constructed using projection operators and (in the case of degeneracies) an orthonormalization procedure.

The irreducible representations for $\mathrm{O}(m)$ and their realizations by hyperspherical harmonics are well-known Hamermesh and Physics (1989); Wybourne (1974), and so we just summarize a few facts here for the readers’ convenience. The subgroup $\mathrm{SO}(m)$ is a Lie group with $m(m-1)/2$ generators in the Lie algebra. Denote these generators as $L_{ij}$ for $i<j$ with $i,j\in\{1,\ldots,m\}$, where $L_{ij}$ generates a rotation in the $ij$-plane. The quadratic Casimir of $\mathrm{SO}(m)$ is the sum of all of these generators squared

[TABLE]

For $\mathrm{SO}(3)$ this is the familiar angular momentum squared operator with eigenvalues $\lambda(\lambda+1)$. The $\mathrm{SO}(3)$ irreps are labeled by $\lambda$ and have degeneracy $2\lambda+1$. In $m>3$ dimensions, the operator $L^{2}$ is the hyperangular momentum squared operator with eigenvalue $\lambda(\lambda+m-2)$ and an irrep labeled by $\lambda$ has degeneracy Avery (1989)

[TABLE]

Once the representation of total inversion is chosen, the irreps of $\mathrm{SO}(m)$ also naturally carry a representation of $\mathrm{O}(m)$ . For example, inversion is represented in the $\lambda$ irrep by multiplication by $(-1)^{\lambda}$ for $\mathrm{O}(3)$.

To reduce an $\mathrm{O}(m)$ irrep $\lambda$ into the irreps of Coxeter group $G_{m}$, the $G$ elements are sorted into conjugacy classes $K_{i}$ with $k_{i}$ elements. Each irrep $W$ of $G_{m}$ has a unique pattern of characters $\chi^{W}(K_{i})$. In particular we are interested in the $G_{m}$ irrep $W=A$ of all anti-invariant states, meaning $\chi^{A}(K_{i})=1$ when $K_{i}$ is a conjugacy class whose elements are an even composition of reflections and $\chi^{A}(K_{i})=-1$ when $K_{i}$ is a class comprised of odd compositions. Further, each conjugacy class has a character $\chi^{\lambda}(K_{i})$ in the $G_{m}$-reducible $\mathrm{O}(m)$-irrep denoted by $\lambda$. When these characters are known, then the number of times the $G_{m}$ irrep $A$ appears in the decomposition of the $\mathrm{O}(m)$ irrep $\lambda$ is Hamermesh and Physics (1989)

[TABLE]

Note that $\chi^{A}(K_{i})=(\chi^{A}(K_{i}))^{*}$ because Coxeter groups are ambivalent. The number $a_{\lambda}$ is an integer that counts how many solutions there are with relative angular momentum $\lambda$. The projection operator onto the anti-invariant irrep $A$ is given by

[TABLE]

where $D^{\lambda}(g)$ is the representation of group element $g$ acting irreducibly on the $d(\lambda)$-dimensional representation space. If $a^{\lambda}=1$, any vector in the $\lambda$ irrep space with a non-zero projection will be proportional to the solution we seek. If $a_{\lambda}>1$, then a set of orthonormal solutions can be found by projecting multiple vectors and then applying Gram-Schmidt orthogonalization.

As an example, consider the Coxeter group $H_{3}$. This group has ten conjugacy classes summarized in Table 2. The $H_{3}$ character $\chi^{A}(K_{i})$ for the anti-invariant irrep is $+1$ for the five even conjugacy classes and $-1$ for the five odd classes. The $\mathrm{O}(m)$ character $\chi^{\lambda}(K_{i})$ for irrep $\lambda$ is

[TABLE]

where $\phi_{i}$ is the angle of rotation and $\pi_{i}$ is the reflection parity for the conjugacy class $K_{i}$. Plugging this into (19), we find the pattern of degeneracies given in eq. [6c] of the main text. The same method is used in Harmer (2007) to find which $\mathrm{O}(m)$ irreps have symmetric irreps of the spherical triangle groups in their reduction.

The construct the actual states $\Upsilon_{\lambda}(\theta,\phi)$, we use the projection operator (20) acting on the spherical harmonics. Instead of explicitly constructing the $(2\lambda+1)\times(2\lambda+1)$ unitary matrices $D^{\lambda}(g)$ that act on the spherical harmonics $Y^{\mu}_{\lambda}(\theta,\phi)$ for each of the 120 elements of $H_{3}$, we choose a slightly different method that takes advantage of two facts: (1) the spherical harmonics can be written as polynomials of the relative coordinates; and (2) we already have the $3\times 3$ matrices $O(g)\in\mathrm{O}(3)$ that represent $H_{3}$ as rotation and reflections.

The first step is to express the spherical harmonics in terms of the relative coordinates $(z_{1},z_{2},z_{3})$. We work with the real form of the spherical harmonics, defined as

[TABLE]

where

[TABLE]

The real spherical harmonics can be written in terms of the relative coordinates:

[TABLE]

Noting the relations

[TABLE]

we can show that $\rho^{\lambda}Y_{\lambda\mu}(\hat{\bf z})$ are homogeneous polynomials of order $\lambda$ in $(z_{1},z_{2},z_{3})$.

The projection (20) is applied using the transformation matrix $O(g)$

[TABLE]

This projection will be zero unless $\lambda$ is in the spectrum given by eq. [6c] in the main text. Note that it may also be zero for any particular $\mu$, but there must be as many linearly-independent polynomials (that also solve the spherical Laplacian) as there are solutions for $n_{1}$ and $n_{2}$ for a given $\lambda$ in eq. [6c]. For $H_{3}$, the first time there are multiple solutions is when $\lambda=45$. Explicit calculation for $\lambda=\lambda_{0}=15$ confirms eq. [5] from the main text for the ground state of $H_{3}$, and we have performed the same calculations for the other $N=4$ Coxeter groups.

Appendix C: Numerical Method for Solving Spherical Triangle

Here we present the numerical method to calculate the energy spectrum for $N=4$ particles with arbitrary masses, which is an extension of the method found in Dehkharghani et al. (2016) to spherical triangles. The general strategy was outlined in the main text: here we provide details about the coordinates, the flattening, and the exact diagonalization used to construct Fig. 4.

After separation of variables, we must solve for the hyperangular wave function $\Upsilon(\Omega)$ within a sector $\Omega_{p}$. This function must satisfy

[TABLE]

with Dirichlet boundary conditions on the three bounding coincidence planes $Z_{p_{1}p_{2}}$, $Z_{p_{2}p_{3}}$, and $Z_{p_{3}p_{4}}$. To solve the problem for arbitrary masses, a particular rotation $J$ in ${\bf z}=JS{\bf x}$ must be specified. For numerical simplicity, we choose the $H$-type four-body coordinates Rittenhouse et al. (2011). Then the rotation $J$ aligns the coordinate plane $Z_{12}$ with the plane $z_{1}=0$ and aligns the coordinate plane $Z_{34}$ with $z_{2}=0$ (see Fig. 1 of main text). The other four coincidence planes are given by the following equations:

[TABLE]

The numerical problem is solved independently in each sector $\Omega_{p_{1}p_{2}p_{3}p_{4}}$ and checked for consistency with the similar sector $\Omega_{p_{4}p_{3}p_{2}p_{1}}$. In the following, as an example we solve the sector $\Omega_{1342}$, which is limited by $Z_{34}$, $Z_{13}$, and $Z_{24}$ coincidence planes. We isolate $z_{1}$ in the above equations and introduce (non-standard) spherical coordinates with $z_{1}=\rho\cos\theta$, $z_{2}=\rho\sin\theta\cos\phi$ and $z_{3}=\rho\sin\theta\sin\phi$ so that $\cos\theta>0$ in sector $\Omega_{1342}$. Then substituting spherical coordinates into the coincidence plane equations and dividing by $z_{1}=\rho\cos\theta>0$, we find

[TABLE]

with

[TABLE]

Introducing new coordinates $u=\tan\theta\cos\phi$ and $v=\tan\theta\sin\phi$, the boundary conditions

[TABLE]

are simple and describe a triangle in flat space. However, the differential operator $\Delta_{\Omega}$ has become more complicated:

[TABLE]

Next we introduce a final coordinate transformation

[TABLE]

This can be inverted as

[TABLE]

Notice that with this choice the coincidence planes are mapped into the nicely symmetric form

[TABLE]

These are the boundaries of a right, isosceles triangle with corners at $(s,t)=(-1,-1)$, $(-1,1)$ and $(1,-1)$. The transformation of $\Delta_{\Omega}$ induced by the coordinate change $(u,v)$ to $(s,t)$ is a lengthy but straightforward and we do not show it here.

A complete, normalized basis for Lebesgue square-integrable functions on this domain bounded by (29) is provided by the functions $h_{n,m}(s,t)$:

[TABLE]

where $n$ and $m$ are positive integers and $n<m$. The matrix elements of the transformed spherical Laplacian can be calculated in this basis and then diagonalized to find the spectrum. In our calculations we set $n,m$ to go up to $N_{max}=80$. With this upper bound the first 300 energies were quite converged up to the second decimal place. We know this because we did convergence analysis from 60-80 and found out that the eigenenergies up to the 2nd decimal place were not changing. The calculation time for $N_{max}=80$ was approximately 10 days on a reasonably powerful desktop computer. If one is interested in higher excited state energies, then one needs to increase $N_{max}$ in order to get a better precision at the higher end of the spectrum.

Additionally, numerical results were compared to the exact algebraic results for the integrable Coxeter sector for several mass families in $A_{3}$, $C_{3}$, and $H_{3}$ to confirm the uncertainty estimates. And as described in the main text, we compared the level density of the spectrum to the prediction of Weyl’s Law in order to establish that all eigenstates were found by this method.