Quantum spin liquid in the semiclassical regime

TL;DR

This paper demonstrates the existence of a $Z_2$ quantum spin liquid ground state in a semiclassical large-$S$ magnet modeled by the spin-$S$ Kitaev honeycomb model, revealing topological order and distinct energy gaps.

Contribution

It shows that quantum spin liquids can occur in the semiclassical, large-$S$ limit, extending the understanding of topological order beyond low-spin quantum regimes.

Findings

Existence of $Z_2$ quantum spin liquid in large-$S$ limit

Presence of topological and symmetry-related degeneracies

Two distinct energy gaps with different scaling behaviors

Abstract

Quantum spin liquids have been at the forefront of correlated electron research ever since their original proposal in 1973, and the realization that they belong to the broader class of intrinsic topological orders, along with the fractional quantum Hall states. According to received wisdom, quantum spin liquids can arise in frustrated magnets with low spin $S$, where strong quantum fluctuations act to destabilize conventional, magnetically ordered states. Here we present a magnet that has a $Z_2$ quantum spin liquid ground state already in the semiclassical, large-$S$ limit. The 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 described by 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. The results apply to both integer and half-integer spins.

Equations

$$\mathcal{H}=K\Big{(}\!\!\sum_{\langle ij\rangle\in\text{`x'}}\!\!\!\!S_{i}^{x}S_{j}^{x}+\!\!\!\!\sum_{\langle ij\rangle\in\text{`y'}}\!\!\!\!S_{i}^{y}S_{j}^{y}+\!\!\!\!\sum_{\langle ij\rangle\in\text{`z'}}\!\!\!\!S_{i}^{z}S_{j}^{z}\Big{)}\,,$$

$$\delta E_{\text{ani}}=-(|K|S/16){\sum }_i({\tilde{a}}_i^4+{\tilde{b}}_i^4+{\tilde{c}}_i^4),$$

$$\begin{array}[]{lll}\mathbf{w}_{1}\!=\!\kappa\eta_{1}{\bf x},&\mathbf{u}_{1}\!=\!-\kappa\eta_{1}{\bf z},&\mathbf{v}_{1}\!=\!{\bf y},\\ \mathbf{w}_{2}\!=\!\eta_{2}{\bf y},&\mathbf{u}_{2}\!=\!-\kappa\eta_{1}{\bf z},&\mathbf{v}_{2}\!=\!-\kappa\eta_{1}\eta_{2}{\bf x},\\ \mathbf{w}_{3}\!=\!\kappa\eta_{3}{\bf z},&\mathbf{u}_{3}\!=\!\eta_{1}\eta_{2}\eta_{3}{\bf y},&\mathbf{v}_{3}\!=\!-\kappa\eta_{1}\eta_{2}{\bf x},\\ \mathbf{w}_{4}\!=\!\eta_{4}{\bf x},&\mathbf{u}_{4}\!=\!\eta_{1}\eta_{2}\eta_{3}{\bf y},&\mathbf{v}_{4}\!=\!\eta_{1}\eta_{2}\eta_{3}\eta_{4}{\bf z},\\ \mathbf{w}_{5}\!=\!\kappa\eta_{5}{\bf y},&\mathbf{u}_{5}\!=\!-\kappa\eta_{1}\eta_{2}\eta_{3}\eta_{4}\eta_{5}{\bf x},&\mathbf{v}_{5}\!=\!\eta_{1}\eta_{2}\eta_{3}\eta_{4}{\bf z},\\ \mathbf{w}_{6}\!=\!\eta_{6}{\bf z},&\mathbf{u}_{6}\!=\!-\kappa\eta_{1}\eta_{2}\eta_{3}\eta_{4}\eta_{5}{\bf x},&\mathbf{v}_{6}\!=\!-\kappa B_{h_{\alpha}}{\bf y},\end{array}$$

$$B_{h_{\alpha }}={\eta }_1{\eta }_2{\eta }_3{\eta }_4{\eta }_5{\eta }_6,$$

$$\displaystyle\begin{array}[]{ll}S_{1}^{z}S_{2}^{z}\mapsto S_{1}^{u}S_{2}^{u},{}&S_{2}^{x}S_{3}^{x}\mapsto S_{2}^{v}S_{3}^{v},\\ S_{3}^{y}S_{4}^{y}\mapsto S_{3}^{u}S_{4}^{u},{}&S_{4}^{z}S_{5}^{z}\mapsto S_{4}^{v}S_{5}^{v},\\ S_{5}^{x}S_{6}^{x}\mapsto S_{5}^{u}S_{6}^{u},{}&S_{6}^{y}S_{1}^{y}\mapsto-\kappa B_{h_{\alpha}}S_{6}^{v}S_{1}^{v}\,.\end{array}$$

$$\delta E(h_{\alpha })=c+J_m{\eta }_1{\eta }_2{\eta }_3{\eta }_4{\eta }_5{\eta }_6=c+J_m{B}_{h_{\alpha }},$$

$$|h_{\beta }⟩=|\kappa {\eta }_8\mathbf{z}{⟩}_9|{\eta }_8\mathbf{z}{⟩}_8|\kappa {\eta }_2\mathbf{y}{⟩}_7|{\eta }_2\mathbf{y}{⟩}_2|\kappa {\eta }_1\mathbf{x}{⟩}_1|{\eta }_1\mathbf{x}{⟩}_{10}.$$

$$⟨{\bar{h}}_{\beta }|{\mathcal{H}}_{\text{eff}}|h_{\beta }⟩\equiv J_e=3\times 2^{5-18S}{S}^{5-6S}[(2S-1)!{]}^{3}K.$$

$${\mathcal{H}}_{\text{eff}}$$

$$|X_{1}\!=\!1,X_{2}\!=\!1\rangle\!=\!\mathcal{N}\!\prod\nolimits_{p}\big{(}1+{\bf B}_{p}\big{)}~{}|\text{FM}_{x}\rangle~{},$$

$${\mathbf{W}}_{\text{BSS}}(h_{\beta })=\exp [i\pi (S_9^x+S_8^y+S_7^z+S_2^x+S_1^y+S_{10}^z)].$$

$$⟨{\bar{h}}_{\beta }|{\mathbf{W}}_{\text{BSS}}(h_{\beta })|h_{\beta }⟩=(-\kappa {)}^{2S}.$$

$$W_{\text{BSS}}(h_{\beta })=-(-\kappa {)}^{2S+1}.$$

$${\mathcal{H}}_{\text{magn}}(\{B_p\})={\sum }_{i=1}^N{\omega }_i(\{B_p\})b_i^{†}{b}_i,$$

$$|t_d|/|K|=3S^2{2}^{-6S}/(1-2^{-12S}).$$

$$\displaystyle\begin{array}[]{c}|1\rangle=|\kappa\eta_{1}{\bf y}\rangle_{1}~{}|\eta_{2}{\bf x}\rangle_{2}~{}|\kappa\eta_{2}{\bf x}\rangle_{3}~{}|\eta_{4}{\bf z}\rangle_{4}~{}|\kappa\eta_{4}{\bf z}\rangle_{5}~{}|\eta_{1}{\bf y}\rangle_{6}\\ |2\rangle=|\kappa\tilde{\eta}_{1}{\bf z}\rangle_{1}~{}|\tilde{\eta}_{1}{\bf z}\rangle_{2}~{}|\kappa\tilde{\eta}_{3}{\bf y}\rangle_{3}~{}|\tilde{\eta}_{3}{\bf y}\rangle_{4}~{}|\kappa\tilde{\eta}_{5}{\bf x}\rangle_{5}~{}|\tilde{\eta}_{5}{\bf x}\rangle_{6}\end{array}\,.$$

$$|\Omega |=|⟨1|2⟩|=2^{-6S},$$

$$\displaystyle\begin{array}[]{c}E_{0}\equiv\langle 1|\mathcal{H}|1\rangle=\langle 2|\mathcal{H}|2\rangle=-3|K|S^{2}\\ v\equiv\langle 1|\mathcal{H}|2\rangle=-6|K|S^{2}\Omega\end{array}$$

