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

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.

The complexity of interacting quantum systems can be partially tamed by extrapolating from solvable models, especially in one dimension [1][2][3]. Prominent examples include the Lieb-Liniger model with zero-range contact interactions in free space [4], the Tonks-Girardeau gas with hard-core contact interactions [5], the Calogero-Moser (CM) model with inverse square interactions either free or in a harmonic trap [6,7], and the extended family of Calogero-Sutherland-Moser (CSM) models [8]. 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 [9,10]. These systems are well-described by a one-dimensional model with contact interactions [11]. Dynamical effects predicted by this model like delayed thermalization due to integrability at the hard-core limit have been observed [12]. 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 [13,14] as well as for applications in quantum technologies [15][16][17].
However, the famous models mentioned above primarily consider equal-mass particles [18]. 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 massimbalanced 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 pos-sess potentially observable signatures in the energy level statistics [19,20], particle correlations [21], and thermalization dynamics [22,23]. Measurements probing the relationship between ergodicity and entanglement in a closed quantum system have recently been performed in superconducting qubits [24], 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 [25][26][27][28][29][30][31][32]. 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 hardcore interactions in free space [33][34][35]. Exact solutions for hard-core interactions form the basis for approximation schemes for strongly-interacting systems, providing a valuable benchmark for testing numerical methods [36][37][38]. This possibility is of special importance for mass imbalanced systems where ab initio calculations for strong interaction are especially challenging [39].
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 arXiv:1704.01433v3 [quant-ph] 28 Aug 2017 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 [40]. 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 [41][42][43]. 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 [8], and so do Gaudin models [3].
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 [44,45]. 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 [2,[46][47][48]. 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 [49]). As an ex-ample, 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 [50]). 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 [7,51,52]. 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 [36,37]. 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 x = (x 1 , x 2 , . . . , x N ), the Hamiltonian has the form In the limit g → ∞ of hard-core contact interactions, the order of the particles is a dynamical invariant [38]. 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 p1 ≤ x p2 ≤ · · · ≤ x p N is labeled by a permutation p = {p 1 , p 2 , . . . , p N }, or more briefly p = p 1 p 2 · · · 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 Ω p bounded by (N−1) hard walls. We establish this equivalence by making a mass-dependent transformation T of the position coordinates z = T x, and then separating out the scaled center-of-mass z N and the scaled  Figure 1. This figure depicts the relative configuration space for four particles in the simplest case when all masses are the same, an example of the one-parameter mass family of the Coxeter group A3. The gray sphere represents an equipotential of the harmonic trap in the mass-normalized coordinates (z1, z2, z3). The six colored disks that intersect the plane represent the coincidence planes Z12 (red), Z13 (yellow), Z14 (green), Z23 (cyan), Z24 (blue), Z34 (magenta). Twelve sectors are visible and are labeled by Ωp, where p is the order of the four particles. For two sectors Ω1234 and Ω1324, the three angles are also labeled. relative hyperradius ρ where M is the total mass. The remaining relative coordinates are the (N − 2) hyperangles {φ, θ 1 , . . . , θ N−3 } that cover the sphere 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γ ij the unit normal vector to the Z ij coincidence hyperplane. The specific angular ordering sector Ω p is bounded by the intersection of the sphere with the (N−1) coincidence hyperplanes Z p1p2 , Z p2p3 , . . . , Z p N−1 p N .Each sector Ω p has (N −2) angles ω ijk of the form corresponding to the intersections of coincidence planes Z ij and Z jk that share a particle, and (N −3)(N −2)/2 angles of ω ij,kl = π/2 for the intersections of coincidence planes Z ij and Z kl that do not share a particle. For four particles each ordering sector Ω p = Ω ijkl is a spherical triangle bounded by three great circles; see Fig. 1. Solving the (hyper)spherical Helmholtz equation on an angular sector Ω p with Dirichlet boundary conditions is an example of quantum billiards. The problem of quantum and classical billiards in planar triangles is wellstudied [53][54][55][56][57][58][59][60][61], 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 [33]). This serves as our guide for the following result for spherical quantum billiards. The dynamics in an angular sector Ω p is integrable and exactly solvable when: • The sector Ω 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 ω ijk and ω ij,kl define a spherical kaleidoscope.
All finite reflection groups (in all dimensions) were classified by Coxeter [40,62]. 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 nonbranching 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 I. 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 ω ij,kl = π/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 [25], then the reducible Coxeter groups like I 2 (q) × A 1 could be employed.
For each group in Table I, there exists a one-parameter family of mass sequences for which there are integrable Table I. This table provides information about the connected, non-branching, finite Coxeter reflection groups. For N particles, each rank m = N − 1 Coxeter group defines a oneparameter family of masses for which the system is exactly solvable. For each group Gm, the following data are provided [63]: the Coxeter bracket [q1, . . . , qm] from which one determines the angles of the integrable sector; the number of reflections λ0 in the group which determines the relative angular momentum of the ground state solution; and the order G of the group which gives the number of integrable sectors required to tile the sphere. Note that there are two series of groups Am and Cm that provide integrable mass families for any number of particles.
sectors. The bracket notation for the Coxeter group [q 1 , q 2 , . . . , q N −2 ] determines the sector angle by eq. 3, where ω i(i+1)(i+2) = π/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 π/q 1 , . . . , π/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 . . . 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 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 ⊂ G m of all reflections R ∈ G m . We denote the number of reflection planes (and the order of K) by λ 0 and denote the normals to these planes byγ(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  Note that for the Coxeter groups in the A-series, Cseries and the exceptional group F 4 , the sector angles are all π/2, π/3 or π/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 countablyinfinite 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 oneparameter family of masses such that the complete spectrum of energy eigenstates can be exactly solved in the ordering sector Ω 1···N and its inverted sector Ω N ···1 . The ground state in each of these 'Coxeter sectors' is nondegenerate and its hyperangular wave function can be expressed as whereẑ = (z 1 , . . . , z N −1 )/ρ is a unit vector expressed in hyperspherical coordinates and N λ0 is a normalizing factor. For all R ∈ K, the function (5) is reflection antisymmetric Υ λ0 (Rẑ) = −Υ λ0 (ẑ) and therefore vanishes on all reflections planes, including the coincidence planes. Note that ρ λ0 Υ λ0 (ẑ) is the lowest degree anti-invariant polynomial of the corresponding group [62]. The function (5) is defined on the entire sphere, but its restriction to the ordering sectors Ω 1···N or Ω N ···1 provides that sector's ground state with energy ω(λ 0 + N/2). Exploiting separability, a tower of states are laddered from the ground state manifold with energies ω(n+2ν +λ 0 +N/2), where n is the center-of-mass excitation and ν 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 [64,65]. 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 [66,67]. 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ẑ i that are eigenstates of the relative angular momentum L 2 rel with eigenvalue λ(λ + N − 3) [68,69]. 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 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 λ > λ 0 for the relative angular momentum support such solutions. For the three groups A 3 , C 3 and H 3 , the allowed spectra of λ are C 3 : λ = 9 + 4n 1 + 6n 2 , and (6b) In each case, the first number in the sum is λ 0 , the relative angular momentum of the ground state and the number of reflections in the Coxeter group. Then the nonnegative 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 λ, 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 [62]. Incorporating the center-of-mass and hyperradial degrees of freedom, all energy eigenstates are uniquely identified by four quantum numbers: {n, ν, n 1 , n 2 }. For N ≥ 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 [30]: 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 [19]) 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 [20].
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.

LIOUVILLE INTEGRABILITY AND SUPERINTEGRABILITY
The separability of the model (1) provides four functionally-independent integrals of the motion for N ≥ 3 particles with any masses: the center-of-mass Hamiltonian H COM , the relative Hamiltonian H rel , the total angular momentum squared L 2 , and the relative angular momentum squared L 2 rel . Three of these integrals {H COM , H rel , L 2 rel } are in involution, but L 2 does not  [19]. The integrable sector Ω1234 is depicted in the top left graph and agrees with the prediction from (6c). The blue lines depict the Poissonian statistics expected for an integrable system; the red lines are Wigner-Dyson distribution derived from random matrix theory expected for quantum ergodic systems with time-reversal symmetric Hamiltonians. The bottom graph shows the quality of the Weyl's Law (7) in the integrable sector Ω1234 as well as the non-integrable sectors.
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 ≥ 3 when in a sector with Coxeter mass ratios. Our evidence is: (1) using the method of images [54,73], 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 [51,52]. 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, functionallyindependent polynomials that remain unchanged under any of the group transformations. They are known and tabulated for all the reflection groups [74]. For 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 where {σ} are the set of vectors describing the six fivefold rotation axes of H 3 . From the three polynomials q 2 , q 6 , and q 10 , we define the three operators J m ≡ 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 2 rel and so it does not give an additional integral of motion.
However, the operator J 6 completes the commuting set {H COM , H rel , L 2 rel } 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 COM , L 2 , H rel , L 2 rel , 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 where a j = (−i∂ zj − iz j )/ √ 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 higherdimensional 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 [51,52,75], 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 [76] can be tuned from its "bare" value to almost zero [77]. In particular the effective mass can be made 3 times greater than its bare counterpart in a lattice of a depth V 0 = 7E R , and 23 times greater for V 0 = 16E R respectively. Here E R = 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 [35], the role of massive particles is played by bosonic solitons in an atom waveguides [78][79][80]. 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 Li atoms with m F = −1 and m F = 0 in a magnetic field of 855 G constitutes an example [81].
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 [82]. The ergodicity of classical flat triangular billiards is conjectured to strongly depend on the rationality of the billiard angles [58,61]. 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 [22,83,84], which is a faithful quantum analogue of classical deviations from ergodicity (c.f. [85] 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 [86]. 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 [87].
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 wellknown [30,89], 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 → ∞ these planes are impenetrable. The angle between coincidence hyperplanes X ij and X jk that share a particle is π/3 (only possible for N ≥ 3); the angle between hyperplanes X ij and X kl that do no share a particle is π/2 (only possible for N ≥ 4).
As a first step, we scale the position coordinates into unitless position variables y i = m i ω/ x i . Then the Hamiltonian (1) becomes where µ ij = m i m j /(m i + m j ) andg ij = g ij µ ij ω/ . This scaling transformation y = Sx 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 ≡ SX ij with normalsβ ij . The contact interaction then has the form i<jg ij δ β ij · y .
The angle ω ijk between coincidence planes Y ij and Y jk with normalsβ ij andβ jk is now Whereas in the equal mass case, ω ijk is always π/3, for three arbitrary masses it can range from 0 (m j much a b c Figure 5. This figure depicts the coordinate transformation for three particles with the A2 masses m1 = 1/13M , m2 = 9/52M , and m = 3/4M , where M is the total mass. Subfigure (a) depicts the configuration space in natural particle positions x; the gray ellipsoid represents an equipotential for the equal-frequency (and therefore not equal-strength); the red, green, and blue disks represent the X12, X23, and X31 coincidence planes, respectively. Subfigure (b) depicts the configuration space in mass-scaled and rotated coordinates z = Jx, with a spherical equipotential and transformed coincidence planes Z12, Z23, and Z31. Subfigure (c) depicts the project of (b) into relative coordinate z1-z2 plane. The sector with angle ω123 = π/3 is the Coxeter sector.
lighter than other two masses) to π/2 (m j much heavier). The angle ω ij,kl between coincidence planes Y ij and Y kl that do not share a particle remains π/2. In the limit g → ∞ the Hamiltonian (11) separates in hyperspherical coordinates with radius R 2 = i y 2 i . The interaction term [12] is proportional to 1/R, so it is not separable for finite values ofg ij , but asg ij → ∞ there is no distinction between 1/R or 1/R 2 times the sum of delta functions. So hyperspherical symmetry of (11) emerges and there is SO(2, 1) dynamical symmetry in the total hyperradial coordinate R [89]. 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 z = Jy so that the component z N ≡ Z is the scaled center-of-mass Z = i y i 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 ≡ JY ij and their normalsγ ij ≡ Jβ ij , but leaves the angles between planes like ω ijk and ω ij,kl invariant. Since all the normal vectorsγ ij have zero Z-components, the Hamiltonian in z-coordinates separates into H = H COM + H rel , where H COM is the Hamiltonian for a one-dimensional harmonic oscillator in the center-of-mass Z coordinate and the relative Hamiltonian is Finally, we go to hyperspherical coordinates in the relative space, where the relative hyperradius ρ is and there are (N − 2) angles charting the sphere S N−2 , conventionally chosen as Ω = {φ, θ 1 , . . . , θ N−3 } with φ ∈ [0, 2π) and θ i ∈ [0, π]. The relative Hamiltonian now becomes where ∆ Ω is the angular part of the Laplacian in relative configuration space. As before, the relative Hamiltonian (16) is ρ-Ω separable in the limitg ij → ∞. A general energy eigenstate can be separated into a product of center-of-mass ζ(Z), relative hyperradial R(ρ) and relative hyperangular Υ(Ω) functions where n is the center-of-mass quantum number, the function ζ n (Z) is the one-dimensional harmonic oscillator wave function, and ν is the relative hyperradial quantum number. At this point, λ is just derived from the judiciously-parametrized relative hyperangular separation constant λ(λ + N − 3) and µ is just an additional label to distinguish any possible degenerate states for a given λ. If there was no angular potentialg ij = 0, then λ would be a non-negative integer, the eigenfunctions on S N−2 would be the hyperspherical harmonics, and µ would be a collective index to label degeneracies [68,69]. However since there is an angular potential in (16), the hyperangular solutions are unknown and we must explicitly solve for λ, including any possible degeneracies. Whatever value λ takes (including non-integer values), the relative hyperradial function R νλ (ρ) is the standard solution for the radial factor of an (N − 1)-dimensional isotropic harmonic oscillator in hyperspherical coordinates [69] with energy ω(2ν + λ + (N −1)/2): where A ν,λ = 2 ν!/Γ(ν +λ+(N −1)/2). We have achieved our desired result: this series of coordinate transformations has reduced solving the Nparticle Hamiltonian (1) with equal frequencies and infinite-strength contact interactions into solving hardwall 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 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 Ω 12···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 The irreducible representations for O(m) and their realizations by hyperspherical harmonics are wellknown [90,91], and so we just summarize a few facts here for the readers' convenience. The subgroup 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 ∈ {1, . . . , m}, where L ij generates a rotation in the ij-plane. The quadratic Casimir of SO(m) is the sum of all of these generators squared For SO(3) this is the familiar angular momentum squared operator with eigenvalues λ(λ + 1). The SO(3) irreps are labeled by λ and have degeneracy 2λ + 1. In m > 3 dimensions, the operator L 2 is the hyperangular momentum squared operator with eigenvalue λ(λ + m − 2) and an irrep labeled by λ has degeneracy [68] Once the representation of total inversion is chosen, the irreps of SO(m) also naturally carry a representation of O(m) . For example, inversion is represented in the λ irrep by multiplication by (−1) λ for O(3).
To reduce an O(m) irrep λ into the irreps of Coxeter group G m , the G elements are sorted into conjugacy Table II. This table provides information about the conjugacy classes of H3. The first column is the Schönflies notation for the elements in the class Ki. The second column gives the angle of rotation φi of the element realized in O(3). The third column is whether it is generated by an even or odd number of reflections πi = ±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 Ki and the number of elements ki in the class.
Elements classes K i with k i elements. Each irrep W of G m has a unique pattern of characters χ W (K i ). In particular we are interested in the G m irrep W = A of all anti-invariant states, meaning χ A (K i ) = 1 when K i is a conjugacy class whose elements are an even composition of reflections and χ A (K i ) = −1 when K i is a class comprised of odd compositions. Further, each conjugacy class has a character χ λ (K i ) in the G m -reducible O(m)-irrep denoted by λ. When these characters are known, then the number of times the G m irrep A appears in the decomposition of the O(m) irrep λ is [90] a λ = 1 G Ki k i χ A (K i )χ λ (K i ).
Note that χ A (K i ) = (χ A (K i )) * because Coxeter groups are ambivalent. The number a λ is an integer that counts how many solutions there are with relative angular momentum λ. The projection operator onto the anti-invariant irrep A is given by where D λ (g) is the representation of group element g acting irreducibly on the d(λ)-dimensional representation space. If a λ = 1, any vector in the λ irrep space with a non-zero projection will be proportional to the solution we seek. If a λ > 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 II. The H 3 character χ A (K i ) for the anti-invariant irrep is +1 for the five even conjugacy classes and −1 for the five odd classes. The O(m) character χ λ (K i ) for irrep λ is where φ i is the angle of rotation and π 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 [92] to find which O(m) irreps have symmetric irreps of the spherical triangle groups in their reduction. The construct the actual states Υ λ (θ, φ), we use the projection operator (20) acting on the spherical harmonics. Instead of explicitly constructing the (2λ+1)×(2λ+ 1) unitary matrices D λ (g) that act on the spherical harmonics Y µ λ (θ, φ) 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 × 3 matrices O(g) ∈ O(3) that represent H 3 as rotation and reflections.
The real spherical harmonics can be written in terms of the relative coordinates: