Quantum spin liquid in the semiclassical regime

Quantum spin liquids (QSLs) have been at the forefront of correlated electron research ever since their proposal in 1973, and the realization that they belong to the broader class of intrinsic topological orders. According to received wisdom, QSLs can arise in frustrated magnets with low spin S, where strong quantum fluctuations act to destabilize conventional, magnetically ordered states. Here, we present a Z2 QSL ground state that appears already in the semiclassical, large-S limit. This state has both topological and symmetry-related ground-state degeneracy, and two types of gaps, a “magnetic flux” gap that scales linearly with S and an “electric charge” gap that drops exponentially in S. The magnet is the spin-S version of the spin-1/2 Kitaev honeycomb model, which has been the subject of intense studies in correlated electron systems with strong spin–orbit coupling, and in optical lattice realizations with ultracold atoms.


Supplementary Notes
Supplementary Note 1. Semiclassical expansion around the states of the star dimer pattern Lattice superstructure & Hamiltonian Here we provide the details of the semiclassical expansion around the states of the star dimer pattern of Fig. 2 of the main text. The analysis is based on the six-sublattice decomposition shown in Fig. 8 of the main text, but we repeat the details here for completeness. The six-sublattice decomposition is shown again in Supplementary Figure 1, with a superlattice defined by the primitive translation vectors T 1 and T 2 . Any given site i of the lattice can be labeled as i = (R, ν), where R is a primitive vector of the superlattice and ν = 1-6 is the sublattice index. In this parametrization, the positions of the empty hexagons h α are labeled by R. The classical state is parametrized in terms of the η-variables, as shown in Fig. 2 of the main text. We will also use the local coordinate frames given in Eq. (3) of the main text, and define for each empty hexagon h R γ R ≡ κB R . (1) With these conventions and definitions, the Hamiltonian reads

Semiclassical expansion
In our semiclassical expansion we will keep up to four boson terms. So it suffices to keep the following terms from the standard [1] Holstein-Primakoff expansion for each site i = (R, ν): where c i , c + i are bosonic operators. We have: • Terms of the type S u i S u j (where i and j belong to the same empty hexagon): • Terms of the type S v i S v j (where i and j belong to the same empty hexagon): • Terms of the type S w R,ν S w R+Tνµ,µ : The terms that couple different empty hexagons are of the form S w R,ν S w R+Tνµ,µ . For simplicity, we will label (R, ν) → i and (R + T νµ , µ) → j. We have: The quartic term decouples as follows: We now repeat the main arguments described in the Methods section of the main paper that simplify this expression. The state around which we expand does not break the local BSS flux operators defined on the empty hexagons: = (−1) λ R S exp − iπ κ η 1 n R,1 + η 3 n R,3 + η 5 n R,5 + η 2 n R,2 + η 4 n R,4 + η 6 n R,6 , where n i = c + i c i is the boson number operator and λ R = κ(η R,1 + η R,3 + η R,5 ) + (η R,2 + η R,4 + η R,6 ), see main text. The invariance of the Hamiltonian and the state around which we expand under this operation translates into the invariance of the parity of the number κ η 1 n R,1 + η 3 n R,3 + η 5 n R,5 + η 2 n R,2 + η 4 n R,4 + η 6 n R,6 . But since κ and η can only take the values +1 and −1, it follows that the parity of this number is the same as the parity of the total number N R of bosons in any given empty hexagon: This means that terms that change the parity of N R are not allowed in the expansion. This excludes terms of the type c i c j or c i c + j , where i and j belong to different empty hexagons (see definition above). Equivalently, the meanfield parameters m ij and δ ij vanish by symmetry, and this is true to all orders in the Holstein-Primakoff expansion. We therefore get: i.e. empty hexagons decouple from each other and the S w i S w j terms give, for each empty hexagon R alone, a contribution where the constant p j = p R+Tνµ,µ refers to a neighboring hexagon and has to be found self-consistently in the general case.
It is useful to add here one more consequence of the BSS flux conservation. In the classical, reference state, where all n i vanish, the BSS fluxes are equal to (−1) λ R S (see main text and [2]). Spin wave fluctuations dress the reference state but cannot change the BSS fluxes, because these are integer numbers. (13) then implies that the dressed ground state contains only terms with an even number of bosons N R .

