Quantum Kagome Ice

Actively shought since the turn of the century, two-dimensional quantum spin liquids (QSLs) are exotic phases of matter where magnetic moments remain disordered even at extremely low temperatures. Despite ongoing searches, QSLs remain elusive, due to a lack of concrete knowledge of the microscopic mechanisms that inhibit magnetic order in real materials. Here, we study a theoretical model for a broad class of frustrated magnetic rare-earth pyrochlore materials called"quantum spin ices". When subject to an external magnetic field along the [111] crystallographic direction, the resulting spin interactions contain a mix of geometric frustration and quantum fluctuations in decoupled two-dimensional kagome planes. Using large-scale quantum Monte Carlo simulations, we identify a simple set of interactions sufficient to promote a groundstate with no magnetic long-range order, and a gap to excitations, consistent with a $Z_2$ spin liquid phase. This suggests a systematic experimental procedure to search for two-dimensional QSLs within the broader class of three-dimensional pyrochlore quantum spin ice materials.


INTRODUCTION
In a two-dimensional (2D) quantum spin liquid (QSL) state, strong quantum fluctuations prevent the ordering of magnetic spins, even at zero temperature. The resulting disordered phase can potentially be a remarkable state of matter, supporting a range of exotic quantum phenomena. Some, such as emergent gauge structures and fractional charges, are implicated in a wide range of future technologies like hightemperature superconductivity 1,2 and topological quantum computing. 3 It is therefore remarkable that, despite extensive examination of the basic theoretical ingredients required to promote a 2D QSL in microscopic models, 4,5 the state remains elusive, with only a few experimental candidates existing today. 6,7 Recently, the search for QSL states has turned to consider quantum fluctuations in the so-called spin ice compounds. 8 In these systems, magnetic ions reside on a pyrochlore lattice-a non-Bravais lattice consisting of corner-sharing tetrahedra. Classical magnetic moments (described by Ising spins) on the pyrochlore lattice can be geometrically frustrated at low temperatures, leading to spin configurations that obey the so-called "ice rules", a mapping to the proton-disorder problem in water ice 9 . The ice rules result in a large set of degenerate ground states -a classical spin liquid with a finite thermodynamic entropy per spin. 10,11 Two canonical materials, Ho 2 Ti 2 O 7 and Dy 2 Ti 2 O 7 , have been demonstrated to manifest spin ice behaviour, and experiments and theory enjoy a healthy dialog due to the existence of classical microscopic models capable of describing a wide range of experimental phenomena. 10 Classical spin ice pyrochlores are conjectured to lead to QSLs in the presence of the inevitable quantum fluctuations at low temperatures. 4,8 The effects of certain types of quantum fluctuations on the spin ice state have been investigated theoretically 12 and numerically, 13,14 where they have been demonstrated to lift the classical degeneracy and promote a three-dimensional (3D) QSL phase with low-energy gapless excitations that behave like photons. 12,13 In several related pyrochlore compounds, particularily Tb 2 Ti 2 O 7 , Yb 2 Ti 2 O 7 , Pr 2 Zr 2 O 7 , and Pr 2 Sn 2 O 7 , quantum effects have been observed, which make them natural candidates to search for such 3D QSLs. [15][16][17] In an attempt to elucidate the microscopic underpinnings of these and related materials, recent theoretical studies have produced a general lowenergy effective spin-1/2 model for magnetism in rare earth pyrochlores. [18][19][20] In an important development, Huang, Chen and Hermele 21 have shown that, on the pyrochlore lattice, strong spin-orbit coupling can lead to Kramers doublets with dipolar-octupolar character in dand f -electron systems. This allows for a specialization of the general effective model to one without the debilitating "sign problem" -amenable to solution through quantum Monte Carlo (QMC) methods -thus admitting a systematic search for QSL phases via large-scale computer simulations.

