Quantum Spin Liquid in the semiclassical regime

Quantum Spin Liquid in the semiclassical regime

Ioannis Rousochatzakis School of Physics and Astronomy, University of Minnesota, Minneapolis, MN 55455, USA    Yuriy Sizyuk School of Physics and Astronomy, University of Minnesota, Minneapolis, MN 55455, USA    Natalia B. Perkins School of Physics and Astronomy, University of Minnesota, Minneapolis, MN 55455, USA
July 5, 2019

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 , where strong quantum fluctuations act to destabilize conventional, magnetically ordered states. Here we present a magnet that has a quantum spin liquid ground state already in the semiclassical, large- 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 and an ‘electric charge’ gap that drops exponentially in . The magnet is described by the spin- 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.

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 , 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- version of the celebrated Kitaev honeycomb model is a topological 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 QSL lives on top of a honeycomb superlattice of ‘frozen’ spin dimers, Baskaran et al. (2008) which take only two possible configurations, instead of . 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 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- 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- description breaks down around . For lower , tunneling processes that shift the dimer positions become quickly relevant and compete with the ‘freezing’ energy scale . 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 in terms of another type of spin liquid will be discussed.

Figure 1: (a) Classical ground states of the Kitaev model. Here , or if belongs to the A or B sublattice, and . (b) The Cartesian states of BSS map to dimer coverings, with (yellow) dimers representing satisfied bonds. The spin orientation of each dimer is described by an Ising-like variable . The shaded hexagon has the shortest loop with no dimers.

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


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

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 or for or B, respectively, and . Then, for every pair of NN sites, and , we can minimize their mutual interaction by requiring that or or , 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, Baskaran et al. (2008) Indeed, the energy from the three bonds emanating from any site add up to . 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 finite, and equal to , where . Modulo these Ising-like variables, the Cartesian states map to dimer coverings of the lattice [Fig. 1 (b)]. There are coverings, Wu (2006); Baxter (1970); Kasteleyn (1963) and 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 . 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 , with along the classical directions, and write . Then we split into a diagonal part , describing fluctuations in the local field , and a perturbation , 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 , and of Fig. (1), give , , ,where and . Using the spin length constraints and disregarding overall constants, gives the anisotropy term


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 , 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 . 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 .

Figure 2: (a) ‘Star’ dimer (yellow) pattern selected from spin waves. Baskaran et al. (2008) The ’s describe the orientation of the two spins per dimer. They sit at the middle of the bonds of a honeycomb superlattice (dashed). (b) The resulting Toric code description of Eq. (9) on the honeycomb superlattice. The three- and six-body operators and of Eq. (9) are defined on vertices and plaquettes of the superlattice. For the torus geometry, C and C (similarly for C and C) are non-contractible loops that wrap the system in different directions.

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 , where is the genus of the system). To see this, we repeat the spin wave analysis using our -parametrization. We begin by rewriting in the local frame. Let us take the empty hexagon of Fig. 2 (a) and choose and in the following way (and similarly for every other empty hexagon):


where the product of the six -variables on empty hexagons,


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 . For example, . On the other hand, the intra-hexagon terms map as follows


