Exactly solvable spin-1/2 XYZ models with highly-degenerate, partially ordered, ground states

Exactly solvable models play a special role in Condensed Matter physics, serving as secure theoretical starting points for investigation of new phenomena. Changlani et al. [Phys. Rev. Lett. 120, 117202 (2018)] have discovered a limit of the XXZ model for $S=1/2$ spins on the kagome lattice, which is not only exactly solvable, but features a huge degeneracy of exact ground states corresponding to solutions of a three-coloring problem. This special point of the model was proposed as a parent for multiple phases in the wider phase diagram, including quantum spin liquids. Here, we show that the construction of Changlani et al. can be extended to more general forms of anisotropic exchange interaction, finding a line of parameter space in an XYZ model which maintains both the macroscopic degeneracy and the three-coloring structure of solutions. We show that the ground states along this line are partially ordered, in the sense that infinite-range correlations of some spin components coexist with a macroscopic number of undetermined degrees of freedom. We therefore propose the exactly solvable limit of the XYZ model on corner-sharing triangle-based lattices as a tractable starting point for discovery of quantum spin systems which mix ordered and spin liquid-like properties.


I. INTRODUCTION
Calculations in condensed matter theory must generally bridge two different gaps. The first is the gap between model and experiment: any model simple enough to be successfully studied cannot capture every aspect of a real many body system, though we hope to capture the most important and interesting ones. The second is the gap between models and the actual calculation of physical quantities. Even when the model in question is simple to write down, it is usually necessary to employ some kind of approximation scheme during calculations.
Sometimes, however, the second gap is absent. There are models which are both physically relevant and for which exact calculations are possible. These have been important in the development of modern Condensed Matter physics, especially where they have been used to establish the theoretical possibility of novel phenomena or phases of matter. Notable examples of this include the Shastry-Sutherland model 1 , which established a featureless, gapped ground state in a many-body quantum spin model, and later found realization in SrCu 2 (BO 3 ) 2 2,3 and the Kitaev honeycomb model 4 , which a gave an example of a Z 2 quantum spin liquid with emergent Majorana fermions, and was later found to be relevant to various spin-orbit coupled magnets [5][6][7][8][9] .
An interesting example of an exactly solvable spin-1/2 model has been pointed out by Changlani et al. 10 . They considered the nearest neighbor XXZ model on the kagome lattice [ Fig. 1] and showed that at the special point J z /J ⊥ = −1/2 (denoted as the XXZ0 point 11 ) the model has a set of exact, 1. A three-color covering of the kagome lattice. Every triangle must include a site with each of three colors. The number of configurations satisfying this local rule grows exponentially with the system size. The exact ground states of the kagome XYZ model [Eq. (2)] map on to this threecoloring problem when the exchange parameters are chosen to obey Eq. (3). The coordinate axesxi,ŷi indicate a choice of local basis for the spins, such that Eq. (2) is consistent with the symmetry of the kagome lattice, with theẑ axis being uniformly perpendicular to the kagome plane. In contrast to the depicted coloring, this local spin basis is uniform across the lattice.
degenerate, ground states, the number of which grows exponentially with system size. These ground states can be written as simple product states, despite the fact that the Hamiltonian is composed of non-commuting terms and that general excited eigenstates are entangled. The ground states correspond with the solutions of the threecoloring problem on the kagome lattice 12 , in which the vertices of the lattice are colored with three different colors such that no triangle has two vertices of the same color [ Fig. 1]. For a non-trivial, quantum many-body model to have such a huge degeneracy of ground states, with such a simple structure, is remarkable. The XXZ0 point is also significant in the wider phase diagram, being a point at which several different phases, including distinct quantum spin liquids, meet. It was suggested that the phase diagram of the model in the surrounding parameter space could be understood starting from this point 10,13 .
Here we show how the XXZ0 point can be generalized to a wider set of exactly solvable anisotropic S = 1/2 exchange models. Specifically, we consider the nearestneighbor XYZ model on the kagome lattice: i.e., the case where each component of the spin has a different associated exchange constant. If the basis for the spin operators S α i is defined with a local coordinate frame, as shown in Fig. 1, then this model is consistent with the symmetry of the lattice and can be considered as a limit of a more general model of anisotropic exchange 14 , as shown in Appendix A.
We show that for any parameter set {J} = (J x , J y , J z ) fulfilling the conditions: there exists a large manifold of exact product ground states, which generally differs from the solutions of the XXZ0 model, but retains the correspondence with the three-coloring problem.
We further show that despite having an extensive number of undetermined degrees of freedom, the ground states possess infinite range correlations of some spin components, coexisting with algebraic correlations. This is qualitatively reminiscent of the phenomenon of "magnetic moment fragmentation" 15-17 studied in both kagome [18][19][20] and pyrochlore 15,[21][22][23] systems, but is particularly remarkable because it occurs in the exact ground states of a quantum model. This is in contrast to most known examples of moment fragmentation, which either occur in settings where the spins can be considered classical 15,18,19,23 or as a feature of the semi-classical dynamics rather than of the ground state 21,22 .
The Article is structured as follows: Section II shows the construction of the family of exactly solvable models and their highly-degenerate ground states; Section III demonstrates the partial order of these ground states; Section IV demonstrates the presence of gapless excitations above these ground states despite the absence of continuous symmetry in the Hamiltonian; Section V discusses how our construction can be generalized to other lattices and spin lengths and gives a criterion for the presence of exact ground states of the type discussed here; lastly, Section VI contains a brief summary of our results and outlook for future work.

