Dynamics of the Kitaev-Heisenberg Model

3 MB·

by  Matthias Gohlke

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´eel 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´eeland
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 ob-
served 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´eel
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 numeri-
cal 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 ener-
gies 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´eelorder.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´eelorder. 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´eel- 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 sponta-
neous symmetry breaking from finite size scaling in the
context of exact diagonalization, one can conclude that
the 2D limit achieves N´eel 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´eel 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´eel order which arises 1.8 α=1.500π
inthe2Dlimitwouldhaveaverysmallspingap. Thisis α=1.530π
different for XXZ-type perturbations, which induce N´eel 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 symme-
try 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´eelorderweseeinFig.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´eel 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) ±