A QUANTUM KAGOME ICE MODEL
While the possibility for 3D QSLs in the above compounds is intriguing, spin ice materials offer a compelling mechanism for dimensional reduction to 2D, since singleion anisotropy constrains magnetic moments to point along the local tetrahedral symmetry axes in the pyrochlore lattice. This mechanism consists of the application of an external magnetic field along the global [111] crystallographic direction that partially lifts the spin ice degeneracy by "pinning" one spin per tetrahedron. As illustrated in Fig. 1a and triangular layers of the original pyrochlore structure. To simplest approximation, the system becomes a two-dimensional system of stacked kagome planes, [22][23][24][25] where spins on the intervening triangular planes align in the direction of the field (becoming energetically removed from the problem), while those in the kagome plane remain partially disordered. These kagome spins retain a fraction of the zero-field spin ice entropy, though still preserving the spin ice rules (two-in, two-out) of the parent pyrochlore system. This leads to classically disordered state, termed "kagome ice", 22,24-26 evidenced to date in several experimental studies on spin ice materials. [27][28][29] The above observations lead to a natural microscopic mechanism to search for 2D QSL behaviour. 30 First, one begins with classical nearest-neighbor spin ice in an applied [111] field, so as to promote the aforementioned "kagome ice" state. This model maps to a projected pseudo-spin Ising model with a symmetry-breaking Zeeman field h, arising from a combination of the physical [111] field and the original pyrochlore spin exchange interaction (Fig. 1c). For moderate h, the classical ground state retains an extensive degeneracy, before becoming a fully-polarized ferromagnetic state for h/J z > 2. Next, to include the effect of quantum fluctuations, one may add exchanges from the recent quantum spin ice models. [18][19][20][21] We consider only those quantum fluctuations discussed by Huang et. al. 21 , to obtain a pseudo-spin Hamiltonian on the kagome lattice, Here, S r are spin-1/2 operators, with a global z axis (S z r = 1/2 = and S z r = −1/2 = in Fig. 1c). This Hamiltonian cannot be solved exactly by analytical techniques; however large-scale QMC simulations are possible in a parameter regime devoid of the prohibitive "sign problem", which occurs for J ± > 0.
One can imagine a 2D QSL state arising conceptually by considering the quantum fluctuations J ± and J ±± as perturbations on the classical kagome ice limit, where only diagonal terms J z > 0 and h J z are present. Previously, large-scale QMC simulations have been performed on the kagome model in the limit J ± > 0 and J ±± = 0, 31,32 (a parameter regime where the Hamiltonian retains U(1) invariance). In that case, quantum fluctuations promote an in-plane ferromagnetic (FM) phase for h = 0, and a "valence bond-solid" (VBS -a conventional symmetry broken phase) for h > 0. Thus, it happens that fluctuations of the form induced by J ± are not sufficient to promote a 2D QSL state.
However, there remains the theoretical possibility of a gapped Z 2 QSL phase promoted by the J ±± quantum fluctuations. As detailed in the Supplemental Information, the local constraints of classical kagome ice can be translated into a charge-free condition on the dual honeycomb lattice. Then, the full Hamiltonian (1) can be re-cast as a system of interacting bosonic spinons coupled to a compact U(1) gauge field on the dual lattice. In the limit of J ± = 0, this theory is expected to exhibit two distinct phases. One is a "confined" phase, corresponding to a conventional spin-ordered state; the other is a "deconfined" Z 2 QSL phase. 20,21,33 From these simple arguments it is conceivable that these two phases exist in the phase diagram of Eq. (1). In the next section,  we set J ± = 0 and explore this possibility for all parameter regimes J ±± /J z and h/J z , using non-perturbative, unbiased QMC simulations.