$$t_d=\frac{v-E_0\Omega }{1-{\Omega }^2}=-\frac{3K{S}^2{2}^{-6S}}{1-2^{-12S}}\times \text{sgn}(\Omega ),V=-\Omega t_d.$$

$${\gamma }_{\mathbf{R}}\equiv \kappa B_{\mathbf{R}}.$$

$$\mathcal{H}$$

$$S_i^w=S-c_i^{+}{c}_i,S_i^{+}\simeq \sqrt{2S}(c_i-\frac{n_i}{4S}{c}_i),S_i^{-}\simeq \sqrt{2S}(c_i^{+}-c_i^{+}\frac{n_i}{4S})$$

$$S_i^u\simeq \frac{\sqrt{S}}{\sqrt{2}}(c_i+c_i^{+}-\frac{n_i}{4S}{c}_i-c_i^{+}\frac{n_i}{4S}),S_i^v\simeq -i\frac{\sqrt{S}}{\sqrt{2}}(c_i-c_i^{+}-\frac{n_i}{4S}{c}_i+c_i^{+}\frac{n_i}{4S})$$

$$\displaystyle S_{i}^{u}S_{j}^{u}\simeq\frac{S}{2}\Big{(}c_{i}c_{j}+c_{i}c_{j}^{+}+h.c.\Big{)}-\frac{1}{8}\Big{(}c_{i}n_{j}c_{j}+c_{i}^{+}n_{j}c_{j}+c_{j}n_{i}c_{i}+c_{j}^{+}n_{i}c_{i}+h.c.\Big{)}$$

$$\displaystyle S_{i}^{v}S_{j}^{v}\simeq-\frac{S}{2}\Big{(}c_{i}c_{j}-c_{i}c_{j}^{+}+h.c.\Big{)}+\frac{1}{8}\Big{(}c_{i}n_{j}c_{j}-c_{i}^{+}n_{j}c_{j}+c_{j}n_{i}c_{i}-c_{j}^{+}n_{i}c_{i}+h.c.\Big{)}$$

$$\boxed{p_i=⟨c_i^{+}{c}_i⟩},\boxed{q_i=⟨c_i{c}_i⟩},\boxed{m_{ij}=⟨c_i{c}_j^{+}⟩},\boxed{{\delta }_{ij}=⟨c_i{c}_j⟩}.$$

$$c_i{n}_j{c}_j\simeq m_{ij}{c}_j{c}_j+2{\delta }_{ij}{n}_j+2p_j{c}_i{c}_j+q_j{c}_i{c}_j^{+}-m_{ij}{q}_j-2{\delta }_{ij}{p}_j$$

$$c_i^{+}{n}_j{c}_j\simeq {\delta }_{ij}{c}_j{c}_j+2m_{ij}{n}_j+2p_j{c}_i^{+}{c}_j+q_j{c}_i^{+}{c}_j^{+}-{\delta }_{ij}{q}_j-2m_{ij}{p}_j$$

$$\displaystyle\boxed{S_{i}^{u}S_{j}^{u}\simeq\tau_{ij}+\Big{(}f_{ji}c_{j}c_{j}+f_{ij}c_{i}c_{i}+g_{ij}c_{i}c_{j}+g_{ij}c_{i}c_{j}^{+}+h.c.\Big{)}+4f_{ij}(n_{j}+n_{i})}$$

$$\boxed{f_{ij}=-\frac{m_{ij}+{\delta }_{ij}}{8},g_{ij}=\frac{S}{2}-\frac{2(p_i+p_j)+q_i+q_j}{8},{\tau }_{ij}=-2f_{ij}[q_i+q_j+2(p_i+p_j)]}$$

$$\displaystyle\boxed{S_{i}^{v}S_{j}^{v}\simeq\tau^{\prime}_{ij}+\Big{(}{f_{ij}^{\prime}}c_{j}c_{j}+{f_{ij}^{\prime}}c_{i}c_{i}+g^{\prime}_{ij}c_{i}c_{j}-g^{\prime}_{ij}c_{i}c_{j}^{+}+h.c.\Big{)}-4f^{\prime}_{ij}(n_{i}+n_{j})}$$

$$\boxed{f_{ij}^{\prime }=\frac{m_{ij}-{\delta }_{ij}}{8},g_{ij}^{\prime }=-\frac{S}{2}+\frac{2(p_i+p_j)-q_i-q_j}{8},{\tau }_{ij}^{\prime }=-2f_{ij}^{\prime }[q_i+q_j-2(p_i+p_j)]}$$

Full text

Quantum spin liquids (QSLs) describe systems that evade magnetic long-range order down to zero temperature, and manifest a number of remarkable phenomena, such as topological ground state degeneracies, emergent gauge fields, and fractional excitations with non-trivial statistics. Anderson (1973); Fazekas and Anderson (1974); Kalmeyer and Laughlin (1987); Wen et al. (1989); Moessner et al. (2001); Kitaev (2003, 2006); Levin and Wen (2005); Balents (2010); Savary and Balents (2017); Zhou et al. (2017) The rich phenomenology of QSLs derives from an intrinsic tendency to form massive quantum superpositions of local, ‘product-like’ wavefunctions. Notable examples are the resonating valence bond (RVB) state, Anderson (1973); Rokhsar and Kivelson (1988); Moessner et al. (2001); Misguich et al. (2002); Hao et al. (2014) the gapped QSL of the Toric code, Kitaev (2003) and the gapless QSL phase of the spin-1/2 Kitaev honeycomb model. Kitaev (2006)

Typically, such massive superpositions arise in frustrated magnets with low spin $S$, which ideally have an infinite number of competing states and a strong tunneling between them. Lacroix et al. (2011) Here we show that the spin-$S$ version of the celebrated Kitaev honeycomb model is a topological $Z_{2}$ QSL already in the semiclassical limit. Specifically, the leading semiclassical fluctuations give rise to an effective low-energy description in terms of a pseudospin-1/2 Toric code. Kitaev (2003) The ‘magnetic flux’ term of the Toric code arises from the zero-point energy of spin waves above the classical ground states, while the ‘electric charge’ term stems from the tunneling between different classical states. The ensuing $Z_{2}$ QSL lives on top of a honeycomb superlattice of ‘frozen’ spin dimers, Baskaran et al. (2008) which take only two possible configurations, instead of $(2S\!+\!1)^{2}$. These two states are the pseudospin-1/2 degrees of freedom of the Toric code. The frozen dimer pattern breaks translational symmetry, so the QSL possesses an extra degeneracy associated to symmetry breaking, besides the topological one.

The $Z_{2}$ gauge structure is not an emergent property of the low-energy sector of the problem, but descends from the gauge structure of the original spin-$S$ model, which was discovered in a seminal study by Baskaran, Sen and Shankar (BSS). Baskaran et al. (2008) As such, the gauge structure is not only present in the low-energy sector, but also in the single-particle, spin-wave excitation channel, which we analyze in detail beyond the quadratic level.

The large-$S$ description breaks down around $S\!\sim\!3/2$. For lower $S$, tunneling processes that shift the dimer positions become quickly relevant and compete with the ‘freezing’ energy scale $\delta E_{\text{f}}$. Including these processes leads to a picture of ‘decorated quantum dimers’, where both the dimer positions and the orientations of the two spins in each dimer are allowed to resonate. The ensuing picture for $S\!=\!1$ in terms of another type of spin liquid will be discussed.

**Results

Model & classical ground states.** The spin-$S$ Kitaev model on the honeycomb lattice is described by the Hamiltonian

[TABLE]

where ‘x’, ‘y’ and ‘z’ denote the three orientations of nearest neighbor (NN) bonds, see Fig. 1 (a), and $K$ is the coupling constant. Note that there is a four-sublattice duality transformation Rousochatzakis et al. (2015a) that maps the positive $K$ to the negative $K$ model, but we shall discuss the general case here for completeness. We shall also define $\kappa\!=\!-\text{sgn}(K)$.

The classical ground states of this model were first analyzed by BSS. Baskaran et al. (2008) There, the authors identified an infinite number of so-called ‘Cartesian’ states, which map to dimer coverings of the honeycomb lattice, modulo a factor of two for the orientation of the two spins per dimer. They further showed that the Cartesian states are connected to each other by continuous valleys of other ground states, leading to a huge ground state degeneracy. Soon after, Chandra, Ramola and Dhar Chandra et al. (2010) showed that the manifold actually consists of infinitely more ground states and possesses an emergent gauge structure that leads to power-law correlations.