II. CONSTRUCTION OF EXACTLY SOLVABLE HAMILTONIAN AND ITS GROUND STATES
In this Section we show how to obtain the result that the ground states of Eq. (2) can be found exactly for all sets of parameters satisfying (3).
Any nearest-neighbor Hamiltonian on the kagome lattice can be written as a sum of single-triangle Hamiltonians: The spectrum of the single triangle Hamiltonian H XYZ, is composed of a quadruplet with energy and of two doublets with energies where Λ = J x J y + J y J z + J z J x . When Λ = 0 and e q < 0 (cf. Eq. (3)), the spectrum is composed of a 6-fold degenerate ground state e q = e d− = e 0 , and of an excited doublet e d+ . In this case it follows that the Hamiltonian can be written as a sum of non-commuting projectors: with P being the operator that projects onto the pair of single-triangle excited states: Since e d+ > e 0 , the coefficient in front of the projection operator in Eq. (7) is positive and any state which is annihilated by all of the P is a ground state.
We can search for states annihilated by P amongst the set of site-product states: As we show below, there are many product states on the lattice which are annihilated by P . Our strategy is to first identify these product states on a single triangle (Section II A) and then generate ground states on the lattice by tiling single triangle ground states across the system.
A. Product ground states on a single triangle For a single triangle, the possible product wave functions [Eq. (9)] are parametrized by 6 variables: the θ j and φ j of each of the three sites. To be ground states, these wavefunctions need to satisfy four constraints: the real and imaginary parts of their overlaps with |d+, ± [Eq. (8)] must vanish. This suggests a 6 − 4 = 2 dimensional surface of exact product state ground states for a single triangle. We have verified that such a surface, that we shall call M, indeed exists for all {J} obeying (3). An exact parameterization of M has also been found, but its derivation is quite involved and is therefore presented in Appendix B.
Nonetheless, the solutions for two limits of (3) are easily found, and the general case can be understood as an interpolation between these two limits [ Fig. 2]. The first limit is the XXZ0 limit 10 J x = J y = 2/3, J z = −1/3 in which case M is composed of two spheres: These two spheres meet at the points θ 1 = 0 and π corresponding to the all-up and all-down states, as illustrated under Fig. 2 a). The other easily solved limit is the Ising limit: J x = J y = 0, J z = 1. In this case, product ground states have θ = 0 (spin fully up) on one site and θ = π (spin fully down) on another, with θ, φ on the remaining site completely free. The product ground state manifold M is thus composed of six spheres, each connected to two other spheres at the points where two spins have θ = 0 or θ = π. This is illustrated in Fig. 2 b). The six spheres correspond to the six ways of assigning an up spin, a down spin, and a free spin to the three sites of a triangle. The surface of each individual sphere corresponds to the direction of the expectation value of the free spin.
The manifold for more general exactly solvable {J} interpolates between these two limits by smoothing out the singular points where the spheres meet so as to produce a manifold with the topology of a torus. This interpolation can be made precise by parameterizing the Fig. 2 a)], κ = π/3 an Ising point [ Fig. 2 b)], and the M for generic κ ∈ 0, π/3 smoothly deforms as we vary κ from one limiting case to the other [ Fig. 2 c)]. In between, M has the topology of a torus, and this torus pinches at two (six) points as κ → 0 (κ → π/3).
We have checked numerically that the surface of single triangle solutions for general parameters has genus g = 1 (see Appendix C and Supplemental Material 24 ). This is itself an unusual situation, which can be contrasted (e.g.) with the exact ground states of a Heisenberg ferromagnet which cover the surface of a sphere.
Having obtained the exact solutions on a single triangle, the question is then how much freedom there is to Schematic representation of the topology of the product ground state manifold of a triangle M for: a) thê z-XXZ0 point (κ = 0, Jx = Jy = 2/3, Jz = −1/3), b) thê x-Ising point (κ = π/3, Jx = 1, Jy = Jz = 0), and c) generic κ ∈ 0, π/3 between these two limits. The M of the XXZ0 point is made of two spheres with different chiralities, denoted + and −, that touch at their poles. The points of contact are marked with red dots, and the respective spin configurations are written on top of these contact points. The Ising point M is made of six spheres that have one ↑ spin, one ↓ spin, and one free spin, which we denoted with a 0. These six spheres touch adjacent spheres at the specified red points. In between, M has the topology of a torus. Although drawn here as (deformed) circles, lines, etc., all of the above sketches should be understood as representing 2D surfaces embedded in the 6D space of all possible θj, φj in a triangle. tile these solutions over the whole lattice. It turns out we can choose any direction for the first S i and still find configurations for the surrounding spins such that the triangle and system are in a ground state. Specifying the first S i will in general remove the continuous freedom identified for the single triangle in II A, leaving only a discrete set of possibilities for the neighboring spins to remain in the ground state. If we fix S i to S 0 and consider one of the two triangles connected to the site i, then the remaining spins j, k on that triangle can take only one of two configurations. These two configurations are related to one another by a permutation symmetry that swaps the remaining pair of sites, with the values S 1 and S 2 being fixed by S 0 .
Propagating this throughout the system, consistency between triangles will force all sites to take one of the three expectation values (S 0 , S 1 , S 2 ). Some examples of allowed triads of expectation values (S 0 , S 1 , S 2 ) for J x /J z = 1/4, J y /J z = −1/5 are shown in Fig. 3. Provided that the first triangle ijk was in a ground state, the permutation symmetry of H XYZ guarantees that any triangle where each of the three expectation values is represented once is also a ground state.
Provided that (S 0 , S 1 , S 2 ) are all distinct from one another (which is true for generic members of the ground state manifold), these conditions are in precise correspondence to the rules of the three coloring model. Thus, the set of product state solutions to H XYZ,Λ=0 is given by the space of three-color configurations, along with two continuous degrees of freedom which were used up by fixing the first spin.
This establishes our main result: that H XYZ with parameters chosen to obey Eq. (3) has a set of exact, product, ground states, the number of which grows exponentially with system size. The wavefunctions described by Eq. (9) are not all linearly independent, so the actual ground state degeneracy is not the same as the number of product state solutions. Nevertheless, if the number of product state solutions grows exponentially with system size, it should be expected that the number which are linearly independent also grows exponentially. This was verified for the case of the XXZ0 model in Ref. 10.
We have checked using exact diagonalization on small clusters that there is indeed a collapse of many excited states towards the ground state approaching the points in parameter space given by (3), consistent with the establishment of a macroscopic ground state degeneracy in the thermodynamic limit. This is shown for a 24-site cluster in Fig. 4.
It is possible that there are also additional ground states, beyond the product ground states found here. Changlani et al. found numerically that there are indeed some additional ground states in the XXZ0 limit for lattices with periodic boundaries, although not with open boundaries 10 .
The set of exactly solvable XYZ models described by (3) includes both the XXZ0 model ({J} ∝ ( 2 3 , 2 3 , −1 3 ) and permutations) and the Ising model ({J} ∝ (0, 0, 1) and permutations), three times each. The Ising model represents a singular limit for the set of product ground states because if two spins on a triangle are fixed to be one up and one down, the third spin is actually completely undetermined, increasing the freedom in the construction of ground states on the lattice.