QUANTUM MONTE CARLO RESULTS
We implement a finite-temperature Stochastic Series Expansion 34-36 (SSE) QMC algorithm with directed loop updates in a 2 + 1 dimensional simulation cell, designed specifically to study the Hamiltonian Eq. (1) with J ± = 0 (for details, see the Methods Summary). Note, this Hamiltonian explicitly breaks U(1) invariance, retaining global Z 2 symmetries. By a canonical transformation, S ± → ±iS ± ; we simulate only J ±± < 0, without loss of generality. 21 Various measurements are possible in this type of QMC simulation. Simplest are the standard SSE estimators for energy, magnetization m z = m = 1 V i S z i , and uniform spin susceptibility The latter two allow us to map out the broad features of the phase diagram. Further, we measure the off-diagonal spin structure factor 37 Here, r i points to the sites of the underlying triangular lattice (containing N s sites) of the kagome lattice (containing V = 3 × N s sites). The vectors α are the position of each site within the unit cell with respect to the vector r i . This quantity allows us to define, for this spin Hamiltonian, the analogue of a condensate fraction in bosonic systems, 38,39 which detects transverse magnetic ordering. We define f 0 = n M V as the ratio of largest eigenvalue n M of the "single-particle" density matrix ρ i,j = S + ri S − rj to the number of sites V . The eigenvalues of ρ i,j coincide with n q = α n α,α q for a translationally invariant system. Figure 2 shows the QMC phase diagram for the J ± = 0 model of Eq. (1), using data for the condensate fraction f 0 . Careful finite-temperature and finite-size scaling, performed up to lattice sizes of V = L × L × 3 = 39 × 39 × 3 and β = J z /T = 96, is detailed in the Supplemental Information. The magnetization curve and the uniform spin susceptibility across the phase boundary at fixed h/J z = 0.833 are presented in Fig. 2. The data clearly indicates the existence of two magnetized "lobes" on the phase diagram for J ±± /J z < 0.5 and h/J z = 0, where the zero-momentum condensate fraction of a surrounding FM phase is destroyed by a phase transition (which appears to be first order). The lobes have magnetizations of m ≈ −1/6 and m ≈ +1/6 for h/J z < 0 and h/J z > 0, respectively. The FM phase has a finite uniform susceptibility χ z , while the lobe phases retain a small but finite χ z that can be understood by the nature of the quantum fluctuation (S + r S + r + S − r S − r ) as a spin pair interaction, which does not conserve the total magnetization S z tot . As discussed above, the phase in these lobes is a candidate for supporting a 2D QSL state.
In order to examine this hypothesis, we perform a detailed search for ordered structures in the lobes. In related models, particularly the spin-1/2 XXZ model on kagome (i.e. J ±± = 0 and J ± > 0), 31,32 the analogous lobes support a conventional VBS phase, which is evident in the diagonal structure factor: S αβ If there there is long-range order then S q = α S αα q will scale with system size for at least one value of q.
We also measure the bond-bond structure factor using a four-point correlation function. where Nearest neighbor sites i aα and j aα belong to bond α in a unit cell located at position r a . Again, if there is pair long-range order then BB q = α BB αα q should scale with system size for at least one value of q, with which we define B q = BB q /V . Figure 3 illustrates the various q-dependent structure factors for spin and bond order. These structure factors display diffuse peaks at various wave vectors, notably q = 0, q = K = (2π/3, 0), and symmetry-related momenta. Such peaks would indicate the presence of longrange order, should they sharpen, and survive in intensity in the infinite-size limit, where S/V would correspond to an order parameter squared. In the insets to Fig. 3, we examine this through a standard finite-size scaling analysis, for several candidate peaks for each of the structure factors. Further scaling analysis, including larger system sizes, is presented in the Supplemental Information. In each case, the QMC data indicates a scaling of each order parameter to zero in the limit V → ∞. Note in particular, the largest value of B q corresponds to q = 0, which remains finite as V → ∞, meaning that the bond expectation values S + iaα S + jaα = 0 is finite in the lobes. This is expected as S + iaα S + jaα represents the "kinetic energy" of quantum fluctuations in the system, thus it should be finite in all phases. More importantly, the data indicates that in the limit of V → ∞ this quantity is the same on all bonds of the unit cell of the kagome lattice, meaning that there is no breaking of space-group symmetry (see Supplemental Information).
Finally, since the above data suggests the existence of a phase that is homogeneous, disordered, and quantummechanically fluctuating at extremely low temperatures, one should also examine whether the energy for excitations out of this ground state is gapped or gapless. Although a direct measurement of the gap is not possible in this type of SSE QMC method, we can indirectly probe its existence by looking at the decay of real-space correlations. In Fig. 4, we compare the decay of single- particle correlations between the m z = ±1/6 magnetization lobes, and the adjacent FM ordered phase. For the system size studied, it is clear that that correlations in the lobe are consistent with exponential decay, and therefore indicative of a gap. In contrast, in the FM phase the correlations quickly reach a finite value, indicating symmetry breaking. Similarly, the diagonal part of the spin correlation function is consistent with exponential decay both in the lobes and in the FM phase (not illustrated). Thus, our QMC results have elucidated a phase diagram for our kagome pseudo-spin XYZ model (with J ± = 0) which contains a predominant FM phase, surrounding "lobes" of an exotic disordered m z = ±1/6 magnetization phase. Since our dual gauge theory (detailed in the Supplemental Information) indicates that these lobes may realize a QSL with an emergent Z 2 gauge symmetry, it is clear that further simulation work should be carried out to address this hypothesis. To confirm the presence of a Z 2 QSL, one requires evidence of either excitations consistent with this gauge structure (e.g. magnetic spinons or non-magnetic visons at non-zero temperature), or a smoking gun such as the topological entanglement entropy. 40,41 Such evidence, although demonstrated in the past with SSE QMC, 42,43 is resource-intensive to obtain, requiring high numerical precision at very low temperatures, and thus outside of the scope of the present manuscript.
However, we also note that, due to the presence of only a discrete symmetry in our kagome XY Z model, an emergent Z 2 structure is not strictly required by the Lieb-Schultz-Mattis-Hastings (LSMH) theorem, [44][45][46] which states that a system with half-odd-integer spin in the unit cell cannot have a gap and a unique ground state. In higher-symmetry Hamiltonians, the requirements of the LSMH theorem are satisfied in a gapped QSL phase by the topological degeneracy, which is a consequence of the emergent discrete gauge symmetry. For our Hamiltonian with a gapped QSL arising in a model with only global discrete symmetries, an emergent gauge structure is not required. Rather, it is possible that the groundstate is a quantum paramagnet. On the other hand, other types of emergent gauge structure, topological order, or other exotic phenomena are theoretically possible. Fortunately, the nature of this Hamiltonian, which is among the first to show a 2D QSL phase with only nearest-neighbor interactions, lends itself exceedingly well to study by sign-problem-free QMC simulations. We therefore expect a large number of studies in the near future will help elucidate the precise nature of this QSL phase.

