Phase diagram of interacting spinless fermions on the honeycomb lattice: A comprehensive exact diagonalization study
PPhase diagram of interacting spinless fermions on the honeycomb lattice:a comprehensive exact diagonalization study
Sylvain Capponi ∗ and Andreas M. Läuchli † Laboratoire de Physique Théorique, CNRS UMR 5152, Université Paul Sabatier, F-31062 Toulouse, France Institut für Theoretische Physik, Universität Innsbruck, A-6020 Innsbruck, Austria (Dated: September 25, 2018)We investigate the phase diagram of spinless fermions with nearest and next-nearest neighbour density-density interactions on the honeycomb lattice at half-filling. Using Exact Diagonalization techniques of the fullHamiltonian and constrained subspaces, combined with a careful choice of finite-size clusters, we determinethe different charge orderings that occur for large interactions. In this regime we find a two-sublattice Néel-likestate, a charge modulated state with a tripling of the unit cell, a zig-zag phase and a novel charge ordered stateswith a 12 site unit cells we call Néel domain wall crystal, as well as a region of phase separation for attractiveinteractions. A sizeable region of the phase diagram is classically degenerate, but it remains unclear whether anorder-by-disorder mechanism will lift the degeneracy. For intermediate repulsion we find evidence for a Kekuléor plaquette bond-order wave phase. We also investigate the possibility of a spontaneous Chern insulator phase(dubbed topological Mott insulator), as previously put forward by several mean-field studies. Although we areunable to detect convincing evidence for this phase based on energy spectra and order parameters, we find anenhancement of current-current correlations with the expected spatial structure compared to the non-interactingsituation. While for the studied t − V − V model the phase transition to the putative topological Mott insulatoris preempted by the phase transitions to the various ordered states, our findings might hint at the possibilityfor a topological Mott insulator in an enlarged Hamiltonian parameter space, where the competing phases aresuppressed. PACS numbers: 71.10.Fd, 71.27.+a, 71.30.+h, 75.40.Mg
I. INTRODUCTION
The seminal discoveries of quantum Hall phases have pro-vided the first examples of topological phases of matter, i.e.phases which cannot be adiabatically connected to conven-tional ones. More recently, several other topological insula-tors or superconductors have been theoretically predicted andexperimentally observed , giving rise to an intense activityof the community. For such materials, a very thorough un-derstanding has been achieved thanks to the fact that a non-interacting starting point is sufficient: indeed, topological in-sulators can be obtained from band theory in the presence ofstrong spin-orbit coupling.On the other hand, given the variety of electronic phasesthat are also generated by strong interactions, a highly topicalquestion is to investigate the appearance of topological insula-tors due to correlations. This new route would not require in-trinsic spin-orbit coupling, and thus could be advantageous tomany applications. Nevertheless, we face the well-known dif-ficulty of strongly correlated systems for which we lack unbi-ased tools, and often have to resort to some hypothesis. In thiscontext, a very promising result has been obtained by Raghu et al. who have shown the emergence of quantum anoma-lous (spin) Hall (QAH) states on the honeycomb lattice withstrong density-density repulsions between next-nearest neigh-bors. These phases can be characterized by the existence ofspontaneously appearing charge (respectively spin) currentsfor spinless (respectively spinful) fermions. This finding hasalso been reproduced in the spinless case using more sys-tematic mean-field approaches . Nevertheless, in the latterstudy, the QAH extension in the phase diagram has shrunksubstantially due to the stability of charge modulation with larger unit cell.The possibility to generate topological phases with longer-range interactions has triggered a large theoretical activityto investigate its experimental feasibility using metallic sub-strate , cold-atomic gases on optical lattices or RKKY inter-action among many others.In the light of several recent approximate and numericalstudies with conflicting conclusions regarding the pres-ence of the putative topological Mott state and competing or-dered phases, we feel that a comprehensive state-of-the-artexact diagonalization study could shed some light on these is-sues. Based on our simulations, we will provide numerical ev-idence that quantum fluctuations at small V > are responsi-ble for an enhanced response similar to the QAH state, but un-fortunately our numerical data (obtained on finite clusters upto N = 42 sites) indicates that these short-range fluctuationspresumably do not turn into a true long-range ordered phase.In addition, we will provide a detailed analysis of the chargeand bond orderings that occur at large interaction strength, al-lowing us to sketch a global phase diagram of the minimalmodel studied so far. In Sec. II, we detail the model and meth-ods we use. We then discuss the large interaction (classical)limit in Sec. III. Our numerical results are given in Sec. IV,and we conclude in Sec. V. II. MODEL, METHOD AND PHASE DIAGRAMA. Hamiltonian
Following the proposal by Raghu et al. a r X i v : . [ c ond - m a t . s t r- e l ] A ug V t V MX K ( a ) ( b ) ⇤ X ⇤ FIG. 1. (Color online) (a) Illustration of the honeycomb lattice with V , V interactions and the hopping t . (b) First (solid line) and second(dashed line) Brillouin zone of the honeycomb lattice including thelocation of a few special points in the Brillouin zone. V interaction, we will consider the following Hamiltonian : H = − t (cid:88) (cid:104) ij (cid:105) ( c † i c j + h.c. ) (1) + V (cid:88) (cid:104) ij (cid:105) ( n i − / n j − / V (cid:88) (cid:104)(cid:104) ij (cid:105)(cid:105) ( n i − / n j − / depicted in Fig. 1(a), where c i and c † i are the spinlessfermionic operators, t = 1 is the nearest-neighbor hoppingamplitude, V and V are the density-density repulsion or at-traction strengths respectively on nearest- and next-nearestneighbors. Note that we will focus on the half-filled casewhere thanks to particle-hole symmetry, the chemical poten-tial is known to be zero exactly.Despite the particle-hole symmetry, it seems that quantumMonte-Carlo (QMC) approaches exhibit a sign problem thatprohibits them. In fact, only in the simpler case V = 0 whichexhibits a direct phase transition from semi-metallic (SM)phase to charge density wave (CDW) as a function of V /t ,a very interesting recent proposal was made allowing to refor-mulate the problem without sign-problem , thus amenable toaccurate, unbiased QMC simulations . Quite remarkably, thecase V ≥ and V ≤ can also be studied with QMC with-out a minus-sign problem using a Majorana representation .In the absence of a QMC algorithm for the frustrated case, weuse Exact Diagonalization (ED) to get unbiased numerical re-sults. Of course, one is limited in the sizes available but, sincewe are investigating in particular the topological QAH phasewhich is a translationally invariant ( q = 0 ) instability (seebelow), we can in principle use any finite-size lattice, eventhough some of them do not have the full space-group sym-metry. On the other hand, when looking at competing chargeinstabilities, it is crucial to consider only clusters compatiblewith such orderings and we will devote some discussion onthis topic in the next subsection. B. Finite lattices
Since we plan to provide a systematic study on several clus-ters, we have implemented lattices sizes ranging from 12 to 42sites, with various shapes and spatial symmetries. In particu-lar, these clusters may or may not have some reflection or ro-tation symmetries (they all have inversion symmetry). More-over, some of these lattices possess the ± K points (see theBrillouin zone depicted in Fig. 1(b)) where the tight-bindingdispersion exhibits two Dirac cones, leading to a 6-fold de-generacy of the half-filled free fermion ground-state (GS) sothat correlations require some care for these clusters. We referto Appendix A for further details on these finite-size clusters.Our large choice of clusters provides a substantial step for-ward with respect to other recent ED studies, where resultswere obtained solely on N = 18 and N = 24 in Ref. 10 or on N = 24 and N = 30 in Ref. 11.In order to illustrate the need for a systematic study of clus-ter geometries we present some data for the ground state en-ergy per site. For instance, at vanishing V /t , as shown inFig. 2(a), clusters with the K points yield the lowest energiesfor large V /t , i.e. are compatible with the appearing order (tobe discussed later). V -2.5-2-1.5-1 e V -2-1.5-1 e (a) V = 0 (b) V / t = 4 FIG. 2. (Color online) GS energy per site e vs V /t for (a) V = 0 and (b) V /t = 4 on various clusters. We refer to the Appendix Afor further details on these finite-size clusters. Lines are guide to theeyes. At V /t = 4 and large V /t , the situation becomes differentas now clusters involving M point(s) seem to have the lowestenergies, see Fig. 2(b). This is already a clear indication ofa phase transition between CDW (known to be stable at small V /t and large V /t ) and at least two other phases that we willcharacterise later.Clearly, the choice of clusters is important depending onthe parameters and the kind of instability. It was also pointedout recently that the choice of boundary conditions might becrucial too. By studying the N = 18 cluster with open bound-ary conditions (OBC), as opposed to standard periodic bound-ary conditions (PBC), the authors of Ref. 12 have found twolevel crossings for V /t (cid:39) . and V /t (cid:39) . respectively(for V = 0 ). This was taken as an indication for a stable QAHphase in this region. We have checked that in this case, theselevel crossings correspond to two low-energy states havingopposite parities with respect to inversion symmetry. How-ever, these crossings disappear when one considers the nextlargest cluster N = 32 with OBC, so that it probably corre-sponds to a short-range feature only. We generally think thatusing PBC is more appropriate to minimize finite-size effects,so this will be the case here.Last, we would like to recall some features of the semi-metallic (SM) phase that exists in the absence of interactions.Since there is a vanishing density of states at the Fermi level,one needs a finite strength of a short-ranged interaction to trig-ger an instability , so that SM phase must have a finite ex-tension in the phase diagram.Regarding possible gap opening mechanism of the SMphase, Refs 18 and 19 have listed all explicit (i.e. external)weak-coupling perturbations which can open a gap. In thespinless case considered here, the three particle-hole relatedgaps are: i) the Néel-like charge density wave, which breaksthe A-B sublattice symmetry, ii) the Kekulé distortion pattern,which breaks translation symmetry by adopting a tripling ofthe unit-cell of modulated bond-strengths (this order parame-ter has a real and an imaginary part, thus corresponds to twomasses), and iii) the integer quantum Hall mass , inducedby breaking the time-reversal invariance and parity symmetryupon adding complex Peierls phases on next-nearest neighborhoppings, without enlarging the size of the unit-cell. In ad-dition to the particle-hole gaps, there is also the possibility toopen gaps by the addition of superconducting order parame-ters , we will however not address these instabilities inthis work.The model Hamiltonian (1) considered here features all theusual symmetries. If the semi-metallic phase is gapped outby interactions, then the gap opening has to happen throughthe interactions by spontaneously breaking some of the sym-metries. The case i) quoted before is a well known instabil-ity, since the Néel CDW state is an obvious strong-couplingground state at large V /t . The other instabilities ii) and iii)currently lack a strong coupling picture, and need to be con-firmed by numerical simulations. We note however that allthree particle-hole instabilities have been reported in mean-field studies. C. Overview of the phase diagram
We start by drawing the global phase diagram that summa-rizes our main findings, see Fig. 3. Its main features are theexistence of several types of charge or bond ordering for inter-mediate to large V and/or V interactions: Néel CDW, chargemodulation (CM), zigzag (ZZ) phase, Néel domain wall crys-tal (NDWC), and plaquette/Kekulé phase (P-K) that we willclarify later. The large orange region (ST*) in the upper rightpart of the phase diagram features a degeneracy at the semi-classical level, and it is presently unclear whether and how anorder-by-disorder mechanism will lift the degeneracy. Whilesome of these phases (CM, P-K, CDW) had been predicted us-ing mean-field studies and confirmed numerically in some re-gions , the others (including NDWC and ST* phase for re- pulsive interactions and the ZZ phase for attractive V ) had notbeen advocated before. Note already that the plaquette/Kekulé(P-K) phase only exists in some bounded region for interme-diate ( V , V ) values, and does not extend to strong coupling.There is also a large region of phase separation, mostly forstrong attractive interactions, in agreement with the results ofRef. 23 for V = 0 . Possibly superconductivity is present inparts of the attractive region of the phase diagram, but we didnot focus on this instability here.Last but not least, we do not have any convincing evidencefor the stability of the QAH phase, as found in similar recentnumerical studies but in contradiction with another nu-merical study using open boundary conditions and entangled-plaquette ansatz. ST*Néel CDWCM P - K Semimetal -4 -2 0 2 4 6 8 10 12 14 V / t -2-101234567 V / t ZZ Phase Separation
N D W C V /t V / t FIG. 3. (Color online) Phase diagram in the ( V /t, V /t ) parameterspace obtained from several exact diagonalization techniques (seetext). Dashed lines represent the classical transition lines, see Fig. 4.The semi-metal, which is the ground-state for non-interacting spin-less fermions, has a finite extension in the phase diagram because ofits vanishing density of states at the Fermi level. We will argue in theremainder of this article that several other phases can be stabilised forintermediate and/or large interactions: Néel CDW, plaquette/Kekulé(P-K), Néel domain wall crystal (NDWC), zigzag (ZZ) phase, andcharge modulation (CM). The region (ST*) is degenerate at the semi-classical level, and it is presently unclear whether and how an order-by-disorder mechanism will lift the degeneracy. Note also the largeregion of phase separation mostly in the attractive quadrant. Filledsymbols correspond to numerical evidence (using level spectroscopyor measurements of correlations, see Sec. IV) obtained mostly on a N = 24 cluster which contains the most important points in its Bril-louin zone and features the full lattice point group symmetry of thehoneycomb lattice. Star symbols denote likely first order transitions,witnessed by level crossings on the same cluster. Our numerical re-sults do not support any region of topological QAH phase. We will now turn to the presentation of various numericaldata and considerations that we have used to come up withthis global phase diagram.
Néel CDWPhase Separation V V t = 0( a ) V /V = 4 V /V = Stripy* Zigzag*
Néel CDWStripy* Zigzag* Phase Separation V V Charge Modulation (CM) (K point) ( b ) Néel Domain Wall Crystal (X point) t/ ( V , V ) V /V = 4 V /V = FIG. 4. (Color online) Classical phase diagram (a) and qualitative first order in t/ ( V , V ) phase diagram (b). The "Zigzag*" and the "Stripy*"phases feature a nontrivial ground state degeneracy (hence the "*" suffix in "Zigzag*" and "Stripy*"). The three points V = 1 , V = ± and V = 1 , V = 0 feature an extensive ground state degeneracy at the classical level. Upon including the first order correction due to thefinite hopping t , two of these points spawn new phases. The V = 1 , V = 0 point develops into the charge modulation (CM) phase, whilethe V = 1 , V = +4 point broadens into a novel "Néel domain wall crystal" (NDWC), which is sketched in Fig. 12. All regions and lines in(b) beyond the CM and NDWC phases have no first order (in t ) quantum corrections. Note that the radial direction in the two panels has nophysical relevance, i.e. it is not parametrizing t/V , . III. LIMIT OF LARGE INTERACTIONS (ISING LIMIT)
The Hamiltonian ( ) reduces to a free fermion problem inthe absence of interactions, yielding a semi-metallic state withtwo Dirac points at half filling. In the opposite limit t = 0 , themodel reduces to a (generally frustrated) classical Ising modelwith competing antiferromagnetic (for V , > ) or ferromag-netic (for V , < ) nearest and next-nearest neighbor inter-actions. Since the phase diagram at finite but small hopping isexpected to be influenced by this (frustrated) Ising limit, wefind it worthwhile to explore the classical phase diagram first.Our results are summarized in the phase diagrams sketchedin Fig. 4 where we have considered the pure classical case t = 0 (left panel), as well as the perturbative effect of asmall t (right panel). We will discuss these two cases andtheir phases in the following. We will exchangeably use the V , V notation or the θ notation: V = cos θ , V = sin θ , with θ ∈ [0 , π ) . V / V / V / V V V FIG. 5. (Color online) Sketch of a Hamiltonian unit which rendersthe classical V − V problem "frustration free". A. Classical limit: t = 0 First we point out that this classical Hamiltonian can bewritten in a "frustration free" form, when the Hamiltonian isdecomposed into parts acting on triangles spanned by three V bonds and their central site (this object can alternativelybe considered as a tripod, i.e. a site plus its three nearest V neighbors), whose three V couplings are counted only withhalf their strength (this is necessary because each such bondbelongs to two different V triangles, while the V bonds be-long to a unique triangle), see Fig. 5. In this decomposi-tion each such tripod acts on 4 sites, and depending on θ theHamiltonian on a tripod has a different number of local groundstates. We have numerically verified that all ground states ofthe full problem are simultaneous ground states of each tripodterm, thus the characterization "frustration free".By finite cluster ground state enumeration techniques (withand without exploiting the frustration free property) we haveestablished the classical phase diagram displayed in the leftpanel of Fig. 4. For V = 0 , V = 1 ( θ = π/ ), thesystem decouples into two interpenetrating triangular latticeswith antiferromagnetic interactions. The ground state man-ifold is thus spanned by the product of the ground spacesof each triangular sublattice, and is macroscopically degen-erate (i.e. featuring an extensive ground state entropy). For arctan(1 / < θ < π/ the ground state manifold is muchreduced, but still grows with system size, cf. Tab. II in Ap-pendix B. On clusters featuring one or several M point(s) inthe Brillouin zone, some of these states can be understood aspristine stripy ordered states (sketched in the appropriate re-gion in Fig. 4), combined with line defects, as illustrated inFig. 7. It is however not clear to us whether all states in theground space manifold on all lattices can be generated alongthese ideas. At θ = arctan(1 / we find again a single pointwith a rapidly growing ground state degeneracy. For − π/ < θ < arctan(1 / we find the two Néel ground states, whichare obvious ground states at θ = 0 , corresponding to nearestneighbor interactions only. For π − arctan(1 / < θ < π/ we find a region of phase separation, where the system wantsto be either empty or completely filled with fermions. At θ = π − arctan(1 /
4) ( V = − , V = 1) we have anotherpoint with an extensive ground state degeneracy, albeit with asmaller number than on the reflected side θ = arctan(1 / .Finally the region π/ < θ < π − arctan(1 / hosts a degen-erate set of states, which seem to be related to zigzag chargeconfigurations (sketched in the appropriate region in Fig. 4),with some defect structures in addition.The explicit ground space counting for a number of dif-ferent clusters is listed for further reference in Tab. II in Ap-pendix B. B. Perturbed classical limit: t = 0 + After having established the classical phase diagram, weinclude the effect of a finite hopping amplitude t > on theclassical ground states. We limit ourselves here to perturbativearguments (mostly first order in t ), and check these predictionsby full scale ED in the following sections.We proceed by discussing the fate of the regions in the clas-sical phase diagram, and then analyze the behavior on the de-generate points.
1. Néel CDW state
We start with the least degenerate region, − π/ < θ < arctan(1 / , where there are only two Néel ground states.These two states are symmetry related (and thus cannot besplit by diagonal terms) and it is not possible to connect thetwo states by a finite order in perturbation theory in the ther-modynamic limit. So based on perturbation theory this phaseis expected to have a finite region of stability when t (cid:54) = 0 .
2. Phase separation region
The phase separation region π − arctan(1 / < θ < π/ is also stable at finite t (cid:54) = 0 because both the empty and thefull system are inert under the action of the hopping term.
3. Zigzag state (ZZ)
Next we consider the extended region π/ < θ < π − arctan(1 / , featuring zigzag ground states and some typesof domain walls. To first order in t/V , it is not possible tohop and stay in the ground state manifold. Exact diagonaliza-tions of the full microscopic models on various clusters (seeSec. IV) with true classical ground states reveal however anorder-by-disorder mechanism, where a diagonal term in high-order perturbation theory seems to lift the degeneracy and se-lect the pristine zigzag state, see Fig. 6. This state is six-fold -2.518067-2.518066-2.518065-2.518064-2.518063-2.518062 E ne r g y pe r s i t e V =-8, V =6 N=24
42 states1+2+3 zigzag states 1+2+3 zigzag states1+1 zigzag states32 3222 3333 66 6
FIG. 6. (Color online) Splitting of the ground state degeneracy at V /t = − , V /t = 6 , yielding the two or six zigzag states as theground states. Black (red) symbols denote the two different particle-hole symmetry sectors. All levels are singly degenerate unless the(spatial symmetry) degeneracy is indicated by a blue number. Theperturbatively generated diagonal terms splits the ground space intosymmetry related families, whereby the two or six zigzag states be-come lowest in energy. Note the very small energy difference, indi-cating a high-order perturbative effect. degenerate. A similar state translated to spin language appearsin the Heisenberg-Kitaev model on the honeycomb lattice.
4. Stripy* region
Finally we consider the extended region arctan(1 / <θ < π/ , featuring stripy ground states and some degree ofdomain walls. To first order in t/V , it is not possible to hopand stay in the ground state manifold. However on a finitecluster there is a finite order at which the moves sketched inFig. 7 become possible. As the order at which this processfirst occurs diverges with the linear extent of the cluster, suchtunneling events seem irrelevant at small t/V , in the ther-modynamic limit. As we currently do not fully understand thestructure of the ground space manifold, we can not rule out theappearance of other perturbative processes connecting groundstates without wrapping around the torus.In order to shed some light on this issue, we diagonalizethe full Hamiltonian at V /t = 8 , V /t = 6 for a numberof samples. The resulting spectra are shown in Fig. 8. Incontrast to the zigzag case studied before, here we do not ob-serve a clearly structured order-by-disorder selection of a newground state. The results for these system sizes can be ra-tionalized as the spectrum of an abstract tight binding model,where the sites are the different ground states in the classi-cal ground space, while the hopping matrix elements are ofthe cluster wrapping type just discussed. As the perturba-tive orders of these processes strongly depend on the shapeof the cluster, the strong correlation of the bandwidth of thesplit ground space with the formal order in perturbation the- J.B. Fouet, P. Sindzingre, C. Lhuillier: An investigation of the quantum J − J − J model on the honeycomb lattice 13 Fig. 22.
The N = 18 , , , , ,
32 samples. but frustrates the two others as well as the non coplanarsolutions. The other samples are frustrating but can allowa collinear order if twisted boundary conditions (TBC) areused. This is the case for the N = 18 sample if a twist of π is applied along t or t .To search for spiral order, we used TBC and sweep thewhole range [0 , π ] of twist angles φ , in the t and t directions. These specific boundary conditions are definedas: S ( R i + t j ) = e iφ j S z ( R i ) S ( R i ) e − iφ j S z ( R i ) . (13)This allows to look for boundary conditions which wouldnot frustrate helical ground-states. This approach was foundeffective for the Heisenberg model on the triangular lat-tice to deal with samples frustrating the three-sublatticeN´eel order [1]. For such samples the ground-state energywas found to reach its minimum for the twists which re-lease the frustration: at that point the spectra recover thecharacteristic features of N´eel order.A VBC with the pattern of Read and Sachdev [31](considered in Sect. 3.3) fits in the N = 30 sample butnot on a N = 32 sample. References
1. B. Bernu, C. Lhuillier, and L. Pierre, Phys. Rev. Lett. ,2590 (1992), B. Bernu, P. Lecheminant, C. Lhuillier, andL. Pierre, Phys. Rev. B , 10048 (1994).2. L. Capriotti, A.E. Trumper and S. Sorella, Phys. Rev. Lett. , 3899 (1999).3. P. Sindzingre et al. , Phys. Rev. Lett. , 2953 (2000) andreference therein.4. L. Capriotti and S. Sorella, Phys. Rev. Lett. , 3173(2000) and reference therein.5. R. R. P. Singh et al. , Phys. Rev. B , 7278 (1999). S I ( A, B )( C, D ) (
A, B, C ) (
A, B ) (
A, B, C, D ) G I t R π/ σ R ! π/ σN el Γ Γ − − Γ − Γ − − Γ − − Table 1.
Character table of the permutation group S . Firstline indicates classes of permutations. Second line gives an el-ement of the space symmetry class corresponding to the classof permutation. These space symmetries are: t the one steptranslation ( A → B ), R π/ ( R ! π/ ) the three-fold rotationaround a site of the D ( B )-sublattice, and σ the axial symme-try keeping invariant C and D . N el is the number of elementsin each class. N = 32 S n Γ ( S ) 2 0 3 1 4 2 4 2 4 n Γ ( S ) 1 0 2 1 2 1 2 1 1 n Γ ( S ) 3 0 5 2 6 3 6 3 5 n Γ ( S ) 0 4 4 7 6 8 7 8 6 n Γ ( S ) 0 4 3 6 5 7 5 6 4 Table 2.
Number of occurrences n Γ i ( S ) of each irreduciblerepresentation Γ i in { ˜ E } as a function of the total spin S . N = 32 S n Γ ( S ) 1 0 1 0 1 0 1 0 1 n Γ ( S ) 1 0 1 0 1 0 1 0 1 n Γ ( S ) 0 1 0 1 0 1 0 1 0 Table 3.
Number of occurrences n Γ i ( S ) of each irreduciblerepresentation Γ i in { ˜ E } as a function of the total spin S .6. V. N. Kotov, J. Oitmaa, O. Sushkov and Z. Weihong, Phys.Rev. B , 14613 (1999), and cond-mat/9912228.7. A. Chubukov and T. Jolicoeur, Phys. Rev. B , 11137(1992).8. S. E. Korshunov, Phys. Rev. B , 6165 (1993).9. P. Lecheminant, B. Bernu, C. Lhuillier, and L. Pierre,Phys. Rev. B , 9162 (1995).10. G. Misguich, B. Bernu, C. Lhuillier, and C. Waldtmann,Phys. Rev. Lett. , 1098 (1998); G. Misguich, C. Lhuil-lier, B. Bernu, and C. Waldtmann, Phys. Rev. B , 1064(1999).11. J. Oitmaa and D.D. Betts, Can. J. Phys.
897 (1978).12. J.D. Reger, J.A. Riera and A.P. Young, J. Phys.:CondensMatter , 1855 (1989).13. J. Oitmaa, C.J. Hamer and Z. Weihong, Phys. Rev. B ,9834 (1992).14. Z. Weihong, J. Oitmaa and C.J. Hamer, Phys. Rev. B ,11689 (1991).15. A. Mattsson, P. Fr¨ojdh and T. Einarson, Phys. Rev. B ,3397 (1994). J.B. Fouet, P. Sindzingre, C. Lhuillier: An investigation of the quantum J − J − J model on the honeycomb lattice 13 Fig. 22.
The N = 18 , , , , ,
32 samples. but frustrates the two others as well as the non coplanarsolutions. The other samples are frustrating but can allowa collinear order if twisted boundary conditions (TBC) areused. This is the case for the N = 18 sample if a twist of π is applied along t or t .To search for spiral order, we used TBC and sweep thewhole range [0 , π ] of twist angles φ , in the t and t directions. These specific boundary conditions are definedas: S ( R i + t j ) = e iφ j S z ( R i ) S ( R i ) e − iφ j S z ( R i ) . (13)This allows to look for boundary conditions which wouldnot frustrate helical ground-states. This approach was foundeffective for the Heisenberg model on the triangular lat-tice to deal with samples frustrating the three-sublatticeN´eel order [1]. For such samples the ground-state energywas found to reach its minimum for the twists which re-lease the frustration: at that point the spectra recover thecharacteristic features of N´eel order.A VBC with the pattern of Read and Sachdev [31](considered in Sect. 3.3) fits in the N = 30 sample butnot on a N = 32 sample. References
1. B. Bernu, C. Lhuillier, and L. Pierre, Phys. Rev. Lett. ,2590 (1992), B. Bernu, P. Lecheminant, C. Lhuillier, andL. Pierre, Phys. Rev. B , 10048 (1994).2. L. Capriotti, A.E. Trumper and S. Sorella, Phys. Rev. Lett. , 3899 (1999).3. P. Sindzingre et al. , Phys. Rev. Lett. , 2953 (2000) andreference therein.4. L. Capriotti and S. Sorella, Phys. Rev. Lett. , 3173(2000) and reference therein.5. R. R. P. Singh et al. , Phys. Rev. B , 7278 (1999). S I ( A, B )( C, D ) (
A, B, C ) (
A, B ) (
A, B, C, D ) G I t R π/ σ R ! π/ σN el Γ Γ − − Γ − Γ − − Γ − − Table 1.
Character table of the permutation group S . Firstline indicates classes of permutations. Second line gives an el-ement of the space symmetry class corresponding to the classof permutation. These space symmetries are: t the one steptranslation ( A → B ), R π/ ( R ! π/ ) the three-fold rotationaround a site of the D ( B )-sublattice, and σ the axial symme-try keeping invariant C and D . N el is the number of elementsin each class. N = 32 S n Γ ( S ) 2 0 3 1 4 2 4 2 4 n Γ ( S ) 1 0 2 1 2 1 2 1 1 n Γ ( S ) 3 0 5 2 6 3 6 3 5 n Γ ( S ) 0 4 4 7 6 8 7 8 6 n Γ ( S ) 0 4 3 6 5 7 5 6 4 Table 2.
Number of occurrences n Γ i ( S ) of each irreduciblerepresentation Γ i in { ˜ E } as a function of the total spin S . N = 32 S n Γ ( S ) 1 0 1 0 1 0 1 0 1 n Γ ( S ) 1 0 1 0 1 0 1 0 1 n Γ ( S ) 0 1 0 1 0 1 0 1 0 Table 3.
Number of occurrences n Γ i ( S ) of each irreduciblerepresentation Γ i in { ˜ E } as a function of the total spin S .6. V. N. Kotov, J. Oitmaa, O. Sushkov and Z. Weihong, Phys.Rev. B , 14613 (1999), and cond-mat/9912228.7. A. Chubukov and T. Jolicoeur, Phys. Rev. B , 11137(1992).8. S. E. Korshunov, Phys. Rev. B , 6165 (1993).9. P. Lecheminant, B. Bernu, C. Lhuillier, and L. Pierre,Phys. Rev. B , 9162 (1995).10. G. Misguich, B. Bernu, C. Lhuillier, and C. Waldtmann,Phys. Rev. Lett. , 1098 (1998); G. Misguich, C. Lhuil-lier, B. Bernu, and C. Waldtmann, Phys. Rev. B , 1064(1999).11. J. Oitmaa and D.D. Betts, Can. J. Phys.
897 (1978).12. J.D. Reger, J.A. Riera and A.P. Young, J. Phys.:CondensMatter , 1855 (1989).13. J. Oitmaa, C.J. Hamer and Z. Weihong, Phys. Rev. B ,9834 (1992).14. Z. Weihong, J. Oitmaa and C.J. Hamer, Phys. Rev. B ,11689 (1991).15. A. Mattsson, P. Fr¨ojdh and T. Einarson, Phys. Rev. B ,3397 (1994). J.B. Fouet, P. Sindzingre, C. Lhuillier: An investigation of the quantum J − J − J model on the honeycomb lattice 13 Fig. 22.
The N = 18 , , , , ,
32 samples. but frustrates the two others as well as the non coplanarsolutions. The other samples are frustrating but can allowa collinear order if twisted boundary conditions (TBC) areused. This is the case for the N = 18 sample if a twist of π is applied along t or t .To search for spiral order, we used TBC and sweep thewhole range [0 , π ] of twist angles φ , in the t and t directions. These specific boundary conditions are definedas: S ( R i + t j ) = e iφ j S z ( R i ) S ( R i ) e − iφ j S z ( R i ) . (13)This allows to look for boundary conditions which wouldnot frustrate helical ground-states. This approach was foundeffective for the Heisenberg model on the triangular lat-tice to deal with samples frustrating the three-sublatticeN´eel order [1]. For such samples the ground-state energywas found to reach its minimum for the twists which re-lease the frustration: at that point the spectra recover thecharacteristic features of N´eel order.A VBC with the pattern of Read and Sachdev [31](considered in Sect. 3.3) fits in the N = 30 sample butnot on a N = 32 sample. References
1. B. Bernu, C. Lhuillier, and L. Pierre, Phys. Rev. Lett. ,2590 (1992), B. Bernu, P. Lecheminant, C. Lhuillier, andL. Pierre, Phys. Rev. B , 10048 (1994).2. L. Capriotti, A.E. Trumper and S. Sorella, Phys. Rev. Lett. , 3899 (1999).3. P. Sindzingre et al. , Phys. Rev. Lett. , 2953 (2000) andreference therein.4. L. Capriotti and S. Sorella, Phys. Rev. Lett. , 3173(2000) and reference therein.5. R. R. P. Singh et al. , Phys. Rev. B , 7278 (1999). S I ( A, B )( C, D ) (
A, B, C ) (
A, B ) (
A, B, C, D ) G I t R π/ σ R ! π/ σN el Γ Γ − − Γ − Γ − − Γ − − Table 1.
Character table of the permutation group S . Firstline indicates classes of permutations. Second line gives an el-ement of the space symmetry class corresponding to the classof permutation. These space symmetries are: t the one steptranslation ( A → B ), R π/ ( R ! π/ ) the three-fold rotationaround a site of the D ( B )-sublattice, and σ the axial symme-try keeping invariant C and D . N el is the number of elementsin each class. N = 32 S n Γ ( S ) 2 0 3 1 4 2 4 2 4 n Γ ( S ) 1 0 2 1 2 1 2 1 1 n Γ ( S ) 3 0 5 2 6 3 6 3 5 n Γ ( S ) 0 4 4 7 6 8 7 8 6 n Γ ( S ) 0 4 3 6 5 7 5 6 4 Table 2.
Number of occurrences n Γ i ( S ) of each irreduciblerepresentation Γ i in { ˜ E } as a function of the total spin S . N = 32 S n Γ ( S ) 1 0 1 0 1 0 1 0 1 n Γ ( S ) 1 0 1 0 1 0 1 0 1 n Γ ( S ) 0 1 0 1 0 1 0 1 0 Table 3.
Number of occurrences n Γ i ( S ) of each irreduciblerepresentation Γ i in { ˜ E } as a function of the total spin S .6. V. N. Kotov, J. Oitmaa, O. Sushkov and Z. Weihong, Phys.Rev. B , 14613 (1999), and cond-mat/9912228.7. A. Chubukov and T. Jolicoeur, Phys. Rev. B , 11137(1992).8. S. E. Korshunov, Phys. Rev. B , 6165 (1993).9. P. Lecheminant, B. Bernu, C. Lhuillier, and L. Pierre,Phys. Rev. B , 9162 (1995).10. G. Misguich, B. Bernu, C. Lhuillier, and C. Waldtmann,Phys. Rev. Lett. , 1098 (1998); G. Misguich, C. Lhuil-lier, B. Bernu, and C. Waldtmann, Phys. Rev. B , 1064(1999).11. J. Oitmaa and D.D. Betts, Can. J. Phys.
897 (1978).12. J.D. Reger, J.A. Riera and A.P. Young, J. Phys.:CondensMatter , 1855 (1989).13. J. Oitmaa, C.J. Hamer and Z. Weihong, Phys. Rev. B ,9834 (1992).14. Z. Weihong, J. Oitmaa and C.J. Hamer, Phys. Rev. B ,11689 (1991).15. A. Mattsson, P. Fr¨ojdh and T. Einarson, Phys. Rev. B ,3397 (1994). J.B. Fouet, P. Sindzingre, C. Lhuillier: An investigation of the quantum J − J − J model on the honeycomb lattice 13 Fig. 22.
The N = 18 , , , , ,
32 samples. but frustrates the two others as well as the non coplanarsolutions. The other samples are frustrating but can allowa collinear order if twisted boundary conditions (TBC) areused. This is the case for the N = 18 sample if a twist of π is applied along t or t .To search for spiral order, we used TBC and sweep thewhole range [0 , π ] of twist angles φ , in the t and t directions. These specific boundary conditions are definedas: S ( R i + t j ) = e iφ j S z ( R i ) S ( R i ) e − iφ j S z ( R i ) . (13)This allows to look for boundary conditions which wouldnot frustrate helical ground-states. This approach was foundeffective for the Heisenberg model on the triangular lat-tice to deal with samples frustrating the three-sublatticeN´eel order [1]. For such samples the ground-state energywas found to reach its minimum for the twists which re-lease the frustration: at that point the spectra recover thecharacteristic features of N´eel order.A VBC with the pattern of Read and Sachdev [31](considered in Sect. 3.3) fits in the N = 30 sample butnot on a N = 32 sample. References
1. B. Bernu, C. Lhuillier, and L. Pierre, Phys. Rev. Lett. ,2590 (1992), B. Bernu, P. Lecheminant, C. Lhuillier, andL. Pierre, Phys. Rev. B , 10048 (1994).2. L. Capriotti, A.E. Trumper and S. Sorella, Phys. Rev. Lett. , 3899 (1999).3. P. Sindzingre et al. , Phys. Rev. Lett. , 2953 (2000) andreference therein.4. L. Capriotti and S. Sorella, Phys. Rev. Lett. , 3173(2000) and reference therein.5. R. R. P. Singh et al. , Phys. Rev. B , 7278 (1999). S I ( A, B )( C, D ) (
A, B, C ) (
A, B ) (
A, B, C, D ) G I t R π/ σ R ! π/ σN el Γ Γ − − Γ − Γ − − Γ − − Table 1.
Character table of the permutation group S . Firstline indicates classes of permutations. Second line gives an el-ement of the space symmetry class corresponding to the classof permutation. These space symmetries are: t the one steptranslation ( A → B ), R π/ ( R ! π/ ) the three-fold rotationaround a site of the D ( B )-sublattice, and σ the axial symme-try keeping invariant C and D . N el is the number of elementsin each class. N = 32 S n Γ ( S ) 2 0 3 1 4 2 4 2 4 n Γ ( S ) 1 0 2 1 2 1 2 1 1 n Γ ( S ) 3 0 5 2 6 3 6 3 5 n Γ ( S ) 0 4 4 7 6 8 7 8 6 n Γ ( S ) 0 4 3 6 5 7 5 6 4 Table 2.
Number of occurrences n Γ i ( S ) of each irreduciblerepresentation Γ i in { ˜ E } as a function of the total spin S . N = 32 S n Γ ( S ) 1 0 1 0 1 0 1 0 1 n Γ ( S ) 1 0 1 0 1 0 1 0 1 n Γ ( S ) 0 1 0 1 0 1 0 1 0 Table 3.
Number of occurrences n Γ i ( S ) of each irreduciblerepresentation Γ i in { ˜ E } as a function of the total spin S .6. V. N. Kotov, J. Oitmaa, O. Sushkov and Z. Weihong, Phys.Rev. B , 14613 (1999), and cond-mat/9912228.7. A. Chubukov and T. Jolicoeur, Phys. Rev. B , 11137(1992).8. S. E. Korshunov, Phys. Rev. B , 6165 (1993).9. P. Lecheminant, B. Bernu, C. Lhuillier, and L. Pierre,Phys. Rev. B , 9162 (1995).10. G. Misguich, B. Bernu, C. Lhuillier, and C. Waldtmann,Phys. Rev. Lett. , 1098 (1998); G. Misguich, C. Lhuil-lier, B. Bernu, and C. Waldtmann, Phys. Rev. B , 1064(1999).11. J. Oitmaa and D.D. Betts, Can. J. Phys.
897 (1978).12. J.D. Reger, J.A. Riera and A.P. Young, J. Phys.:CondensMatter , 1855 (1989).13. J. Oitmaa, C.J. Hamer and Z. Weihong, Phys. Rev. B ,9834 (1992).14. Z. Weihong, J. Oitmaa and C.J. Hamer, Phys. Rev. B ,11689 (1991).15. A. Mattsson, P. Fr¨ojdh and T. Einarson, Phys. Rev. B ,3397 (1994). N = 24 N = 32 FIG. 7. (Color online) Illustration of the stripy charge ordering pattern and the line defects for the N = 24 (left panel) and N = 32 (rightpanel) clusters. The left part in each panel indicates a pristine stripy configuration, while the right part shows how another ground state can bereached by a one-dimensional collective transposition of particles. For the N = 24 cluster the line move wraps completely around the sampleand connects two ground states without defect lines, while for the N = 32 the indicated shifts of particles generate a ground state with a defectline, where the region along the lines resembles a rotated pristine stripy state. For the N = 24 cluster this move is generated in th orderperturbation theory, while the move on the N = 32 cluster is effective in th order. -2.62-2.6-2.58-2.56-2.54 E ne r g y pe r s i t e V =8, V =6 -2.586-2.585-2.584-2.583-2.582 E ne r g y pe r s i t e V =8, V =6 N=24
N=28
N=32
42 states
N=24b
12 states
N=12
FIG. 8. (Color online) Splitting of the classical ground state degen-eracy at V /t = 8 , V /t = 6 for clusters with N = 12 (left panel)and N = 24 , b, and sites (right panel). Black (red) symbolsdenote the two different particle-hole symmetry sectors. All levelsare singly degenerate unless the (spatial symmetry) degeneracy is in-dicated by a blue number. The order at which processes of the typedescribed in Fig. 7 happen on the clusters are indicated at the bot-tom of the frame for each system size. In contrast to the case withattractive V there is no obvious system-size independent order-by-disorder mechanism at work on the clusters considered. ory seems to validate our picture. It is thus presently not clearwhether the ground state degeneracy will be lifted at a finiteorder in perturbation theory in the thermodynamic limit ornot. The qualitative presence of a high-order diagonal termin the zigzag case discussed before could also carry over tothe present stripy* case, but is masked on moderate finite sizesamples by the leading cluster-wrapping splitting terms.Now we turn to the points in the classical phase diagramwhich feature a massive ground state degeneracy. These arethe points θ = π/ ( V = 0 ), θ = arctan(1 /
4) ( V =4 , V = 1) and θ = π − arctan(1 /
4) ( V = − , V = 1) .
5. Charge modulated (CM) state At V = 0 we have two independent triangular sublatticeswith antiferromagnetic interactions. It is known that each statewhich has two or one charges on each V -triangle is a validground state. When working at half filling the ground spaceof the full problem can be obtained by enumerating all groundstates on one triangular sublattice and then summing up theproduct of the two sublattice Hilbert spaces over all particlebipartitions such that the total number of particles amounts tohalf filling (i.e. N/ ).The effect of a finite t/V is to hop particles from one trian-gular sublattice to the other one. There are indeed many statesin the ground state manifold where such first order hoppingsare possible while staying in the ground state manifold. Theproblem thus maps to a hopping problem on an abstract latticeliving in the constrained Hilbert space, where the sites are theallowed configurations. This abstract lattice is not homoge-neous and features sites with different coordination numbers.It is not known to us how to solve this hopping problem an-alytically. As a starting point it is instructive to identify themaximally flippable states (i.e. abstract sites with the largestcoordination number) of this problem. By explicit enumera-tion we have identified the 18 uud - ddu configurations as themaximally flippable states (the number is independent of N ,given that the sample geometry is compatible with uud - ddu states). These states feature a threefold degenerate charge ana-log of the magnetic uud states on one sublattice, while theother sublattice has three possible charge analogs of the mag-netic ddu state. Another factor of two is obtained by swappingthe role of the two sublattices.The question is now whether the ground state of the full ef-fective model is localized on those maximally flippable states,or whether the hopping network favors delocalized wave func-tions. We have numerically solved the effective model forsamples with N = 18 , , and sites and show the low-lying energy spectra in Fig. 9. We note the presence of 18 low-lying states (indicated by the full arrow in Fig. 9), separatedby a gap (dashed arrow) from the higher lying states. Whilethis picture is not yet very sharp for N = 18 and N = 24 , for N = 42 and N = 54 however the bandwidth of the 18 low-lying states is quite small compared to the gap to the higherstates. We interpret this spectrum as the localization of thelow-lying eigenstates on the 18 maximally flippable states.This is further corroborated by the fact that these eigenfunc-tions have their maximal amplitudes on the 18 uud - ddu states.We have also obtained the static charge structure factor S ( Q ) = 1 N (cid:88) k,l e − i Q · ( R k − R l ) ( n k − / n l − / (2)in the ground state of the effective model. In Fig. 10 we dis-play this quantity for N = 24 and N = 42 , as well as foran equal superposition of the 18 uud - ddu states on N = 24 sites. Note that we plot the data in an extended Brillouinzone scheme. There are obvious Bragg peaks at the momenta K , K (cid:48) = − K stemming from the tripling of the unit cell ineach triangular sublattice. Furthermore in the plain superposi-tion of the 18 uud - ddu states a sublattice imbalance is present,as the densities on the two triangular sublattices are differ-ent (1/3 versus 2/3), leading to another set of Bragg peaksat the Γ ∗ momentum. In the actual ground state of the ef-fective model this peak is renormalized substantially. If ouradvocated picture of a localization of the ground state wavefunctions on the 18 uud - ddu states holds true, then we alsoexpect the Bragg peaks at Γ ∗ to survive, although with a siz-able reduction of the amplitude. In a recent paper Daghoferand Hohenadler also investigated this issue and concludedan absence of the sublattice imbalance based on exact diag-onalizations of the original Hamiltonian (1) on N = 18 and N = 24 . Our results plotted in the extended zone schemehowever reveal at least the presence of a local maximum in S ( Q ) at Γ ∗ , consistent with our claim of a weak sublatticeimbalance in the thermodynamic limit. Note that this statewould be insulating, i.e. unrelated to the presumed metallicpinball liquid phase on the triangular lattice. We close the discussion of the V = 0 case by noting thatthis point enlarges to a finite region with respect to the addi-tion of V when t > . The reason is that the ground state ofthe hopping problem in the degenerate ground state manifoldis able to gain energy proportional to t in the order-by-disorderframework. When switching on V , one needs a finite amountof V to destabilize the CM phase, because the stripy* regionand the zigzag phase do not gain energy from the hoppingterm at first order in t .
6. Néel domain wall crystal (NDWC)
Now let us discuss the behavior at θ = arctan(1 / . Asnoted before, this point features a sizable ground state degen-eracy at the classical level. Since - to the best of our knowl-edge - this manifold has not been studied so far, we proceedby numerically diagonalizing the first order effective Hamil-tonian obtained by projecting the fermionic hopping Hamilto-nian into the classical ground state manifold. In Fig. 11 wepresent the energy per site for various clusters. It appears thatamong the studied clusters the samples with N = 24 and
18 24 42 54 N E - E G S [t] s p li tt i ng & gap [t] splittinggap s t a t e s G ap FIG. 9. (Color online) Low lying energy spectrum in units of thehopping t for the CM phase at V = 0 , V = 1 , t/V (cid:28) . Datafor system sizes N = 18 , , , are shown. The full arrowdenotes the bandwidth of the lowest 18 levels. Note the collapse ofthe bandwidth for N = 42 and N = 54 , while the (particle-hole) gapabove the 18 states (dashed arrows) remains of order t , as illustratedin the inset. -4 -2 0 2 4 k x -4-2024 k y KK’ Γ * Γ * Γ * Γ * Γ * Γ * K’K’KK KKK’K’ FIG. 10. (Color online) Charge structure factor of the effective modelat V = 0 , V = 1 , t/V (cid:28) plotted in the enlarged Brillouin zone.The filled symbols are data for the ground state at N = 24 whilethe empty symbols are for N = 42 . The dashed symbols denote thestructure factor in an equal superposition of the 18 uud - ddu statesfor N = 24 . Note the presence of small peaks at momentum Γ ∗ ,indicating a possibly weak sublattice imbalance. sites have a particularly low energy per site ( N = 28 is close,but its energy per site is higher than N = 36 , so we discardit). It seems that this low energy correlates with a high max-imal flippability of configurations on these clusters , com-pared to the other clusters. In Fig. 12 we display one suchmaximally flippable state for N = 24 . Our understanding ofthis state is that it is formed of alternating strips of the twoNéel CDW states. Along the domain walls between the two
24 26 28 30 32 34 36 38 40 42-0.21-0.2-0.19-0.18-0.17-0.16-0.15-0.14-0.13-0.12 E / N [t]
24 26 28 30 32 34 36 38 40 420123456 G ap [t] FIG. 11. (Color online) Energetics of the first order effective model(1st order in t/V , ) at V /V = 4 . (Top) Energy per site in units of t for various cluster sizes. The clusters with and sites have aparticularly low energy per site. (Bottom) Single-particle charge gap ∆ /t . This gap amounts to ∼ t for the system sizes with a lowenergy per site. CDW states there are many bonds which are flippable to firstorder in t/V , . This state is 18 fold degenerate. There are3 directions in which the domain walls can run, and for fixedorientation, there are 3 distinct bond choices for the domainwalls. The remaining factor of two is due to the assignmentof the two Néel CDW states. We have further checked a sam-ple with N = 72 and again found 18 maximally flippablestates, in agreement with our analysis. In Fig. 13 we displaythe low lying eigenstates of the effective first-order Hamilto-nian in the degenerate subspace for system sizes N = 24 , and sites. The collapse of the lowest 18 states is clearlyvisible, providing strong evidence for symmetry breaking ofthe NDWC type. The results of the effective model simula-tions for the charge structure factor are shown in Fig. 14 andfurther corroborate this picture.We also expect the point V /V = 4 to enlarge to an actualphase at small but finite t/V , , because the two neighbor-ing phases in the Ising limit (the stripy* region and the NéelCDW) both do not gain energy to first order in t/V , . V = − , V = 1 The ground state configurations at V = − , V = 1 donot seem to be flippable to first order in t/V , V , and we havetherefore refrained from a detailed study of this line.To conclude this section we show the resulting semiclassi-cal qualitative phase diagram in the right part of Fig. 4. IV. NUMERICAL RESULTS
In this section, we will present our systematic numericalstudy using ED and typical observables that were used to com-pute our global phase diagram in Fig. 3. We will first present
Néel ANéel BNéel B
FIG. 12. (Color online) The Néel Domain Wall Crystal (NDWC):sketch of a classical ground state at V = 4 , V = 1 on the N =24 sample, which is maximally flippable within the classical groundstate manifold with respect to the hopping t . The shaded regionsdenote the two Néel domains, and the orange circled bonds along thedomain walls are flippable to first order in t . The green box indicatesa twelve-site unit cell.
24 36 72 N E - E G S [t] s p li tt i ng & gap [t] splittinggap s t a t e s s t a t e s s t a t e s FIG. 13. (Color online) Low lying energy spectrum in units of thehopping t for V = 4 , V = 1 , t/V (cid:28) . Data for system sizes N = 24 , , are shown. The full arrow denotes the bandwidthof the lowest 18 levels ( N = 36 does not have rotational invariance,in that case we only expect 6 possible ground states). Note the col-lapse of the bandwidth for larger clusters, while the (particle-hole)gap above the 18 states (dashed arrows) remains approximately con-stant, as illustrated in the inset. some useful tools to detect phase transitions as well as to char-acterize phases. Then we will provide numerical evidence forthe finite extension of the phases that were discussed in thelimit of large interactions in Sec. III. Last but not least, we willdiscuss the stability of phases (with a stronger quantum char-acter) at intermediate coupling, namely the semi-metal (SM),the plaquette/Kekulé phase and the absence of the conjecturedtopological Mott insulator phase. -4 -2 0 2 4 k x -4-2024 k y KK’ Γ * Γ * Γ * Γ * Γ * Γ * K’K’KK KKK’K’ FIG. 14. (Color online) Charge structure factor of the effective modelat V = 4 , V = 1 , t/V (cid:28) plotted in the enlarged Brillouin zone.The filled (empty) symbols are data for the ground state at N = 24 ( N = 72 ). The Bragg peaks appear at the X ∗ point, c.f. Fig. 1(b). A. Numerical tools
1. Energetics
Let us start by showing some full low-energy spectra, wherewe label each eigenstate according to its quantum numbers. In Fig. 15-16, we present these spectra for a fixed V = 0 and V /t = 4 respectively on the N = 24 cluster which has thefull point-group symmetry C v and also possesses the relevant K and M points in its Brillouin zone (see Appendix A).When V /t = 0 and V /t is small, we expect a finite regionof SM phase with low-energy excitations at the Dirac point( ± K ). Increasing V /t , there is a sudden level crossing when V /t (cid:39) . . Beyond this value, if we count the lowest energylevels well separated from the others, we get 18 states, whichmay indicate some ordering of some kind in the thermody-namic limit (see discussion in Sec. III B).When V /t = 4 , starting at small V /t , there are twoalmost-degenerate states at the Γ point with opposite particle-hole symmetries, compatible with the Néel CDW state. In-creasing V /t leads to some intermediate phase < V /t < with low-energy states at momentum K compatible with aKekulé-like (or plaquette) pattern, and finally for large V /t > the emergence of 6 low-energy states with momenta Γ and M . B. Stability of the classical Ising-like phases
Let us now move to a more precise characterization of thepossible patterns occurring for large interactions. Quite gener-ically, large density interactions are expected to lead to somekind of charge ordering. In the equivalent spin language, theseterms are of the Ising type and the resulting Ising model would V / t e n e r gy l e v e l s Γ A1 Γ A2 Γ E1 Γ E2 Γ D1 Γ D2XeXoKA1KA2KEMe1Me2Mo1Mo2
FIG. 15. (Color online) Full spectrum vs V /t at fixed V /t = 0 labeled with symmetry sectors on N = 24 cluster. More precisely,states with momentum Γ can be labeled with C v standard point-group notations; with momentum X only reflection can be used todetect even (e) vs odd (o) states; at the K point, we use C v point-group notations; at the M point, we can use two reflections to la-bel states which are even/even (e1), even/odd (e2), odd/even (o1) orodd/odd (o2). Open (filled) symbols denote respectively states whichare even (odd) with respect to particle-hole symmetry. V / t e n e r gy l e v e l s Γ A1 Γ A2 Γ E1 Γ E2 Γ D1 Γ D2XeXoKA1KA2KEMe1Me2Mo1Mo2
FIG. 16. (Color online) Same as Fig. 15 at fixed V /t = 4 on N =24 cluster. For intermediate values V /t (cid:39) . , low-energy states(indicated with an ellipse) are compatible with a plaquette or Kekuléphase. correspondingly exhibit some spin ordering. For model (1),they were discussed in Sec. III and are denoted under thenames: CDW, NDWC, stripy* region, CM or zigzag.
1. CDW phase
When V = 0 , since the honeycomb lattice is bipartite, asimple charge-density wave (CDW) where particles are local-ized on one sublattice is expected for V > , so that de-0 FIG. 17. (Color online) Connected density correlation functions be-tween a reference site (open circle) and other sites for V = 1 and V = 2 (at fixed V = 0 ) on N = 32 cluster. Periodic boundaryconditions are used. Blue and red correspond to positive/negativecorrelations. generacy is two. Since this phase does not break translationsymmetry but only inversion, symmetry analysis predicts thaton a finite cluster there should be a collapse of the lowest Γ A and Γ D states. This is indeed what is observed in Fig. 16at small V /t for instance. When V ≥ , our numerical dataare compatible with a vanishing gap at zero momentum (datanot shown) because of the inversion symmetry breaking in theCDW phase. Note however that the system remains an insu-lator.A clearer picture of this CDW state is provided by density-density correlations. On Fig. 17, we plot these connected cor-relations obtained with a N = 32 cluster showing a clear tran-sition from a semi-metallic phase to the Néel CDW for large V /t ≥ (at V = 0 ). About the critical value, let us mentionthat a recent unbiased QMC approach has concluded that ( V /t ) c = 1 . .Based on our previous semiclassical analysis, we expectthat the CDW phase cannot extend beyond the line V = 4 V ,at least towards strong coupling.
2. Charge modulated phase
As we have discussed in Sec. III, the large V /t > limit isless trivial since the classical Ising model possesses an ex-tensive ground-state degeneracy. However, one expects ongeneral grounds that quantum fluctuations will provide someselection mechanism within these low-energy states. Our dis-cussion in Sec. III B 5 has led us to conclude in favor of somecharge-modulation (CM) state which is 18-fold degenerateand exhibits charge imbalance between the two sublattices,in agreement with Refs. 6 and 10 and slightly different fromthe similar CM state without charge imbalance .Our ground-state energy plot in Fig. 2(a) has indeed shownthat the emergent phase is compatible with clusters having the K point, as is the case for the proposed tripled unit cell. More-over, the low-energy spectra on these clusters, such as N = 24 in Fig. 15 (or similarly N = 36 , data not shown), do confirmthe presence of 18 low-energy states at large V in agreementwith the expected degeneracy. The question whether all those states will collapse or if there will remain a finite gap (andthus a smaller degeneracy) is left for future studies, but oureffective model diagonalizations in the strong coupling limitindicate a complete collapse of the 18 levels, c.f. Fig. 9.
3. Stripy* region
For large V ∼ V > , the lowest energies are found onclusters having the M point (see Fig. 2(b)), which indicates adifferent instability. Moreover, our symmetry resolved spec-tra (see Fig. 16 for instance) predict a 6-fold degeneracy forlarge V /t (with 3 states at momentum Γ and 3 at equivalentmomenta M) in disagreement with both the CM phase or theKekulé one (see later). As discussed in Sec. III B 4 based alsoon ED, the fate of the stripy* region at finite t/V , V is notclear. Results for N = 32 close to the phase transitions to theneighboring phases still do not shed light on a possible selec-tion process among the classically degenerate ground states.
4. Néel domain wall crystal
Since the ratio V /V = 4 is a special point in the classicalphase diagram we investigate here the phase diagram alongthe line V /t , with the ratio V /V = 4 fixed. In Fig. 18 wedisplay the low energy spectrum for the N = 24 site sam-ple along this line. At very small V /t we expect the sta-ble semimetallic phase emerging out of V /t = 0 . We thendetect an intermediate Kekulé or plaquette phase for values (cid:46) V /t (cid:46) . These values are confirmed by directly com-puting kinetic correlations and checking if they match the ex-pected pattern for a plaquette phase. However for larger val-ues V /t (cid:38) the system enters the realm of the NDWC phase,involving the physics discussed in subsection III B 6. The oc-currence of two very distinct phases beyond the semimetallicphase is also visible in the ground state correlations shown inFig. 19.
5. Zizag phase
As already discussed in Sec. III B, ED with large negative V /t and positive V /t provide spectral evidence for a quan-tum order-by-disorder selection of a pristine zigzag state, asshown for instance in Fig. 6. For the clusters having the 3 M points in their Brillouin zone (see Appendix A), we doobserve 6-fold quasi degenerate ground-states. By comput-ing connected density correlations on any of these states, weobtain a pattern compatible with zigzag state. Indeed, a sim-ple direct computation of connected density correlations onthe zigzag state leads to values ± / or ± / dependingon distances between sites, which agrees to a very high accu-racy with the exact results computed on N = 24 or N = 32 ground-states in this region, see Fig. 20.Having clarified the phases at large interaction strengths,which are connected to the classical limit, we now con-1 V / t e n e r gy l e v e l s Γ A1 Γ A2 Γ E1 Γ E2 Γ D1 Γ D2XeXoKA1KA2KEMe1Me2Mo1Mo2
FIG. 18. (Color online) Same as Fig. 15 as a function of V /t at fixed V /V = 4 on N = 24 cluster. The orange boxes denote the twoparameter values for which the correlations shown in Fig. 19 werecalculated. V /t = 4 V /t = 1 V /t = 40 V /t = 10 kinetic energycorrelations densitycorrelations Kekule/Plaquette phaseNéel Domain Wall Crystal FIG. 19. (Color online) Kinetic energy and density correlationscomputed on N = 24 cluster along the V /V = 4 line, i.e. for ( V /t, V /t ) respectively equal to (4 , and (40 , . sider the intermediate interaction regime where new quantumphases are expected to emerge. C. Emergence of a resonating quantum plaquette phase orKekulé one
In Refs. 6 and 10, it has been argued using mean-field or N = 18 ED respectively that a Kekulé pattern can be stabi-lized when V and V are close. This would correspond to a3-fold degeneracy
30 31 of the ground-state (more precisely onestate at Γ and two states at ± K ).However, while the Kekulé phase was seen to remain stableup to large V ∼ V in Refs. 6 and 10, our numerical data FIG. 20. (Color online) Connected density correlations computed on N = 24 (left panel) and N = 32 (right panel) cluster for V /t = − and V /t = 6 . Blue and red correspond to positive/negativecorrelations. on large clusters do not confirm these findings (see Sec. IV B)since it is replaced either by the stripy* region or the NDWCphase. However, we do confirm its existence, but for interme-diate values only.A first evidence is given by the low-energy spectrum of N = 24 at fixed V /t = 4 (see Fig. 16) where in some inter-mediate V /t range, the lowest energies are compatible withthis plaquette/Kekulé phase. A second one is provided viaa direct computation of the (connected) kinetic energy corre-lations on a larger N = 42 cluster in Fig. 21, which are inperfect agreement with the expected pattern. FIG. 21. (Color online) Kinetic energy correlations at fixed V /t = 8 and V /t = 2 on N = 42 cluster. Periodic boundary conditionsare used. Blue and red correspond to positive/negative correlations.Width is proportional to the correlation. There is a small caveat here, since it is rather difficult to dis-tinguish a static Kekulé pattern from a plaquette crystal, wherethe three strong bonds on a hexagon resonate. Indeed bothphases break the same symmetries and have the same degen-eracy so that a clear identification is often difficult.
Letus mention that for S = 1 / Heisenberg model on the hon-2eycomb lattice with competing J , J (and J ) interactions,there exists an intermediate region with a clear plaquette char-acter. Moreover, in this case, the plaquette phase seemsconnected to the ordered phase through a possibly continu-ous phase transition.
In our spinless fermionic model, wehave not fully clarified the nature of this intermediate phase,but this would be an interesting project as well as the nature ofthe transition both to the semi-metallic and the Néel phases,which is unfortunately beyond what can be done using EDtechnique, but which might be fermionic examples of phasetransitions beyond the Landau paradigm.
D. Stability of the topological phase
The proposal made by Raghu et al. is that, when V > dominates, the system undergoes a transition to a topolog-ically distinct phase that spontaneously break time-reversalsymmetry. A simple picture of this phase can be obtained bydrawing circulating currents on a given sublattice. In order tocheck this picture, we have computed current-current corre-lations on a given sublattice. Note that these correlations arestrictly zero for all distances when V = V = 0 . They areplotted for different V in Fig. 22 in the case of N = 24 lat-tice and at fixed V = 0 where the instability is expected to bethe strongest. We have used the expected q = 0 orientationpattern of the QAH phase, i.e. where all currents would bein the clockwise direction on each hexagon for instance. Wesee on Fig. 22 that when we increase V /t , the signal starts byincreasing in amplitude, and then some defects (in red) start toappear on the cluster, in particular in the case where QAH issupposed to be the strongest around V /t (cid:39) . Let us empha-size again that QAH is perfectly compatible with all clusters(since it is a q = 0 instability) and as a result no frustrationis expected. Therefore, the appearance of defects reveal somecompeting phases. Still the overall pattern is almost perfect,which is remarkable and does not occur in the V only casefor instance. For even larger V , the signal suddenly disap-pears (see Figs. 23 and 24 for several clusters, other data notshown), signalling the entrance in the CM phase (see above).Given that the semimetallic phase should have a finite exten-sion for small V /t , this constrains the putative QAH phase ina small region around V /t (cid:39) .In order to estimate the finite-size effects, we have per-formed extensive simulations on a large number of samples(data not shown). We can clearly identify some clusters onwhich the patterns are far from being perfect (i.e. have manydefects), which can be related either to lower symmetry, or topoorly two-dimensional aspect ratio. Indeed, when looking atclusters having C symmetry at least, we do observe nice cur-rent correlations with no or few defects at small V /t . There-fore, there are definitely correlations compatible with QAHphase, but one needs to investigate whether these are short-range features only or not.In order to do so, we have to measure some order parame-ter. We have chosen the current structure factor, i.e. the sumof current correlations with the signs taken from the QAH pat- tern: J q = (cid:88) (cid:104)(cid:104) ij (cid:105)(cid:105) ε ij (cid:104)J ij J i j (cid:105) (3)where the sum runs over all N b bonds on a given sublattice and the ε ij = ± correspond to the expected QAH orienta-tion (all currents clockwise or anticlockwise depending on thechoice of reference bond ( i j ) ). When using clusters withless symmetry, we average over the inequivalent directions ofthe reference bond. This is indeed a valid order parameterand long-range order requires that it diverges in the thermo-dynamic limit as N b ∼ N . Overall, we see in Fig. 25 a clearincrease of current correlations for moderate V , indicatingthat such an interaction indeed favours current formations, inagreement with mean-field prediction. There are of coursenon-trivial finite-size effects due to the fact that not all clus-ters have the same symmetry (see Appendix A).In the inset of Fig. 25, we plot this structure factor dividedby the number of terms N b for ≤ V /t ≤ where the QAHseems to be the strongest from our data, and also from themean-field estimate of the QAH region : . (cid:47) V (cid:47) . .For clusters without the K point, the behaviour of this cur-rent structure factor has some non-monotonic size dependencedue to the fact that clusters have different shapes, but its ab-solute value is quite small, and finite-size extrapolations arecompatible with a vanishing value. We also plot separatelythe data points for clusters having the K point (which have alarger signal due to the degenerate ground state at half fillingat V = 0 ) but if we focus on N = 24 and N = 42 which arethe best two-dimensional clusters of this kind (good aspect ra-tio and C symmetry at least), data are also compatible withthe absence of long-range order in the thermodynamic limit.Our results are in agreement with other recent numerical stud-ies but in contrast with another work. We hope that fu-ture progress in simulations of two-dimensional strongly cor-related systems will allow to confirm the absence of a QAHphase for this microscopic model, but will be able to stabiliseit in a larger parameter space. Indeed, it is obvious that thesemi-metallic state has a strong response in the QAH channel,so that if other instabilities could be prevented (using addi-tional interactions to frustrate them), there remains an inter-esting possibility to stabilise QAH phase in a truly long-rangeordered fashion.
V. CONCLUSION
We have investigated a model of spinless fermions on hon-eycomb lattice with (generally competing) density-density in-teractions using extensive numerical exact diagonalizations onvarious clusters. Based on a careful analysis, our main resultis summarized in the global phase diagram in Fig. 3. On theone hand, we have clarified the phase diagram in the strongcoupling regime ( | V , | (cid:29) t ) using some non-trivial resultson the classical model and some quantum order-by-disorderarguments. On the other hand, we have used numerical simu-lations of the full quantum model on finite clusters. As a re-sult, for intermediate or large interactions, we have provided3 FIG. 22. (Color online) Current-current correlations on a given sublattice between a reference bond (black) and other bonds for V /t = 1 , and (from left to right) with fixed V = 0 on N = 24 sample. Periodic boundary conditions are used. Blue and red correspond topositive/negative correlations according to the orientations expected in the QAH phase (taken here to be clockwise on all hexagons). Width isproportional to the correlation.FIG. 23. (Color online) Same as Fig. 22 for V = 0 and V /t = 1 on various clusters having C symmetry. From left to right: N = 26 , , and . FIG. 24. (Color online) Same as Fig. 23 for V = 0 and V /t = 2 . evidence for several crystalline phases: Néel CDW, CM, pla-quette/Kekulé, Néel domain wall crystal, the stripy* regionand the zigzag phase, some of which had not been consideredyet. When considering the role of V > interaction alone, ithad been suggested in the literature that it could stabilizea topological phase, which has triggered an intense activity.We have measured the corresponding (charge) current corre-4 V / t -0.1-0.0500.050.10.150.2 c u rr e n t s t r u c t u r e f ac t o r J q N=26N=28N=30bN=32N=34N=38N=40 J q / N b V =1V =1.5V =2 FIG. 25. (Color online) Current structure factor J q as a function of V /t for several lattices. Inset: scaling of this structure factor dividedby the number of bonds N b = 3 N/ − vs /N for several V /t .Open symbols are used for clusters without the K point. Hatchedsymbols are used for the clusters N = 24 and N = 42 which havethe K point as well as C symmetry at least. All data points arecompatible with a vanishing signal in the thermodynamic limit. lations, and while there is definitely some increase due to the V interaction, our finite-size analysis indicate that the topo-logical phase is not present in the thermodynamic limit, andthere is a direct transition between the semi-metallic phase anda charge-modulated one. However, it is very interesting to no-tice that V interaction is able to create short-range currentcorrelations, so that some additional ingredients could possi-bly stabilize it. Future work along these lines should also con-sider the spinful case where a topological spontaneous quan-tum spin hall phase with spin currents has been proposed. As a last remark, let us remind that other models are knownto exhibit Dirac cones and thus could potentially also hostsuch a topological phase: for instance, on the square lattice π -flux model at half-filling, a topological insulating phaseis predicted using mean-field technique but numerical studyhave not confirmed its stability. It is also noteworthy that ageneralized tight-binding model on the kagomé lattice alwayscontain Dirac cones at some particular filling , which leavesroom to investigate the exciting possibility to stabilize a topo-logical phase with interactions. Note added:
While finalizing this manuscript we becameaware of concurrent iDMRG work with consistent resultsregarding the presence and absence of the various phases inthe repulsive regime. ACKNOWLEDGMENTS
We thank M. Daghofer, A. Grushin, M. Hohenadler, J.Motruk and F. Pollmann for discussions. We thank F. Poll-mann and collaborators for sharing their consistent iDMRGresults prior to publication. This work was performed us-ing HPC resources from CALMIP and GENCI-IDRIS (Grant2014-050225) as well as LEO2/3 and VSC2/3. SC thanksInstitut Universitaire de France (IUF) for financial sup- port. AML was supported by the FWF (I-1310-N27/DFGFOR1807).
Appendix A: Lattice geometries
Since we are using several kinds of lattices, we give theirdefinitions in Table I using a and b translations to define thetorus with periodic boundary conditions. Name a b sym. group K (2) M (3) X (6)12 (3 , √
3) (3 , −√ Z no yes (1) no18 (6 ,
0) (3 , √ C v yes no no24 (6 , √
3) (0 , √ C v yes yes yes (4 , √
3) ( − , √ Z no yes (1) no26 (7 , √
3) (2 , √ C no no no28 (7 , √
3) (0 , √ Z no yes (1) no30a (3 , √
3) ( − , √ C × C yes no no30b (5 , √
3) ( − , √ C × C no no no32 (8 ,
0) (4 , √ C v no yes no34 (9 , √
3) (2 , √ Z no no no36 (6 ,
0) (6 , √ C × C yes yes (1) yes (2)38 (8 , √
3) (1 , √ C no no no42 (9 , √
3) (3 , √ C yes no no54 (9 , √
3) (0 , √ C v yes no no72 (8 ,
0) (8 , √ C v yes yes yes TABLE I. Finite lattices studied in this work. Listed are the numberof spins N ; the basis vectors a , b in the plane; the symmetry pointgroup; presence of the K point; M point; X point. Appendix B: Counting of Ising Groundstates
In this section we provide a table with the dimensions of theclassical ground spaces at half filling for the clusters describedin Appendix A.5 N s − / − / 96 6 5’526 6 446 b - 12 - 12 - - - 3’042 ∅ - - 8’964 2 858 a - - 30’618 ∅ b - - 11’300 ∅ 630 42 1’764 42 1’520 - - 230’014 14 3’454 - - - - 5’208 - - 1’211’004 - 8’880 - - 59’711’598 - - - 186 - 186 8’705’390TABLE II. Ground state degeneracy in the classical limit t = 0 at half filling . ∅ denotes a ground state energy on this cluster that ishigher than on other clusters. - indicates that the ground space hasnot been explored for this instance. The column labels V and we usethe parametrization V = 1 − V for V > and V = 1 − | V | for V < . ∗ [email protected] † [email protected] R. E. Prange and S. M. Girvin, The Quantum Hall Effect, 2nd ed. (New York: Springer-Verlag, 1990). M. Z. Hasan and C. L. Kane, Rev. Mod. Phys. , 3045 (2010). X.-L. Qi and S.-C. Zhang, Rev. Mod. Phys. , 1057 (2011). 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). A. G. Grushin, E. V. Castro, A. Cortijo, F. de Juan, M. A. H.Vozmediano, and B. Valenzuela, Phys. Rev. B , 085136 (2013). T. Pereg-Barnea and G. Refael, Phys. Rev. B , 075127 (2012). A. Dauphin, M. Müller, and M. A. Martin-Delgado, Phys. Rev. A , 053618 (2012). T. Liu, B. Douçot, and K. Le Hur, ArXiv:1409.6237. 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´c, N. Chancellor, and I. F. Herbut, Phys. Rev. B , 165123(2014). E. F. Huffman and S. Chandrasekharan, Phys. Rev. B ,111101(R) (2014). L. Wang, P. Corboz, and M. Troyer, New Journal of Physics ,103008 (2014). Z.-X. Li, Y.-F. Jiang, and H. Yao, Phys. Rev. B , 241117(R)(2015). R. Shankar, Rev. Mod. Phys. , 129 (1994). V. N. Kotov, B. Uchoa, V. M. Pereira, F. Guinea, and A. H. Cas-tro Neto, Rev. Mod. Phys. , 1067 (2012). S. Ryu, C. Mudry, C.-Y. Hou, and C. Chamon, Phys. Rev. B ,205319 (2009). I. F. Herbut, V. Juriˇci´c, and B. Roy, Phys. Rev. B , 085116(2009). F. D. M. Haldane, Phys. Rev. Lett. , 2015 (1988). B. Roy and I. F. Herbut, Phys. Rev. B , 035429 (2010). S.-K. Jian, Y.-F. Jiang, and H. Yao, Phys. Rev. Lett. , 237001(2015). P. Corboz, S. Capponi, A. M. Läuchli, B. Bauer, and R. Orús,EPL , 27005 (2012). J. Chaloupka, G. Jackeli, and G. Khaliullin, Phys. Rev. Lett. ,097204 (2013). G. H. Wannier, Phys. Rev. , 357 (1950). C. Hotta and N. Furukawa, Phys. Rev. B , 193107 (2006). N = 24 : 18 states with 8 off-diagonal matrix elements each; N = 28 : 14 states with 10 off-diagonal matrix elements each; N = 36 : 6 states with 12 off-diagonal matrix elements each; N = 72 : 18 states with 24 off-diagonal matrix elements each. We use the full lattice symmetry, as well as particle-hole symme-try. For space symmetries, we label sectors using momentum andirreducible representation corresponding to the point-group com-patible with this momentum. For instance, on N = 24 clusterwhich possesses all symmetries, we have used C v group at the Γ point, Z group at the X point, C v group at the K point and Z × Z group at the M point. Note that the choice of the real space Fermi normal ordering canhave an impact on the absolute symmetry sectors of the finite sizeclusters, relative quantum numbers should however be invariant. A. F. Albuquerque, D. Schwandt, B. Hetényi, S. Capponi,M. Mambrini, and A. M. Läuchli, Phys. Rev. B , 024406(2011). We disagree with the advocated general 4-fold degeneracy men-tioned in Ref. 10. The four-fold degeneracy is actually an artifactof the N = 18 sample . N. Read and S. Sachdev, Phys. Rev. B , 4568 (1990). R. Moessner, S. L. Sondhi, and P. Chandra, Physical Review B , 144416 (2001). Z. Zhu, D. A. Huse, and S. R. White, Phys. Rev. Lett. , 127205(2013). R. Ganesh, J. van den Brink, and S. Nishimoto, Phys. Rev. Lett. , 127203 (2013). S.-S. Gong, D. N. Sheng, O. I. Motrunich, and M. P. A. Fisher,Phys. Rev. B , 165138 (2013). Our numerical current structure factor (see its definition below)are indeed the strongest when V = 0 on clusters N = 24 and N = 42 for instance. We sum over all bonds which do not share any common site withthe reference bond, i.e there are N b = 3 N/ − terms. Y. Jia, H. Guo, Z. Chen, S.-Q. Shen, and S. Feng, Phys. Rev. B , 075101 (2013). K. Asano and C. Hotta, Phys. Rev. B , 245125 (2011). J. Motruk, A. G. Grushin, F. de Juan, and F. Pollmann, Phys. Rev.B , 085147 (2015). P. Corboz, M. Lajkó, K. Penc, F. Mila, and A. M. Läuchli, Phys.Rev. B87