Mean field parameters: General relations
Let us define the six eigenvectors of the matrix g · M that correspond to non-negative eigenvalues by X ν , ν = 1-6. Using: we get the following expressions for the mean-field parameters: where X ν denotes the ν-th eigenvector of g · M. Note that the last expressions in each line do not depend on the arbitrary phase for the eigenvectors X ν , which come out arbitrary when we diagonalize the matrix g · M numerically.

Mean field parameters: Symmetry constraints
We have already mentioned that all mean-field parameters defined above are real quantities. Here we give a list of symmetry operations (of the Hamiltonian and of the classical state around which we expand) which reduce strongly the number of independent mean-field parameters.
• Symmetry Σ 1 . This is a π-rotation in real space around the center of the hexagon, followed by π 2 -rotations around the local w-axes in spin space: These relations are equivalent with • Symmetry Σ 2 . This is a reflection through the bonds (3,6) in real space, followed by π 2 -rotations around the local w-axes in spin space: These relations are equivalent with: • Symmetry Σ 3 . This is a reflection through the middle of the bonds (1,2) and (4,5) in real space, followed by zero or π-rotations around the local-w axes in spin space: These relations are equivalent with: • Symmetry Σ 4 . This is a reflection through the bonds (1,4) in real space, followed by π 2 -rotations around the local w-axes in spin space: Combining Σ 1 -Σ 5 gives the following constraints for the mean-field parameters: The mean field parameter m The numerical, self-consistent treatment of the decoupled spin-wave Hamiltonian gives a vanishing mean-field parameter m. This result does not arise from symmetry and is true only in the asymptotic large-S limit. For general S, m is a very small number. To see this we consider the self-consistent mean-field Hamiltonian for a single hexagon, that corresponds to the decoupled semiclassical problem that we are dealing with: where h loc is the self-consistent field exerted from neighboring hexagons and we have taken K = 1 without loss of generality. In what follows we shall use the Néel operator L defined as and the relations • For S = 1/2, the numerical, self-consistent solution gives h loc = 0.37888 and m = 0. However, this relation is special to S = 1/2 because the self-consistent ground state |g of H MF has the special property L|g = 0. And according to the above relations, this implies that g|S + 1 S − 2 |g = 0, which is equivalent with m = 0.
• For S = 1 and higher, the ground state does not obey the property L|g = 0 and m is therefore finite. The numerical solution for S = 1 gives h loc = 0.83643 and m = 0.0011412, which is a very small number.
• In the large-S limit, the parameter m must eventually vanish (consistent with the numerical results from the decoupled, large-S spin-wave Hamiltonian). The reason behind this is that as we increase S, the ground state |g comes closer and closer to the classical vacuum |0 (with spins fully polarized along their local w-axes), which has the property L|0 = 0 (because |0 is an eigenstate of each S w ν individually). In fact, this relation remains true when we include the leading effect of semiclassical corrections coming from V. At this leading level, the ground state wavefunction is given by [4] where R = 1−|0 0| E 0 −H 0 is the usual resolvent operator. To show that L|g 1 = 0 we use the fact that L commutes with H 0 (and therefore with R as well) and furthermore L|0 = 0. These properties give: We further have: S + 1 S + 2 − S + 2 S + 3 + S + 3 S + 4 − S + 4 S + 5 + S + 5 S + 6 + γS + 6 S + 1 + h.c.
At higher orders n > 1, the ground state |g n does not satisfy this property (i.e. L|g n = 0), and a finite m is therefore expected (as found explicitly for S = 1 above, by the exact treatment of the equivalent spin Hamiltonian H MF ). Nevertheless, the important point is that m vanishes asymptotically for large S, and it is generally a very small number otherwise (m = 0.0011412 at S = 1).
Two-fold degeneracy structure of the spin-wave spectrum Fig. 5 of the main text shows that the six spin-wave energies organize into three degenerate pairs. The symmetry origin of this degeneracy can be seen by considering the effect of the operation Σ 1 discussed above. We have: does not depend on the configuration of η's altogether. The Hamiltonian for the terms along the string becomes H = K(S u 1 S u 2 + S v 2 S v 3 + S u 3 S u 4 + · · · ) − |K|(S w 1 S w 1 + S w 2 S w 2 + · · · ) (51)