III. PARTIAL ORDERING IN THE GROUND STATES
In this Section we argue that the exact ground states identified in Section II possess infinite range correlations, coexisting with the algebraic correlations implied by the three-color mapping and are in this sense partially ordered.
The allowed ground state configurations of spin expectation values on a single triangle generally have a finite value of m = S 0 + S 1 + S 2 . This is illustrated in Fig. 5, where the minimum and maximum value of |m| over the possible single triangle ground states is plotted as a function of J y /J z . Note that m would not be the same as the magnetisation in a real system, because of the local basis used to define the spins [ Fig. 1]. Maximum and minimum possible ground state values of |m| = | S0 + S1 + S2 | within the set of exact product ground states on a single triangle as a function of Jy/Jz, for sets of parameters obeying condition (3). Apart from at the XXZ0 point (Jy/Jz = −1/2), all allowed product ground states have a finite m. Once m is fixed on the first triangle it will be the same on all triangles throughout the lattice, apart from in the Ising limit, (Jy/Jz = 0) where some individual spins can be rotated freely. This establishes that, for general parameter sets obeying (3), the ground states carry infinite range spin correlations, coexisting with the algebraic correlations implied by the three-color mapping.
Since the remaining triangles in the lattice must only feature permutations of the spin expectation values on the first triangle (see Section II B), and since m is invariant under those permutations, we conclude that all triangles must have the same value of m. This implies that M = m has a macroscopic value in any given product state and that there are infinite range correlations of the spin component parallel to m. These infinite-range correlations must coexist with the algebraic correlations that are also present due to the mapping between the spin configurations and the three-color model.