The crucial aspect of the present study is the use of a convenient parametrization of the classical ground state manifold, which reveals the topological terms arising from quantum fluctuations in an explicit way. This parametrization is shown in Fig. 1 (a). We denote the two sublattices of the honeycomb by A and B. Next, we parametrize each spin as $\mathbf{S}_{i}\!=\!(a_{i},b_{i},c_{i})$ or $\kappa(a_{i},b_{i},c_{i})$ for $i\!\in\!\text{A}$ or B, respectively, and $a_{i}^{2}\!+\!b_{i}^{2}\!+\!c_{i}^{2}\!=\!S^{2}$. Then, for every pair of NN sites, ${\bf S}_{i}$ and ${\bf S}_{j}$, we can minimize their mutual interaction by requiring that $a_{i}\!=\!a_{j}$ or $b_{i}\!=\!b_{j}$ or $c_{i}\!=\!c_{j}$, if the two sites share, respectively, an ‘x’ or ‘y’ or ‘z’ type of bond. To see if the ensuing states are ground states we check that they saturate the lower bound of the energy per site, $E_{\text{min}}/N\!=\!-|K|S^{2}/2$. Baskaran et al. (2008) Indeed, the energy from the three bonds emanating from any site $i$ add up to $-|K|(a_{i}^{2}\!+\!b_{i}^{2}\!+\!c_{i}^{2})\!=\!-|K|S^{2}$. And since each bond is shared by two sites, these configurations saturate the lower bound and are therefore ground states. The Cartesian states of BSS arise by keeping only one component of $(a_{i},b_{i},c_{i})$ finite, and equal to $\eta_{i}S$, where $\eta_{i}\!=\!\pm 1$. Modulo these Ising-like variables, the Cartesian states map to dimer coverings of the lattice [Fig. 1 (b)]. There are $(1.381)^{N/2}$ coverings, Wu (2006); Baxter (1970); Kasteleyn (1963) and $(1.662)^{N}$ Cartesian states in total. Baskaran et al. (2008)

The semiclassical analysis leading to the Toric code proceeds in three steps. The first is to show that fluctuations select the Cartesian over the non-Cartesian states, which identifies the positions and spin orientations of the dimers as the relevant degrees of freedom. In the second step, which was carried out by BSS, Baskaran et al. (2008) fluctuations freeze the positions of the dimers to a given pattern, leaving their spin orientation as the only relevant degrees of freedom below the associated freezing energy scale $\delta E_{\text{f}}$. At this point, our parametrization reveals, in addition, a topological structure that was not noticed previously. The third step is to include quantum-mechanical tunneling between states with different orientations of the dimers.

Order-by-disorder I: Selection of Cartesian states. The first crucial ingredient of the effective description in terms of dimers is to show that fluctuations select the Cartesian over the non-Cartesian states. BSS made this hypothesis based on an analogy to a related 1D problem. Here we prove it by real space perturbation theory (RSPT). Lindgård (1988); Long (1989); Heinilä and Oja (1993); Chernyshev and Zhitomirsky (2014) We introduce local frames $({\bf u}_{i},{\bf v}_{i},{\bf w}_{i})$, with $\mathbf{w}_{i}$ along the classical directions, and write ${\bf S}_{i}\!=\!S_{i}^{w}{\bf w}_{i}\!+\!S_{i}^{u}{\bf u}_{i}\!+\!S_{i}^{v}{\bf v}_{i}$. Then we split $\mathcal{H}$ into a diagonal part $\mathcal{H}_{0}\!=\!h\sum_{i}(S\!-\!S_{i}^{w})$, describing fluctuations in the local field $h\!=\!KS$, and a perturbation $\mathcal{V}\!=\!\mathcal{H}\!-\!\mathcal{H}_{0}$, which couples fluctuations on different sites. The essential physics is captured by the leading, short-wavelength corrections from second-order RSPT. The three types of bonds, say $(1\text{-}10)$, $(1\text{-}6)$ and $(1\text{-}2)$ of Fig. (1), give $\delta E_{1,10}\!=\!\xi(1\!-\!\widetilde{a}_{1}^{2})^{2}$, $\delta E_{1,6}\!=\!\xi(1\!-\!\widetilde{b}_{1}^{2})^{2}$, $\delta E_{1,2}\!=\!\xi(1\!-\!\widetilde{c}_{1}^{2})^{2}$,where $\xi\!=\!-|K|S/8$ and $(\widetilde{a}_{i},\widetilde{b}_{i},\widetilde{c}_{i})\!=\!(a_{i},b_{i},c_{i})/S$. Using the spin length constraints and disregarding overall constants, gives the anisotropy term

[TABLE]

similar to the one found in Rousochatzakis and Perkins (2017); Jackeli and Avella (2015). This anisotropy selects the Cartesian states, confirming the hypothesis of BSS. Baskaran et al. (2008)

Order-by-disorder II: Dimer freezing. Next, we discuss the lifting of the degeneracy within the infinite sub-manifold of Cartesian states, starting with the corrections from spin waves. As shown by BSS, i) the linear spin wave Hamiltonian splits into non-interacting modes propagating along loops without dimers, and ii) the minimum zero-point energy arises by maximizing the number of the shortest such ‘empty’ loops, like the shaded hexagon of Fig. 1 (b). This gives the ‘star’ or ‘columnar’ dimer pattern of Fig. 2 (a), which is known from the context of the quantum dimer model and the frustrated Heisenberg model on the honeycomb lattice. Moessner et al. (2001); Albuquerque et al. (2011); Schlittler et al. (2015) In this pattern, the only dynamical degrees of freedom remaining are the Ising-like variables $\eta\!=\!\pm 1$, which specify the direction of the two spins shared by each given dimer.

The physics of the dimer freezing is actually more involved from what is predicted from the linear spin-wave theory, but let us postpone this discussion for later and focus on the spin states associated to the ‘star’ pattern. There are three ways to place this pattern in the lattice and each dimer has two configurations, so at first sight, the number of selected spin states is $3\!\times\!2^{N/2}$. BSS showed, however, that the minimum zero-point energy is associated with spin-wave modes that have antiperiodic boundary conditions (ABC) around the empty hexagons, which reduces the number of states to $3\!\times\!2^{N/3}$.

However, this is not the full story yet. It turns out that the boundary condition on the spin wave modes actually endows the selected manifold with a topological magnetic flux term (and, in particular, the above number of states has to be multiplied $2^{2g-1}$, where $g$ is the genus of the system). To see this, we repeat the spin wave analysis using our $\eta$-parametrization. We begin by rewriting $\mathcal{H}$ in the local frame. Let us take the empty hexagon $h_{\alpha}$ of Fig. 2 (a) and choose ${\bf u}_{i}$ and ${\bf v}_{i}$ in the following way (and similarly for every other empty hexagon):

[TABLE]

where the product of the six $\eta$-variables on empty hexagons,

[TABLE]

is the magnetic flux that plays a central role in the following. With the above choice of the local frames, the couplings between empty hexagons map to terms of the type $\kappa S_{i}^{w}S_{j}^{w}$. For example, $S_{1}^{x}S_{10}^{x}\!\mapsto\!\kappa S_{1}^{w}S_{10}^{w}$. On the other hand, the intra-hexagon terms map as follows

[TABLE]

Thus, in the rotated frame, the only dependence of the Hamiltonian on $\eta$’s is via the products $\{B_{h_{\alpha}}\}$ on the empty hexagons $\{h_{\alpha}\}$. And since the choice of the local frame does not alter the physics, it follows that classical states that belong to the ‘star’ pattern of Fig. 2 (a) and have the same $\{B_{h_{\alpha}}\}$ share the same semiclassical spin wave spectrum, at all orders in $1/S$. (We shall see below that this property reflects a local gauge symmetry of the model. Baskaran et al. (2008).) The same is true for the renormalization of the ground state energy and therefore the order-by-disorder effect. Let us show the latter explicitly and we shall return to the spin-wave modes further below.

