Interaction driven phases in the half-filled honeycomb lattice: an infinite density matrix renormalization group study
Johannes Motruk, Adolfo G. Grushin, Fernando de Juan, Frank Pollmann
IInteraction driven phases in the half-filled honeycomb lattice: an infinite densitymatrix renormalization group study
Johannes Motruk, Adolfo G. Grushin, Fernando de Juan,
2, 3 and Frank Pollmann Max-Planck-Institut für Physik komplexer Systeme, Nöthnitzer Str. 38, 01187 Dresden, Germany Materials Science Division, Lawrence Berkeley National Laboratories, Berkeley, CA 94720, USA Department of Physics, University of California, Berkeley, CA 94720, USA (Dated: August 21, 2018)The emergence of the Haldane Chern insulator state due to strong short range repulsive interac-tions in the half-filled fermionic spinless honeycomb lattice model has been proposed and challengedwith different methods and yet it still remains controversial. In this work we revisit the problemusing the infinite density matrix renormalization group method and report numerical evidence sup-porting i) the absence of the Chern insulator state, ii) two previously unnoticed charge orderedphases and iii) the existence and stability of all the non-topological competing orders that werefound previously within mean field. In addition, we discuss the nature of the corresponding phasetransitions based on our numerical data. Our work establishes the phase diagram of the half-filledhoneycomb lattice model tilting the balance towards the absence of a Chern insulator phase for thismodel.
I. INTRODUCTION
Topological phases of matter are remarkably robuststates; their responses to external fields are governed bytopological invariants and thus many of their most impor-tant properties are insensitive to local perturbations.
Topologically protected responses result in intrinsicallynovel phenomena including fractionalization of quantumnumbers or transport governed by quantum anomalies and can lead to diverse applications, ranging from faulttolerant quantum computation to spintronics. Addedto the remarkable experimental discoveries of new topo-logical phases in both two and three dimensions, theseideas have boosted a sustained and voluminous scientificeffort in the last decade that attempts to classify themand determine when and how can they emerge. This work addresses a concrete question regarding theemergence of topological phases that has so far remainedcontroversial: the existence of the Chern insulator phasetriggered by short range repulsive interactions in thefermionic spinless half-filled honeycomb tight bindingmodel. The Chern insulator state, first identified byHaldane in the particular case of the honeycomb lattice,is a zero field analogue of the integer quantum Halleffect; the Hall conductivity contribution of each bandis quantized in integer units of e /h and determined bya topological invariant, the Chern number. Realising aChern insulator in nature is not a trivial task; time-reversal must be broken with the resulting magnetic fluxaveraging to zero over the unit cell, and hence over theentire sample.In a proof of principle, a set of mean field studies discussed how the Chern insulator can emerge fromrepulsive short range interactions. It was shown thatspinless fermions hopping in the half-filled honeycomblattice with nearest and next-to-nearest neighbor inter-actions, V and V respectively [see Fig. 1], displayed aChern insulator phase as described by Haldane. The non-interacting realization of this Chern insulator phasehas complex next-to-nearest neighbor hoppings [seeFig. 2(g)]. Within the mean field paradigm it occurs in aregion with V > V where V spontaneously breaks timereversal symmetry generating complex next-to-nearestneighbor hopping strengths. The latter condition wasshown not to be a generic requirement; Chern insulatorphases can be realized in mean field with only V atthe expense of enlarging the unit cell and doping thesystem. In addition, analogous topological phaseshave been obtained in other lattices such as the π − fluxmodel and the Kagome lattice. The important question that still remains to be an-swered is whether the interaction induced Chern statesurvives after the effect of quantum fluctuations is in-cluded. The Haldane Chern insulator phase competeswith more conventional but also interesting orders, thatcan jeopardize its emergence and are depicted schemat-ically in Fig. 2. In the V (cid:29) V limit, a charge orderedstate depicted in Fig. 2(b) is expected, where the A andB sublattices are populated differently. Being con-nected to the classical ground state at V /t → ∞ withonly one occupied sublattice, this state (CDW I) hasbeen found to be very robust both in mean field andexact diagonalization. For V (cid:29) V on the otherhand it was shown under the mean field paradigm that the phase space originally attributed to the Hal-dane Chern insulator was severely reduced by the pres-ence of a sublattice charge modulated state (or CMs)depicted in Fig. 2(c). The CMs phase was corrobo-rated to survive quantum fluctuations in exact diag-onalization and within a variational Monte Carloapproach. The exact nature of the phase was furtherdiscussed in Ref. 19 by computing the charge structurefactor; its suppressed Fourier component at Γ suggestedthe absence of charge imbalance between the sublatticesand was therefore termed CM phase. In addition, at in-termediate V (cid:38) V a Kekulé bond order depicted a r X i v : . [ c ond - m a t . s t r- e l ] S e p FIG. 1. The top left illustration shows a schematic represen-tation of the interacting tight binding model considered inthis work: fermions hopping in a half-filled honeycomb latticewith nearest-neighbor hopping t (solid lines) and interactingvia nearest- and next-to-nearest neighbor interactions, V and V respectively (curved lines) as defined by (1). The right il-lustration shows the iDMRG unit cell used in this work with × unit cells yielding a cylinder with a circumference of L y = 12 sites depicted in the bottom left illustration. in Fig. 2(d) was found to be stable within mean fieldtheory. This three-fold degenerate bond order triples theoriginal honeycomb two atom unit cell breaking itstranslational symmetry. Although numerical evidenceconsistent with such bond order was found in exact
FIG. 2. Different orders allowed in the half-filled spinless hon-eycomb lattice considered in this work: (a) semimetal, (b)charge density wave I (CDW I), (c) sublattice charge mod-ulated (CMs), (d) Kekulé, (e) CDW II, (f) CDW III and(g) Haldane Chern insulator phases. Phases (a)-(d) and (g)were found within mean field theory by Refs. 9, 10, and 12.The survival of the Haldane Chern insulator phase (g) waschallenged within exact diagonalization with periodic bound-ary conditions by Refs. 18 and 19 but found in Ref. 20 foropen boundary conditions. The generic charge imbalance be-tween the sublattices in phase (c) was not found within exactdiagonalization. The two sublattices are depicted in blueand orange, < ∆ < (cid:37) < / with (cid:37) + ∆ < / describe thedeviations from half-filling per site, i.e. (cid:104) n (cid:105) = 1 / ± α with α ∈ { (cid:37), (cid:37) ± ∆ } . diagonalization, further insights are needed to cor-roborate its existence in the thermodynamic limit.Finally, the Haldane Chern insulator phase of Fig. 2(g),stable within mean field, was found however to beabsent in exact diagonalization with periodic boundaryconditions and cluster perturbation theory sug-gesting that quantum fluctuations indeed can destabilizeits emergence. Nonetheless, this interpretation waschallenged by Ref. 20 that observed hints of this phasein exact diagonalization with open boundary conditionsand variational Monte Carlo.The contradictory numerical evidence regarding theexistence of the Chern insulator phase in particular,and other competing phases in general, needs of analternative approach that includes quantum fluctua-tions while minimizing finite size effects. In this workwe therefore revisit the controversies left unsolvedby previous studies using the infinite density matrixrenormalization group method (iDMRG) . Thisvariational method determines the ground state of sys-tems of size L x × L y where L x is in the thermodynamiclimit and L y goes beyond what is achievable in exactdiagonalization. Traditionally a method for finding theground state of one-dimensional systems, iDMRG hasrecently been successfully applied to two-dimensionalsystems. The infinite and finite DMRG method wereindeed shown to be able to characterize the propertiesof fractional quantum Hall, Z quantum spin liquid, chiral spin liquid and bosonic and fermionic fractionalChern insulating states. As long as entanglementremains low and the state has short correlation lengthsthe ground state can be represented faithfully by aproduct of matrices –termed matrix product state(MPS)– of a computationally affordable dimension χ ,known as the bond dimension. This requirementis commonly met by gapped systems and in particularby the orders argued above to be expected instabilitiesin the half-filled honeycomb lattice of interest here.Moreover, the iDMRG method deals with an infinitesystem and thus, given a suitable size of the unit cellused for the simulations, it can directly probe groundstates that spontaneously break the original symmetriesof the hamiltonian and the phase transitions amongthem, while allowing for quantum fluctuations to play arole.Motivated by these advantages we have mapped outthe { V , V } phase diagram for the half-filled honeycomblattice in the infinite cylinder geometry with the iDMRGmethod [see Fig. 1]. Our main findings are summarizednext. First we show compelling numerical evidence thatsupports the absence of the Chern insulator state in awide region of { V , V } phase space. Second, we find twonew phases not reported previously neither within meanfield nor the first exact diagonalization results and analyze their semiclassical features. Third, weprovide numerical support for the existence and stabilityof all the competing orders that were found previouslywithin mean field. In particular we characterize thetwo more controversial states, the CMs and Kekuléstates, by computing bond and charge expectationvalues and present further arguments of the existenceof the former from the semiclassical large interactionlimit. Fourth, we provide numerical evidence to assessthe first or second order character of the relevant phasetransitions. Concretely we address this issue by identi-fying the features that the correlation length ξ and theentanglement entropy S present as a function of { V , V } .This work is structured as follows. After describingthe method and model in section II we characterize thedifferent phases that we find in section III. The corre-sponding phase transitions among them are analyzed insection IV and we end with a discussion and prospect ofour results in section V. Appendix A includes details ofour semiclassical analysis as well a discussion of the CMsstate order parameter. II. MODEL AND METHOD
We investigate a system of spinless fermions hoppingon a honeycomb lattice with real nearest neighbor hop-ping t ≥ interacting via nearest and next-to-nearestneighbor interactions { V , V } ≥ { , } respectively [seeFig. 1]. The Hamiltonian for this system can be writtenas H = − t (cid:88) (cid:104) i,j (cid:105) ( c † i c j +h . c . )+ V (cid:88) (cid:104) i,j (cid:105) n i n j + V (cid:88) (cid:104)(cid:104) i,j (cid:105)(cid:105) n i n j , (1)where c i ( c † i ) annihilates (creates) an electron at the i -thsite of the honeycomb lattice.In order to find the ground state of the system inthe { V , V } phase space we employ the iDMRG algo-rithm on an infinite cylinder geometry [see Fig. 1].As discussed above, such infinite geometry allows toprobe ground states with spontaneous symmetry break-ing while taking into account quantum fluctuations. Fea-sible system sizes for our purposes are cylinders of cir-cumference L y = 6 , , and , given the exponentialgrowth of computational cost with the circumference.For our calculations we choose L y = 12 and a unit cell of36 sites depicted in Fig. 1.We pick this geometry for the following two reasons.First, we have to keep in mind the structure of the recip-rocal space. Since we work in the thermodynamic limitin the x -direction along the cylinder, the momentum in x -direction is a continuous quantity. However in the fi-nite y -direction around the cylinder, the momentum is adiscrete variable. In the non-interacting case, the modelforms a Dirac semimetal and the Fermi surface is locatedat the K and K (cid:48) points. In order to capture the lowenergy physics correctly, it is of key importance that K and K (cid:48) are allowed momenta in our unit cell which im-plies that only L y = 6 and are suitable circumferences of our cylinder. Second, we have the freedom of fixingthe length of our unit cell in the x -direction since thecomputational for that only scales linearly. By choosingthree rings of 12 sites each, all orders are commensuratewith our unit cell and a further enlargement would notlead to different results in the parameter region we inves-tigated. All data we show in the remainder of the textis computed for L y = 12 and bond dimension χ = 1600 unless otherwise stated. III. PHASE DIAGRAM
We have mapped the phase diagram as a function of { V , V } with the method described above. Our mainresults are summarized in Fig. 3. We find six differentphases: two new charge order phases labeled CDW II andCDW III and four previously reported mean field orders,the Kekulé, CMs, CDW I and semimetal phases.We have characterized each phase through their chargeand bond ground state expectation values. At each sitethe charge expectation value is defined by n i = (cid:68) c † i c i (cid:69) . (2)The bond ground state expectation value on the otherhand is defined as t ij = (cid:68) c † i c j + h . c . (cid:69) . (3)Next we discuss the features of the different phases in thephase diagram in Fig. 3. A. Semimetal phase
We start by discussing the semimetal phase in thephase diagram shown in Fig. 3. At V /t = V /t = 0 the honeycomb lattice with nearest neighbor hopping itsknown to be described by a low energy theory in terms oftwo massless Dirac fermions. Short range interactionsare irrelevant in the renormalization group sense andtherefore they can only drive a transition to an orderedstate when they have a magnitude comparable to thenearest neighbor hopping strength t . Such perturbativeanalysis guarantees that the semimetal is stable withinthe region { V , V } (cid:46) t , only allowing for a uniformrenormalization of the hopping strength t by interactions.For the semimetal phase in Fig. 3, and to numericalaccuracy, we find by computing the charge expectationvalue (2) that n i = 1 / for all sites, indicating that thisphase is not charge ordered.From (3) and choosing i, j to be nearest neighbors wefind a small asymmetry between bonds pointing alongand around the cylinder axis, with a relative difference of ∼ − . Such asymmetry should –up to a very small ef-fect due to the cylinder geometry which is always presentfor finite L y – vanish in a perfect semimetallic phase and FIG. 3. Left: Phase diagram obtained with iDMRG calculations on an infinite cylinder of circumference L y = 12 , χ = 1600 .The phase boundaries, especially at second order transitions (see Sec. IV) are to be taken with some error which can beestimated from Figs. 4 to 8. We draw sharp lines here for better clarity. Right: Representative charge and bond strengthpatterns for the four phases with largest unit cell discussed in the main text. The area of the blue circles is proportional tothe particle number expectation value on the respective site given by Eq. (2). The thickness of the ellipsoids on the bonds isproportional to the amplitude t ij between nearest neighbors defined by Eq. (3). The unit cells for each phase are depicted bythe red polygons. These correspond to (a) CMs phase with V /t = 0 . , V /t = 3 . , (b) Kekulé phase with V /t = 5 . , V /t = 1 . ,(c) CDW II phase with V /t = 5 . , V /t = 3 . , (d) CDW III phase with V /t = 9 . , V /t = 2 . . For (c) dashed lines indicatedefect lines of rotated hexagons that have zero energetic cost in the classical limit (see main text). indeed it is severely reduced as the bond dimension χ isincreased. This suggests that the cylinder geometry im-plemented in iDMRG artificially differentiates bonds inits two perpendicular directions. The asymmetry inducedby the finite bond dimension vanishes as χ is increasedand thus it is not a physical effect. This is consistentwith a semimetal state; the logarithmic divergence of en-tanglement of a metallic state requires a matrix productstate with χ → ∞ . The bond asymmetry reduces theentanglement by shifting the Dirac cones away from the K and K (cid:48) points. Together, the previous numerical evidence are consistentwith the semimetallic state. We note that numericallythe semimetal state extends beyond { V , V } (cid:46) t towardshigher interaction strengths through a narrow semimetaltrench in the phase diagram. The larger size of this re-gion at L y = 6 (not shown) suggest that it is likely toshrink as the circumference of the cylinder is increasedbut a definitive statement requires going beyond the nu-merically accessible sizes. B. Charge density wave I (CDW I)
Upon increasing V we identify a charge density wavestate labeled CDW I in the phase diagram of Fig. 3. Thisstate has a two site unit cell and is characterized by afinite order parameter of symmetry B under the sym-metry group C (cid:48)(cid:48) v B = (cid:104) n A (cid:105) − (cid:104) n B (cid:105) , (4)that can be calculated using (2), where n A and n B arethe fermion densities in the A and B sublattice sites ofthe two site unit cell respectively. The resulting chargeorder is depicted schematically in Fig. 2(b). For instance, at V /t = 4 and V /t = 0 . we find that B = 0 . . Themagnitude of B increases with V and to numerical ac-curacy this phase has no appreciable bond order.For large V (cid:29) t such a state is a natural insta-bility since the energy is minimized by a charge im-balance between the two sublattices. Indeed, it hasbeen found in a mean field approximation, ex-act diagonalization, and quantum Monte Carlosimulations. C. Sublattice charge modulated phase (CMs)
For small t (cid:28) V the ground state is classically de-generate. Within mean field it was shown that as longas V > V the system chooses a charge ordered patternwith charge modulation of wavevector K or K (cid:48) , termedthe CMs phase and depicted in Fig. 2(c). The orderparameters for any such modulation can be taken as thecorresponding Fourier components of the charge expecta-tion values of Eq. (2), namely (cid:104) n A,B ( K ) (cid:105) = (cid:104) n A,B ( K (cid:48) ) (cid:105) ∗ ,where A and B indicate the two sublattices. It is con-venient to arrange these order parameters in a basisthat has well defined transformation properties under thesymmetry of the lattice. Under the action of the symme-try group C (cid:48)(cid:48) v , these order parameters transform as thefour dimensional G (cid:48) representation G (cid:48) x = Re[ n A ( K ) − n B ( K )] , (5) G (cid:48) y = Im[ n A ( K ) − n B ( K )] , (6) G (cid:48) x = Im[ n A ( K ) + n B ( K )] , (7) G (cid:48) y = Re[ n A ( K ) + n B ( K )] . (8)The CMs state corresponds to a finite expectation valueof both the B order parameter and G (cid:48) = (1 , , , andthe states generated from these by operations of the sym-metry group. A scalar order parameter for the CMsphase can be defined as B Re[( G (cid:48) ) − G (cid:48) ( G (cid:48) ) ] , where G (cid:48) i = G (cid:48) ix + iG (cid:48) iy (see appendix A). The CMs structureis schematically shown in Fig. 2(c). It minimizes thelarge cost attributed to V by reversing two nearest neigh-bor dimers at the expense of paying the small energeticcost determined by V . Within exact diagonalization ithas been argued that the CMs state survives quantumfluctuations. However, the sublattice charge imbal-ance predicted in mean field was suggested to vanish inRef. 19 due to the absence of an enhanced Fourier com-ponent of the charge structure factor at the Γ point. Wenote that the CMs scalar order parameter defined abovequantifies the charge imbalance since B has to be non-zero for it not to vanish.With the iDMRG method we find that the CMs appearsdirectly above the semimetallic phase [see Fig. 3]. A sam-ple of the six site unit cell charge and bond pattern ob-tained from (2) for V /t = 0 . , V /t = 3 . is shown inFig. 3(a). The circles at each site have an area propor-tional to the strength of the charge at the given site,while the thickness of the links between the sites rep-resent the bond strengths. This bond and charge pat-tern survives in the region labeled CMs of the phase dia-gram and coincides with that obtained within mean fieldtheory depicted schematically in Fig. 2(c). The cor-responding charge values for V /t = 0 . , V /t = 3 . aregiven by (cid:37) = 0 . and ∆ = 0 . . Moreover, we find aclear sublattice imbalance that grows as V is increasedthat justifies the label CMs rather than the simpler CM.We note also that the bond ordering of this phase is as-sociated to the charge order: two neighboring sites withsimilar high (or low) charge densities that deviate fromhalf-filling suppress the hopping between them due to thelarge (small) number of filled (empty) states at that site.We have also accounted for the existence of the CMs statefrom a strong coupling perturbation theory analysis. Byexactly diagonalizing the interaction part of H in Eq. (1)and then finding the first order hopping corrections t/V , in perturbation theory, we determined the ground state | GS (cid:105) of a × cluster. We then computed different scalarcorrelation functions of the charge order parameters, andfound that only (cid:10) GS | B Re[( G (cid:48) ) − G (cid:48) ( G (cid:48) ) ] | GS (cid:11) is fi-nite. This confirms that in the strong coupling limit, theground state is indeed of the CMs form. The details ofthe procedure are presented in appendix A. D. Kekulé bond order
The next phase we identify is the Kekulé bond order.Like the CMs, it has a six site unit cell depicted schemati-cally in Fig. 2(d) with uniform charge order and two typesof bonds, strong and weak. Under the mean field approx-imation the state arises between the CMs and the CDW Istate. Although in previous exact diagonalization hintsof this state are also observed, the evidence supporting its occurrence is not entirely conclusive and its tripledunit cell turns the analysis of larger clusters challenging.Within iDMRG we find that this state is stable in afinite but smaller region compared to both exact diago-nalization and mean field. In contrast to the anisotropyfound in the semimetal phase, the Kekulé order remainsstable as the bond dimension χ is increased and staysfinite upon extrapolation to infinite χ . We show a rep-resentative of the state’s numerically obtained chargeand bond order patterns in Fig. 3 (b) calculated withEqs. (2) and (3), respectively. The bond thickness is pro-portional to its strength. From Fig. 3(b) it is apparentthat the state has indeed two bond strengths arrangedas in Fig. 2(d) and no charge order. Therefore the sys-tem chooses to break the original two site unit cell trans-lational symmetry falling into one of the three distinctKekulé ground states. At V /t = 5 . and V /t = 1 . thetwo types bonds are found numerically to be t = 0 . and t = 0 . . E. Charge density wave II (CDW II)
All phases that we have described so far do not differfrom those predicted by mean field and the firstexact diagonalization studies.
At finite V andsufficiently large V we observe a different charge orderstate (CDW II), see Fig. 3 left panel.To gain a first insight on this phase before analyzing thenumerical data it is instructive to discuss its classicallimit. A strong coupling analysis (see appendix A fordetails) reveals quite generically that at t = 0 there aretwo classical regimes separated by the V = 4 V line.For V > V the state CDW I is favored with the samecharge pattern as in Fig. 2(b) but classical occupations0 or 1. At V = 4 V there exists a classically degenerateground state manifold, to be discussed in the next sec-tion in the context of the CDW III phase. The analysisfor V < V with t = 0 reveals that there are essentiallytwo types of degenerate states according to their trans-lational symmetry. The first is a stripe-like phase witha four atom unit cell depicted in Fig. 2(e) (see also 41)that is six fold degenerate in the thermodynamic limit.The second is a set of states that can be understood bytaking the stripe-like phase and rotating ◦ hexagonsalong a defect line that should wrap around the chosencluster. Classically, the energy per site of these two typesof states does not have first order hopping correctionsand it is given by E CDW II = V / V / O ( t /V , ) (see appendix A for details). The first nontrivial orderin perturbation theory depends on the chosen cluster;eventually the degeneracy between the six-fold degen-erate stripe-like phase and the phase with defects willbe lifted, making one of them more favorable at finite t . The question then becomes which type of groundstate does the system choose in the thermodynamic limit.Although we cannot attempt to fully answer this ques-tion some light can be shed through our iDMRG numer-ical data. We indeed find that both configurations, withand without defects are ground states with very simi-lar energy, even close to the transition to the semimetalphase. A sample of one of this charge ordered patternsis shown in Fig. 3(c) for V /t = 5 . and V /t = 3 . . Thisparticular sample state has one defect line that wrapstwice around the cluster (indicated by dashed lines); ro-tating the hexagons along these lines by ◦ in an anti-clockwise direction we recover the more symmetric stripe-like configuration in Fig. 2(e). As with the CMs state, wefind that its bond order follows directly from its chargeorder, suppressed for neighboring sites that have equalfilling.By initializing the algorithm with each inequivalent clas-sical ground state configuration it is possible to calculateto high precision each energy to then find, by compari-son, the lowest energy state at that particular { V , V } .We have determined that within the CDW II region thelowest energy is often achieved by a superposition of twoclassical ground states and their particle-hole conjugates.As expected from general arguments this state showsa degenerate low energy entanglement spectrum. Inspec-tion of the two point density-density correlation functionsindicates that the more favored superposition of classicalground states are those which favor a uniform chargedensity pattern that at the same time distinguishes thebonds around and along the cylinder geometry. F. Charge density wave III (CDW III)
At high V and between the previously discussedCDW II and CDW I states we find a third type ofcharge order with a twelve site unit cell labeled CDW III.A representative pattern of its charge and bond orderas obtained numerically with (2) and (3) respectivelyis shown in Fig. 3(d). As schematically represented inFig. 2(f), this state has three different occupations: half-filling and / ± (cid:37) . For the particular case of V /t = 9 . and V /t = 2 . we find that (cid:37) = 0 . . Its bond or-der is determined by the charge occupations similar tothe CMs and CDW II phases. As with the CDW IIphase, this phase escaped identification in the originalmean field and the first subsequent exact diagonal-ization studies due to its large unit cell.The emergence of the CDW III can be understood froma semiclassical point of view (see appendix A). At t = 0 ,the classical ground state of lowest energy for a clustercommensurate with this phase is either in the CDW IIor CDW I class, and the two phases are separated bythe line V = 4 V . Exactly at this line, however, a thirdclassical state is degenerate in energy with these but, un-like the previous two, is affected by first order quantumcorrections. This implies that around the line V ∼ V there is a finite strip of a new phase in the phase di-agram once finite t is included. We have checked thatthe two point correlation function calculated numerically from the iDMRG data and semiclassically agree qualita-tively, resulting in the charge distribution of the CDW IIIshown in Fig. 3(d). We find furthermore that these areconsistent with those reported in the latest exact diago-nalization study. IV. PHASE TRANSITIONS
In order to unravel the character of the correspondingphase transitions we now study two quantities that aresensitive to a phase change, the entanglement entropy S and the correlation length ξ . The entanglement entropyis defined by S = − Tr ρ L ln ρ L , (9)where ρ L = Tr R | ψ (cid:105)(cid:104) ψ | is the reduced density matrix ofthe left semi-infinite half L of the cylinder after tracingout the right half R . The correlation length ξ can becomputed from the resulting state of the iDMRG cal-culation according to Ref. 27. Both quantities are ex-pected to show a finite discontinuity when crossing afirst order transition while the correlation length divergeswith increasing bond dimension at a second order criticalpoint. In what follows, we focus on the lines labelledby cuts A–E in the phase diagram in Fig. 3.
A. Cut A: semimetal–CDW I
The first horizontal cut, labeled A addresses the char-acter of the phase transition between the semimetallicphase and the simplest charge density wave (CDW I) byfixing V = 0 . This phase transition has been previouslyaddressed with the quantum Monte-Carlo method and was determined to be of second order character. Thetransition point with divergent correlation length was de-termined to be at V /t = 1 . via finite size scaling. InFig. 4 we show the correlation length as a function of V for different values of the bond dimension χ .Firstly, for V (cid:46) . t we observe that the correlationlength drops as a function of V and diverges as thebond dimension χ is increased. This behaviour is ex-pected for a critical state such as the semimetal phase;the logarithmic divergence of entanglement of a metal-lic state requires an matrix product state with χ → ∞ .For V (cid:38) . t the correlation length continues to dropbut has no longer a significant dependence on χ . Thisis characteristic of a gapped phase such as the chargedensity wave. In the inset of Fig. 4, we plot the entangle-ment entropy as a function of the logarithm of the bonddimension for two values of V representative for the twophases. In the semimetal phase for V = 0 . , we observea finite entanglement scaling of S ∝ ln χ characteristicof a critical phase, whereas the entanglement entropy at V = 2 . in the gapped CDW I phase is independent ofthe bond dimension. FIG. 4. Correlation length at V = 0 , labeled cut A in Fig. 3probing the transition between the semimetal and CDW Iphases. The smooth behavior of ξ indicates a second or-der transition. The dashed grey line at V /t = 1 . showswhere the transition has been detected in quantum MonteCarlo simulations. Inset: Scaling of the entanglement en-tropy S with the bond dimension χ for two values of V inthe semimetal and CDW I phase. Unlike the CDW I phase,the semimetal phase presents the typical critical entanglementscaling S ∝ ln χ . The crossover between the two phases is smooth, sig-naling a second order phase transition, in agreement withquantum Monte Carlo studies.
With iDMRG it ishowever not possible to pin point the exact value of V where the transition happens via finite size scaling dueto the few cylinder sizes we have available, as discussedin section II. However, from Fig. 4 it is possible to de-fine a crossover region of . (cid:46) V (cid:46) . , consistent withthe quantum Monte Carlo data, where the transitionoccurs. B. Cut B: CMs–CDW II
In cut B we fix V = 3 . t and study the phase transitionbetween the CMs and the CDW II phases. The entangle-ment entropy is shown in Fig. 5 as a function of V . Thisquantity shows a clear discontinuity at V /t ≈ . sig-naling a direct first order phase transition between thesetwo phases. We observe the same behavior when lookingat the correlation length not presented here. Addition-ally, the energy of the ground state as a function of V shows a kink at the same critical V which further sup-ports the presence of a level crossing at that point. Asexpected for gapped phases, the entanglement entropyonly depends very weakly on the bond dimension, in theCDW II phase even a very low χ of 400 is sufficient tofaithfully represent the state.A first order phase transition is also expected from thestrong coupling approach. As detailed in appendix A,the energy per site of the CM phase in a × cluster is FIG. 5. Entanglement entropy at V /t = 3 . , labeled cut Bin Fig. 3. The finite discontinuity of S at V /t ≈ . clearlysignals a first oder transition between the CMs and CDW IIphases. E CMs = V / V / − . t for V /V < / , while theenergy per site of the CDW II phase is E CDW II = V / V / . At t = 0 , V = 0 both states have the same energy.At finite t , the CM phase has lower energy due to the firstorder quantum correction, which is absent for CDW II.Since E CMs grows faster with V than E CDW II , theremust be a crossing at a critical value of the interactionsignaling a first order phase transition into the CDW IIstate.
C. Cut C: Semimetal–CMs
The cut at V /t = 0 . , labeled cut C in Fig. 3, probesthe transition between the semimetal phase and the CMsphase. In Fig. 6 we present the entanglement entropy S as a function of V , the interaction that drives the phasetransition.As explained in Sec. III A, the entanglement entropyof the semimetal phase depends strongly on the bonddimension χ . At V /t ≈ . we observe a sharp transi-tion to a decaying entanglement entropy that does notdepend strongly on χ as V is increased. However, thechange in entanglement entropy is not as abrupt as in theCMs–CDW II case of the previous subsection as wouldbe expected from a level crossing in a first order transi-tion. We can see a kink in S , but for decreasing V itsmoothly approaches the value of the semimetal phase.In addition, the energy is a smooth function V whichtogether with the previously mentioned behaviour of S suggests that the transition is of second order. The cor-relation length not depicted here shows a similar smoothbehavior as in the transition between the semimetal andCDW I phases in Fig. 4 which provides further evidencein favor of a second order transition. FIG. 6. Entanglement entropy at V /t = 0 . , labeled cut Cin Fig. 3, probing the phase transition between the semimetaland CMs phases. Going towards the phase boundary fromhigh V , the entanglement entropy smoothly approaches thevalue in the semimetal indicating a second order transition. D. Cut D: CDW I–Kekulé–Semimetal–CDW II
The next cut we focus on is labeled cut D in Fig. 3 at V /t = 6 . The correlation length as a function of V isshown in Fig. 7. Starting from low V we first cross thephase boundary between the charge density wave CDW Iand the bond-ordered Kekulé phase. As expected for agapped phase, the correlation length deep in the CDW Iphase depends only weakly on the bond dimension. Aswe approach the phase transition, this behaviour changessince the many body gap is closing. The strongly in-creasing correlation length which peaks at V /t ≈ . suggests a second order transition between the CDW Iand Kekulé phases. Moreover, we can see the presence ofa remnant charge order on top of the bond modulation inthe Kekulé phase close to the critical point. This signalsa slow onset of charge ordering as a finite size effect aswe approach the phase boundary from the Kekulé sidewhich would be in contradiction to a level crossing forwhich the charge order should change suddenly. A sim-ilar second order transition was reported in the relatedantiferromagnetic Heisenberg model between a Néel anda plaquette phase. This transition was conjecturedto show deconfined quantum critical behavior.
The fact that the correlation has not fully convergedwith χ even in the center of the the Kekulé region in-dicates that this phase is strongly entangled. This canbe attributed to the fact that it is not a strong couplingphase in which the gap is determined by the interactionbut rather by the hopping energy scale. As V is increasedwe cross the boundary between the Kekulé and thesemimetal phase at V /t ≈ . . From the data presentedin Fig. 7 it is not possible to reach a conclusion about thenature of the transition. However, from the Landau the-ory perspective the semimetal–Kekulé transition has tobe first order for the following reason. The Kekulé bond FIG. 7. Correlation length at V /t = 6 , labeled cut D inFig. 3. In this plot, several phase transitions are made vis-ible. With increasing V , the CDW I phase turns into theKekulé phase which itself goes into the semimetal phase. Thetransition from semimetal to CDW II can be identified bythe peak in the correlation length. Inset: The vicinity of theKekulé–semimetal transition showing the discontinuity of theentanglement entropy as a function of V order is in the twofold representation E (cid:48) = ( E (cid:48) , E (cid:48) ) and the lowest order non-trivial term in the Landau freeenergy is the scalar ( E (cid:48) ) − E (cid:48) ( E (cid:48) ) which is of thirdorder, therefore signaling a first order transition. A firstorder transition is also consistent with the entanglemententropy which displays a small discontinuity when cross-ing the phase boundary (see inset of Fig. 7).In the semimetal phase, we observe an astonishinglysmall correlation length for a critical state. As explainedin Sec. III A the hopping strengths around and along thecylinder are strongly renormalized. This shifts the Diracnodes away form the K and K (cid:48) points thus leading to aneffectively gapped state for a finite cylinder circumferencewith small correlation length.Finally, at V /t ≈ . the semimetal turns into theCDW II. The diverging correlation length on the CDW IIside of the transition points towards a second order phasetransition. Moreover, the strong dependence of ξ on thebond dimension suggests that the ground state is becom-ing critical as we approach the transition point, espe-cially compared to the behaviour of ξ deep in the CDW IIphase. On the other hand, from the semimetal side thecorrelation length is discontinuous, which could signal afirst order transition. However, the latter observation hasto be taken with care since the semimetal phase, beinga critical state, cannot be represented faithfully with afinite bond dimension. E. Cut E: CDW I–Kekulé–CDW III–CDW II
The last cut we address is labeled cut E in Fig. 3 andfixes V /t = 9 . . The entanglement entropy is presented FIG. 8. Entanglement entropy at V /t = 9 . , labeled cut Ein Fig. 3. The smooth behavior of S at the phase bound-ary between CDW I and Kekulé phase indicates a second or-der transition in accordance with the data of the correlationlength in Fig. 7. At V /t ≈ . , the entanglement entropy isdiscontinous, which clearly indicates a first order transitionbetween the CDW III and CDW II phases. in Fig. 8. With increasing V we cross the second orderphase transition from the CDW I to the Kekulé phaseat V /t ≈ . discussed above. The entanglement datashown in Fig. 8 supports this statement since S smoothlyincreases in the CDW I phase as we approach the phaseboundary.At V /t ≈ . the system undergoes the transition tothe CDW III phase. Unfortunately, a conclusive state-ment about the order of the critical point cannot bedrawn from the present data; especially since the en-tanglement entropy in the CDW III phase still shows asignificant dependence on the bond dimension.The last phase boundary in this cut lies between theCDW III and CDW II phases at V /t ≈ . . Here, wecan see a clear discontinuity of the entanglement entropy,similar to that in cut B (Fig. 5). Moreover, we see akink in the energy when crossing the phase boundary.Together these features indicate a first order phase tran-sition. V. DISCUSSION AND CONCLUSIONS
The spontaneous emergence of the Haldane Chern in-sulator phase due to interactions in the half-filled spine-less honeycomb lattice model has been questioned by dif-ferent methods since its proposal. In this workwe have used the iDMRG method and have mapped outthe phase diagram of this model. The algorithm allowsfor quantum fluctuations to play a role but minimizesthe finite size effects due to its intrinsic infinite cylindergeometry.One of our main results is the absence of the HaldaneChern insulator state in this model. Even when initial- izing the iDMRG calculations with a time reversal sym-metry broken chiral wave function, the state did not re-main stable and evolved into the respective competingphases upon applying the algorithm. It seems thereforethat quantum fluctuations indeed jeopardize emergenceof the CI phase and work against the naive mean field ex-pectations. Although there is no reason to believe thatthe situation is different for the π -flux model, which alsohosts a mean field Chern insulating phase that is absentin exact diagonalization, there is more hope for othermodels with quadratic band point touchings. In particu-lar, the interacting kagome lattice hosts a Chern insulatorstate within mean field theory when the filling is tuned tothe quadratic band point touching between the flat bandand one of the two dispersive bands. Its presence ispredicted within the renormalization group approach which guarantees its robustness at sufficiently low inter-action strength. Adequate substrate engineering can alsolead to such quadratic band point touching, potentiallyfavoring the Chern insulator state. In addition, we have theoretically accounted for two novelcompeting phases, the CDW II and CDW III. Their po-tentially large unit cells turns their identification withinexact diagonalization or mean field studies challenging.On the one hand the CDW II phase stems from a de-generate strong coupling phase that has no first ordercorrection in powers of t/V , . Although the degener-ate manifold includes a distinguishable stripe-like phase,our iDMRG calculations find that a superposition statebetween classical ground states has lower energy. Theparticular superposition that is realized seems to be de-termined by the cylinder geometry where bond strengthalong and around the cylinder are inequivalent. TheCDW III on the other hand, occurs in a finite regionaround the classical transition between the CDW I andII states defined by the line V = 4 V , and has a cleartwelve site unit cell. In this region, the first order cor-rections t/V , lift the classical degeneracy and stabilizeCDW III state. Both the semiclassical limit and our nu-merical iDMRG data show consistent charge structurefactors confirming the semiclassical nature of the state.The CDW III phase presents as well a distinctive bondordering that was understood from its charge order.From the evidence presented in this work we haveestablished that the CMs state occurs from first orderperturbation theory in powers of t/V of an otherwiseclassically degenerate ground state. This part of ourwork establishes the existence and robustness of thissingle particle charge order beyond numerical approxi-mations. Our semiclassical analysis, corroborated by theiDMRG numerical evidence, provides an explanation ofwhy this phase has a many-body gap of the order of thehopping strength t/V , as was previously noticed, rather than determined by V . Both the semiclassicaltreatment and the numerical data we obtain establishthat the CMs state has a finite sublattice imbalance, afeature that was still under dispute. Moreover, wehave reported a characteristic bond order not addressed0in previous studies. As in the CDW II and CDW IIIcases such bond order is determined by the CMs chargeorder pattern. Our results show that the CMs is stablefor small V and finite t . Increasing the former orreducing the latter leads to the CDW II ground statethrough a first order phase transition.Furthermore, we have established the existence of theKekulé bond order which occurs in a region of the phasediagram that appears to be smaller than that predictedby mean field and exact diagonalization.Finally, we have used the fact that the iDMRG methodtreats an infinite system as opposed to a finite clusterto analyse the different phase transitions using inparticular the entanglement entropy S and correlationlength ξ . The conclusions drawn from analyzing thesequantities have been complemented by the continuous ordiscontinuous dependence of the ground state energy onthe { V , V } interactions. However, conclusions aboutthe order of phase transitions from and to the semimetalphase should be taken with care due to the gaplessnature of the semimetal ground state.To conclude we have established the phase diagram ofspinless fermions hopping on the half-filled honeycomblattice with nearest- and next-to-nearest neighbor inter-actions and characterized the phase transitions amongthe different phases with iDMRG. Our results providesolid evidence that Chern insulating phases are far moreelusive than previously thought and so alternative routesare necessary to drive these kind of topological statesfrom strong electronic correlations in general. VI. ACKNOWLEDGEMENTS
We thank A. Läuchli for discussions, useful insightsregarding the nature of the CDW II and CDW III phasesand sharing consistent exact diagonalization results priorto publication. Appendix A: Strong coupling perturbation theory
In this appendix we discuss the exact diagonalizationmethod to determine the ground state in the strong cou-pling limit t (cid:28) V , . Consider splitting the Hamiltonianin Eq. (1) into H = H V + H t with H t = − t (cid:88) (cid:104) i,j (cid:105) ( c † i c j + h . c . ) (A1) H V = V (cid:88) (cid:104) i,j (cid:105) n i n j + V (cid:88) (cid:104)(cid:104) i,j (cid:105)(cid:105) n i n j . (A2)In the strong coupling limit, we can obtain the groundstate of H by diagonalizing H V first and considering H t as a perturbation. The eigenstates of H V H V ψ n,m = E n ψ n,m , (A3) where m accounts for degeneracies, are simple to computebecause charge is conserved at every site in the absence ofhopping. H V is thus already diagonal in the occupationbasis, i.e. | ψ n,m (cid:105) = (cid:89) C n,mi =1 c † i | (cid:105) , (A4)with C n,mi = 0 , the occupation coefficients for the n, m eigenstate. The classical ground state manifold isspanned by (cid:12)(cid:12) ψ ,m (cid:11) , with m = 1 , . . . , M .For small t/V , , we can disregard the classical stateswith n > , and project the H t Hamiltonian into theclassical ground state manifold. Dropping the label n = 0 from now on ( ˜ H t ) m m = (cid:104) ψ m | H t | ψ m (cid:105) . (A5)We can now diagonalize ˜ H t and select the eigenstates oflowest energy (cid:15) ˜ H t v α = (cid:15) v α , (A6)where α = 1 , . . . , D and D is the true ground state de-generacy. The final ground state of H is spanned by | GS (cid:105) α = M (cid:88) m =1 v mα | ψ m (cid:105) . (A7)Obtaining the coefficients v mα is relatively simple becausethe effective size of the Hilbert space M is much smallerthan the full size of the Hilbert space of H . To distinguishdifferent phases, we recall that in finite size exact diag-onalization there is no spontaneous symmetry breaking,because all states related by symmetry are degenerateand will be present in the ground state manifold. In thesimpler case when first order quantum corrections arezero, a phase can be characterized by inspection of theclassical charge patterns of every eigenstate. With quan-tum corrections, however, the ground state can only becharacterized by computing correlation functions evalu-ated in the ground state, from which order parameterscan be obtained. Correlation functions for charge orderparameters can be expressed in general as ρ ij... = tr (cid:104) GS | c † i c i c † j c j . . . | GS (cid:105) , (A8)which is given explicitly by ρ ij... = (cid:88) m,α ( v αm ) ∗ v αm C mi C mj . . . (A9)The behavior of correlation functions can then be usedto identify the different phases. We have used thismethod to diagonalize three different clusters: two 12site clusters with the periodicities of the CDW II andCDW III phases, and a cluster of × unit cells or 18sites, i.e. the L y = 6 version of the cluster in Fig. 1. Theresults can be summarized as follows.1In the case of the CDW II cluster, we find only twopossible classical ground states with degeneracies 8 and2 for generic values of V /V . The projection of H t inthese subspaces is zero for both cases, which means thatboth are stable strong coupling phases. The first can beidentified with the classical version of the CDW II state,where the two values of the occupancies are 0 or 1. The8 states correspond to the 8 equivalent ways to ensemblethe CDW II pattern in the given unit cell: two of themcorrespond to the stripe-like phase shown schematicallyin Fig. 3(c) while the other six correspond to those stateswith defect lines determined by rotated hexagons, as de-scribed in sec. III E. The second is the CDW I state,where the two states have a single sublattice (A or B)fully occupied. The energy per site of these states isgiven by E CDW II = 12 V + 14 V V V < (A10) E CDW I = 3 V < V V (A11)In the case of the CDW III cluster, for generic values of V /V we find two possible classical ground states withdegeneracies 2 and 2. The first correspond to the stripe-like phase mentioned before, which is also commensu-rate with this cluster, and is therefore assigned the labelCDW II. The second corresponds to CDW I. The ener-gies per site remain the same as those in Eqs. A10-A11.In addition, in this cluster there is another classical statewith degeneracy 18 and energy E ∗ = V + V whichrequires consideration. This energy is always higher than E CDW II or E CDW I , except at the special point V = 4 V where the three cross, E CDW II = E CDW I = E ∗ . Thereason why this state must be considered is that quan-tum corrections split it to first order in t , thus loweringits energy. Since neither CDW II or CDW I is affected by t to first order, there is a small region around V = 4 V where the energy of this state is lowest. Inspection of thephase diagram in Fig. 3 suggests that this is the CDW IIIphase. The ground state of this phase has degeneracy oneand energy E CDW III = 5 V V − . t V V ∼ (A12)To confirm the nature of this state, we have computedthe Fourier transform of the correlation function C ( q ) = (cid:80) ij e i ( x i − x j ) q ρ ij , and confirmed that it compares favor-ably with the Fourier transform of the charge pattern ofCDW III shown in Fig. 3(d).In the case of the × cluster, at t = 0 we find fourpossible classical ground states as a function of increas-ing V , with degeneracies 234, 108, 36, 2. The projectionof H t into these classical ground states is finite only forthe first two, and the pattern of degeneracies at low en-ergies becomes (4 , , , . . . ) for both. This low energy pattern is the same as the one that is obtained for theCMs phase in exact diagonalization, pointing to thefact that these two states represent the CMs phase (fur-ther evidence to support this claim is also shown below.).The third state is an intermediate state with the sameenergy as the CDW III classical states, but which is notcorrected by quantum fluctuations at any V /V . Thefourth is again the CDW I. The energies per site of thethree phases in the × cluster are given by E CMs = (cid:40) V + V − . t V V < V + V − . t < V V < (A13) E ∗ = 5 V V < V V < (A14) E CDW I = 3 V < V V (A15)As explained in the main text, the energy per site ofthe CDW II phase E CDW II is lower than the classicalenergies of any state of the × cluster. Since theintermediate state does not have quantum corrections,its energy satisfies E int > E CDW II and it is never re-alized. The CMs state, however, has E CMs < E
CDW II for small enough V once the quantum corrections are in-cluded. Further evidence that this is the CMs state canbe obtained by computing correlation functions. Fouriertransforming the real space indices ij . . . and taking thecombinations defined in the main text for the representa-tions B and G , we can compute any correlation functionof order parameters (cid:104) GS | O O . . . | GS (cid:105) with O = G, B .The simplest scalars that signal the presence of chargeorder are simply the two point functions S = B , (A16) S = | G | + | G | . (A17)We have computed these correlation functions in theground states of the × cluster with V < V andfinite t and found that they are both finite. However,this is not enough information to distinguish the CMsphase. In order to detect its particular order we havealso computed the correlation functions of higher orderfor this ground state S = Im( G − G G ) , (A18) S = B Im( G G ) , (A19) S = [Im( G G )] , (A20) S = B Re( G − G G ) , (A21)and found that only S is finite, therefore confirming thepresence of the CMs phase with the strong coupling ap-proach for V < V . As for the CDW III we have alsocompared the semicalssical and numerical structure fac-tors and found good agreement, supporting the semiclas-sical interpretation of the CMs phase.2 F. Pollmann, S. Mukerjee, A. M. Turner, and J. E. Moore,Phys. Rev. Lett. , 255701 (2009). M. Z. Hasan and C. L. Kane, Rev. Mod. Phys. , 3045(2010). X.-L. Qi and S.-C. Zhang, Rev. Mod. Phys. , 1057(2011). C. Nayak, S. H. Simon, A. Stern, M. Freedman, andS. Das Sarma, Rev. Mod. Phys. , 1083 (2008). G. E. Volovik,
The Universe in Helium Droplet (Claren-don, Oxford, UK, 2003). T. Senthil, arXiv:1405.4015 (2014). F. D. M. Haldane, Phys. Rev. Lett. , 2015 (1988). C.-Z. Chang, J. Zhang, X. Feng, J. Shen, Z. Zhang,M. Guo, K. Li, Y. Ou, P. Wei, L.-L. Wang, Z.-Q. Ji,Y. Feng, S. Ji, X. Chen, J. Jia, X. Dai, Z. Fang, S.-C.Zhang, K. He, Y. Wang, L. Lu, X.-C. Ma, and Q.-K. Xue,Science , 167 (2013). S. Raghu, X.-L. Qi, C. Honerkamp, and S.-C. Zhang, Phys.Rev. Lett. , 156401 (2008). C. Weeks and M. Franz, Phys. Rev. B , 085105 (2010). T. Pereg-Barnea and G. Refael, Phys. Rev. B , 075127(2012). A. G. Grushin, E. V. Castro, A. Cortijo, F. de Juan,M. A. H. Vozmediano, and B. Valenzuela, Phys. Rev. B , 085136 (2013). E. V. Castro, A. G. Grushin, B. Valenzuela, M. A. H.Vozmediano, A. Cortijo, and F. de Juan, Phys. Rev. Lett. , 106402 (2011). Y. Jia, H. Guo, Z. Chen, S.-Q. Shen, and S. Feng, Phys.Rev. B , 075101 (2013). J. Wen, A. Rüegg, C.-C. J. Wang, and G. A. Fiete, Phys.Rev. B , 075125 (2010). L. Wang, P. Corboz, and M. Troyer, New Journal ofPhysics , 103008 (2014). Z.-X. Li, Y.-F. Jiang, and H. Yao, New Journal of Physics , 085003 (2015). N. A. García-Martínez, A. G. Grushin, T. Neupert,B. Valenzuela, and E. V. Castro, Phys. Rev. B , 245123(2013). M. Daghofer and M. Hohenadler, Phys. Rev. B , 035103(2014). T. Ðurić, N. Chancellor, and I. F. Herbut, Phys. Rev. B , 165123 (2014). C. Chamon, Phys. Rev. B , 2806 (2000). C.-Y. Hou, C. Chamon, and C. Mudry, Phys. Rev. Lett. , 186809 (2007). B. Roy and I. F. Herbut, Phys. Rev. B , 035429 (2010). B. Roy, V. Juricic, and I. F. Herbut, Phys. Rev. B , 041401 (2013). I. P. McCulloch, arXiv:0804.2509 (2008). S. R. White, Phys. Rev. Lett. , 2863 (1992). J. A. Kjall, M. P. Zaletel, R. S. K. Mong, J. H. Bardarson,and F. Pollmann, Phys. Rev. B , 235106 (2013). M. P. Zaletel, R. S. K. Mong, and F. Pollmann, Phys.Rev. Lett. , 236801 (2013). Y.-C. He, D. N. Sheng, and Y. Chen, Phys. Rev. B ,075110 (2014). Y.-C. He, N. Sheng, D., and Y. Chen, Phys. Rev. Lett. , 137202 (2014). H.-C. Jiang, Z. Wang, and L. Balents, Nature Physics ,902 (2012), 1205.4289. L. Cincio and G. Vidal, Phys. Rev. Lett. , 067208(2013). Z. Liu, D. L. Kovrizhin, and E. J. Bergholtz, Phys. Rev.B , 081106 (2013). A. G. Grushin, J. Motruk, M. P. Zaletel, and F. Pollmann,Phys. Rev. B , 035136 (2015). A. H. Castro Neto, F. Guinea, N. M. R. Peres, K. S.Novoselov, and A. K. Geim, Rev. Mod. Phys. , 109(2009). R. Shankar, Rev. Mod. Phys. , 129 (1994). V. N. Kotov, B. Uchoa, V. M. Pereira, F. Guinea, andA. H. Castro Neto, Rev. Mod. Phys. , 1067 (2012). B. Amorim, A. Cortijo, F. de Juan, A. G. Grushin,F. Guinea, A. Gutiérrez-Rubio, H. Ochoa, V. Parente,R. Roldán, P. San-José, J. Schiefele, M. Sturla, andM. A. H. Vozmediano, arXiv:1503.00747 (2015). D. M. Basko, Phys. Rev. B , 125418 (2008). F. de Juan, Phys. Rev. B , 125419 (2013). S. Capponi and A. M. Läuchli, Phys. Rev. B , 085146(2015). Z.-X. Li, Y.-F. Jiang, and H. Yao, arXiv:1408.2269 (2014). R. Ganesh, J. van den Brink, and S. Nishimoto, Phys.Rev. Lett. , 127203 (2013). A. F. Albuquerque, D. Schwandt, B. Hetényi, S. Capponi,M. Mambrini, and A. M. Läuchli, Phys. Rev. B , 024406(2011). T. Senthil, A. Vishwanath, L. Balents, S. Sachdev, andM. P. A. Fisher, Science , 1490 (2004). We thank Hong Yao for bringing this point to our atten-tion. K. Sun and E. Fradkin, Phys. Rev. B , 245122 (2008). K. Sun, H. Yao, E. Fradkin, and S. A. Kivelson, Phys.Rev. Lett. , 046811 (2009).49