CONCLUSION AND OUTLOOK
Through extensive quantum Monte Carlo (QMC) simulations, we have studied a sign-problem-free model of frustrated quantum spins interacting on a twodimensional (2D) kagome lattice. This model is dece-dent from a more general quantum XYZ Hamiltonian discussed by Huang, Chen and Hermele, 21 derived for the three-dimensional pyrochlore lattice, when subject to a magnetic field along the [111] crystalographic direction. For a large range of Hamiltonian parameters, the QMC data uncovers an exotic disordered phase which breaks no symmetries, has strong quantum mechanical fluctuations and exponentially decaying correlations -a gapped quantum spin liquid (QSL) phase. This discovery is consistent with an analytical dual gauge theory (detailed in the Supplemental Information), which indicates that, in the limit of small quantum fluctuations, the phase could be a 2D QSL with an emergent Z 2 gauge symmetry.
Our work suggests a new experimental avenue to search for the highly-coveted QSL phase in two dimensions. Previous efforts have focussed largely on SU(2) Hamiltonains on kagome or triangular lattice materials. 6,7 In contrast, we propose to concentrate the search on the quantum spin ice pyrochlore materials, subject to an external field along the [111] direction. Such kagome ice phases have been identified in various materials in the past. A closer look at several quantum spin ice candidates is warranted, particularly in materials where strong quantum fluctuations are known to exist, such as Tb 2 Ti 2 O 7 , Yb 2 Ti 2 O 7 , Pr 2 Zr 2 O 7 , and Pr 2 Sn 2 O 7 . Also, in light of recent experiments 47 which suggest that the oft-studied classical spin ice state is only metastable in Dy 2 Ti 2 O 7 , it would seem prudent to re-examine the kagome ice state of this material using similar long-timescale techniques, to ascertain whether evidence of a QSL state may be present yet dynamically inhibited in the short-timescale studies performed to date.