IV. GAPLESS EXCITATIONS
In addition to the many ground states, this exactly solvable model also possesses gapless excitations.
We demonstrate this below by making use of the continuous degeneracy of the ground states identified in Section II A, and a trial wave function for the excitations based on the single mode approximation. Starting from a translationally invariant member of the set of ground states one can construct excitations using the generators of rotations within the single triangle ground state manifold, and show that these are gapless in the long wavelength limit.
Let us emphasize that the considerations of this section depend only on the continuous degeneracy of the ground states, and not on the discrete degeneracy of how one can tile the one triangle solution across the whole lattice. Thus the argument of this section would work equally well for (e.g.) exactly solvable XYZ models on the triangular lattice, which would retain the continuous degeneracy of the kagome case considered here, but for which only 6 consistent three-colorings are possible.
The gapless excitations are similar to Nambu-Goldstone modes, except that unlike ordinary Nambu-Goldstone modes the associated symmetry generator G does not commute with the Hamiltonian but rather the commutator annihilates a ground state: One way of motivating the above is to consider a oneparameter family of ground states |Ψ t that we know exists because of the continuous degeneracy. Then it is always possible to find an unitary operator U t that maps |Ψ 0 → |Ψ t . We may then identify the quasisym- (13) follows. Qualitatively, we may say that G becomes an exact symmetry only in the zero-energy limit.
Of course, it may be the case that some of the ground state quasi-symmetry generators are also genuine symmetry generators. The XXZ0 point, considered in Section II A, is a good example. From (10) it is evident that a global rotation around theẑ axis maps a ground state to a ground state. In fact, this is an exact symmetry of the XXZ Hamiltonian. On the other hand, rotating the spins by θ i → θ i + δθ also maps a ground state to a ground state, but is not a genuine symmetry. Note how in the case of the θ i → θ i + δθ transformation, the generator G depends on the starting ground state |Ψ 0 .
To begin with, we choose a member of the set of product ground states to construct excitations around. In such a ground state, each spin expectation value takes one of three values S r , S g , and S b , according to whether the site is red, green, or blue in the three-color representation of the state. |Ψ 0 is chosen to be the state corresponding to a translationally invariant threecoloring of the lattice. Later we comment on more general red-green-blue patterns.
The existence of a continuous two-dimensional manifold of exact product ground states on the single triangle implies the existence of infinitesimal global rotations which keep every triangle, and therefore the system as a whole, in a ground state.
There are two independent generators of these rotations for every starting ground state (spin triad). We will pick one of them, which we label G = i g i . The site generator g i depends on the color of site i in the three-color tiling of the ground state.
We now do a site-dependent coordinate transformation on the spin basis {x, y, z} → {u i , v i , w i }, with the new coordinate axes chosen such thatŵ i aligns with the expectation value of spin i in the ground state and so that the generators g i are proportional to rotations aroundû i .
The spin operators in this basis then satisfy: where U i is a real scaling factor that depends only on the coloring of the site i. (In general, U i also depend on the three spin vectors S r , S g , and S b , but these are the same across |Ψ 0 .) It is normalized according to Since the localû i axes have been chosen to correspond with the rotation axes that keep the system in a ground state, it must be true that: We then consider the following variational wavefunction for the excitations, based on the single mode approximation, where N is the number of unit cells. The wavefunction |ex, q describes a spin wave-like excitation. |ex, q is orthogonal to |Ψ 0 because Ψ 0 |S u i |Ψ 0 vanishes everywhere. At finite q and in the thermodynamic limit N → ∞, it is also orthogonal to the other members of the ground state manifold. Because of this, the expectation value of the energy in |ex, q is an upper bound on the energy of excitations with momentum q: The denominator Ψ 0 |G(−q)G(q)|Ψ 0 = 1 4 , and we can set E 0 = 0 (i.e., measure all energies relative to the ground state), so that Due to Eq. (17), the q → 0 limit of Eq. (21) vanishes: implying gapless excitations. The part of Eq. (21) which is linear in q also vanishes: which also follows from Eq. (17). This leads us to the conclusion that the dispersion of excitations has a quadratic upper-bound at small q: This also agrees with a linear spin wave analysis which finds quadratically dispersing excitations around the translationally invariant exact ground states.
Although not obvious, a detailed analysis of the properties of the coupling matrix J ij in the new basis {u i , v i , w i } shows that the second generator of symmetry has g i = U i S v i with the same scaling factors U i from Eq. (15). Thus one finds that the commutator [G, G ] has a non-vanishing expectation value in the ground state |Ψ 0 , implying that the two generators represent only one degree of freedom (cf. [x, p] = i ) 25 . In the spin wave analysis this is reflected in the fact that there is only one gapless mode, despite the two broken quasi-symmetry generators.
For more general coloring patterns (that are not translationally symmetric), G(q) may still be used to probe the low-lying excitations, but this time q cannot be identified with the momentum. The above argument thus suggests that low-lying excitations still exist even for more general RGB patterns.
With that said, one should keep in mind that the above argument may fail if the upper bound from (20) is discontinuous at q = 0, (see e.g. Supplementary Material of [26]).
To verify that our upper bound is continuous, we evaluate it by noting that only averages of the form Ψ 0 |S u i S µ i S w i+δ |Ψ 0 for µ = u, v are non-vanishing, giving: From the fact that |Ψ 0 is an exact eigenstate, it follows that in the new basis the exchange coefficients satisfy The exchange coefficients J ij µν also depend only on the coloring of the sites i and j and satisfy J ij µν = J ji νµ . Moreover, from the fact that the change in energy within one triangle is to second order equal to zero when we vary the spins with G, it follows that: Applying to Eq. (25), we obtain: which is indeed continuous and vanishing in the limit q → 0. We therefore conclude that the partially ordered exact ground states of the XYZ model have gapless excitations, despite the absence of continuous symmetry in the original Hamiltonian.