We introduce the usual Holstein-Primakoff bosons $a_{i}$ via the transformation, Holstein and Primakoff (1940) $S_{i}^{\dagger}\!=\!S_{i}^{u}+iS_{i}^{v}\!=\!(2S\!-\!a_{i}^{\dagger}a_{i})^{1/2}a_{i}$ and $S_{i}^{w}\!=\!S\!-\!a_{i}^{\dagger}a_{i}$. To order $\mathcal{O}(S)$, empty hexagons decouple, leading to a quadratic, six-site boson problem, with two sublattices and periodic (PBC) or antiperiodic (ABC) boundary conditions, for $\kappa B_{h_{\alpha}}\!=\!-1$ or $1$, respectively. So the BSS result that ABC give the lowest zero-point energy amounts to imposing $\kappa B_{h_{\alpha}}\!=\!1$ for all empty hexagons $h_{\alpha}$. More explicitly, by combining the energies Baskaran et al. (2008) $\delta E_{\text{PBC}}$ and $\delta E_{\text{ABC}}$ for PBC and ABC, respectively, we get the contribution to the zero-point energy from $h_{\alpha}$,

[TABLE]

where $c\!=\!\frac{\delta E_{\text{PBC}}\!+\!\delta E_{\text{ABC}}}{2}$ and $J_{m}\!=\!\frac{\delta E_{\text{PBC}}\!-\!\delta E_{\text{ABC}}}{2}$. This shows that the corrections to the ground state energy depend explicitly on the fluxes $\{B_{h_{\alpha}}\}$, and that states with the same set of fluxes have the same zero-point energy.

The linear spin-wave theory of BSS Baskaran et al. (2008) gives $\delta E_{\text{PBC}}\!=\!2|K|S$ and $\delta E_{\text{ABC}}\!=\!\sqrt{3}|K|S$, and so $J_{m}\!=\!\frac{2\!-\!\sqrt{3}}{2}KS$. However, as shown in Fig. 3 and emphasized below, the linear theory overestimates $|J_{m}|$ strongly due to the presence of a large percentage (four out of six) of ‘spurious’ zero modes.

We are now ready to identify the first crucial ingredient of the Toric code description announced above. The $\eta_{i}$ variables live on the midpoints of the bonds of a honeycomb superlattice (Fig. 2 (a)), and Eq. (9) tells us that promoting these variables to Pauli matrices ${\bm{\eta}}_{i}^{z}$ leads to the magnetic flux term of the Toric code Kitaev (2003) on this superlattice.

Order-by-disorder III: Tunneling. The second ingredient of the Toric code, the electric charge term, stems from processes that flip the three $\eta$’s around a vertex of the superlattice. Let us take, e.g., the spin coherent state of the $h_{\beta}$ hexagon of Fig. 2 (a),

[TABLE]

The leading processes that transform this state to its time reversed state $|\bar{h}_{\beta}\rangle$, with $\eta_{1}$, $\eta_{2}$ and $\eta_{8}$ flipped, appear in $(6S)$-th order of RSPT, with $\mathcal{V}\!=\!K(S_{7}^{x}S_{8}^{x}\!+\!S_{9}^{y}S_{10}^{y}\!+\!S_{1}^{z}S_{2}^{z})$. The corresponding off-diagonal matrix element $J_{e}$ of the resulting effective Hamiltonian $\mathcal{H}_{\text{eff}}$ depends, unlike $J_{m}$, on the choice of the local axes $({\bf u}_{i},{\bf v}_{i})$. Here we fix $J_{e}$ to be a real number by choosing the local axes such that $\mathcal{V}\!\mapsto\!K(S_{7}^{u}S_{8}^{u}\!+\!S_{9}^{u}S_{10}^{u}\!+\!S_{1}^{u}S_{2}^{u})$. Following e.g., the steps of the Supplemental Material of Rousochatzakis and Perkins (2017), we get

[TABLE]

In the language of the ${\bm{\eta}}$ operators, this matrix element is represented by $J_{e}{\bm{\eta}}_{1}^{x}{\bm{\eta}}_{2}^{x}{\bm{\eta}}_{8}^{x}$, which involve the three $\eta$’s around the vertex that sits at the center of $h_{\beta}$ (see Fig. 2 (a)).

Toric code model. Collecting the potential energy (disregarding $c$) and the tunneling terms above gives [see Fig. 2 (b)]:

[TABLE]

where $v$ and $p$ label, respectively, the vertices and the plaquettes of the honeycomb superlattice. In terms of the original lattice, the former sit at the centers of non-empty hexagons of type $h_{\beta}$, while the latter enclose the empty hexagons of type $h_{\alpha}$. Essentially then, $v$ and $p$ label $h_{\beta}$ and $h_{\alpha}$, respectively.

The remarkable properties of the model (12) stem from the relations ${\bf A}_{v}^{2}\!=\!{\bf B}_{p}^{2}\!=\!1$ and the fact that $\{{\bf A}_{v},{\bf B}_{p},\mathcal{H}_{\text{eff}}\}$ is a set of mutually commuting operators. Kitaev (2003) This model is a $Z_{2}$ lattice gauge theory, Wegner (1971); Kogut (1979) with the local gauge transformations generated by ${\bf A}_{v}$. In the following, we discuss the most important properties Kitaev (2003); Levin and Wen (2005); Savary and Balents (2017) of the Toric code. Without loss of generality, we will consider the $K\!<\!0$ case, where both $J_{m}$ and $J_{e}$ are negative.

Topological sectors. On a torus, $\prod_{v}{\bf A}_{v}\!=\!\prod_{p}{\bf B}_{p}\!=\!1$ and so there are $N_{v}\!=\!2^{N/3-1}$ and $N_{p}\!=\!2^{N/6-1}$ independent choices of ${\bf A}_{v}$ and ${\bf B}_{p}$, respectively, leading to $2^{N/2-2}$ states. So the quantum numbers $\{A_{v},B_{p}\}$ do not exhaust all $2^{N/2}$ states of ${\bm{\eta}}$’s. The missing quantum numbers are provided by the nonlocal operators ${\bf X}_{1}\!=\!\prod_{\text{C}_{\text{X1}}}\!\!{\bm{\eta}}^{x}$ and ${\bf X}_{2}\!=\!\prod_{\text{C}_{\text{X2}}}\!\!{\bm{\eta}}^{x}$, defined on the non-contractible loops C${}_{\text{X1}}$ and C${}_{\text{X2}}$ of Fig. 2 (b). These operators commute with ${\bf A}_{v}$ and ${\bf B}_{p}$, and with each other, and in addition ${\bf X}_{1}^{2}\!=\!{\bf X}_{2}^{2}\!=\!1$. The quantum numbers $\{A_{v},B_{p},X_{1},X_{2}\}$ then exhaust the Hilbert space of ${\bm{\eta}}$’s.

Ground states. The ground states have $A_{v}\!=\!B_{p}\!=\!1$, $\forall v,p$. On a torus, there are four such states, which correspond to the choices of the winding numbers $X_{1}$ and $X_{2}$. One of them is

[TABLE]

where $\mathcal{N}$ is a normalization factor, and $|\text{FM}_{x}\rangle\!=\!|\!\!\rightarrow\cdots\rightarrow\rangle$ is the fully polarized state along ${\bf x}$, which has $A_{v}\!=\!1$, $\forall v$. Expanding the product over $(1\!+\!{\bf B}_{p})$ shows that this state is the equal amplitude superposition of all possible loops of overturned spins (spins pointing along $-{\bf x}$, which correspond to electric flux lines) on top of the FM background, see Fig. 4 and Kitaev (2006); Levin and Wen (2005). The remaining three ground states of the Toric code, $|X_{1},X_{2}\rangle\!=\!|\text{-}1,1\rangle$, $|1,\text{-}1\rangle$ and $|\text{-}1,\text{-}1\rangle$, arise by replacing the reference state $|\text{FM}_{x}\rangle$ in (13) with ${\bf Z}_{2}|\text{FM}_{x}\rangle$, ${\bf Z}_{1}|\text{FM}_{x}\rangle$ and ${\bf Z}_{1}{\bf Z}_{2}|\text{FM}_{x}\rangle$, respectively, where ${\bf Z}_{1}\!=\!\prod_{C_{Z_{1}}}\!\!{\bm{\eta}}^{z}$ and ${\bf Z}_{2}\!=\!\prod_{C_{Z_{2}}}\!\!{\bm{\eta}}^{z}$, defined along $C_{Z_{1}}$ and $C_{Z_{2}}$ of Fig. 2 (b). These operators flip $X_{2}$ and $X_{1}$, respectively, because of the anti-commutation relations $\{{\bf Z}_{1},{\bf X}_{2}\}\!=\!0$ and $\{{\bf Z}_{2},{\bf X}_{1}\}\!=\!0$.