METHODS SUMMARY
We developed a Stochastic Series Expansion 34 (SSE) QMC algorithm in the global S z basis designed to study the Hamiltonian Eq. (1) with J ± = 0 at finite temperature, using a 2 + 1-dimensional simulation cell. Within the SSE, the Hamiltonian was implemented with a triangular plaquette breakup, 36 which helps ergodicity in the regime where J z /J ±± is large. Using this Hamiltonian breakup, the standard SSE directed loop equations 35 were modified to include sampling of off-diagonal operators of the type S + r S + r + h.c. The resulting algorithm is highly efficient, scaling linearly in the number of lattice sites V and the inverse temperature β. This scaling is modified to V 2 β in the cases where a full q-dependent structure factor measurement is required.
The program was implemented in Fortran and verified by comparing results for small clusters with exact diagonalization data. For each set of parameters in Eq. (1), the simulation typically requires 10 7 QMC steps, with approximately 10% additional equilibration steps. The data presented in this paper required computational resources equivalent to 100 CPU core-years, run on a highperformance computing (HPC) cluster with Intel Xeon

ACKNOWLEDGMENTS
We would like to thank F. Becca, A. Burkov, L. Cincio, T. Senthil, and M. Stoudenmire for enlightening discussion, and M. Gingras for a critical reading of the manuscript. We are particularly indebted to Gang Chen for bringing the models discussed in Ref. 21  To understand the possible phases of the Hamiltonian (1), we concentrate on its classical limit first: At h = 0, the energy is minimized by any spin configurations with two spins up (down) one spin down (up) on every triangle. The number of such spin configurations increases exponentially with the system size. A finite magnetic field h > 0 partially lifts the extensive degeneracy: the energy of (5) is minimized by any spin configuration with two spin pointing up and one spin pointing down on each triangle 48 . The ground state degeneracy remains extensive. It has been well established that that the manifold of such states can be described by electric flux configurations with no charges present 49 .
To make the assertion explicit, we use x to label the centers of triangles forming a honeycomb lattice. r labels sites on the kagome lattice. A nearest neighbor bond xx connecting sites x and x on the honeycomb lattice goes through site r on the kagome lattice. We translate S z r to electric flux E x,x using the following relation (see Fig.1(c) and Fig.5): where η x = ±1 if site x belongs to A or B sublattice of the honeycomb lattice.μ (µ = 1, 2, 3) are the three vectors connecting site x with its nearest neighbors.1 = (0, 1), 2 = ( √ 3/2, −1/2) and3 = (− √ 3/2, −1/2). The total charge on site x satisfies the Gauss's law: Use these relations, we rewrite (5) as: By construction Q x is constrained to take only integer values. For h > 0, it is now explicit that the ground states of the classical model (5) correspond to divergence-free electric field configurations, or Q x = 0. We now turn our attention to transverse Hamiltonian terms. A transverse operator, S + r for example, creates a pair of positive and negative charges. We follow Savary and Balents 19 to introduce conjugate operator φ x : [φ x , Q x ] = iδ xx . S + r is mapped to: FIG. 5. Divergence-free gauge field configuration on the honeycomb lattice corresponding to the pseudo-spins configuration in Fig.1(c). Blue (orange) arrows indicate electric flux of 1/3 (2/3).
Here ψ ± x increases or decreases Q x by 1 and s + xµ = s + x+μ,µ . Using (8), and (9), we write (1) as: The Hamiltonian is invariant under a local U (1) gauge transformation: The azimuthal angle of pseudo spin s r is the vector gauge field 0 ≤ A xµ < 2π: the gauge theory is compact. Equation (10) is our main result. The low energy physics of the spin model (1) is thus described by bosonic matter interacting with a U (1) compact gauge field. The theory possesses two types of ground states. 33 The first class consists of confined phases where either matter fields are confined or charge 1 particles are condensed. In terms of spins, such phases correspond to long-range ordered phases where some symmetry of the Hamiltonian is broken. For the other category of ground states, charge n (n > 1) excitations are condensed. The resulting phase, the so-called charge-n Higgs phase, does not need to break any symmetry of the Hamiltonian, possesses local Z n gauge structure and intrinsic topological order. These are the Z n spin liquids. The simplest example is the Z 2 spin liquid with n = 2. In such a state, the spinon pairing introduces a mass to the gauge field, much like the condensation of cooper pairs in a superconductor generates a mass to the gauge field. This means that, at low energies, the local gauge redundancy reduces from U (1) to Z 2 through the Anderson-Higgs mechanism. Both the gauge field and the spinons should be gapped, which is corroborated in our simulations in the sense that both diagonal (gauge field) and off-diagonal (spinon) spin correlation functions are consistent with an exponential decay as a function of distance.
A simple mean-field decoupling of the spinon interaction term in Eq. (10) for J ± = 0 provides further insight into the properties of possible phases. We decompose the four-rotor term in both the hopping channel and the pairing channel 20,21 through the following mean-field parameters: A finite g indicates semi-classical long-range order since S + = s + ψ † x ψ x+ηxμ ∼ g; this is the FM phase obtained in our numerical simulations. On the other hand, finite n and f with g = 0 represents a state with condensed double charged matter fields. It is a gapped Z 2 spin liquid state. 20,21,33 We find that in the lobes S + r ∼ g vanishes while S + r S + r ∼ nf remains finite, in agreement with a Z 2 spin liquid phase.