V. RECIPE FOR CONSTRUCTING FURTHER EXACTLY SOLVABLE MODELS
Our work exposes a simple recipe for the construction of further highly degenerate exactly solvable models, of the same kind as those discussed here.
First, one must define a Hamiltonian for quantum spins of length S on a set of corner sharing units (for instance, triangles or tetrahedra), with the Hamiltonian on each unit being symmetric under permutations of the sites. Then one must tune the parameters such that the degeneracy d 0 of the ground state of the single unit Hamiltonian is large enough that product wavefunctions on the single unit have enough free parameters to be made orthogonal to all of the excited states.
If n is the number of sites in each unit then, the total number of states in the single unit spectrum is (2S + 1) n and a single unit product state has 4Sn degrees of freedom. Requiring that the real and imaginary parts of the overlap with every excited state in the single unit spectrum vanishes gives 2((2S + 1) n − d 0 ) constraints, so to be able to find product-like solutions we need: If the single unit Hamiltonian can be tuned to have a sufficiently large d 0 then one can search for product-like ground states. If, in these ground states, all n sites in the unit are distinguishable from one another (e.g., if the spin expectation values differ on the sites), and the singleunit Hamiltonian has a permutation symmetry under the swapping of sites, then the problem of tiling solutions over the whole lattice becomes an n-coloring problem. We anticipate that this recipe can be used to construct further exactly solvable models on other lattices. In the case of S = 1/2 spins on the kagome lattice, the minimal d 0 is 5 (cf. d 0 = 6 for the models considered in this manuscript). In a model with d 0 = 5, there would be no continuous degrees of freedom left in the ground state, but only discrete degrees of freedom from the threecoloring pattern.
The exactly solvable "two-coloring" models identified in Ref. 27 can also be seen as an example of the construction described above, with the corner-sharing units being single bonds (n = 2) and the requirement for exact solvability d 0 ≥ 2. The "two-coloring" models in Ref. 27 have d 0 = 3.

VI. SUMMARY AND OUTLOOK
We have demonstrated that the exactly solvable model pointed out in Ref. 10 is a member of a wider family of exactly solvable models, and that the ground states of these models possess coexisting infinite-range and algebraic correlations. The combination of the large ground state degeneracy, and the coexistence of infinite-range and algebraic correlations suggest an analogy with the phenomenon of "magnetic moment fragmentation" [15][16][17][18][19][21][22][23] but the case here is distinguished by the fact that it occurs in the exact ground states of a quantum model.
Our conclusions generalize straightforwardly to XYZ models on other triangle-based lattices, such as the hyperkagome lattice.
The models described here, fall into a wider category of frustrated systems possessing a large number of low energy states, defined by local constraints. In such cases, a description of the low energy physics in terms of emergent gauge fields and charges is frequently useful [28][29][30] , and this may be an interesting avenue to investigate further for the present case.
The exact ground states we have defined are not themselves quantum spin liquids, since they lack quantum entanglement and have a large non-topological degeneracy. Nevertheless, highly degenerate points of a model are often a good starting point for discovery of spin liquids because small perturbations can stabilize a variety of non-trivial superpositions of the degenerate states. In this case, given the partial ordering of the ground states, perturbing the models identified here may be a way to stabilize phases which combine the interesting features of a quantum spin liquid (entanglement, fractional excitations) with spontaneous symmetry breaking. With this in mind, we suggest that a numerical study of the ground states of H XYZ for parameter sets close to the exactly solvable limit may be a rewarding subject for future work.
In view of the large number of ground states, and low lying excited states, the physics of the exactly solvable models at finite temperature is also an interesting topic for the future. In particular, it is possible that at T > 0 thermal order-by-disorder will lead to a breaking of the ground state degeneracy. Quantum (T = 0) order-bydisorder is conclusively ruled out for the models described in this work, because the exact ground states and their energies are known. It may nevertheless be the case that thermal fluctuations can distinguish between the states, leading to an entropic selection at finite temperature. The low energy gapless modes discussed in Section IV are likely to play a key role in this selection, if it occurs.
A further unusual feature of the XYZ model presented here is that even the solutions to the single triangle problem are topologically non-trivial, in that the continuous manifold of single-triangle solutions has the topology of a torus. Whether interesting phenomena can be derived from this by, for example, adiabatically moving the system around this manifold is something which remains to be seen.
Finally we note that these exact ground states may also be of interest from the point of view of non-equilibrium physics. With a simple alteration of the Hamiltonian 31 the exact ground states can be turned into many-body quantum scars: highly excited states which violate the eigenstate thermalization hypothesis. The XYZ models proposed here offer a chance to explore this in a setting which lacks continuous spin rotation symmetry.
Appendix A: Relationship between XYZ model and more general anisotropic exchange models on the kagome lattice Here we elaborate on the relationship between the XYZ model studied in the main text and the model dicussed in Ref. 14. The model discussed in Ref. 14 gives the most general set of nearest neighbour interactions consistent with the symmetries of the lattice, including a mirror symmetry in the kagome plane itself. In total, the symmetries considered comprise the inversion and translation symmetries of the kagome lattice along with the three point group symmetries shown in Fig. 6.
The Hamiltonian of this model is: The three coupling matrices on the triangle, when written in the global basis (x (g) ,ŷ (g) ,ẑ) [ Fig. 7], have the form: with the numbering convention for sites on a triangle as shown in Fig. 7.
Transforming to a local coordinate system, with local axes in the xy-planex i ,ŷ i on each site of the triangle as shown in Fig. 7, while maintaining the same globalẑ axis, the coupling matrices become uniform: In the special case D z = √ 3J ⊥ , this becomes an XYZ model with This establishes that the XYZ model discussed in the main text occurs as a limit of the general symmetryallowed anisotropic exchange Hamiltonian for kagome magnets.

Appendix B: Derivation of the exact ground states
Here we show how to obtain analytic expressions for the ground states in the exactly solvable limits of the XYZ model.
The ground state wave functions are product states of the form Eq. (9) Such a wave function is completely determined (up to a global phase) by the expectation values of the spin components on each site: The vector S i = ( S x i , S y i , S z i ) is constrained to lie on a sphere: In the following, we find analytic expressions for S i for all product state ground states of the single triangle Hamiltonian, at the exactly solvable points of the XYZ model. The ground states on the lattice follow from this by tiling the lattice with single triangle solutions in the way described in the main text. The exactly solvable parameter space is given by J x J y + J y J z + J z J x = 0 and J x + J y + J z > 0. Setting J x + J y + J z = 1 as the unit of energy we can then write the exchange parameters as a function of a single variable κ: The points κ = 0, 2π/3, and 4π/3 correspond to the XXZ0 points along theẑ,ŷ, andx axis, respectively, whereas κ = π/3, π, and 5π/3 correspond to the Ising points along thex,ẑ, andŷ axis, respectively. If we attribute β = x → 1, y → 2, and z → 3, then the above can be compactly written as: The above parameterisation has the nice property that it allows us to focus on the parameter range κ ∈ [0, π/3], with the rest related to this range by a permutation of the {J} coefficients and spin components. In particular: This permutation symmetry (of the x, y, z components) relates systems with different exchange coefficients and should not be conflated with the dynamical permutation symmetry of exchanging the spins on a triangle, whilst keeping {J} the same. The latter we used in Section II B to derive the tiling rules.
Since the XXZ0 and Ising points were analyzed in the main text under Section II A, below we focus on the case of generic κ ∈ 0, π/3 .