Importantly, the ground state sector of the original Kitaev spin model is 12-fold and not 4-fold degenerate, because there are three ways to place the dimer pattern of Fig. 2 (a) into the lattice and each sector has its own Toric code description.

Excitations of Toric code (12). The elementary excitations are pairs of static charges (vertices with $A_{v}\!=\!-1$), or pairs of static fluxes (plaquettes with $B_{p}\!=\!-1$). Their energy is $\Delta_{e}\!=\!4|J_{e}|$ and $\Delta_{m}\!=\!4|J_{m}|$, respectively. So $\Delta_{m}$ scales roughly linearly with $S$ (see Fig. 3), whereas $\Delta_{e}$ is exponentially small in $S$, as follows from Eq. (11), and practically vanishes for $S\!\geq\!1$ and realistic values of $K$. These excitations describe deconfined particles (the energies do not depend on the distance between the charges or fluxes) and they also carry nontrivial mutual statistics. Kitaev (2003)

Origin of gauge structure & BSS fluxes. The local $Z_{2}$ gauge symmetry of (12) is not an emergent property, but descends from the $Z_{2}$ gauge structure of the original spin-$S$ model, discovered by BSS. Baskaran et al. (2008) This structure stems from the presence of local conserved operators defined on the hexagons of the original lattice, which are called BSS fluxes in the following. For the $h_{\beta}$ hexagon of Fig. 2 (a), the BSS flux operator reads:

[TABLE]

Now, the BSS fluxes on non-empty hexagons have the same effect as the ${\bf A}_{v}$ operators, e.g.​ ${\bf W}_{\text{BSS}}(h_{\beta})|h_{\beta}\rangle\!\to\!|\bar{h}_{\beta}\rangle$ (modulo some prefactor, see (15) below). So the local gauge symmetry of (12) indeed descends from that of the full model.

Let us now examine the ground state BSS flux pattern. Unlike the original classical states associated with the ‘star’ pattern, where only the empty hexagons have well defined ${\bf W}_{\text{BSS}}$, Baskaran et al. (2008) the QSL ground states of (12) have well defined ${\bf W}_{\text{BSS}}$ on all hexagons. Indeed, using the same choice of local axes as the ones used above for the tunneling we find:

[TABLE]

Now, the resonating QSL state $|1,1\rangle$ of Eq. (13) satisfies $J_{e}A_{v}|1,1\rangle\!=\!-|1,1\rangle$, and therefore contains the combination $\frac{1}{\sqrt{2}}\left(|h_{\beta}\rangle\!-\!\text{sgn}(J_{e})|\bar{h}_{\beta}\rangle\right)$. So the ground state expectation value $W_{\text{BSS}}(h_{\beta})$ of the operator ${\bf W}_{\text{BSS}}(h_{\beta})$ is equal to

[TABLE]

For half-integer $S$, in particular, $W_{\text{BSS}}(h_{\beta})\!=\!-1$, irrespective of $\kappa$. For the empty hexagons, such as $h_{\alpha}$, a well-defined flux is already fixed by the zero-point energy, as shown by BSS. Baskaran et al. (2008) Specifically, $W_{\text{BSS}}(h_{\alpha})\!=\!(-1)^{\lambda S}$, where $\lambda\!=\!\kappa(\eta_{1}\!+\!\eta_{3}\!+\!\eta_{5})\!+\!\eta_{2}\!+\!\eta_{4}\!+\!\eta_{6}$, which is even. So, for integer $S$, $W_{\text{BSS}}(h_{\alpha})$ is always equal to $1$, while for half-integer $S$, $W_{\text{BSS}}(h_{\alpha})\!=\!-\kappa B_{h_{\alpha}}\!=\!-1$, because of the ABC condition on spin waves.

The BSS fluxes are in fact well defined in all eigenstates of (12), not just in the ground states. An excited state with an electric charge sitting on $h_{\beta}$ has $W_{\text{BSS}}(h_{\beta})\!=\!(-\kappa)^{2S+1}$, opposite to the one in the ground state. On the other hand, an excited state with a magnetic charge on $h_{\alpha}$ has again $W_{\text{BSS}}(h_{\alpha})\!=\!1$ for integer $S$ as in the ground state, but $+1$ for half-integer $S$. These results also mean that i) magnetic fluxes are related to BSS fluxes on empty hexagons for all $S$, and ii) electric charges are related to BSS fluxes on non-empty hexagons for half-integer $S$.

More generally, the fact that the BSS fluxes are well defined on all hexagons is consistent with Elitzur’s theorem Elitzur (1975); Fradkin (2013); Batista and Nussinov (2005) that local gauge symmetries cannot be broken spontaneously. Following the works of Baskaran et al. (2007, 2008), this also necessitates that static and dynamic two-spin correlation functions are identically zero beyond NN separation, consistent with the Toric code description.

Spin wave modes. In the frozen dimer pattern of Fig. 2 (a), the local Hilbert space for each spin-$S$ dimer has dimension $(2S\!+\!1)^{2}$, and Eq. (12) describes the dynamics inside the subspace of $|m_{1},m_{2}\rangle\!=\!|S,\kappa S\rangle$ and $|\text{-}S,\text{-}\kappa S\rangle$, where the projections $m_{1}$ and $m_{2}$ are defined along the local quantization axes. To this Hamiltonian (12), we should also add the terms that describe the coherent spin-wave bosonic modes

[TABLE]

describing the elementary, single-particle excursions outside this $2\!\times\!2$ manifold, with $\Delta m\!=\!\pm 1$. Note that the important constants arising from the spin wave theory have been assigned to $J_{m}$ already, and that the $b_{i}$ bosons are the eigenmodes of the spin-wave Hamiltonian, either at the quadratic or the self-consistent quartic order [see Supplementing material, Eq. (A22)]. Also, as mentioned above, the spin-wave frequencies $\omega_{i}$ depend on the set $\{B_{p}\}$ only, and are therefore the same for all states with the same $\{B_{p}\}$ but different $\{A_{v}\}$. This entails a huge, $2^{\frac{N}{3}+2g-1}$-fold degeneracy in the spin-wave branches, for each given set of $\{B_{p}\}$. We emphasize that the magnons discussed here do not describe the elementary excitations above some magnetically ordered state. Instead, they describe coherent excitations that are present in the spectrum independently of the elementary flux and charge excitations.

We now examine the actual structure of the magnon spectrum. At the quadratic level, BSS have shown Baskaran et al. (2008) that the spectrum consists of six flat bands, with $\omega_{1\text{-}4}({\bf k})\!=\!0$ and $\omega_{5,6}({\bf k})\!=\!\sqrt{3}|K|S$, where the momentum ${\bf k}$ belongs to the magnetic Brillouin zone. However, the problem with the quadratic theory is that the modes $1$-$4$ are not true zero modes, i.e. they will be gapped out by interactions. Such spurious zero modes are typical Chandra and Doucot (1988); Harris et al. (1992); Chubukov (1992); Khaliullin (2001); Dorier et al. (2005); Mulder et al. (2010); Chalker (2011); Rousochatzakis et al. (2015b, a) artifacts of the harmonic theory and reflect the modes that connect different classical minima. As commented above, the large number of such spurious zero modes in the present model leads to unreliable estimates for the relevant energy scales of the problem. This necessitates that we push the semiclassical expansion to quartic order, and treat the problem via a standard self-consistent decoupling scheme (see Supplementing material).

A key finding of this analysis is that spin waves remain localized inside the empty hexagons even at the interacting spin wave level, because of the local conservation laws associated with the BSS fluxes. In the language of the Holstein-Primakoff bosons, $a_{i}$, the conservation of BSS fluxes on empty hexagons (which remain well defined in the classical states of the ‘star’ pattern) translates into the conservation of the parity of the total number of magnons inside the empty hexagons (see Supplementing Material). As a result, individual hopping of magnons from one empty hexagon to the next is forbidden by symmetry. Pair hopping also does not occur because, as discussed above, the inter-hexagon couplings take the form $S_{i}^{w}S_{j}^{w}$, which gives rise to a term of the type $a_{i}^{\dagger}a_{i}a_{j}^{\dagger}a_{j}$, that leaves no room for pair hopping upon decoupling. Altogether then, the $2^{\frac{N}{3}+2g-1}$-fold degenerate branches corresponding to a given flux sector $\{B_{p}\}$ are flat in momentum space.

Fig. 5 (a) shows the magnon frequencies for the ground state flux sector, where all $B_{p}\!=\!1$, along with the corresponding results from the quadratic theory. All spurious modes are gapped out, and the spectrum organizes into three degenerate pairs due to symmetry (see Supplementing material). This figure also tells us that all modes sit far above the energy scales $|J_{m}|$ and $|J_{e}|$ of Eq. (12). In addition, the spin length corrections $\delta S$ of Fig. 5 (b) shows that spin waves do not reduce the spin length appreciably (at maximum it is about 15% for $S\!=\!1/2$), so the ${\bm{\eta}}$ variables are well defined objects.

Physics at low $S$. We now turn our discussion to what can go wrong with the above semiclassical picture as we lower $S$. The dimer freezing in the star pattern of Fig. 2 (a) stems from the zero-point energy of spin waves. However, this analysis disregards the quantum tunneling between different dimer patterns. The leading process is the one around a hexagon, see Fig. 6. The states associated with different dimer patterns are not orthonormal, but we can estimate the relevant tunneling amplitude $t_{d}$ using the truncation method of Rousochatzakis et al. (2014) (see Methods):

[TABLE]

At large $S$, $t_{d}$ is extremely small, and the spin-wave analysis of the dimer freezing has solid ground. This would in fact remain true down to $S\!=\!1$, if we were to use linear spin wave theory. However, this theory overestimates strongly the freezing energy scale (like $|J_{m}|$) due to the spurious zero modes mentioned above. As a result, $t_{d}$ becomes relevant below $S\!\sim\!3/2$. To see this, let us take as a representative freezing energy scale, the energy difference $\delta E_{\text{f}}^{(6,\infty)}$ between the star pattern and the ‘staggered’ pattern of Fig. 7, where the empty loops have infinite length. At the level of interacting spin wave theory, these energies are shown in Fig. 8 along with $|t_{d}|$ (where we divide by $N$ and by $6$, respectively, so that we compare energies per site). The results show clearly that dimers become mobile below $S\!\sim\!3/2$. (By contrast, linear spin-wave theory gives $\delta E_{\text{f}}^{(6,\infty)}\!/(NK)\!=\!(\frac{\sqrt{3}}{6}\!-\!\frac{1}{\pi})S$, Baskaran et al. (2008) which is much larger than $|t_{d}|/6$ down to $S\!=\!1$.)

It follows that in order to understand the physics of the $S\!=\!3/2$ and $S\!=\!1$ cases, we need to return to the Cartesian basis, and allow both the position of the dimers and their spin orientation to resonate. Such a ‘decorated quantum dimer’ description may appear quite more involved, but it may actually not be the case for the particular $S\!=\!1$ case. The reason is that $t_{d}/6$ is more than ten times larger than $\delta E_{\text{f}}/N$ for $S\!=\!1$ (see Fig. 8) and, from the standard quantum dimer model on the honeycomb lattice, Moessner et al. (2001); Schlittler et al. (2015) we know that $t_{d}$ stabilizes a resonating ‘plaquette’ dimer pattern, known also from the context of the frustrated Heisenberg model Albuquerque et al. (2011); Ganesh et al. (2013); Zhu et al. (2013). Including the much smaller $J_{e}$ term will include the resonances with the dimers of the opposite spin orientations. It would be interesting to check numerically this generalized semiclassical picture for $S\!=\!1$, and moreover whether certain features of this picture carry over to the exactly solvable $S\!=\!1/2$ case.

Discussion. It is shown that the low-energy sector of the large-$S$ Kitaev honeycomb model is described by a Toric code on a honeycomb superlattice. This should be contrasted with the effective square-lattice Toric code that arises in the spin-1/2 model when one of the three types of bonds has much stronger coupling than the other two. Kitaev (2006) Here, the magnetic and electric flux terms of the effective description arise respectively from the zero-point energy of spin waves and quantum-mechanical tunneling between different orientations of frozen dimers. This picture breaks down for $S\!\lesssim\!3/2$ where tunneling between different dimer patterns becomes relevant.

The prospects for realizing $S\!>\!1/2$ Kitaev magnets remain at present limited, although there are reports for nearly perfect honeycomb magnets with Co2+ ions, such as Na2Co2TeO6 and Na3Co2SbO6, Viciu et al. (2007) with peculiar spatial magnetic correlations. Lefrançois et al. (2016) These systems show single-ion anisotropy, but it is worth checking via ab initio methods if a strong Kitaev term is also present, as in the layered spin-1/2 iridates and ruthenates. Jackeli and Khaliullin (2009); Chaloupka et al. (2010); Witczak-Krempa et al. (2014); Trebst (2017) In parallel, there are proposals for emulating the model with trapped ions Schmied et al. (2011), superconducting quantum circuits, You et al. (2010) coupled cavity arrays, Xiang et al. (2012) and ultracold atoms in optical lattices, Duan et al. (2003); Micheli et al. (2006); Gorshkov et al. (2013); Manmana et al. (2013) which in particular offer the possibility for $S\!>\!1/2$ extensions of the model. Micheli et al. (2006); Gorshkov et al. (2013); Manmana et al. (2013)

Finally, we point out that the uniform Misguich et al. (2002) or staggered Wan and Tchernyshyov (2013); Hwang et al. (2015) charge sectors of Eq. (12) describe another well known $Z_{2}$ spin liquid, the RVB state of the spin-1/2 Heisenberg kagome antiferromagnet. Sachdev (1992); Yan et al. (2011); Depenbrock et al. (2012); Jiang, Hong-Chen and Wang, Zhenghan and Balents, Leon (2012); Rousochatzakis et al. (2014); Ralko et al. This highlights the universal topological features of QSLs arising from very different settings, across both isotropic and highly anisotropic magnets.

Acknowledgements. We thank G. Baskaran, A. Ralko, A. Tsirlin, Y. Wan and M. D. Schulz for fruitful discussions. Part of this work was done at the Perimeter Institute in Waterloo, which is supported by the Government of Canada through Industry Canada and by the Province of Ontario through the Ministry of Economic Development and Innovation. We also acknowledge the support from NSF Grant No. DMR-1511768.

**Methods

Derivation of Eq. (18).** To calculate the tunneling $t_{d}$ around a single hexagon we consider the simplest $2\!\times\!2$ truncation approach described in Rousochatzakis et al. (2014) (see also Schwandt et al. (2010)). Namely, we take a hexagon cluster and project the Hamiltonian into the $2\!\times\!2$ basis of dimer states shown in Fig. 6:

[TABLE]

The magnitude of the overlap $\Omega$ between the two states is

[TABLE]

and the matrix elements of the cluster Hamiltonian are

[TABLE]

Orthonormalizing the basis leads to the effective $2\!\times\!2$ Hamiltonian \left(\!\!\!\begin{array}[]{cc}E_{0}\!+\!v\!\!&\!\!t_{d}\\ t_{d}\!\!&\!\!E_{0}\!+\!v\end{array}\!\!\!\right), where the tunneling amplitude $t_{d}$ and the potential energy $V$ are given by Rousochatzakis et al. (2014); Schwandt et al. (2010)

[TABLE]

The latter is much smaller than $|t_{d}|$ and can be ignored.