Thus, in the rotated frame, the only dependence of the Hamiltonian on ’s is via the products on the empty hexagons . 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 share the same semiclassical spin wave spectrum, at all orders in . (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.

Figure 3: The coupling as a function of , extracted from linear (dashed) and non-linear (solid) spin wave theory.
Figure 4: The quantum spin liquid ground states of the Toric code (9) correspond to massive, equal-amplitude superpositions of all possible loops of spins (red solid lines) pointing along (red arrows), on top of a FM background of spins pointing along (blue arrows).

We introduce the usual Holstein-Primakoff bosons via the transformation, Holstein and Primakoff (1940) and . To order , empty hexagons decouple, leading to a quadratic, six-site boson problem, with two sublattices and periodic (PBC) or antiperiodic (ABC) boundary conditions, for or , respectively. So the BSS result that ABC give the lowest zero-point energy amounts to imposing for all empty hexagons . More explicitly, by combining the energies Baskaran et al. (2008) and for PBC and ABC, respectively, we get the contribution to the zero-point energy from ,


where and . This shows that the corrections to the ground state energy depend explicitly on the fluxes , 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 and , and so . However, as shown in Fig. 3 and emphasized below, the linear theory overestimates 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 variables live on the midpoints of the bonds of a honeycomb superlattice (Fig. 2 (a)), and Eq. (6) tells us that promoting these variables to Pauli matrices 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 ’s around a vertex of the superlattice. Let us take, e.g., the spin coherent state of the hexagon of Fig. 2 (a),


The leading processes that transform this state to its time reversed state , with , and flipped, appear in -th order of RSPT, with . The corresponding off-diagonal matrix element of the resulting effective Hamiltonian depends, unlike , on the choice of the local axes . Here we fix to be a real number by choosing the local axes such that . Following e.g., the steps of the Supplemental Material of Rousochatzakis and Perkins (2017), we get


In the language of the operators, this matrix element is represented by , which involve the three ’s around the vertex that sits at the center of (see Fig. 2 (a)).

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


where and 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 , while the latter enclose the empty hexagons of type . Essentially then, and label and , respectively.

The remarkable properties of the model (9) stem from the relations and the fact that is a set of mutually commuting operators. Kitaev (2003) This model is a lattice gauge theory, Wegner (1971); Kogut (1979) with the local gauge transformations generated by . 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 case, where both and are negative.

Topological sectors. On a torus, and so there are and independent choices of and , respectively, leading to states. So the quantum numbers do not exhaust all states of ’s. The missing quantum numbers are provided by the nonlocal operators and , defined on the non-contractible loops C and C of Fig. 2 (b). These operators commute with and , and with each other, and in addition . The quantum numbers then exhaust the Hilbert space of ’s.

Ground states. The ground states have , . On a torus, there are four such states, which correspond to the choices of the winding numbers and . One of them is


where is a normalization factor, and is the fully polarized state along , which has , . Expanding the product over shows that this state is the equal amplitude superposition of all possible loops of overturned spins (spins pointing along , 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, , and , arise by replacing the reference state in (10) with , and , respectively, where and , defined along and of Fig. 2 (b). These operators flip and , respectively, because of the anti-commutation relations and .

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 (9). The elementary excitations are pairs of static charges (vertices with ), or pairs of static fluxes (plaquettes with ). Their energy is and , respectively. So scales roughly linearly with (see Fig. 3), whereas is exponentially small in , as follows from Eq. (8), and practically vanishes for and realistic values of . 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 gauge symmetry of (9) is not an emergent property, but descends from the gauge structure of the original spin- 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 hexagon of Fig. 2 (a), the BSS flux operator reads:


Now, the BSS fluxes on non-empty hexagons have the same effect as the operators, e.g. (modulo some prefactor, see (12) below). So the local gauge symmetry of (9) 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 Baskaran et al. (2008) the QSL ground states of (9) have well defined on all hexagons. Indeed, using the same choice of local axes as the ones used above for the tunneling we find:


Now, the resonating QSL state of Eq. (10) satisfies , and therefore contains the combination . So the ground state expectation value of the operator is equal to


For half-integer , in particular, , irrespective of . For the empty hexagons, such as , a well-defined flux is already fixed by the zero-point energy, as shown by BSS. Baskaran et al. (2008) Specifically, , where , which is even. So, for integer , is always equal to , while for half-integer , , because of the ABC condition on spin waves.

The BSS fluxes are in fact well defined in all eigenstates of (9), not just in the ground states. An excited state with an electric charge sitting on has , opposite to the one in the ground state. On the other hand, an excited state with a magnetic charge on has again for integer as in the ground state, but for half-integer . These results also mean that i) magnetic fluxes are related to BSS fluxes on empty hexagons for all , and ii) electric charges are related to BSS fluxes on non-empty hexagons for half-integer .

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- dimer has dimension , and Eq. (9) describes the dynamics inside the subspace of and , where the projections and are defined along the local quantization axes. To this Hamiltonian (9), we should also add the terms that describe the coherent spin-wave bosonic modes