Method of Lagrange multipliers
For generic values of κ, not equal to Zπ/3, we find the ground state by minimising the expectation value of the energy, using Lagrange multipliers to enforce Eq. (B2). If this expectation value Ψ|H XYZ |Ψ coincides with the known ground state energy e 0 = − 1 4 (J x + J y + J z ) = −1/4, then |Ψ is a ground state.
The energy of the product state can be written directly in terms of the spin expectation values: To simplify the notation, we will from this point forward denote S β i as S β i and treat them as classical variables. The problem is then to minimize E Ψ with respect to the nine variables S β i (three components β on each of three sites of a triangle i), while respecting the three constraints Eq. (B2). This can be achieved with the method of Lagrange multipliers.

(B21)
Since we're using polar coordinates to express µ j in Eq. (B17), (r, u) and (−r, u + π) represent the same µ j and the n = 1 and n = 2 radial solutions are equivalent.
(vi) Lastly, for a fixed n, r n (u, κ) and µ j (u, n, κ) as functions of (u, κ) are smooth everywhere on their domains of definitions, excluding the cusps listed under (v). This is a stronger statement than being smooth in only u or κ, and it follows from the structure of Eqs. (B21) and (B17). With that said, in the arguments that follow we shall only be needing the piecewise continuity and differentiability in u of these functions.
With the values of µ j determined, the vector of spin componentsS [Eq. (B12)] is given by a linear combination of the null vectors of M. Row reduction of M β gives where in the second step we have used Eq. (B14) and assumed that at least two i out of three satisfy J β + µ i = 0. It turns out that for all κ, there are special values of u, denoted u * , for which two J β + µ i (with the same β) vanish, invalidating the second step above. For generic κ = Zπ/3 (i.e., not Ising or XXZ0 points), however, these special u * are restricted to a discrete number of points. Although not obvious, allS solutions of these special u * can be obtained by a limiting procedure of the generic u solutions.
Below we first determine when J β + µ i = 0 and find that these special u * are, for generic κ, given by u * = Zπ/3. Then we find the spinsS at these special u * explicitly to confirm that there are no points isolated from the generic u solutions. Lastly, we find theS for generic u.