**Author contributions

**All authors contributed to the analysis and interpretation of the results, and the preparation of the manuscript.

**Competing financial interests

**The authors declare no competing financial interests.

Appendix A Semiclassical expansion around the states associated with the star dimer pattern

A.1 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. To this end, we shall use the six-sublattice decomposition of Fig. 9, with a superlattice defined by the primitive translation vectors ${\bf T}_{1}$ and ${\bf T}_{2}$. Any given site $i$ of the lattice can be labeled as $i=(\mathbf{R},\nu)$, where ${\bf R}$ is a primitive vector of the superlattice and $\nu=1$-$6$ is the sublattice index. In this parametrization, the positions of the empty hexagons $h_{\alpha}$ are labeled by ${\bf R}$. The classical state is parametrized in terms of the $\eta$-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_{{\bf R}}$

[TABLE]

With these conventions and definitions, the Hamiltonian reads

[TABLE]

A.2 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 Holstein and Primakoff (1940) Holstein-Primakoff expansion for each site $i=({\bf R},\nu)$:

[TABLE]

where $c_{i},c_{i}^{+}$ are bosonic operators. We have:

[TABLE]

Below we shall make use of the following mean-field parameters

[TABLE]

These parameters are all real numbers because when written in the local coordinate frames above, the Hamiltonian has real matrix elements, and in addition the states around which we expand are real. This also implies the relations $m_{ij}=m_{ji}$ and $\delta_{ij}=\delta_{ji}$. Next, we can decouple the quartic terms as follows:

[TABLE]

Let us now write down the resulting expressions for each type of interaction that appears in the Hamiltonian.

• •

Terms of the type $S_{i}^{u}S_{j}^{u}$ (where $i$ and $j$ belong to the same empty hexagon):

[TABLE]

where

[TABLE]

• •

Terms of the type $S_{i}^{v}S_{j}^{v}$ (where $i$ and $j$ belong to the same empty hexagon):

[TABLE]

where

[TABLE]

• •

Terms of the type $S_{{\bf R},\nu}^{w}S_{{\bf R}+{\bf T}_{\nu\mu},\mu}^{w}$: The terms that couple different empty hexagons are of the form $S_{{\bf R},\nu}^{w}S_{{\bf R}+{\bf T}_{\nu\mu},\mu}^{w}$. For simplicity, we will label $({\bf R},\nu)\to i$ and $({\bf R}+{\bf T}_{\nu\mu},\mu)\to j$. We have:

[TABLE]

The quartic term decouples as follows:

[TABLE]

Now, the state around which we expand does not break the local BSS flux operators defined on the empty hexagons:

[TABLE]

where $n_{i}=c_{i}^{+}c_{i}$ is the boson number operator and $\lambda_{\bf R}=\kappa(\eta_{{\bf R},1}+\eta_{{\bf R},3}+\eta_{{\bf R},5})+(\eta_{{\bf R},2}+\eta_{{\bf R},4}+\eta_{{\bf 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 \kappa\big{(}\eta_{1}n_{\bf R,1}+\eta_{3}n_{\bf R,3}+\eta_{5}n_{\bf R,5}\big{)}+\big{(}\eta_{2}n_{\bf R,2}+\eta_{4}n_{\bf R,4}+\eta_{6}n_{\bf R,6}\big{)}. But since $\kappa$ and $\eta$ 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_{\bf R}$ of bosons in any given empty hexagon:

[TABLE]

This means that terms that change the parity of $N_{\bf 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 mean-field parameters $m_{ij}$ and $\delta_{ij}$ vanish by symmetry, and this is true to all orders in the Holstein-Primakoff expansion. We therefore get:

[TABLE]

i.e. empty hexagons decouple from each other and the $S_{i}^{w}S_{j}^{w}$ terms give, for each empty hexagon ${\bf R}$ alone, a contribution

[TABLE]

where the constant $p_{j}=p_{{\bf R}+{\bf T}_{\nu\mu},\mu}$ 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)^{\lambda_{\bf R}S}$ (see main text and Baskaran et al. (2008)). Spin wave fluctuations dress the reference state but cannot change the BSS fluxes, because these are integer numbers. Equation (38) then implies that the dressed ground state contains only terms with an even number of bosons $N_{\bf R}$.

A.3 Semiclassical expansion around the states of the star dimer pattern with uniform $\gamma_{\bf R}$

In the following we shall focus on the classical states that have the same $\gamma_{\bf R}$ on all empty hexagons. This means that $p_{{\bf R},\nu}$ is independent of ${\bf R}$, and we can therefore replace $p_{{\bf R}+{\bf T}_{\nu\mu},\mu}\to p_{\mu}$ in the above contribution from the $S_{i}^{w}S_{j}^{w}$ terms. Collecting all terms referring to a given hexagon and dropping the index ${\bf R}$ we get:

[TABLE]

where we have defined:

[TABLE]

Next we define $\mathbf{C}^{+}=\left(c_{1}^{+},\cdots c_{6}^{+},c_{1},\cdots c_{6}\right)$ and write:

[TABLE]

where the nonzero matrix elements of the matrix $\mathbf{M}$ are as follows:

[TABLE]

We next define the commutator matrix

[TABLE]

where $1\!\text{l}_{6}$ stands for the identity $6\!\times\!6$ matrix, and then perform a standard Blaizot and Ripka (1986) Bogoliubov transformation $\mathbf{C}\!=\!\mathbf{S}\!\cdot\!\mathbf{B}$, which must conserve the commutation relations $\tilde{\mathbf{g}}\!=\!\mathbf{g}$. This relation gives

[TABLE]

The matrix $\mathbf{S}$ must also satisfy the relation

[TABLE]

and at the same time diagonalize the Hamiltonian:

[TABLE]

where $\mathbf{\Omega}_{M}$ is diagonal and can be found from the eigenvalue equation $(\mathbf{g}\mathbf{M})\cdot\mathbf{S}=\mathbf{S}\cdot(\mathbf{g}\mathbf{\Omega}_{M})\equiv\mathbf{S}\cdot\mathbf{\Omega}_{gM}$. It can be shown Blaizot and Ripka (1986) that the eigenvalues of $\mathbf{g}\cdot\mathbf{M}$ come in pairs $(\omega_{\nu},-\omega_{\nu})$, where $\nu=1$-$6$. We finally get

[TABLE]

The ground state energy is, in particular, given by $E_{0}=f_{0}+\frac{1}{2}\sum_{\nu}\omega_{\nu}$.

A.4 Mean field parameters: General relations

Let us define the six eigenvectors of the matrix ${\bf g}\cdot{\bf M}$ that correspond to non-negative eigenvalues by ${\bf X}_{\nu}$, $\nu=1$-$6$. Using:

[TABLE]

we get the following expressions for the mean-field parameters:

[TABLE]

where ${\bf X}_{\nu}$ denotes the $\nu$-th eigenvector of ${\bf g}\cdot{\bf M}$. Note that the last expressions in each line do not depend on the arbitrary phase for the eigenvectors ${\bf X}_{\nu}$, which come out arbitrary when we diagonalize the matrix ${\bf g}\cdot{\bf M}$ numerically.

A.5 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 $\bm{\Sigma}_{1}$. This is a $\pi$-rotation in real space around the center of the hexagon, followed by $\pi/2$-rotations around the local $\bf{w}$-axes in spin space:

[TABLE]

These relations are equivalent with $\boxed{c_{1}\to-ic_{4},~{}~{}c_{2}\to-ic_{5},~{}~{}c_{3}\to-ic_{6},~{}~{}c_{4}\to+i\gamma c_{1},~{}~{}c_{5}\to+i\gamma c_{2},~{}~{}c_{6}\to+i\gamma c_{3}}$

• •

Symmetry $\bm{\Sigma}_{2}$. This is a reflection through the bonds (3,6) in real space, followed by $\pi/2$-rotations around the local $\bf{w}$-axes in spin space:

[TABLE]

These relations are equivalent with: $\boxed{c_{1}\to-i\gamma c_{5},~{}~{}c_{2}\to-i\gamma c_{4},~{}~{}c_{3}\to-i\gamma c_{3},~{}~{}c_{4}\to-i\gamma c_{2},~{}~{}c_{5}\to-i\gamma c_{1},~{}~{}c_{6}\to+ic_{6}}$

• •

Symmetry $\bm{\Sigma}_{3}$. This is a reflection through the middle of the bonds (1,2) and (4,5) in real space, followed by zero or $\pi$-rotations around the local-$\bf{w}$ axes in spin space:

[TABLE]

These relations are equivalent with: $\boxed{c_{1}\to-\gamma c_{2},~{}~{}c_{2}\to-\gamma c_{1},~{}~{}c_{3}\to c_{6},~{}~{}c_{4}\to c_{5},~{}~{}c_{5}\to c_{4},~{}~{}c_{6}\to c_{3}}$

• •

Symmetry $\bm{\Sigma}_{4}$. This is a reflection through the bonds (1,4) in real space, followed by $\pi/2$-rotations around the local $\bf{w}$-axes in spin space:

[TABLE]

These relations are equivalent with: $\boxed{c_{1}\to+i\gamma c_{1},~{}~{}c_{2}\to-ic_{6},~{}~{}c_{3}\to-ic_{5},~{}~{}c_{4}\to-ic_{4},~{}~{}c_{5}\to-ic_{3},~{}~{}c_{6}\to-ic_{2}}$

• •

Symmetry $\bm{\Sigma}_{5}$. This is a $\pi/6$-rotation in real space, followed by a $\pi/2$-rotation around the local ${\bf w}$-axes in spin space:

[TABLE]

These relations are equivalent with: $\boxed{c_{1}\to-ic_{2},~{}~{}c_{2}\to-ic_{3},~{}~{}c_{3}\to-ic_{4},~{}~{}c_{4}\to-ic_{5},~{}~{}c_{5}\to-ic_{6},~{}~{}c_{6}\to+i\gamma c_{1}}$

Combining $\bm{\Sigma}_{1}$-$\bm{\Sigma}_{5}$ gives the following constraints for the mean-field parameters:

[TABLE]

A.6 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:

[TABLE]

where $h_{\text{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 $\mathcal{L}$ defined as

[TABLE]

and the relations

[TABLE]

• •

For $S=1/2$, the numerical, self-consistent solution gives $h_{\text{loc}}=0.37888$ and $m=0$. However, this relation is special to $S=1/2$ because the numerical, self-consistent ground state $|g\rangle$ of $\mathcal{H}_{\text{MF}}$ has the special property $\mathcal{L}|g\rangle=0$. And according to the above relations, this implies that $\langle g|S_{1}^{+}S_{2}^{-}|g\rangle=0$, which is equivalent with $m=0$.

• •

For $S=1$ and higher, the ground state does not obey the property $\mathcal{L}|g\rangle=0$ and $m$ is therefore finite. The numerical solution for $S=1$ gives $h_{\text{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\rangle$ comes closer and closer to the classical vacuum $|0\rangle$ (with spins fully polarized along their local ${\bf w}$-axes), which has the property $\mathcal{L}|0\rangle=0$ (because $|0\rangle$ is an eigenstate of each $S_{\nu}^{w}$ individually). In fact, this relation remains true when we include the leading effect of semiclassical corrections coming from $\mathcal{V}$. At this leading level, the ground state wavefunction is given by Lindgren (1974)

[TABLE]

where $\mathcal{R}=\frac{1-|0\rangle\langle 0|}{E_{0}-\mathcal{H}_{0}}$ is the usual resolvent operator. To show that $\mathcal{L}|g_{1}\rangle=0$ we use the fact that $\mathcal{L}$ commutes with $\mathcal{H}_{0}$ (and therefore with $\mathcal{R}$ as well) and furthermore $\mathcal{L}|0\rangle=0$. These properties give:

[TABLE]

We further have:

[TABLE]

Using the standard spin commutation relations we find

[TABLE]

from which it follows that $[\mathcal{L},\mathcal{V}]|0\rangle=0$ and therefore $\mathcal{L}|g_{1}\rangle=0$.

At higher orders $n>1$, the ground state $|g_{n}\rangle$ does not satisfy this property (i.e. $\mathcal{L}|g_{n}\rangle\neq 0$), and a finite $m$ is therefore expected (as found explicitly for $S=1$ above, by the exact treatment of the equivalent spin Hamiltonian $\mathcal{H}_{\text{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$).

A.7 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 $\bm{\Sigma}_{1}$ discussed above. We have:

[TABLE]

Let us take the first row of this matrix equation:

[TABLE]

Suppose further that ${\bf X}_{1}=(a_{1},a_{2},a_{3},a_{4},a_{5},a_{6},a_{1}^{\prime},a_{2}^{\prime},a_{3}^{\prime},a_{4}^{\prime},a_{5}^{\prime},a_{6}^{\prime})$. Now, if $b_{1}$ describes an eigenmode, then $\bm{\Sigma}_{1}\cdot b_{1}$ is also an eigenmode with the same energy:

[TABLE]

This means that the eigenvectors corresponding to the positive (or the negative) eigenvalues of ${\bf g}\cdot{\bf M}$ come in pairs:

[TABLE]

If these two modes are linearly independent they belong to a 2-dimensional irreducible representation of the symmetry group generated by $\bm{\Sigma}_{1}$-$\bm{\Sigma}_{5}$. The numerical results show that this is the case for the whole spectrum of the spin-wave Hamiltonian.

Appendix B Semiclassical expansion around the states associated with the staggered dimer pattern

For the numerical data shown in Fig. 8 of the main text we have also performed a non-linear semiclassical expansion around the classical state associated with the staggered pattern of Fig. 7 of the main text. In this pattern, the dimers occupy the horizontal, ‘zz’ bonds, while the empty bonds form infinite strings. Similarly to the above, the strings decouple from each other and it suffices to consider one string only. To this end we use the labeling scheme of Fig. 10 and the following local frames:

[TABLE]

and so on. With this choice of local axes we move all the dependence on $\eta$’s on the last bond at infinity. And since the string is infinite, the energy contribution from that last bond does not matter, and therefore the spin wave expansion does not depend on the configuration of $\eta$’s altogether. The Hamiltonian for the terms along the string becomes

[TABLE]

The Hamiltonian along the string describe a system with a unit cell of two sites, and we can relabel the sites as follows:

[TABLE]

and so on. Keeping only the terms pertaining to the given string and going to momentum space (along the string) gives, in matrix notation:

[TABLE]

where

[TABLE]

and

[TABLE]

where the constants $f_{12}$, $f^{\prime}_{12}$, $g_{12}$, $g^{\prime}_{12}$, $\tau_{12}$ and $\tau^{\prime}_{12}$ are defined again as in Eqs. (33) and (35) above, and $N_{s}$ is the number of sites along the string. Here, the matrix $\mathbf{S}_{k}$ must satisfy:

[TABLE]

Note that the second relation replaces the relation $\bm{L}\cdot{\bf S}\cdot\bm{L}={\bf S}^{\ast}$ that we had in Eq. (60) above.

B.1 Symmetry constraints

• •

Symmetry $\bm{\Sigma}^{\prime}_{1}$. This is a translation by one lattice spacing, followed by a $\pi/2$-rotation around the local ${\bf w}$-axes:

[TABLE]

which is equivalent with $\boxed{c_{R,1}\to-ic_{R,2},~{}~{}~{}c_{R,2}\to-ic_{R+1,1}}$.

• •

Symmetry $\bm{\Sigma}^{\prime}_{2}$. This is a reflection though the bond (2,2’) (see Fig. 10), followed by a $\pi/2$-rotation around the local ${\bf w}$-axes:

[TABLE]

which is equivalent with $\boxed{c_{R,1}\to-ic_{R+1,1},~{}~{}~{}c_{R,2}\to-ic_{R,2}}$.

• •

Symmetry $\bm{\Sigma}^{\prime}_{3}$. This is an inversion through the middle of the bond (1,2) (see Fig. 10), which maps $\boxed{c_{R,1}\to c_{R,2}}$.

Combining the symmetries $\bm{\Sigma}^{\prime}_{1}$-$\bm{\Sigma}^{\prime}_{3}$ gives:

[TABLE]