Absence of breaking of rotational invariance
In Fig. 6 we show a finite-size scaling of both the order parameter f 0 , as well as the zero-momentum bond occupation B 0 , across the transition between the FM phase and the disordered phase in the lobes for T = J z /24. Outside the lobes both f 0 and B 0 remain finite as V → ∞, whereas in the lobes f 0 vanishes (see Fig. 6(c)) but B 0 remains finite (see Fig. 6(d)). Notice that because B 0 is finite, the state in the lobes can still break the rotational symmetry of the lattice by having different expectation values of B α ra over the 6 independent bonds per unit cell in the kagome lattice. This potential symmetry breaking can be quantified by analyzing the size scaling of the different sublattice occupations B αβ 0 = BB αβ 0 /V . In the thermodynamic limit, they are all expected to be equal if the system does not break rotational symmetry. In Fig. 7 we show the finite-size scaling of several matrix elements B αβ 0 in the FM phase and in the lobes ( Fig. 7(a) and Fig. 7(b), respectively). The finite-size scalings clearly suggest that the different B αβ 0 converge to the same value as V → ∞. Using a linear fit to the finite-sized data, we evaluate the extrapolated matrix elements lim V →∞ BB αβ 0 /V in the FM phase and in the lobes, as shown in Fig. 8(a) and (b). These are consistent with a picture where there is no breaking of rotational invariance in both the FM and in the lobe phases.   Finite-size scaling inside the lobe for larger system sizes and T = Jz/16.
In Fig. 9 we re-examine the finite-size scaling of the most relevant candidate peaks in the structure factors presented in Fig. 3  39 × 39, which happens to accomodate the two relevant reciprocal lattice vectors 0 and K. The continuous lines represent linear fits to the numerical data and the extrapolations to the limit of V → ∞ are stable within the error bars upon removal of the smaller system sizes. These results obtained on larger clusters support the conclusion the phases in the lobes correspond to magnetically disordered phases.
Spin structure factors in the lobe at T = Jz/96.
In figure 10 we present the structure factors inside the lobe for for a system with h/J z = 0.8333, J ±± /J z = 0.495, and T = J z /96, on a cluster with N s = 24 × 24. The data confirm that at the lowest temperature reached in our simulations, the conclusions about the existence of a magnetically disordered phase in the lobes still holds.
Details of the FM phase and absence of a phase transition at zero field In figure 11 we present results for spin correlation functions just outside the upper 1/3-magnetization lobe. Both n q and BB q exhibit diverging peaks at zero momentum (which have been removed in the figures), while the diagonal structure factor S q does not. The fact that f 0 remains finite in outside the lobe together with the The "condensate fraction" f0 and the uniform magnetic susceptibility χz, ((a) and (b),respectively), as a function of J±±/Jz and system size at zero magnetic field h and T = Jz/96.