Finding the special u *
In the coming two sections, the functions J β + µ i shall play a prominent role, so we introduce a βi (u, n, κ) = 1 2 J β (κ) + µ i (u, n, κ) .
Using ∂ u µ 3 (u * ) = 0 one can find the special u * without solving any complicated algebraic equations that include the expression (B21). First, we differentiate the cubic equation (B18) and solve for ∂ u r, Next, we differentiate the definition (B17) of µ j and use ∂ u µ 3 (u * ) = 0 to get By combining the above two equations, one obtains the desired equation (r(u * ) cos(u * ) + 1) sin(u * ) = 0 .
From the previous section we know that |a β1i | = |a β1j | = 0, |a β2i | = |a β2j | = 0, and |a β1k |, |a β2k |, |a γk | = 0 are all non-zero. In light of Eq. (B22), the null vectors of M β1 and M β2 are given by whereas the null vector of M γ has the components SinceS is a null vector of M [Eq. (B11)], we may write it as:S where the coefficients X β are fixed by the spin normalization constraints [Eq. (B2)]. Because of the relations among {a βi }, only two of the three spin normalization constraints are independent, to wit Eq. (B31) defines an ellipsoid and Eq. (B32) a cylinder.
These two intersect at: and the spins are given by where the rows represent the β 1 , β 2 , and γ components of the spins. What these equations represent depends on the ratios (a β1k /a β1i ) 2 and (a β2k /a β2i ) 2 . We have three cases to consider: • When (a β1k /a β1i ) 2 > 1 and (a β2k /a β2i ) 2 > 1, the X γ are imaginary and there are noS that satisfy the spin normalization constraints. Said differently, the ellipsoid [Eq. (B31)] and cylinder [Eq. (B32)] do not intersect. On 0 < κ < π/3, this is the case for the n = 0 solutions on all u * = Zπ/3.
Lastly, let us comment on how the i, j, k, β 1 , β 2 , γ are determined. Since from previous study we know that |a βi | are continuous functions of κ and that only at the special Ising and XXZ0 points may additional |a βi | coincide or vanish, it follows that it is sufficient to determine the i, j, k, β 1 , β 2 , γ for a given u * ∈ Zπ/3 at one κ ∈ 0, π/3 to know these indices across this whole range. That is, because of discreteness and continuity, these indices may only change at κ = Zπ/3.
The allowed values of X β are determined by the intersection of 3 ellipsoids that follow from the spin normalization constraints: with semi-axis lengths |a βi | [Eq. (B23)], known to all be non-zero for u = Zπ/3. One might expect Eqs. (B40) to give eight solutions for X β for each κ, u, n given by However, the matrix A actually has a vanishing determinant, signifying that one of the three constraints in Eq. (B40) is linearly dependent on the other two.
In the previous section, we have already seen this linear dependence make the i, j constraints of Eq. (B31) identical. Although we are presently unable to give an analytic proof of the statement det A = 0 for generic u, we have verified it numerically for the full spectrum of possible values of κ, u, n and it holds for all cases.
Given that only two of the constraints in Eq. (B40) are linearly independent, there exist one-parameter families of solutions for X β . We will parametrize these families with a continuous parameter v. Combined with the fact that u is also a continuous parameter, this means that for all sets of generic exchange parameters {J} [Eq. (B4) with κ = Zπ/3], the set of ground states has two continuous free parameters, u and v.
As previously established, |a βi | are piecewise continuous functions of u and κ that can (for generic κ) intersect only at the special u * = Zπ/3 points. Thus if we consider, for instance, u ∈ 0, π/3 and κ ∈ 0, π/3 , then for this whole range all |a βi | differ and, moreover, the corresponding ellipsoids must intersect, or not, in the same way. When one considers the n = 0 solution, one must also pay attention to the divergences of |a βi | that happen at u div = π/6 + Zπ/3.
After analysing the various cases for 0 < κ < π/3, one finds that the n = 0 ellipsoids with semi-axes |a βi (u, n = 0, κ)| never intersect. This agrees with the previous section where the n = 0 case also could not satisfy the spin normalization constraints (the cylinder and ellipsoid did not intersect).
If we consider theS that the above solutions for X β , when combined with Eqs. (B1), (B12), and (B30), give, then one finds that for a fixed sign of X x these solutions are discontinuous at u = Z2π/3. The cause is the fact that X x always has the same sign, whereas a xi , and therefore S x i , change their sign at u = Z2π/3. If one tries to find a sign convention that preserves continuity, then one finds that the six patches of length 2π/3 that Eqs. (B48)-(B50) yield combine in to one periodic patch of length 4π.
In addition to the above, if one considersS for a fixed v as a function of u, the + sign convention in Eq. (B49) again gives discontinuities, but this time at u = π/3 + Z2π/3 and in the S y i components of the spins. A simple sign change eliminates these.
In summary, the X β , now named X, Y, Z, that give a continuous parameterisation of the ground state manifold are X(u, v, κ) = sgn [sin(3u/2)] a x1 (u, κ) × The above expressions are valid even at the special u points if one takes a careful limiting procedure u → u * ∈ Zπ/3. What one obtains through such a limiting procedure agrees with the solutions derived in the previous section explicitly. This confirms that the above solutions exhaust the manifold of all product state ground states for κ = Zπ/3.
In the Supplementary Material 24 , the reader may also find animations of the parameterisation from Eqs. (B54)-(B56) for κ = 0.1, π/6, 0.9. "In the animations, two rows are drawn with three spin spheres and one M torus in each. In the upper row, lines of constant u are drawn whose coloring indicates how v is varied along them. In the lower row, v is held constant and u is varied." • Then we find the best fitting plane of these closest points through a singular-value decomposition, and project all the ∼ 32 closest point on to this plane.
• Next, we find the Delaunay triangulation of this 2D sample of points and add the edges and triangles of only those points that are adjacent to i to the manifold triangulation.
• Lastly, we add these points adjacent to i to a queue, and repeat the above three steps for all points in the queue. The only difference is that now we have to carry out a constrained Delaunay triangulation to ensure that we respect the edges of the previous local-plane triangulations.
A demonstration of our numerical routine, together with an application to an already computed mesh of the ground state manifold M, is given in triangulation demonstration and verification.ipynb. Our routine detects various cusps or irregularities in the manifold or mesh by calculating the standard deviations orthogonal to the best fitting plane. If they are large, a warning is given.
In addition, it is also a good idea to use a mesh that is dense and uniform. Uniformity can be achieved by relaxing the mesh under the influence of short-ranged repulsive forces.
When the triangulation finder.py routine is applied to a sufficiently dense mesh of the ground state manifold M, for κ not too close to the Ising or XXZ0 points, one finds that M is orientable and has genus 1. Thus M is a torus.