Quantum Kagome Ice Juan Carrasquilla,1, ∗ Zhihao Hao,2 and Roger G. Melko1, 2 arXiv:1407.0037v1 [cond-mat.str-el] 30 Jun 2014 2 1 Perimeter Institute for Theoretical Physics, Waterloo, Ontario, N2L 2Y5, Canada Department of Physics and Astronomy, University of Waterloo, Ontario, N2L 3G1, Canada (Dated: July 2, 2014) 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 Z2 spin liquid phase. This suggests a systematic experimental procedure to search for twodimensional 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 superconductivity1,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 ice9 . 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, Ho2 Ti2 O7 and Dy2 Ti2 O7 , 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 theoretically12 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 Tb2 Ti2 O7 , Yb2 Ti2 O7 , Pr2 Zr2 O7 , and Pr2 Sn2 O7 , quantum effects have been observed, which make them natural candidates to search for such 3D QSLs.15–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–20 In an important development, Huang, Chen and Hermele21 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, this [111] magnetic field effectively decouples spins between the alternating kagome 2 [111] (a) (c) (b) a2 a1 FIG. 1. From pyrochlore spin ice to kagome ice. (a) The pyrochlore lattice viewed as a set of alternating kagome (blue) and triangular (green) layers along the [111] direction. Spins on the pyrochlore lattice satisfy the “ice rules”: two-in, two-out of each tetraheron. (b) Two-dimensional projection of the pyrochlore spin configuration onto a kagome plane. At the center of each triangle is a representation of the out-of-plane spin: dots (crosses) refer to a spin pointing inward (outward) of each tetrahedron in (a). (c) The associated pseudo-spin S z configuration of Eq. (1) where, filled (empty) circles represents a pseudo-spin up (down). Our QMC simulations of the pseudo-spin Hamiltonian are defined on periodic tori spanned by the primitive vectors a1 and a2 , where ka1 k = ka2 k = 2. and triangular layers of the original pyrochlore structure. To simplest approximation, the system becomes a two-dimensional system of stacked kagome planes,22–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–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/Jz > 2. Next, to include the effect of quantum fluctuations, one may add exchanges from the recent quantum spin ice models.18–21 We consider only those quantum fluctuations discussed by Huang et. al.21 , to obtain a pseudo-spin Hamiltonian on the kagome lattice, Xh J± + − HXYZ = Jz Srz Srz0 − (S S 0 + h.c.) 2 r r hrr 0 i i X J±± + + + (Sr Sr0 + h.c.) − h Srz . (1) 2 r Here, Sr are spin-1/2 operators, with a global z axis and Srz = −1/2 = # in Fig. 1c). This (Srz = 1/2 = 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 Jz > 0 and h Jz 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 Z2 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” Z2 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, 3 0.0 −0.5 Z2 SL −1.0 hSr+ S + i = r+δ 6 0 −1.5 hSr+ i = 0 −2.0 0.2 lobes 0.4 FM 6 0 hSr+i = 0.4 Ns Ns Ns Ns 0.3 (b) =6×6 = 12 × 12 = 18 × 18 = 24 × 24 0.13 0.12 4 0.2 0.8 0.15 0.14 3 0.1 0.6 0.16 mz Z2 SL + 1.0 hSr+ Sr+δ i= 6 0 + 0.5 hSr i = 0 0.5 f0 h/Jz 1.5 (a) (c) 0.44 0.48 0.52 2 χz 2.0 1 0 0.56 J±± /Jz J±± /Jz FIG. 2. Phase diagram of the model. (a) Phases of the model in Eq. 1 for J± = 0 as a function of function of h/Jz and J±± /Jz . The color scale represents the zero momentum occupation f0 obtained from a Ns = 6 × 6 lattice and temperature T = Jz /24. Magnetization (b) and uniform susceptibility (c) as a function of J±± /Jz entering the Z2 QSL lobe at fixed h/Jz = 0.833. Note that a fully spin-polarized phase occurs for h/Jz ≥ 2, which is not illustrated in this phase diagram. we set J± = 0 and explore this possibility for all parameter regimes J±± /Jz and h/Jz , using non-perturbative, unbiased QMC simulations. QUANTUM MONTE CARLO RESULTS We implement a finite-temperature Stochastic Series Expansion34–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 Z2 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 Pestimators for energy, magnetization mz = hmi ˆ = h V1 i Siz i, and uniform spin susceptibility ˆ 2 i − hmi ˆ 2 . The latter two allow us to map χz = VT hm out the broad features of the phase diagram. Further, we measure the off-diagonal spin structure factor37 nαβ q = 1 X iq[(ri +α)−(rj +β)] + e hSri +α Sr−j +β i. Ns r r (2) i j Here, ri points to the sites of the underlying triangular lattice (containing Ns sites) of the kagome lattice (containing V = 3 × Ns sites). The vectors α are the position of each site within the unit cell with respect to the vector ri . 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 f0 = nVM as the ratio of largest eigenvalue nM of the “single-particle” density matrix ρi,j = hSr+i Sr−j i to the number of sites V . The eigenvalues of ρi,j coincide with nq = P α for a translationally invariant system. nα,α q Figure 2 shows the QMC phase diagram for the J± = 0 model of Eq. (1), using data for the condensate fraction f0 . Careful finite-temperature and finite-size scaling, performed up to lattice sizes of V = L × L × 3 = 39 × 39 × 3 and β = Jz /T = 96, is detailed in the Supplemental Information. The magnetization curve and the uniform spin susceptibility across the phase boundary at fixed h/Jz = 0.833 are presented in Fig. 2. The data clearly indicates the existence of two magnetized “lobes” on the phase diagram for J±± /Jz < 0.5 and h/Jz 6= 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/Jz < 0 and h/Jz > 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 (Sr+ Sr+0 + Sr− Sr−0 ) as a spin pair interaction, z which does not conserve the total magnetization Stot . 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: Sqαβ /Ns = β β hSqα S−q i − hSqα ihS−q i, where Sqα = 1 X iq(ri +α ) z e Sri +α . Ns r (3) i P αα If there there is long-range order then Sq = α Sq will scale with system size for at least one value of q. We also measure the bond-bond structure factor using a 4 Sq 0.4 0.2 0.2 0.2 0.0 qy /2π 0.4 0.0 0.0 −0.2 −0.2 −0.2 −0.4 −0.4 −0.4 −0.4 −0.2 0.0 0.2 −0.4 −0.2 0.4 qx /2π 0.0 0.2 −0.4 −0.2 0.4 qx /2π SK , S0 0.03 0.015 0.010 SK S0 0.2 0.4 0.015 BK 0.020 0.0 qx /2π 0.020 0.04 0.025 n0 /V BBq 0.4 qy /2π qy /2π nq 0.02 0.010 0.005 0.01 0.005 0.000 0.000 0.002 0.004 0.006 0.008 0.010 0.00 0.000 0.002 0.004 0.006 0.008 0.010 0.000 0.000 0.002 0.004 0.006 0.008 0.010 1/V 1/V 1/V FIG. 3. Structure factors and absence of order in the lobe. Top: Off-diagonal nq , diagonal Sq , and bond BBq structure factors inside the lobe for a system with Ns = 24 × 24, h/Jz = 0.8333, J±± /Jz = 0.495, and T = Jz /48. Bottom: The corresponding finite-size scaling of candidate peaks at q values where local maxima occur in the structure factor. The zero-momentum peak of BBq has been removed. BBqαβ = 1 X iq(ra −rb ) α β e hBra Brb i, Ns r r (4) a b where Brαa = Si+aα Sj+aα + Si−aα Sj−aα . Nearest neighbor sites iaα and jaα belong to bond α in a unit cell located at position P ra . Again, if there is pair long-range order then BBq = α BBqαα should scale with system size for at least one value of q, with which we define Bq = BBq /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 Bq corresponds to q = 0, which remains finite as V → ∞, meaning that the bond expectation values hSi+aα Sj+aα i 6= 0 is finite in the lobes. This is expected as hSi+aα Sj+aα i 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 single0 loghSr+i Sr−i+dˆxi four-point correlation function. −5 Fit Lobe FM −10 −15 −20 0 2 4 6 8 10 12 d FIG. 4. Exponential decay of correlation functions. The off-diagonal spin correlation function as a function of ˆ direction for a system with Ns = 24 × 24 distance along the x and T = Jz /48. The dashed line corresponds to a fit of the numerical data to the function f (d) = c − dξ − α ln d 5 particle correlations between the mz = ±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 mz = ±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 Z2 gauge symmetry, it is clear that further simulation work should be carried out to address this hypothesis. To confirm the presence of a Z2 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 Z2 structure is not strictly required by the Lieb-Schultz-Mattis-Hastings (LSMH) theorem,44–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 Z2 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 Tb2 Ti2 O7 , Yb2 Ti2 O7 , Pr2 Zr2 O7 , and Pr2 Sn2 O7 . Also, in light of recent experiments47 which suggest that the oft-studied classical spin ice state is only metastable in Dy2 Ti2 O7 , 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 Expansion34 (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 Jz /J±± is large. Using this Hamiltonian breakup, the standard SSE directed loop equations35 were modified to include sampling of off-diagonal operators of the type Sr+ Sr+0 + 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 107 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 6 CPUs running at 2.83 GHz clock speed. ∗ 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 [email protected] Anderson, P. W. The Resonating Valence Bond State in La2 CuO4 and Superconductivity. Science 235, 1196–1198 (1987). Lee, P. A., Nagaosa, N. & Wen, X.-G. Doping a Mott insulator: Physics of high-temperature superconductivity. Rev. Mod. Phys. 78, 17–85 (2006). Ioffe, L. B. et al. Topologically protected quantum bits using Josephson junction arrays. Nature 415, 503–506 (2002). Balents, L. Spin liquids in frustrated magnets. Nature 464, 199–208 (2010). Yan, S., Huse, D. A. & White, S. R. Spin-Liquid Ground State of the S = 1/2 Kagome Heisenberg Antiferromagnet. Science 332, 1173–1176 (2011). Pratt, F. L. et al. Magnetic and non-magnetic phases of a quantum spin liquid. Nature 471, 612–616 (2011). Han, T.-H. et al. Fractionalized excitations in the spinliquid state of a kagome-lattice antiferromagnet. Nature 492, 406–410 (2012). Gingras, M. J. P. & McClarty, P. A. Quantum spin ice: a search for gapless quantum spin liquids in pyrochlore magnets. Reports on Progress in Physics 77, 056501 (2014). Pauling, L. The Structure and Entropy of Ice and of Other Crystals with Some Randomness of Atomic Arrangement. Journal of the American Chemical Society 57, 2680–2684 (1935). Bramwell, S. T. & Gingras, M. J. P. Spin Ice State in Frustrated Magnetic Pyrochlore Materials. Science 294, 1495–1501 (2001). Gingras, M. Spin Ice. In Introduction to Frustrated Magnetism, vol. 164 of Springer Series in Solid-State Sciences, 293–329 (Springer Berlin Heidelberg, 2011). Hermele, M., Fisher, M. P. A. & Balents, L. Pyrochlore photons: The U (1) spin liquid in a S = 1/2 threedimensional frustrated magnet. Phys. Rev. B 69, 064404 (2004). Banerjee, A., Isakov, S. V., Damle, K. & Kim, Y. B. Unusual Liquid State of Hard-Core Bosons on the Pyrochlore Lattice. Phys. Rev. Lett. 100, 047208 (2008). Shannon, N., Sikora, O., Pollmann, F., Penc, K. & Fulde, P. Quantum Ice: A Quantum Monte Carlo Study. Phys. Rev. Lett. 108, 067204 (2012). Molavian, H. R., Gingras, M. J. P. & Canals, B. Dynamically Induced Frustration as a Route to a Quantum Spin Ice State in Tb2 Ti2 O7 via Virtual Crystal Field Excitations and Quantum Many-Body Effects. Phys. Rev. Lett. 98, 157204 (2007). Kimura, K. et al. Quantum fluctuations in spin-ice-like Pr2 Zr2 O7 . Nat Commun 4 (2013). Article. Fennell, T. et al. Magnetoelastic Excitations in the Pyrochlore Spin Liquid Tb2 Ti2 O7 . Phys. Rev. Lett. 112, 017203 (2014). Onoda, S. Effective quantum pseudospin-1/2 model for Yb pyrochlore oxides. Journal of Physics: Conference Series 320, 012065 (2011). Savary, L. & Balents, L. Coulombic Quantum Liquids in Spin-1/2 Pyrochlores. Phys. Rev. Lett. 108, 037202 (2012). 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 Lee, S., Onoda, S. & Balents, L. Generic quantum spin ice. Phys. Rev. B 86, 104412 (2012). Huang, Y.-P., Chen, G. & Hermele, M. Quantum Spin Ices and Topological Phases from Dipolar-Octupolar Doublets on the Pyrochlore Lattice. Phys. Rev. Lett. 112, 167203 (2014). Matsuhira, K., Hiroi, Z., Tayama, T., Takagi, S. & Sakakibara, T. A new macroscopically degenerate ground state in the spin ice compound Dy2 Ti2 O7 under a magnetic field. Journal of Physics: Condensed Matter 14, L559 (2002). Moessner, R. & Sondhi, S. L. Theory of the [111] magnetization plateau in spin ice. Phys. Rev. B 68, 064411 (2003). Isakov, S. V., Raman, K. S., Moessner, R. & Sondhi, S. L. Magnetization curve of spin ice in a [111] magnetic field. Phys. Rev. B 70, 104418 (2004). Macdonald, A. J., Holdsworth, P. C. W. & Melko, R. G. Classical topological order in kagome ice. Journal of Physics: Condensed Matter 23, 164208 (2011). Wills, A. S., Ballou, R. & Lacroix, C. Model of localized highly frustrated ferromagnetism: The kagom´e spin ice. Phys. Rev. B 66, 144407 (2002). Hiroi, Z., Matsuhira, K., Takagi, S., Tayama, T. & Sakakibara, T. Specific Heat of Kagom´e Ice in the Pyrochlore Oxide Dy2Ti2O7. Journal of the Physical Society of Japan 72, 411–418 (2003). Sakakibara, T., Tayama, T., Hiroi, Z., Matsuhira, K. & Takagi, S. Observation of a Liquid-Gas-Type Transition in the Pyrochlore Spin Ice Compound Dy2 Ti2 O7 in a Magnetic Field. Phys. Rev. Lett. 90, 207205 (2003). Tabata, Y. et al. Kagom´e Ice State in the Dipolar Spin Ice Dy2 Ti2 O7 . Phys. Rev. Lett. 97, 257205 (2006). Molavian, H. R. & Gingras, M. J. P. Proposal for a [111] magnetization plateau in the spin liquid state of Tb2 Ti2 O7 . Journal of Physics: Condensed Matter 21, 172201 (2009). Isakov, S. V., Wessel, S., Melko, R. G., Sengupta, K. & Kim, Y. B. Hard-Core Bosons on the Kagome Lattice: Valence-Bond Solids and Their Quantum Melting. Phys. Rev. Lett. 97, 147202 (2006). Damle, K. & Senthil, T. Spin Nematics and Magnetization Plateau Transition in Anisotropic Kagome Magnets. Phys. Rev. Lett. 97, 067202 (2006). Fradkin, E. & Shenker, S. H. Phase diagrams of lattice gauge theories with Higgs fields. Phys. Rev. D 19, 3682– 3697 (1979). Sandvik, A. W. Stochastic series expansion method with operator-loop update. Phys. Rev. B 59, R14157–R14160 (1999). Sylju˚ asen, O. F. & Sandvik, A. W. Quantum Monte Carlo with directed loops. Phys. Rev. E 66, 046701 (2002). Melko, R. G. Simulations of quantum XXZ models on two-dimensional frustrated lattices. Journal of Physics: Condensed Matter 19, 145203 (2007). Dorneich, A. & Troyer, M. Accessing the dynamics of large many-particle systems using the stochastic series expansion. Phys. Rev. E 64, 066701 (2001). Penrose, O. & Onsager, L. Bose-Einstein Condensation and Liquid Helium. Phys. Rev. 104, 576–584 (1956). 7 39 40 41 42 43 44 45 46 47 48 49 50 Giamarchi, T., Ruegg, C. & Tchernyshyov, O. BoseEinstein condensation in magnetic insulators. Nat Phys 4, 198–204 (2008). Kitaev, A. & Preskill, J. Topological Entanglement Entropy. Phys. Rev. Lett. 96, 110404 (2006). Levin, M. & Wen, X.-G. Detecting Topological Order in a Ground State Wave Function. Phys. Rev. Lett. 96, 110405 (2006). Isakov, S. V., Hastings, M. B. & Melko, R. G. Topological entanglement entropy of a Bose-Hubbard spin liquid. Nat Phys 7, 772–775 (2011). Tang, Y. & Sandvik, A. W. Method to Characterize Spinons as Emergent Elementary Particles. Phys. Rev. Lett. 107, 157201 (2011). Oshikawa, M. Commensurability, Excitation Gap, and Topology in Quantum Many-Particle Systems on a Periodic Lattice. Phys. Rev. Lett. 84, 1535–1538 (2000). Hastings, M. B. Lieb-Schultz-Mattis in higher dimensions. Phys. Rev. B 69, 104431 (2004). Nachtergaele, B. & Sims, R. A Multi-Dimensional LiebSchultz-Mattis Theorem. Communications in Mathematical Physics 276, 437–472 (2007). Pomaranski, D. et al. Absence of Pauling’s residual entropy in thermally equilibrated Dy2 Ti2 O7 . Nat Phys 9, 353–356 (2013). Letter. Moessner, R., Sondhi, S. L. & Chandra, P. TwoDimensional Periodic Frustrated Ising Models in a Transverse Field. Phys. Rev. Lett. 84, 4457–4460 (2000). Damle, K. & Senthil, T. Spin Nematics and Magnetization Plateau Transition in Anisotropic Kagome Magnets. Phys. Rev. Lett. 97, 067202 (2006). Sengupta, K., Isakov, S. V. & Kim, Y. B. Superfluidinsulator transitions of bosons on the kagome lattice at noninteger fillings. Phys. Rev. B 73, 245103 (2006). SUPPLEMENTAL INFORMATION Effective gauge theory To understand the possible phases of the Hamiltonian (1), we concentrate on its classical limit first: H = Jz X hrr 0 i Srz Srz0 − h X Srz . (5) r 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 triangle48 . 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 present49 . 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 hxx0 i connecting sites x and x0 on the honeycomb lattice goes through site r on the kagome lattice. We translate Srz to electric flux Ex,x0 using the following relation (see Fig.1(c) and Fig.5): 1 z − Sr (6) Ex,x+ηx µˆ = ηx 6 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) and ˆ3 = (− 3/2, −1/2). The total charge on site x satisfies the Gauss’s law: Qx = 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 our attention and for stimulating discussions motivating this study. This research was supported by NSERC of Canada, the Perimeter Institute for Theoretical Physics, and the John Templeton Foundation. R.G.M. acknowledges support from a Canada Research Chair. Research at Perimeter Institute is supported through Industry Canada and by the Province of Ontario through the Ministry of Research & Innovation. Numerical simulations were carried out on the Shared Hierarchical Academic Research Computing Network (SHARCNET). 3 X Ex,x+ηx µˆ . (7) µ=1 Use these relations, we rewrite (5) as: X Jz Jz − h H= Q2x − ηx Qx . 2 2 x (8) By construction Qx 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 Qx = 0. We now turn our attention to transverse Hamiltonian terms. A transverse operator, Sr+ for example, creates a pair of positive and negative charges. We follow Savary and Balents19 to introduce conjugate operator φx : [φx , Qx0 ] = iδxx0 . Sr+ is mapped to: −ηx ηx + Sr+ = eiηx (φx −φx+µˆ ) s+ xµ ≡ ψx sxµ ψx+ηx µ ˆ. (9) 8 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 channel20,21 through the following mean-field parameters: 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 Qx by 1 and s+ xµ = s+ . Using (8), and (9), we write (1) as: x+ˆ µ,µ X Jz Jz − h Q2x − ηx Qx − H= 2 2 x J ± X ηx + − −ηx ψx sxµ sx+ηx µˆν ψx+η + h.c − (10) (ˆ µ −ˆ ν ) x 2 µ<ν J±± ηx ηx + + −ηx −ηx ψx ψx sxµ sxν ψx+ηx µˆ ψx+η + h.c . ˆ xν 2 The Hamiltonian is invariant under a local U (1) gauge transformation: ψx± → ψx± e∓iαx , s+ xµ → iηx (−αx +αx+ηx µ ˆ) s+ . xµ e (11a) (11b) The azimuthal angle of pseudo spin sr is the vector gauge field 0 ≤ Axµ < 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 Zn gauge structure and intrinsic topological order. These are the Zn spin liquids. The simplest example is the Z2 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 Z2 through the Anderson-Higgs mechanism. Both the gauge field and the spinons should be gapped, − g ≡ hψx+ ψx+η ˆ i, xµ n ≡ hψx+ ψx+ i = hψx− ψx− i, + + − − f ≡ hψx+η ˆ ψx+ηx ν ˆ i = hψx+ηx µ ˆ ψx+ηx ν ˆ i. xµ (12a) (12b) (12c) A finite g indicates semi-classical long-range order since hS + i = hs+ ihψx† ψx+ηx µˆ i ∼ 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 Z2 spin liquid state.20,21,33 We find that in the lobes hSr+ i ∼ g vanishes while hSr+ Sr+0 i ∼ nf remains finite, in agreement with a Z2 spin liquid phase. Absence of breaking of rotational invariance In Fig. 6 we show a finite-size scaling of both the order parameter f0 , as well as the zero-momentum bond occupation B0 , across the transition between the FM phase and the disordered phase in the lobes for T = Jz /24. Outside the lobes both f0 and B0 remain finite as V → ∞, whereas in the lobes f0 vanishes (see Fig. 6(c)) but B0 remains finite (see Fig. 6(d)). Notice that because B0 is finite, the state in the lobes can still break the rotational symmetry of the lattice by having different expectation values of hBrαa i 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 B0αβ = BB0αβ /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 B0αβ 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 B0αβ converge to the same value as V → ∞. Using a linear fit to the finite-sized data, we evaluate the extrapolated matrix elements lim BB0αβ /V in the FM phase and in the lobes, V →∞ 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. 9 0.12 0.4 0.0 0.025 0.06 0.04 (a) (b) 0.02 0.00 0.010 0.020 f0 B0 0.015 0.010 0.006 FM, T = Jz 48 0.006 0.009 0.002 11 12 13 14 15 16 22 33 44 55 66 0.000 0.000 0.003 0.006 0.009 1/V 1/V FIG. 6. (a) Finite-size scaling of the “single-particle” zeromomentum occupation f0 across the phase transition between the ferromagnet and the disordered phases in the lobes. (b) Same as in (a) but for the bond zero-momentum bond occupation B0 . Figures (c) and (d) zoom in the size scalings in the lobes for f0 and B0 , respectively. The temperature is set to T = Jz /24. 0.0172 0.0170 0.0168 =11 =12 =13 =14 =15 =16 αβ αβ αβ αβ αβ =22 =33 =44 =55 =66 0.0166 0.0164 (a) 0.0162 0.000 0.002 0.004 0.006 0.008 0.010 1/V 0.0018 0.0016 0.0014 0.0012 αβ αβ αβ αβ αβ αβ =11 =12 =13 =14 =15 =16 αβ αβ αβ αβ αβ 0.000640 0.000635 0.000630 0.000625 0.000620 0.000615 Lobe, T = Jz 48 (b) 11 12 13 14 15 16 22 33 44 55 66 αβ FIG. 8. The limit of V → ∞ of several matrix elements BB0αβ as a function of the matrix indices αβ in the FM phase and in the lobes ((a) and (b), respectively). 0.0174 αβ αβ αβ αβ αβ αβ 0.000645 (d) V →∞ 0.003 (a) αβ lim BB0αβ /V (c) 0.000 0.000 BB0αβ /V 0.0000005 0.004 0.005 BB0αβ /V 0.0000010 0.0000000 0.008 −2 0.0000015 V →∞ 0.1 0.08 B0 0.2 =0.54 =0.53 =0.52 =0.51 =0.50 =0.49 =0.48 lim BB0αβ /V 0.10 J±±/Jz J±±/Jz J±±/Jz J±±/Jz J±±/Jz J±±/Jz J±±/Jz 0.3 f0 0.0000020 +1.6302×10 =22 =33 =44 =55 =66 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. 0.0010 (b) 0.0008 FIG. 7. The finite-size data for different BB0αβ as a function of 1/V in the FM and in the lobe ((a) and (b), respectively). In figure 10 we present the structure factors inside the lobe for for a system with h/Jz = 0.8333, J±± /Jz = 0.495, and T = Jz /96, on a cluster with Ns = 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. Finite-size scaling inside the lobe for larger system sizes and T = Jz /16. Details of the FM phase and absence of a phase transition at zero field In Fig. 9 we re-examine the finite-size scaling of the most relevant candidate peaks in the structure factors presented in Fig. 3 using larger system sizes and higher temperature. The parameters of the simulations are h/Jz = 0.8333, J±± /Jz = 0.495, and higher T = Jz /16. The largest cluster in the Fig. 9 corresponds to Ns = In figure 11 we present results for spin correlation functions just outside the upper 1/3-magnetization lobe. Both nq and BBq exhibit diverging peaks at zero momentum (which have been removed in the figures), while the diagonal structure factor Sq does not. The fact that f0 remains finite in outside the lobe together with the 0.0006 0.000 0.002 0.004 0.006 0.008 0.010 1/V 10 0.010 0.010 S0 0.008 S0 SK 0.006 0.004 (a) 0.002 0.000 0.000 SK 0.008 0.001 0.006 0.004 (b) 0.002 0.000 0.000 0.002 0.001 1/V 0.006 0.004 0.003 0.002 (c) 0.001 0.000 0.000 BK B0 0.005 B0, BK n0/V 0.006 n0 /V 0.005 0.002 1/V 0.001 0.004 0.003 0.002 (d) 0.001 0.000 0.000 0.002 0.001 0.002 1/V 1/V FIG. 9. Finite-size scaling of the structure factors S0 (a), SK (b), n0 /V (c), BBK and BB0 (d). The maximum system size in the figures corresponds to V = 39 × 39 × 3. The temperature is set to T = Jz /16, while h/Jz = 0.8333 and J±± /Jz = 0.495. Sq BBq 0.4 0.4 0.2 0.2 0.2 0.0 qy /2π 0.4 qy /2π qy /2π nq 0.0 0.0 −0.2 −0.2 −0.2 −0.4 −0.4 −0.4 −0.4 −0.2 0.0 0.2 −0.4 −0.2 0.4 qx /2π 0.0 0.2 −0.4 −0.2 0.4 qx /2π 0.0 0.2 0.4 qx /2π FIG. 10. Off-diagonal nq , diagonal Sq , and bond BBq structure factors inside the lobe for a system with Ns = 24 × 24, h/Jz = 0.8333, J±± /Jz = 0.495, and T = Jz /96. Sq BBq 0.4 0.4 0.2 0.2 0.2 0.0 qy /2π 0.4 qy /2π qy /2π nq 0.0 0.0 −0.2 −0.2 −0.2 −0.4 −0.4 −0.4 −0.4 −0.2 0.0 qx /2π 0.2 0.4 −0.4 −0.2 0.0 qx /2π 0.2 0.4 −0.4 −0.2 0.0 0.2 0.4 qx /2π FIG. 11. Off-diagonal nq , diagonal Sq , and bond BBq structure factors in the FM phase for a system with Ns = 24 × 24, T = Jz /48, and h/Jz = 0.8333. The zero-momentum peaks of nq and BBq have been removed. 11 0.428 nαβ 0 /Ns 0.426 αβ αβ αβ αβ αβ αβ absence of diverging peaks in Sq means that the spins order in the XY plane. To clarify the nature of this ordered phase, we analyze the sublattice zero-momentum occupation nαβ through 0 finite-size scaling shown in figure 12. In the limit of V → ∞, we expect that in an ordered phase the matrix def − + elements f0αβ = lim nαβ 0 /Ns = hSα ihSβ i, i.e., they =11 =12 =13 =22 =23 =33 0.424 Ns →∞ coincide with products of expectation values of the spin operators, which we extract from linear extrapolations of our data in figure 12. The resulting matrix is 0.422 0.01 0.02 0.03 1/Ns FIG. 12. The scaled sublattice zero-momentum occupation nαβ 0 /Ns as a function of 1/Ns . (a) 1.6 (b) 10 1.4 8 1.2 6 0.8 χ f0 1.0 4 0.6 0.4 0.2 6 2 2 0.42060 0.42053 0.42057 0.42059 0.42053 0.42057 ± 10−5 6 2 2 (13) 0.42059 0.42053 0.42057 6 2 2 0.420 0.00 Ns Ns Ns Ns =6×6 =12×12 =18×18 =24×24 0.0 0.0 0.5 1.0 1.5 2.0 J±±/Jz 2 Ns Ns Ns Ns =6×6 =12×12 =18×18 =24×24 0 0.0 0.5 1.0 1.5 2.0 J±±/Jz FIG. 13. 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. The fact that all matrix elements are real, positive, and have same magnitude means that the three spins in the unit cell are either aligned or antialigned along the x direction (consistent with the two possibilities given the Z2 symmetry of the Hamiltonian), thus the ground state is a ferromagnet. This is consistent with spin-wave theory calculations and with the expectation that in the limit of large J±± Jz , when the terms − Sr+ Sr+0 + Sr− Sr−0 = 2 (−Srx Srx0 + Sry Sry0 ) dominate, the spins gain energy by aligning solely along the x axis, where the spin interaction is unfrustrated. Finally, the zero-field results for the order parameter f0 and the uniform spin susceptibility are shown in Fig. 13(a) and (b). The results clearly point towards the absence of a phase transition as J±± /Jz is decreased. The FM phase remains stable down to J±± /Jz = 0 where the extrapolated order parameter f0 is about 11% of its maximum value, which occurs exactly at the U(1) symmetric point J±± /Jz = 1. Similarly, the uniform susceptibility appears to be finite in the limit J±± /Jz → 0. This scenario is similar to what happens in the spin-1/2 XXZ model on kagome lattice, where numerical evidence shows that a FM phase persists down to J± /Jz → 0,31 which was understood as the absence of vortex condensation in a duality analysis of the model.31,50

© Copyright 2019 ExploreDoc