Table Of Content

Dynamics of the Kitaev-Heisenberg Model Matthias Gohlke,1,∗ Ruben Verresen,1,2,∗ Roderich Moessner,1 and Frank Pollmann1,2 1Max-Planck-Institut fu¨r Physik komplexer Systeme, 01187 Dresden, Germany 2Technische Universita¨t Mu¨nchen, 85747 Garching, Germany (Dated: January 18, 2017) Weintroduceamatrix-productstatebasedmethodtoefficientlyobtaindynamicalresponsefunc- tions for two-dimensional microscopic Hamiltonians, which we apply to different phases of the Kitaev-Heisenberg model. We find significant broad high energy features beyond spin-wave the- oryevenintheorderedphasesproximatetospinliquids. Thisincludesthephasewithzig-zagorder of the type observed in α-RuCl , where we find high energy features like those seen in inelastic 3 neutron scattering experiments. Our results provide an example of a natural path for proximate 7 spin liquid features to arise at high energies above a conventionally ordered state, as the diffuse 1 remnants of spin-wave bands intersect to yield a broad peak at the Brillouin zone center. 0 2 n Introduction. The interplay of strong interactions a) b) N2 N1 a and quantum fluctuations in spin systems can give rise J to new and exciting physics. A prominent example are K2 K1 7 quantum spin liquids (QSL), as fascinating as they are N N 2 1 1 hard to detect: they lack local order parameters and are instead characterized in terms of emergent gauge fields. ] z l On the experimental side, spectroscopic measurements e - provide particularly useful insights into such systems, in x y r particular by probing the fractionalised excitations (e.g. t s deconfined spinons) accompanying the gauge field. Such t. measurementscanberelatedtodynamicalresponsefunc- a FIG. 1. (a) Green, red and blue edges correspond to Ki- tions, e.g. inelastic neutron scattering to the dynamical m taev exchange couplings SγSγ with γ = x,y,z. (b) Allowed structurefactor. Onthetheoreticalside,determiningthe i j k-vectors (red lines) for an infinite long cylinder with cir- - ground state properties of such quantum spin models is d cumference L = 6 and periodic boundary condition along 2 already a hard problem, and it is even more challenging n N . BlacknodespicturethepositionofthegaplessMajorana 2 o to understand the dynamics of local excitations. cones. c Here we present a combination of the density-matrix [ renormalization (DMRG) ground state method and a 1 matrix-product states (MPS) based dynamical algo- and one hosting gapless Majorana and gapped flux ex- v rithm to obtain the response functions for generic two- citations (“B phase”)3. If not stated otherwise, we use 8 dimensionalspinsystems. Withthisweareabletoaccess 7 thedynamicsofexoticphasesthatcanoccurinfrustrated the parametrization J = cosα and Kγ = K = 2sinα. 6 systems. Moreover it is also very useful for regular or- If J = 0 and Kγ bond-independent, the Kitaev model 4 is in the B phase, which is stable under time-reversal deredphaseswhereonewouldconventionallyuselarge-S 0 symmetric perturbations as pointed out by Kitaev. Nu- approximations,whichinsomecasescannotqualitatively . 1 explain certain high energy features1,2. merical studies of the ground state phase diagram of the 0 KHM have shown an extended QSL phase for small J We demonstrate our method by applying it to the 7 and four symmetry broken phases for larger J4. currentlymuch-studiedKitaev-Heisenberg model (KHM) 1 model on the honeycomb lattice The dynamical response functions of the pure Ki- : v taev model are known exactly and reveal characteristic Xi H = (cid:88) KγSiγSjγ +J (cid:88)Si·Sj. (1) features6,7, such as a spectral gap due to a spin flip not only creating gapless Majorana but also gapped flux ex- r (cid:104)i,j(cid:105)γ (cid:104)i,j(cid:105) a citations. This feature is perturbatively stable to small The first term is the pure Kitaev model exhibiting J8, but the influence of J on high-energy features (or strongly anisotropic spin exchange coupling3. Neigh- non-perturbatively at low energies) is unclear and of on- boring spins couple depending on the direction of their goinginterest9. Morepressingly,thereappeartobeprox- bond γ with SxSx, SySy or SzSz (Fig. 1). The sec- imate spin liquids10,11, such as possibly the currently ondistheSU(2)-symmetricHeisenbergterm. TheKHM much-studiedα-RuCl 2,5,11–19,whoselow-energyphysics 3 serves as a putative minimal model for several mate- is consistent with spin waves on an ordered background, rials including Na IrO , Li IrO 4, and α-RuCl 5. The but whose broad high-energy features resemble those of 2 3 2 3 3 pure model is an exactly solvable spin-1/2 model stabi- aKSL.Inparticular,forintermediateenergyscalesthere lizing two different Kitaev quantum spin liquids (KSL): are star-like features2 apparently arising from a combi- a gapped Z one with abelian excitations (“A phase”) nation of spin wave and QSL physics. 2 2 ED lattice transformation that maps zigzag to AF and iPEPS 8 stripy to FM22. Plotted are the ground state energy 7 0.2 and the entanglement or von-Neumann entropy S = − 6 Trρredlogρred of the reduced density matrix ρred for − 0.4 5 a bipartitioning of the cylinder by cutting along a ring. GS − 4 S Both the cusps in the energy density and the discon- E 0.6 3 tinuities of the entanglement entropy indicate first or- − der transitions. A careful finite size scaling is diffi- 2 0.8 cult because of the large bond dimension needed and − 1 thus it is not possible to make definite statements about 0 whether the transitions remain first order in the limit 0.0 0.5 1.0 1.5 2.0 L . The symmetry broken phases can be iden- α/π 2 → ∞ tified by measuring the local magnetization. We iden- FIG. 2. Phase diagram for an infinite cylinder with circum- tify a Néel phase ( 0.185 < α/π < 0.487) that extends − ference L = 12 obtained using iDMRG. The black line cor- around the pure anti-ferromagnetic Heisenberg26 point, 2 respondstothegroundstateenergydensityandtheblueline the corresponding zigzag phase (0.513 < α/π < 0.894), to the entanglement entropy for a bipartition of the cylinder a ferromagnetic phase around the pure FM Heisenberg into a left and right half. The insets illustrate the ordering point (0.894 < α/π < 1.427), and its stripy phase pattern of the magnetic phases. Two spin liquid phases ex- (1.559<α/π <1.815). ThetwoKSLsbetweenNéeland ist around the pure Kitaev model (α=0.5π and 1.5π). The zigzagaswellasbetweenFMandstripyareconfirmedto results of ED4 and iPEPS21 are illustrated on top. begapless. Inparticular, ifL isamultipleofsixweuse 2 thefiniteentanglementscalingapproach27–29andextract the expected chiral central charge c = 1 for both KSLs, Inthisarticle,wefirstrevisitthegroundstatephasedi- eachofthetwoMajoranaconescontributingc=1/2. See agram and confirm the previously found phases. The in- also appendix B. Note that when a gapless spin liquid is finitecylindergeometryallowsustonumericallyconfirm placed on a cylinder, the gauge field generically adjusts thatthegaplessnessoftheKSLisrobustthroughoutthe to open a gap30. In order to see gapless behaviour, we entirephase. SecondlyweusearecentlyintroducedMPS have to initiate the iDMRG simulations in the gapless based time evolution algorithm20 to obtain the dynam- sector to access a metastable state (see appendix C for ical spin structure factor. We benchmark our method additional details). The gapped ground state having a by comparing to exact results for the Kitaev model and non-zerofluxthroughthecylinderoverestimatesthesta- find a good agreement. We calculate the spectra of dif- bilityoftheQSLphases. Itisnotablehowwellthephase ferent (non-soluble) phases of the KHM. Most notably, boundaries agree with those from the infinite projected we identify broad high energy continua even in ordered entangled pair state (iPEPS) simulations21. phases that are reminiscent of the broad features observed in recent experiments on α-RuCl and which are 3 moreover similar to the high energy features in the spin liquid phase, thus providing a concrete realisation of the Dynamical structure factor (k,ω). Startingfrom S concept of a proximate spin liquid. a ground state obtained using iDMRG, we calculate Ground state phase diagram. We use the iDMRG (k,ω) by Fourier transforming the dynamical correla- algorithm on the KHM on infinite cylinders to map out Stion function Cγγ(r,t) = Sγ(t)Sγ(0) . The real-time (cid:104) r 0 (cid:105) the phase diagram. We choose cylinder geometries such correlations can be efficiently obtained using a recently that the correspondingmomentumcuts contain thegap- introducedmatrix-productoperatorbasedtimeevolution less Majorana modes of the Kitaev spin liquid. For the method20. This allows for long range interactions result- pure isotropic Kitaev model, there are gapless Majorana ing from unraveling the cylinder to a one-dimensional cones on the corners of the first Brilluoin zone, Fig. 1b. system which render standard methods like the time- The full KHM has a C symmetry which means that evolving block decimation inefficient. Following the gen- 6 in the 2D limit these cones cannot shift. The iDMRG eral strategy laid out in Refs. [32–34], we perform the method determines the ground state of systems of size simulationsforaninfinitecylinderwithafixedcircumfer- L L where L is in the thermodynamic limit and ence. Notethattheentanglementgrowthandtheresult- 1 2 1 × L a finite circumference of up to 12 sites beyond what inggrowthoftherequirednumberofstatesisgenerically 2 is achievable in exact diagonalization. While tradition- slowasweonlylocallyperturbthegroundstateandthus ally iDMRG is used for finding the ground state of one- long times can be reached even in the cylinder geometry. dimensional systems, it has become a fairly unbiased We show results obtained for 0 t T and to avoid ≤ ≤ methodforstudyingtwo-dimensionalfrustratedsystems. Gibbs oscillations we multiply our real-time data with a The resulting phase diagram for L = 12 is shown Gaussian (σ 0.43T). This corresponds to a broad- 2 t in Fig. 2 (for the iDMRG simulations we keep χ = ening in ω-spa≈ce (σ 2.3). We use linear prediction ω ≈ T 1200 states), which agrees with previous studies4,21–25. to allow room for the tail of the Gaussian in real-time, For this L , the system is compatible with the sub- but confirm that the final results do not depend on its 2 3 a) 4 ence L = 6 (red), we see qualitative similarities, such 2 2D 1.5 Re Sz(t)Sz(0) as a spectral gap (dashed lines; slightly obscured by our h i i i 3 cylinder 1.0 numerics finite-time window), a dip where the fluxes suppress the ) exact,cylinder vanHovesingularityoftheMajoranaspectrum6,compa- ω 0.5 exact,2D 0, rable bandwidth and strong low-energy weight. To bet- =2 0.0 ter resolve the spectral gap, we rely slightly on linear k ( 0.5 prediction35 byusingareal-timeGaussianenvelopewith zz − S1 1.0 σt = 0.56T, corresponding to σω 0.045. Two striking − 0 5 10 15 20 25 30 35 40 quantitativedifferencesare(i)the≈spectralgapwhichfor t this circumference is approximately half that of the 2D 00.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 limit, and (ii) the presence of a delta-peak on this gap b) 4 ω ( 4% of total spectral weight). The latter, present for 2D 1.0 Re Sz(t)Sz(0) a≈ny cylinder, vanishes as L2 . The inset compares cylinder h i i i exactreal-timeresultsonthe→cyl∞inder31 withournumer- 3 numerics 0.5 ω) exact,cylinder ics. Despite the true ground state on this cylinder be- 0, ing gapless and MPS only being able to capture gapped =2 0.0 ground states exactly, we still find good agreement for k ( appreciable times. z z S1 0.5 After this benchmarking, we explore (k,ω) in dif- − 0 5 10 15 20 25 30 35 40 S t ferent phases of the KHM shown in Fig. 4, all with σ 0.06. The pure Heisenberg FM (α = π) can be 0 ω ≈ 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 solved in terms of linear spin wave theory (LSWT) and ω numerically captured with bond dimension χ = 2. In- stead of this special point, in Fig. 4a we show results for FIG. 3. Dynamical structure factor Szz(k =0,ω) from our α = 1.1π (corresponding to K = 0.65J) where we still numerical approach compared with exact result (insets show find excellent agreement with LSWT. Note that there is real time data). Exact results were obtained following [6], an extremely small gap ( 0.05J ) despite the presence except for the blue curve in (b)[31]. (a) Gapped KSL on a ≈ | | of anisotropic couplings, as the entire KHM is SU(2)- cylinder with L2 = 10 and anisotropic couplings Kx = −2 aanndd αKy==3πK.z =−13. (b) Gapless isotropic KSL with L2 =6 sdyemr emffeetcrtiscoinnLthSeWdTy.nWameidcso,nwohtiochbsiesrvpereasnuymsatbrolyngreclyaltiend- 2 to the short correlation length and the absence of frac- tionalexcitations. ThepureHeisenbergAFM(withsmall details35. Thence, XXZanisotropy)inFig.4bshowsappreciabledeviations from LSWT, with second order SWT36 giving better 1 (cid:88)(cid:90) ∞ agreement. Moreover,theweightinthespinwavesisap- γγ(k,ω)= ei(ωt−k·r)Cγγ(r,t)dt S 2π proximately halved, indicating the importance of higher r −∞ order magnon contributions. Staying within the Néel normalized as (cid:82) γγ(k,ω) dkdω = (cid:82) dk. If not stated phasebutapproachingtheQSL,spinwavetheorycannot otherwise,wepreSsentresultsfor (k,ω)=(cid:80) γγ(k,ω). even qualitatively describe Fig. 4c, with much weight in S γS very broad high energy features unaccounted for. We benchmark the method by comparing our numerical approach to exact results for the pure Kitaev model. Lastly we focus on a parameter regime producing zig- Figure 3a shows a comparison for the gapped Kitaev zag ordering like that found in α-RuCl 2,11,12. Fig. 5 3 model in the A phase with K /K = 6, the exact so- shows (k,ω) for four different choices of α: the first x y,z S lution for zz(k = 0,ω) shown in black. Our numerics row contains the exact solution for the pure AFM Ki- S (with resolution σ 0.06 in units shown) for an infi- taev model, and the subsequent rows are all numerical ω ≈ nite cylinder with L = 10 (red) agrees well with such results within the zigzag phase with increasing α. For 2 features as gap, bandwidth and total spectral weight. each α we show (k,ω) at fixed ω: the columns dis- S In the real-time data (inset), whilst the numerics agrees play representative low-, mid- and high-energy features, withtheexactsolutionforthecylindergeometry,itover- with parameters L = 12 and time cut-off T = 10 cor- 2 laps with the 2D result only until a characteristic time responding to σ 0.23. We average over the different ω ≈ scalecorrespondingtotheperturbationtravelingaround symmetry broken directions. In appendix D, we show thecylinderandthenfeelingthestaticfluxesinsertedby results for L = 6 and T = 40, revealing that even at 2 the spin-flip. More generally we expect such timescales this resolution the high-energy features stay very broad. (after which 2D physics becomes 1D) to be particularly Thefirstcolumnshowsthelow-energyphysicsoftheKi- significantforsystemswithfractionalization. ForFig.3b taev model being reconstructed into spin wave bands, we take K = K = K = 2 being in the gapless with minima on the edges of the first Brillouin zone. For x y z KSL phase at α = 3π. Com−paring the exact 2D re- α = 0.7π,0.8π these obey the C -symmetry, indicating 2 6 sult (black) to our numerics for a cylinder of circumfer- that the cylinder geometry locally looks like 2D. Inter- 4 S(k ,k = 0,ω) Energy a) 4 1 2 1.4 8.0 0.5 0.6 0.6 1.1 α=1.1π 12 1 3 9 π 0 / k1 ω2 6 1 − 1 3 0.6 3.6 0.3 1.0 0.5 1.9 1 0 0 b) 4 π α=0,∆XXZ =1.1 6.0 0 k/1 3 1 4.5 − 1.1 6.3 0.9 1.9 0.7 1.7 ω2 3.0 1 1 1.5 π 0 / k1 0 0.0 c) 4 1 − α=0.455π 2.0 1.1 5.9 1.5 3.4 0.4 0.7 3 1.5 1 ω2 π 1.0 0 / k1 1 0.5 1 − 1 0 1 1 0 1 1 0 1 -01.0 -0.5 0.0 0.5 1.0 0.0 α − k1⊥/π − k1⊥/π − k1⊥/π k /π 1 FIG. 5. S(k,ω) at three different energies for four models: KSL at α = 0.5π (analytic result, 2D) and zigzag order at FIG. 4. Dynamical structure factor S(k,ω) for cuts k = α=0.55π,0.7π,0.8π (with L2 =12) (kx,0) in different phases of the KHM with the ω-resolution σω ≈0.06. DashedlinesshowresultsfromLSWT.Insetsshow the data for all allowed cuts. (a) Ferromagnetic phase for a is a six-pointed-star whose arms point towards the edges cylinder with L = 12. (b) Antiferromagnet with small spin 2 of the first Brillouin zone. It is interesting to note that anisotropy without Kitaev term (L = 8). Blue line shows 2 if we do not average over different symmetry broken di- second order spin wave calculations. (c) Antiferromagnetic phase in proximity of the KSL (L2 =6). rections, the low-energy physics strongly breaks the C6 symmetry yet the six-pointed star at intermediate energies persists: thus even if we interpret these high energy features as the overlap of broad spin waves, at this point estingly, the high-energy physics of the ordered phases the effect of symmetry breaking has disappeared. Un- is very similar to that of the pure Kitaev model: we der what conditions such a symmetry restoration occurs have broad features centered around k = 0 which are more generally is an interesting question. diffuse w.r.t. ω, with its characteristic energy and width simultaneously decreasing as α increases. The interplay Conlusion. We have presented a new method for ob- between these low- and high- energy features then gives taining the dynamical properties of generic lattice spin rise to different mid-energy shapes. In fact the six spin modelsin(quasi-)twodimensions,whichweexpecttobe wave bands start on the edges of the first Brillouin zone. useful for many future studies. In the KHM, our study Astheenergyincreases,thesebandsbecomeincreasingly reveals several features beyond spin-wave theory even in diffuse,eventuallyoverlappinginaverybroadblobabove theorderedphases,providingamoredetailedpicturefor the symmetric Γ point k = 0. Both spin waves and the concept of a proximate spin liquid as potentially re- blob sharpen as one moves away from the nearby QSL. alised in α-RuCl3. Comparing with inelastic neutron data for α-RuCl 2, Acknowledgements. We are grateful to Roser 3 we find the best qualitative agreement in Fig. 5 around Valenti, Mike Zaletel and Johannes Knolle for stimulat- α = 0.7π. In particular at intermediate energies there ingdiscussions. InparticularwethankJohannesforpro- 5 viding unpublished data for the dynamical correlations was supported in part by DFG via SFB 1143 and Re- of the isotropic Kitaev model on the cylinder. This work search Unit FOR 1807 through grants no. PO 1370/2-1. ∗ These authors contributed equally to this work. Lett. 105, 027204 (2010). 1 B. Dalla Piazza, M. Mourigal, N. B. Christensen, G. J. 23 H.-C.Jiang,Z.-C.Gu,X.-L.Qi, andS.Trebst,Phys.Rev. Nilsen, P. Tregenna-Piggott, T. G. Perring, M. Enderle, B 83, 245104 (2011). D. F. McMorrow, D. A. Ivanov, and H. M. Rønnow, Nat 24 E.Sela,H.-C.Jiang,M.H.Gerlach, andS.Trebst,Phys. Phys 11, 62 (2015). Rev. B 90, 035113 (2014). 2 A.Banerjee,J.Yan,J.Knolle,C.A.Bridges,M.B.Stone, 25 K. Shinjo, S. Sota, and T. Tohyama, Phys. Rev. B 91, M.D.Lumsden,D.G.Mandrus,D.A.Tennant,R.Moess- 054401 (2015). ner, and S. E. Nagler, arXiv:1609.00103 (2016). 26 Note that due to the mapping of the 2D lattice onto a 1D 3 A. Kitaev, Annals of Physics 321, 2 (2006). chain,Mermin-Wagner-Coleman37 appliesatthepureAF- 4 J. Chaloupka, G. Jackeli, and G. Khaliullin, Phys. Rev. HeisenbergpointandsuppresseslongrangeNéelorder.In Lett. 110, 097204 (2013). factitisreplacedbyagappedsymmetry-preservingstate, 5 K. W. Plumb, J. P. Clancy, L. J. Sandilands, V. V. whichextendsoverafiniteregioninparameterspace.For Shankar, Y. F. Hu, K. S. Burch, H.-Y. Kee, and Y.-J. further discussion, see appendix A. Kim, Phys. Rev. B 90, 041112 (2014). 27 P. Calabrese and J. Cardy, J. Stat. Mech. Theor. Exp. 6 J. Knolle, D. L. Kovrizhin, J. T. Chalker, and R. Moess- 2004, P06002 (2004). ner, Phys. Rev. Lett. 112, 207203 (2014). 28 L. Tagliacozzo, T. R. de Oliveira, S. Iblisdir, and J. I. 7 J. Knolle, G.-W. Chern, D. L. Kovrizhin, R. Moessner, Latorre, Phys. Rev. B 78, 024410 (2008). and N. B. Perkins, Phys. Rev. Lett. 113, 187201 (2014). 29 F.Pollmann,S.Mukerjee,A.M.Turner, andJ.E.Moore, 8 X.-Y. Song, Y.-Z. You, and L. Balents, arXiv:1604.04365 Phys. Rev. Lett. 102, 255701 (2009). [cond-mat] (2016). 30 Y.-C. He, M. P. Zaletel, M. Oshikawa, and F. Pollmann, 9 D. Gotfryd, J. Rusnaˇcko, K. Wohlfeld, G. Jackeli, arXiv:1611.06238 [cond-mat] (2016). J. Chaloupka, and A. M. Ole´s, arXiv:1608.05333 [cond- 31 J. Knolle, (2016), unpublished data. mat] (2016). 32 J. A. Kja¨ll, F. Pollmann, and J. E. Moore, Phys. Rev. B 10 T. Suzuki, T. Yamada, Y. Yamaji, and S.-i. Suga, Phys. 83, 020407 (2011). Rev. B 92, 184411 (2015). 33 H. N. Phien, G. Vidal, and I. P. McCulloch, Phys. Rev. 11 A. Banerjee, C. A. Bridges, J.-Q. Yan, A. A. Aczel, L. Li, B 86, 245107 (2012). M. B. Stone, G. E. Granroth, M. D. Lumsden, Y. Yiu, 34 V. Zauner, M. Ganahl, H. G. Evertz, and T. Nishino, J. J.Knolle,S.Bhattacharjee,D.L.Kovrizhin,R.Moessner, Phys.: Condens. Matter 27, 425602 (2015). D. A. Tennant, D. G. Mandrus, and S. E. Nagler, Nat 35 S.R.WhiteandI.Affleck,Phys.Rev.B77,134437(2008). Mater 15, 733 (2016). 36 Z. Weihong, J. Oitmaa, and C. J. Hamer, Phys. Rev. B 12 R.D.Johnson,S.C.Williams,A.A.Haghighirad,J.Sin- 44, 11869 (1991). gleton, V. Zapf, P. Manuel, I. I. Mazin, Y. Li, H. O. 37 S. Coleman, Commun.Math. Phys. 31, 259 (1973). Jeschke, R. Valenti, and R. Coldea, Phys. Rev. B 92, 38 J.D.Reger,J.A.Riera, andA.P.Young,J.Phys.: Con- 235119 (2015). dens. Matter 1, 1855 (1989). 13 H.-S.KimandH.-Y.Kee,Phys.Rev.B93,155143(2016). 39 S. R. White, R. M. Noack, and D. J. Scalapino, Phys. 14 L. J. Sandilands, Y. Tian, A. A. Reijnders, H.-S. Kim, Rev. Lett. 73, 886 (1994). K. W. Plumb, Y.-J. Kim, H.-Y. Kee, and K. S. Burch, Phys. Rev. B 93, 075144 (2016). 15 L.J.Sandilands,C.H.Sohn,H.J.Park,S.Y.Kim,K.W. Kim, J. A. Sears, Y.-J. Kim, and T. W. Noh, Phys. Rev. B 94, 195156 (2016). 16 J. Chaloupka and G. Khaliullin, Phys. Rev. B 94, 064435 (2016). 17 F. Lang, P. J. Baker, A. A. Haghighirad, Y. Li, D. Prab- hakaran,R.Valenti, andS.J.Blundell,Phys.Rev.B94, 020407 (2016). 18 G. B. Hala´sz, N. B. Perkins, and J. van den Brink, Phys. Rev. Lett. 117, 127203 (2016). 19 A. Koitzsch, C. Habenicht, E. Mu¨ller, M. Knupfer, B. Bu¨chner, H. C. Kandpal, J. van den Brink, D. Nowak, A. Isaeva, and T. Doert, Phys. Rev. Lett. 117, 126403 (2016). 20 M. P. Zaletel, R. S. K. Mong, C. Karrasch, J. E. Moore, and F. Pollmann, Phys. Rev. B 91, 165112 (2015). 21 J.OsorioIregui,P.Corboz, andM.Troyer,Phys.Rev.B 90, 195102 (2014). 22 J. Chaloupka, G. Jackeli, and G. Khaliullin, Phys. Rev. 6 Appendix A: 1D vs 2D physics: symmetry breaking L2=6 0.4 L2=8 From Monte-Carlo studies38 it is known that the L2=10 ground state of the Heisenberg antiferromagnet (AFM) L2=12 0.3 ∆=1.0 on the honeycomb lattice displays symmetry breaking i S| ∆=1.1 Néelorder. However,whenweplacetheHeisenbergAFM h|0.2 on an infinitely long cylinder of finite circumference, it is in principle a 1D system and the presence of a con- 0.1 tinuous symmetry in fact forbids spontaneous symmetry breaking37. Instead we numerically find a gapped state 0.0 which preserves both spin rotation and translation sym- 101 102 103 metry. This is analogous to the results for stacking an χ evennumberofcoupledspin-1 Heisenbergchains39. The 2 transition from 1D to 2D can be understood by noting FIG. 6. The absolute on-site magnetization for the pure that this symmetry-preserving state is effectively Néel- Heisenberg AFM (solid) and for the AFM XXZ model with likewithinacorrelationlengthξ,thelattergrowingwith ∆=1.1 anisotropy (dashed) for different circumferences circumference. Similarly to how one determines spontaneous symmetry breaking from finite size scaling in the context of exact diagonalization, one can conclude that the 2D limit achieves Néel order by scaling with respect to circumference. 2.2 The presence of a gap implies this symmetry- preserving state is stable under SU(2)-breaking pertur- 2.1 bations. ForexampleforL =6itextendsover 0.2π 2 − ≤ α 0.43π, with a Néel order arising for larger α until 2.0 ≤ we hit the spin liquid. The stability of this symmetry- S 1.9 preserving state under Kitaev perturbations is presum- ably related to the fact that the Néel order which arises 1.8 α=1.500π inthe2Dlimitwouldhaveaverysmallspingap. Thisis α=1.530π different for XXZ-type perturbations, which induce Néel 1.7 α=1.540π order for relatively small anisotropies as shown in Fig. 6 α=1.546π 1.6 (with∆=1.1),whereourstateisnumericallyconverged 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 (for large χ) and the physics quickly becomes indepen- log(ξ)/6 dent of circumference. TheDMRGsimulationsuseaparameterχwhichgives FIG.7. EntanglemententropyS andlogarithmofcorrelation an upper bound on the entanglement. By limiting χ length ξ for different bond dimensions. The lines correspond we can find a variational state with ξ < L . Locally to a central charge of c=1. 2 this state then looks 2D and hence we can have symmetry breaking even for the SU(2)-symmetric Heisenberg model, as confirmed in Fig. 6. As we increase χ, eventu- allyξ becomesoftheorderofL ,whichsignalsthetran- 2 sition of 2D to 1D physics and the symmetry-preserving state arises. For L=12 the necessary ξ is already out of the entanglement entropy S scales logarithmically with reach,explainingtheeffectiveNéelorderweseeinFig.2. the correlation length ξ. In the MPS formalism, this Similarly,inthezigzagphasethereisanextendedregion is known as Finite-Entanglement Scaling with S = χ with a gapped symmetry-restored ground state. This c/6logξ , where χ is the bond dimension of the MPS χ is in keeping with the sublattice transformation, which and c is the chiral central charge28,29. maps the zigzag to the Néel phase (in particular α= 3π 4 maps onto α=0). Fig.7showsS andlogξ forvariousMPSbonddimen- sions χ of up to 1024. The lines serve as a guide to the Appendix B: Entanglement scaling of the gapless eye corresponding to a slope with c = 1. We observe a KSL goodmatchofthescalingforthepureKitaevspinliquid at α = 3/2π. This reflects the fact, that the KSL can Matrix-productstates(MPS)cannotcapturealgebraic bemappedtoafreefermionproblemwithtwoMajorana ground state correlations. However, increasing the bond cones in the first Brillouin zone, each contributing 1/2 dimension gives an increasingly accurate estimate of the to the central charge. The gapless nature persists within wave function. Calabrese and Cardy27 have shown that the whole KSL phase and the scaling suggests c=1. 7 ED iPEPS DMRG thedesiredsector. TableIcontainsthephasetransitions L =6 L =12 for the gapped and the gapless sector and compares it 2 2 gapped gapless gapped gapless to exact diagonalization (ED) and infinite Projected En- AF/KSL 0.488 0.487 0.484 0.494 0.485 0.487 tangled Pair States (iPEPS). As the gapped sector has a lower energy, its stability is enhanced and widens the KSL/ZZ 0.510 0.513 0.523 0.513 0.514 0.512 KSL phase. This effect is more pronounced for a small FM/KSL 1.399 1.432 1.405 1.44 1.421 1.428 circumference L =6. 2 KSL/ST 1.577 1.557 1.573 1.548 1.562 1.558 TABLE I. Transition points α/π for different circumferences Appendix D: Dynamics of L =6 cylinder sectorscomparedtoexactdiagonalization(ED)9 andinfinite 2 Projected Entangled Pair States (iPEPS)21. In Fig. 8 we show (k = 0,ω) for the same choices S of α as in Fig. 5, but now with a sharper ω-resolution (corresponding to T = 40) which is possible due to a Appendix C: Ground sectors of the KSL on the smaller circumference (L = 6). The finer features are cylinder 2 mostlikelydiscretizationeffectsduetothefinitecircum- ference, but the main points are that the broadness in Similar to the plaquette operators Wp =(cid:81)j∈ σjγj we ω-space persists despite a finer resolution, and that the define a loop operator around the cylinder as high-energy feature gets squeezed downward as we get (cid:55) further away from the nearby spin liquid. Note that the W = (cid:89) σγj , (C1) latter is a meaningful statement and not just due to an l j overall α-dependent scaling of the Hamiltonian since the j∈loop minimaofthespinbands(asshowninthefirstcolumnof Fig.5)donot comedowninenergy(allatapproximately where γ = x,yz corresponds to the bond that is not i part of the l{oop at}site i. Following Kitaev3, W can be ω =0.4). l expressed in terms of Z gauge field variables u 2 jk 8 (cid:89) W˜ = u . (C2) α=0.5π l jk (j,k)∈loop 6 ) ω 0, α=0.55π For our choice of lattice periodicity, both loop opera- =4 tors are related by a minus sign. Thus, W˜l +1 (pe- k → ( α=0.7π riodic boundary condition of the fermions) translates to S W 1, which corresponds to the gapless sector if the 2 l →− cylinder is chosen such that cuts in reciprocal space go α=0.8π through the nodes of the Majorana cones. The second 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 sector (antiperiodic boundary condition of the fermions) ω is always gapped and has a lower ground state energy than the gapless sector. Regardingthecomputationofthegroundstate,wecan FIG.8. ComplementingFig.5: S(k=0,ω)forα=0.5π (2D analytic result) and α=0.55π,0.7π,0.8π (with L =6). now make use of the loop operator and initialize DMRG 2 with a state ψ that has ψ W ψ = 1 depending on l | (cid:105) (cid:104) | | (cid:105) ±