describing the elementary, single-particle excursions outside this manifold, with . Note that the important constants arising from the spin wave theory have been assigned to already, and that the 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 depend on the set only, and are therefore the same for all states with the same but different . This entails a huge, -fold degeneracy in the spin-wave branches, for each given set of . 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.

Figure 5: (a) Magnon frequencies (in units of ) from linear (dashed) Baskaran et al. (2008) and non-linear spin wave theory (solid), above the classical states associated with the dimer pattern of Fig. 2 (a). (b) Relative spin length correction from interacting spin waves.

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 and , where the momentum belongs to the magnetic Brillouin zone. However, the problem with the quadratic theory is that the modes - 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, , 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 , which gives rise to a term of the type , that leaves no room for pair hopping upon decoupling. Altogether then, the -fold degenerate branches corresponding to a given flux sector are flat in momentum space.

Fig. 5 (a) shows the magnon frequencies for the ground state flux sector, where all , 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 and of Eq. (9). In addition, the spin length corrections of Fig. 5 (b) shows that spin waves do not reduce the spin length appreciably (at maximum it is about 15% for ), so the variables are well defined objects.

Figure 6: Tunneling process that shifts the dimers around a hexagon.
Figure 7: Dimer pattern where empty bonds form infinite strings (red).

Physics at low . We now turn our discussion to what can go wrong with the above semiclassical picture as we lower . 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 using the truncation method of Rousochatzakis et al. (2014) (see Methods):


At large , is extremely small, and the spin-wave analysis of the dimer freezing has solid ground. This would in fact remain true down to , if we were to use linear spin wave theory. However, this theory overestimates strongly the freezing energy scale (like ) due to the spurious zero modes mentioned above. As a result, becomes relevant below . To see this, let us take as a representative freezing energy scale, the energy difference 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 (where we divide by and by , respectively, so that we compare energies per site). The results show clearly that dimers become mobile below . (By contrast, linear spin-wave theory gives Baskaran et al. (2008) which is much larger than down to .)

It follows that in order to understand the physics of the and 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 case. The reason is that is more than ten times larger than for (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 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 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 , and moreover whether certain features of this picture carry over to the exactly solvable case.

Discussion. It is shown that the low-energy sector of the large- 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 where tunneling between different dimer patterns becomes relevant.

The prospects for realizing Kitaev magnets remain at present limited, although there are reports for nearly perfect honeycomb magnets with Co ions, such as NaCoTeO and NaCoSbOViciu 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 extensions of the model. Micheli et al. (2006); Gorshkov et al. (2013); Manmana et al. (2013)

Figure 8: Competition between kinetic and potential (dimer freezing) energy of dimers.

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. (9) describe another well known 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.

Derivation of Eq. (15). To calculate the tunneling around a single hexagon we consider the simplest 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 basis of dimer states shown in Fig. 6:


The magnitude of the overlap between the two states is


and the matrix elements of the cluster Hamiltonian are


Orthonormalizing the basis leads to the effective Hamiltonian , where the tunneling amplitude and the potential energy are given by Rousochatzakis et al. (2014); Schwandt et al. (2010)


The latter is much smaller than 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.


Supplemental material

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 and . Any given site of the lattice can be labeled as , where is a primitive vector of the superlattice and - is the sublattice index. In this parametrization, the positions of the empty hexagons are labeled by . 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


With these conventions and definitions, the Hamiltonian reads

Figure 9: The ‘star’ dimer pattern of Fig. 2 of the main text and the six-sublattice decomposition used here.

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 :


where are bosonic operators. We have:

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


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 and . Next, we can decouple the quartic terms as follows:

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

  • Terms of the type (where and belong to the same empty hexagon):



  • Terms of the type (where and belong to the same empty hexagon):



  • Terms of the type : The terms that couple different empty hexagons are of the form . For simplicity, we will label and . We have: