Interaction effects on topological phase transitions via numerically exact quantum Monte Carlo calculations
aa r X i v : . [ c ond - m a t . s t r- e l ] J un Interaction effects on topological phase transitions via numerically exact quantumMonte Carlo calculations
Hsiang-Hsuan Hung, Victor Chua,
1, 2
Lei Wang, and Gregory A. Fiete Department of Physics, The University of Texas at Austin, Austin, Texas, 78712, USA Department of Physics, The University of Illinois at Urbana-Champaign, Illinois, 61801, USA Theoretische Physik, ETH Zurich, 8093 Zurich, Switzerland
We theoretically study topological phase transitions in four generalized versions of the Kane-Mele-Hubbard model with up to 2 × sites. All models are free of the fermion-sign problem allowingnumerically exact quantum Monte Carlo (QMC) calculations to be performed to extremely low tem-peratures. We numerically compute the Z invariant and spin Chern number C σ directly from thezero-frequency single-particle Green’s functions, and study the topological phase transitions drivenby the tight-binding parameters at different on-site interaction strengths. The Z invariant and spinChern number, which are complementary to each another, characterize the topological phases andidentify the critical points of topological phase transitions. Although the numerically determinedphase boundaries are nearly identical for different system sizes, we find strong system-size depen-dence of the spin Chern number, where quantized values are only expected upon approaching thethermodynamic limit. For the Hubbard models we considered, the QMC results show that corre-lation effects lead to shifts in the phase boundaries relative to those in the non-interacting limit,without any spontaneously symmetry breaking. The interaction-induced shift is non-perturbativein the interactions and cannot be captured within a “simple” self-consistent calculation either, suchas Hartree-Fock. Furthermore, our QMC calculations suggest that quantum fluctuations from inter-actions stabilize topological phases in systems where the one-body terms preserve the D symmetryof the lattice, and destabilize topological phases when the one-body terms break the D symmetry. PACS numbers: 71.10.Fd,71.70.Ej
I. INTRODUCTION
Electron interactions in topological insulators havebeen a topic of intense research in recent years.
It iscrucial to go beyond the mean-field level to captureimportant fluctuation effects originating in the electroniccorrelations, as this can be decisive in determining thephase.
Exact diagonalization studies are inhibitedby significant finite-size effects, though they are un-biased by any particular ansatz in the way mean-fieldtheories are. In this work, we study correlation effectsin topological insulators by considering the Kane-Mele-Hubbard model and several variants with numericallyexact projective quantum Monte Carlo (QMC) calcula-tions. Due to their particle-hole symmetry, these modelsare free of the fermion minus sign problem, and QMCsimulations provide a great opportunity to study cor-relation effects in topological matter with an unbiasedtheoretical approach. We are able to accurately treatcorrelations , compute interacting topological invari-ants such as the Z invariant and spin Chern number,and identify topological phase transitions through thezero-frequency single-particle Green’s function. Wethus avoid complications associated with ground-stateevolution under twisted boundary conditions, wherenumerical computations of a manifold ground states isrequired, and potential subtleties regarding energy gapclosures must be addressed. Strictly speaking, the ground state is altered whentwisted boundary conditions are introduced. The exis-tence of a family of ground states smoothly connected to one and another, and a finite spectral gap are requiredfor the use of twisted boundary conditions. Meetingthese conditions can be especially challenging when ap-proaching a phase transition where excitation gaps canbecome very small. Moreover, the use of the twistedboundary conditions is not practical in large-scale simu-lations, such as QMC; hence current implementations oftwisted boundary conditions in interacting models havemainly focused on small sizes where exact diagonaliza-tion techniques have been used.
In addition, the ini-tial use of twisted boundary conditions for defining thespin Chern number introduced edge effects, which ini-tially cast doubt on its robustness as a bulk topologicalinvariant.
Therefore, it is worth revisiting topological phase tran-sitions from the point-of-view of spin Chern numbers,particularly in the context of systems with interactionsand finite-size effects present. We observe in our nu-merical QMC results a dichotomy in the role of on-siteHubbard interactions: Depending on the underlying lat-tice symmetry, they favor either a topological or trivialphase. Although our results are only limited to the classof models with particle-hole symmetry and S z conser-vation, they could be hints of a more general principleregarding the interplay between point-group symmetryand interactions.The remainder of this paper is organized as follows. InSection II, we introduce the Kane-Mele model and thefour variants of it that we study. We compare and con-trast the particular spatial symmetries exhibited by thesetoy models. In Section III, we follow up with a discussionon the time-reversal invariant topological Z index andthe spin Chern number with a focus on their numericalimplementation in the presence of interactions. Next inSection IV which is the main part of our work, we presentcomputations of topological indices in the presence inter-mediately strong interactions for the models introducedin Section II. This is followed up with discussions, in-terpretation, and speculation regarding these results inSection V. Then in Section VI we conclude with a sum-mary and conclusions. Also included in the appendicesare details regarding our quantum Monte Carlo method-ology and supporting numerical results on the spin Chernnumber. II. THE KANE-MELE MODEL AND SEVERALVARIANTS
The Kane-Mele (KM) model, an early model support-ing a Z topological insulator (TI) on the honeycomblattice, remains central to the study of interactioneffects in TI. The honeycomb lattice is a Bravais (trian-gular) lattice with a two-point basis (labeled as A and B ). The vectors connecting two neighboring sites are a , = ± √ a ˆ x + a ˆ y and a = − a ˆ y , where a is the latticeconstant between two nearest-neighbor sites as shown inFig. 1 (a); we set a = 1 hereafter. The Hamiltonianreads as H KM = − t X h i,j i X σ c † iσ c jσ + iλ SO X hh i,j ii X σ σc † iσ ν ij c jσ , where c † iσ ( c iσ ) creates (annihilates) a spin σ fermion onsite i and σ runs over ↑ and ↓ . Here, hh· · · ii denotessecond-neighbor terms given by vectors b = a − a , b = a − a and b = a − a describing the spin-orbit coupling λ SO . i = 1 , ,
3, and ν ij = 1 for counter-clockwise hopping and ν ij = 1 otherwise. For our numerical study, we consider four time-reversalsymmetric model Hamiltonians which are KM model-like: (i) the generalized Kane-Mele (GKM) model, theKM model with a spin-independent real-valued third-neighbor hopping term, (ii) the dimerized Kane-Mele(DKM), the KM model with a biased nearest-neighborhopping along the a direction, (iii) the t L -KM model,the KM model with 5-th neighbor hopping, and (iv) the t N -dimerized KM model, which is a hybrid of model (i)and (ii).All of the models are generalized versions of the KMmodels, and at half-filling, they preserve the particle-hole symmetry. In the non-interacting limit, they host atopological-insulator/trivial insulator phase transition bytuning tight-binding parameters. However, there existscrucial differences among the models. The GKM model(i) and the t L -KM model (iii) preserve the six-fold rota-tion or C symmetry of the honeycomb lattice, whereasthe DKM model (ii) and the t N -dimerized model (iv)explicitly break it down to C with the bias in the a M M K K M (b) a a (a) a FIG. 1. (Color online) (a) The honeycomb lattice and theunderlying vectors a , , . (b) The first Brillouin zone of thehoneycomb lattice. The Dirac points are K , = ( ± π √ a , ,
0) and M , = ( ± π √ a , π a ) and M =(0 , π a ). direction. In the following we will mainly focus on theGKM and DKM models.The wallpaper or 2D space group of the honeycomblattice is p6m which is symmorphic and has D as itspoint group. The different models considered are meantto represent different modifications of the bare KM modelsuch that the D symmetry is either preserved or bro-ken but which nevertheless exhibits a topological insu-lator phase transition. However, time-reversal, spin- S z ,inversion, and particle hole symmetry are preserved inall these models. More precisely these models exhibitthe Quantum Spin Hall (QSH) phase when topologicallynon-trivial under the Z classification of 2D time-reversalsymmetric topological insulators. It is well known thatthe KM model includes spin-orbit coupling in the formof spin-dependent second nearest neighbor hopping thatfavors the topological insulator phase. We shall see inmodels (i-iv) that hoppings which are additions to thestandard KM model will, at sufficient strengths, over-come this tendency of the spin-orbit coupling and stabi-lize the trivial phase without breaking any symmetries.Resuming our discussion on symmetry, recall that D ∼ = D × Z ( i )2 , where the ( i ) superscript in Z ( i )2 denotes2D inversion about the center of the hexagon. Further-more D ∼ = C ⋊ Z ( m )2 where the ( m ) denotes reflectionabout a vertical mirror plane i.e. Z ( m )2 ≡ σ v in Sch¨onfliesnotation. Moreover, in two dimensions inversion is iso-morphic to rotation by 180 ◦ or Z ( i )2 ∼ = C . Models (i-iv) are selected to maintain Z ( i )2 inversion symmetry butmay either preserve or break D down to Z ( m )2 or com-pletely. A further essential property that is common tothese models is the absence of QMC sign problems intheir respective Hubbard model incarnations, which areobtained by the inclusion of an on-site Hubbard interac-tion.Specializing to two dimensions, the seminal works ofKane and Mele , and later Bernevig et. al. , Schny-der et. al . , Kitaev, and Qi et. al. showed that withonly time-reversal symmetry the non-interacting topolog-ical phases are classified by a Z invariant, which is alsogeneralized to three dimensions by Fu and Kane , andMoore and Balents . The physical content of this binarytopological index is that it enumerates the number parityof Kramers pairs of gapless edge modes at a boundary ofthe system with the vacuum–at least for non-interactinggapped band insulators. With the absence of S z mixing,the non-interacting occupied bands may be further cat-egorized by their S z polarization, and each spin speciesis topologically non-trivial carrying a non-zero integralChern number; i.e . C σ = 0 , σ = ↑ , ↓ . However, due totime-reversal symmetry, C ↑ + C ↓ = 0, and hence we donot expect an Integer Quantum Hall effect. Nevertheless,if the spin Chern number defined by C spin = C ↑ − C ↓ C ↑ = − C ↓ (1)is odd and non-zero, then the ground state of filled bandsare in the non-trivial topological insulator phase or the Z odd phase. The non-mixing of S z sectors, designatesthis non-trivial phase as being the QSH phase where onthe edge, an odd number of helical edge states carry-ing S z -current persists so long as time-reversal symmetryand the bulk band gap are preserved. As was stated, allthe models we have considered will either be trivial orin the QSH phase in the non-interacting limits. In thiswork, we will demonstrate that the classification by C spin will not only be applicable in the non-interacting limit,but can also be extended to finite interaction where ourmain interests lie. It is also clear that under this classifi-cation, a topological phase transition between even andodd values of C spin must proceed by an odd variation∆ C spin ∈ Z + 1.It will be highlighted in the upcoming sections thatcrystal symmetry will play a crucial role in the nature ofsuch topological transitions. (i) Generalized Kane-Mele model We start with the GKM model previously introducedin Ref.[33] whose Hamiltonian is given by H GKM = − t X h i,j i X σ c † iσ c jσ + iλ SO X hh i,j ii X σ σc † iσ ν ij c jσ − t N X hhh i,j iii X σ c † iσ c jσ , (2)where hhh· · · iii third-neighbor terms, and the vectorsconnecting third-neighbor terms are given by c i = a i + b i .At t N = 0 and finite spin-orbit coupling λ SO , themodel Hamiltonian Eq. (2) is reduced to the Kane-Mele model, which is a two-dimensional Z topolog-ical insulator. Like the KM model, the GKM modelis invariant under both the time-reversal symmetry andthe honeycomb space group p6m symmetry with its pointgroup D . t -4-2024 M M K K E /t M (b) t c3N = t/3 trivial insulator (a) Z TI FIG. 2. (Color online) (a) Schematic phase diagram of theGKM model. (b) The noninteracting band structure of theGKM model at t c N = t (here using λ SO = 0 . t ). The pre-sented momenta are chosen along the path depicted as bluearrows in Fig. 1 (b). The gap closes at three TRIM points: M , , , instead of the Dirac point K , . In the large t N limit of Eq. (2), the system isa trivial insulator, implying the GKM model under-goes a symmetry-preserving topological phase transi-tion as a function of t N . The GKM model can berecast as H GKM = P k Φ † k H GKM k Φ k , where Φ k =( c A ↑ k , c B ↑ k , c A ↓ k , c B ↓ k ) is a 4-component spinor and H k reads H GKM k = f ( k ) h ( k ) h ∗ ( k ) − f ( k ) − f ( k ) h ( k ) h ∗ ( k ) f ( k ) , where h ( k ) = g ( k ) − t N P i e i k · c i ; g ( k ) = − t P i e i k · a i ,and f ( k ) = 2 λ SO P i sin ( k · b i ); note that a i , b i and c i are real-space vectors to describe nearest, second andthird-neighbor hoppings. For most t N values, the GKMmodel is gapped. However a simple analysis of the dis-persion will show a gap closure at t c N = t independentof λ SO , and permits a change in the topological Z in-dex. The schematic phase diagram is shown in Fig. 2(a). For t N < t c N , the system is a Z TI, whereas for t N > t c N , the system is a trivial insulator. Thus, thereexists a topological phase transition at t N = t .The non-interacting band structure of the GKM modelis depicted in Fig. 2 (b). The chosen momenta are alongthe high symmetry momentum lines as indicated arrowsin Fig. 1 (b). Of particular note is the C symmetryof the dispersion that relates the three high-symmetry M , , points which are also inversion symmetric points.These become three Dirac points at the critical topolog-ical phase transition; gaps are opened at K , due tothe spin-orbit coupling λ SO . One should contrast thetopological phase transition with the KM model, whichinvolves two Dirac nodes at the K and K points. Notethat at t N = t/ λ = 0, the bands still touch at -1.0-0.50.00.51.0 /a E /t k x (a) t = 0.2t /a -1.0-0.50.00.51.0 /a /a (b) t = 0.45t k x FIG. 3. (Color online) The edge spectra for the noninter-acting GKM model at λ SO = 0 . t and (a) t N = 0 . t , a Z topological insulator and (b) t N = 0 . t , a trivial insulator.A ribbon geometry is used with an armchair edge for periodicboundary conditions (length, x -direction) and zigzag edge foropen boundary conditions (width). the M , , points as well as at the K , points, yielding5 Dirac points.Displayed in Fig. 3 are results of a band structure com-putation in a strip geometry of the GKM model demon-strating the existence of edge states with energies thattraverses the bulk energy gap. By counting the number ofKramers pairs of edge states per edge, we note thatthe GKM model undergoes a topological phase transi-tion as a function of t N , whereby an odd number ofKramers pairs characteristic of a Z topological insula-tor phase turns into an even number characteristic of atopologically trivial phase. Although the two pairs ofKramers helical states in Fig. 3(b) shows that the GKMmodel at t > t c N is a Z trivial insulator, for each spinflavor the spin Chern number, C σ is even and nonzero.Namely, | C σ | = 2 = 0 implying nontrivial edge states, albeit ones not protected by time-reversal symmetry. Inaddition, since the bulk gap of the GKM model closesat the three time-reversal invariant momenta: M , M and M [in Fig. 2 (b)], we expect that the spin Chernnumber will suffer an odd variation | ∆ C σ | = 3 signal-ing a topological transition from the C σ = ± C σ = ∓ C crystal symmetry, since M , M , M transformamongst themselves in a non-trivial irreducible represen-tation of the unbroken C symmetry. It must be men-tioned that a topological transition with ∆ C σ odd mayalso involve an even number of Dirac cones as is the caseas in the KM model where a staggered AB -site potentialcompetes with spin-orbit coupling. In this instance, how-ever, the inversion symmetry is broken from the pristinegraphene band structure allowing the transfer of an oddamount of Chern number. This is because under bro-ken inversion symmetry, the K and K are not required -4-2024 M M K K E /t M (b) t d t cd = 2t trivial insulator (a) Z TI FIG. 4. (Color online) (a) Schematic phase diagram of theDKM model. (b) The noninteracting band structure of theDKM model at t cd = 2 t (here using λ SO = 0 . t ). The pre-sented momenta are chosen along the path depicted as bluearrows in Fig. 1 (b). The gap closes at one TRIM point: M ,instead of the Dirac point K , . Compare with Fig.2. to contribute equally in the transfer of Chern number.The GKM model, however, differs by always maintain-ing inversion symmetry and in fact the M points areindividually inversion symmetric points. Lastly, fromthe perspective of the spin Chern numbers, both phasesare nontrivial–they exhibit edge robust states so long as S z symmetry is preserved, and are classified by the spinChern number. (ii) Dimerized Kane-Mele model The second model we consider, the DKM model is ex-pressed by the following Hamiltonian H DKM = − X h i,j i X σ t ij c † iσ c jσ + iλ SO X hh i,j ii X σ σc † iσ ν ij c jσ , (3)where t ij = t d when r j = r i + a , whereas t ij = t oth-erwise. One can recast the Hamiltonian as H DKM = P k Φ † k H DKM k Φ k , where H k is H DKM k = f ( k ) h ′ ( k ) h ′ ( k ) ∗ − f ( k ) − f ( k ) h ′ ( k ) h ′ ( k ) ∗ f ( k ) , where h ′ ( k ) = − t d e i k · a − t P i =2 , e i k · a i . At t d = t , Eq.(3) is reduced to the KM model. This kinetic Hamilto-nian explicitly breaks the C subgroup of D resultingin the point group Z ( m )2 × Z ( i )2 which is a mirror reflec-tion perpendicular to a and inversion or 180 ◦ rotation.A schematic of its phase diagram and band structure atthe critical point are shown in Fig. 4.A topological phase transition will occur by tuning t d to twice the nearest-neighbor hopping. In this instance -2-1012 /a E /t k x (a) t d = 1.8t /a -2-1012 /a /a (b) t d = 2.2t k x FIG. 5. (Color online) The edge spectra for the noninteract-ing DKM model at (a) t d = 1 . t and (b) t d = 2 . t , for a Z topological insulator and a trivial insulator, respectively. λ SO = 0 . t is used. The anisotropic hopping t d is introducedalong the zigzag direction. the conduction and valence bands touch at a single Diraccone at the M point when t cd = 2 t (again independent ofthe value of λ SO ). This critical point separates a trivialand the topological insulator phase, as shown in Fig.4. In Fig. 5 we show the band structure in a strip geometryfor the DKM model. The topological phase transition inthe DKM model, however, is different from the one in theGKM model as noted by the absence of any helical edgemodes on the trivial insulator side [shown in Fig. 5(b)].Thus, the trivial insulator phase ( t d > t cd ) has zero spinChern number and its variation is | ∆ C σ | = 1 during thetopological phase transition.From the symmetry perspective, the transition in theDKM model greatly differs from the GKM model since C is completely broken leaving only mirror and inversion Z ( m )2 × Z ( i )2 symmetries of the original D point group.Besides the trivial Γ point, the M point – where the single critical Dirac cone appears – is the only inversionsymmetric point which also respects the residual mirrorsymmetry. Thus, the topological phase transition pro-ceeds as a unit change of spin Chern number and hencethe topological Z index. In summary, we see that atleast in the non-interacting limit, point group symmetrycan greatly influence the form of the electronic structureof the critical point straddling a QSH phase and trivialphase. (iii) t L Kane-Mele model
The next model on our list is the t L -KM model whichsupplements the KM model with a four-lattice-constant-range hopping of strength t L . Similar to the t N term inthe GKM model, the tight-binding parameter, t L , in the t L -KM model preserves the D point group symmetry of (a) t d t (b) FIG. 6. (Color online) (a) The t L -KM model and (b) The t N -dimerized KM model. the honeycomb lattice. The model Hamiltonian reads as H t L = − t X h i,j i X σ c † iσ c jσ + iλ SO X hh i,j ii X σ σc † iσ ν ij c jσ − t L X { i,j } X σ c † iσ c jσ , (4)where the first two terms describe the KM model, and inthe third term { i, j } denotes the real-valued hopping withthe distance of 4 a . The lattice structure is shown in Fig.6 (a). Similar to the GKM model, in the non-interactinglimit, there exists a topological phase transition from the Z topological insulator to the trivial insulator state. Inthis instance, the boundary is located at t L = t withthree M , M , M Dirac cones. For simplicity, we do notdiscuss the properties of the edge dispersion as they arequalitatively similar to the GKM model. (iv) t N -Dimerized Kane-Mele model The final model we consider is the t N -dimerized KMmodel which is constructed by the combination of onethird-neighbor hopping (instead of three) and the bonddimerization in the KM model. As shown in Fig. 6 (b),the solid blue lines denote the dimerized bonds with t d strength and the purple dotted lines denote the diagonal t N hopping. The Hamiltonian reads as H t N d = − X h i,j i X σ t ij c † iσ c jσ + iλ SO X hh i,j ii X σ σc † iσ ν ij c jσ − t N X { i,j } = c X σ c † iσ c jσ , (5)where t ij = t d if r j = r i + a ; otherwise t ij = t . Thefirst two terms give the DKM model. The real-valueddiagonal t N hopping is selected along c = √ a ˆ x − a ˆ y .The simultaneous presence of the dimerized bonds and t N bonds breaks the Z ( m )2 mirror reflection and C rota-tional symmetry of D ∼ = ( C ⋊Z ( m )2 ) × Z ( i )2 . Thus the D subgroup is completely broken, however the Z ( i )2 inversionsymmetry is still respected. There also exists a topologi-cal phase transition between the Z topological insulatorstate and the trivial state, and the non-interacting criti-cal condition can be determined to be t N + t d = 2 t whichis again independent of the value of λ SO . At the topo-logical phase boundary, t cd = 2 t − t N , the bands form asingle Dirac cone at M , which hints that the spin Chernnumber has changed by | ∆ C σ | = 1 during the topologicalphase transition. III. NUMERICAL EVALUATION OFTOPOLOGICAL INDICES
For each generalization or variant of the Kane-Melemodel, an on-site Hubbard interaction will be added,and interacting phase diagrams containing the trivial andtopological phases are obtained via QMC simulations.But first we review and discuss the quality of numeri-cally computed topological indices of the finite clustersin the non-interacting limit.The first and most important topological index is the Z invariant of a two-dimensional non-interacting topo-logical insulator. When inversion symmetry is present,the noninteracting Z invariant is determined as ( − ν = Y k i ∈ TRIM Y m ξ m ( k i ) , (6)where ξ m ( k i ) is the parity of 2 m -th occupied Hamilto-nian eigenstate at the time reversal invariant momentum(TRIM); in the KM models, they are Γ and M , , as de-picted in Fig. 1 (b). (2 m − m -th states sharethe same parity and are a Kramers pair, and thereforeshould only be counted once in the determination of thetopological invariant. The time-reversal invariant topo-logical insulator phase is stable in the weakly interactinglimit and the Z index is also well defined in the case ofweak-interactions. It may be obtained conveniently fromthe zero-frequency single-particle Green’s function. Specialized to the presence of inversion symmetry, theFu-Kane expression Eq. (6) with interactions general-izes to ( − ν = Y k i ∈ TRIM ˜ η ( k i ) , (7)where ˜ η ( k i ) are the parity eigenvalues (one per Kramer’spair) of the R-zero eigenstates of the zero-frequencyGreen’s functions at TRIM. R- and L-zeros are termsused to refer to eigenfunctions of the zero-frequency sin-gle particle Green’s function G σ ( iω = 0 , k ). Eigenfunc-tions | v nkσ i with band index n , spin σ , crystal momen-tum k and eigenvalue such that G σ ( iω = 0 , k ) | v nkσ i = λ nkσ | v nkσ i (8)are called R-zeros when λ nkσ >
0, and L-zeros when λ nkσ < In the non-interacting limit, R-zeros corre-spond to occupied states below the Fermi-energy.(Notethat Eq.(8) is often expressed in terms of the inverse
Green’s function. Since our system consists of 2 x 2 matricies for each spin value, we can equivalently expressthe formula directly in terms of the Green’s function.One need only exercise care in the meaning of L-zeros,R-zeros, and singularities of the Greens functions.) Thesingular case G σ ( iω, k ) ∼ /ω as ω → G σ ( iω = 0 , k ) = 0 or λ nkσ = 0 isan indication of the onset of an interaction driven metal-insulator transition in the Brinkman-Rice sense. Eithera pole singularity or zero of G σ ( iω = 0 , k ) may inducea change in the topological index. This expression forthe ( − ν index is immensely useful and convenient indetermining the topological phase of an interacting time-reversal invariant system, but is however limited to theinversion symmetric situations.The second topological index that will concern us isthe spin Chern number defined by Eq. (1) in the QSHcontext. The Chern numbers C σ of the S z projectedbands are expressed in terms of one particle spectralprojectors as C σ = i π Z d k ǫ µν Tr [ P σ ( k ) ∂ µ P σ ( k ) ∂ ν P σ ( k )] , (9)where P σ ( k ) = P n | v nkσ ih v nkσ | is the single particlespectral projector onto R-zero states. Here we haveused the Berry curvature in k − space interpretation of thespin-Chern number as opposed to the original formu-lation in terms of twisted boundary conditions. Wewill often refer to C σ as the spin Chern number as well,since in the case of S z conservation–which applies to allcases considered in this work–the spin Chern number isproportional C σ , up to a sign determined by convention.More importantly, it is the parity–even or oddness–of C σ ,not its sign, that determines the time-reversal topological Z index. When S z is not a good quantum number, theexpression Eq. (9) for the spin Chern number definedin the thermodynamic limit – that is without twistedboundary conditions – may be generalized to the casewithout S z conservation rigorously. Even though wewill not consider these situations in this work, we wouldlike to point out that it is certainly possible to generalizeour numerical methods for the computation of the spinChern number and hence the Z index for interactingsystems where S z is not conserved.The inversion symmetric invariant of Eq. (7) and the S z conserving spin Chern number of Eq. (9) exhibit acomplementary relationship. The former is only appli-cable to inversion symmetric Hamiltonians, but does notrequire the S z conservation. The latter, however, doesnot require inversion symmetry but is nevertheless con-veniently computed only for S z conserving Hamiltonians( e.g . the staggered potential cases ). Moreover, the spinChern number, which may be any integer value in thethermodynamic limit, carries more information and thusa finer topological classification than the Z index, andcan remain quantized even when time-reversal symme-try is broken. However there is an obvious bias towardsfavoring Eq. (7) because by construction it is always inte-gral in finite sized systems. Whilst Eq. (9) will in generalyield non-integral values in finite-size systems where theBerry curvature over the BZ is no longer smooth. Thepracticalities of numerically computing the spin Chernnumber and its sensitivities to finite system size will bethe subject of our next discussion.For an interacting system, we compute the zero-frequency single-particle Green’s function with QMC andthen determine its eigenvectors and eigenvalues. Thedetermination of the topological response of a systemby the zero-frequency Green’s function has been demon-strated in both the non-interacting and interacting limitin Ref.[41]. Both expressions (7) and (9) sidestep difficul-ties associated with using twisted boundary conditions, which requires multiple numerically expensive calcula-tions of a non-degenerate ground state. It is also inappli-cable when artificial edge degeneracies are encounteredand is usually only practical with exact diagonalization. The fact that both of the expressions and their inter-acting generalizations only rely on the zero-frequencysingle-particle Green’s functions, is very convenient sincemore sophisticated numerical simulation methods likeQMC and Dynamical Mean-field Theory (which cannotprovide ground state wave functions) can be implementedin determining the topological phases with interactions.The zero-frequency property also implies that numericalanalytical continuation does not need to be employed.The computation of Eq. (7) for finite-size interactingsystems has been previously performed in Refs. [32,33],and is straightforward. There is, however, a requirementthat only cluster shapes with BZs containing a TRIMpoints may be studied with this method.By contrast computing, Eq. (9) for interacting sys-tems is a relatively new enterprise and we describe ournumerical method for its computation in finite sizes inAppendix B. Our results for the non-interacting GKMmodel are shown in Fig. 7. It is evident that, even inthe non-interacting limit, this method of evaluating thespin Chern number suffers from strong finite-size effects,though it is free of any physical effects associated withtwisted boundary conditions. For small sizes, such as6 × ×
12, the resulting C σ are not well approx-imated by integer values, but the discontinuous jump atthe transition can still be detected by inspection. Onlyupon increasing the system size do the spin Chern num-bers and their discontinuous jumps converge to integers.A finite size scaling analysis is needed to extrapolate tothe thermodynamically limit.The finite-size effects present in these non-interactingcases will be important for interpreting the fully inter-acting system that we will turn to shortly. However, theprecise critical value of t c N = t is clearly seen from thedata for all system sizes to identify the topological phasetransition. In particular, for the 400 ×
400 cluster, thevariation is | ∆ C σ | = 3 across the topological phase tran-sition. As we mentioned earlier, this is consistent withthe gap closing at the M , , points. s p i n C he r n nu m be r t /t spin-up L= 6 L= 12 L= 18 L= 60 L=400 spin-down
L= 6 L= 12 L= 18 L= 60 L=400
FIG. 7. (Color online) The spin Chern numbers C σ vs t N forthe non-interacting GKM model at λ SO = 0 . t for differentcluster sizes. For t N < t , the system is a Z topologicalinsulator with | C σ | = 1. For t N > t , the system is a trivialinsulator but with | C σ | = 2. As was alluded to earlier, the source of the “non-integerness” of the Chern number is associated with theneed to approximate the k -space gradients of projector P σ from a finite set of points in the BZ, cf Eq. (B4). Thisalso implies that – rigorously speaking – an exact inte-ger value is only ever attainable in the thermodynamiclimit. This is an important implication since it meansthat topological classification as captured by the Chernnumber and its myriad generalization is an effect that isonly rigorously stably protected in the thermodynamiclimit. This is intuitively clear since, only in the ther-modynamic limit do the energy gaps between smoothlyconnected Bloch states collapse. The remaining finiteenergy gaps are the band gaps that are the source of thetopological protection of a ground state.The finite-size computations of Eq. (9) shown in Fig.7 with non-integral results are an honest reflection of thelimitation of working with finite-size clusters. We notethat an alternative method by Fukui et. al. sidestepsthis with a construction which always yields an integerresult. However this can be misleading since the accu-racy of the results requires a critical mesh size, whichFukui et. al . have estimated. Moreover, the integer val-ues obtained by their method excludes the possibility ofusing a finite-size scaling analysis to judge the conver-gence of their results and is a weakness in their method.These considerations also apply to the integral inversion Z invariant ( − ν which should and does fluctuate withcluster size: There is a shift in boundaries based on thisinvariant with changing cluster size and shape.The need for large cluster sizes, however, is com-pensated by the QMC method which provides accessto ground state correlators of cluster sizes significantlylarger than those manageable by exact diagonalization.Furthermore, it is not necessary to have a manifold ofground states (which also has to be of a sufficient den-sity) as is required with the twisted boundary conditionsmethod. Another important point to note from Fig. 7is that at a fixed cluster size, the tendency to an integervalue improves, the further away the tuning parameteris from the critical point. Furthermore the many-bodyexcitation gap remains open through that portion of pa-rameter space and the single particle Green’s function atzero frequency develops no poles or singularities, permit-ting us to invoke the principle of adiabatic continuity andinfer the thermodynamic value of the spin-Chern numberof the entire portion of phase space from the finite-sizescaling in the large t N limit and when t N = 0.With this information the sudden discontinuity in thenumerical Chern number can then be used to pinpointthe critical point. This is the general strategy that weemploy in mapping out a phase diagram of both non-interacting and interacting models. In the case of an in-teracting model phase diagram, we have one more tuningparameter which is the interaction strength itself. Thefree model can then be trivially classified and when ro-bust excitation gaps persists above the numerical groundstate, the principle of adiabatic continuity can be usedto reliably map out a phase diagram from sudden jumpsin the numerical Chern number. As a consistency check,we also compare the spin Chern number with the Z in-variant using Eq. (7) at various tight-binding parametersand interaction strength, as shown in Fig.8. IV. EFFECTS OF INTERACTIONS INHUBBARD MODEL EXTENSIONS OF THEKANE-MELE VARIANTS
We now come to the main part of the paper where wediscuss the effects of the Hubbard interaction on interact-ing topological insulator models. To obtain ground statecorrelators and capture correlation effects, we use pro-jective quantum Monte Carlo to study interactingvariants of the KM model, Eq. (2)-(5) to which an on-site Hubbard term is added, H → H + U P i ( n i − where U > n i is the number operator onsite i . In our QMC calculations, the number of sites is N = 2 × L , where L takes the values 6, 12 and 18. Thelargest system sizes are far beyond current capabilities forexact diagonalization studies, rendering our “unbiased”calculations on interaction effects in topological systemsimportant for going beyond mean-field approaches andthe severe finite-size limitations of exact diagonalizationstudies. The QMC methodology is described in detail inAppendix A. (i) Generalized Kane-Mele Hubbard model We first turn our attention to interaction effects in theGKM-Hubbard model, i.e. H GKM + U . Previously inRef. [33], the correlation effects were discovered to resultin a shift of the phase boundary that can be accurately -2-1012341 (- ) C (a) U=3t -2-1012341 (- ) t /t C (b) U=4t FIG. 8. (Color online) The Z invariant ( − ν (upper panels,for L = 12 only) and spin Chern number C σ (lower panels) forthe GKM-Hubbard model as a function of t N /t at (a) U = 3 t and (b) U = 4 t . The spin-up Chern numbers C σ are denotedby solid symbols. The spin-orbit coupling is λ SO = 0 . t andthe systems sizes are chosen as 6 × × ×
18 (green triangles). t c N = t (verticalblue line) is the critical point for the non-interacting limit. computed with QMC simulations: with increasing U , t c N shifts to larger values (compared to the vertical blue linein Fig.8). This behavior was identified by evaluating the Z invariant from exploiting the inversion symmetry ofthe single-particle Green’s function and using Eq. (7).The QMC results showed that at U = 4 t the topo-logical phase transition boundary moves into the trivialinsulator phase by roughly 10%. Thus, correlation stabi-lizes the topological phase in the GKM-Hubbard model.Here we demonstrate that the topological phase tran-sition can also be clearly identified by computing the spinChern numbers, as shown in the lower panels of Fig. 8.We chose intermediate interaction strengths U = 3 t and U = 4 t , which are below the threshold required to inducemagnetic ordering or any other symmetry breaking. Forcomparison, we also depict the Z invariant vs t N forthe 12 ×
12 cluster. For both U values, Figs.8 (a) and (b)show marked changes in the spin Chern number at thesame locations, as the Z invariant varies for the 12 × C ↑ = − C ↓ within tiny error bars, so long as t N is far from the phase boundary. When the value of t N is close to the topological phase transition, one needs toincrease the sampling to recover the relation. Similarto the non-interacting limit, the spin Chern numbers inFigs. 8 converge to integers only as t N is far away fromthe critical point. The spin-up Chern numbers in the Z regime is C ↑ ≃ +1 ( t N = 0 . t ) and turns to C ↑ ≃ − t N = 0 . t ), indi-cated in the 12 ×
12 and 18 ×
18 clusters. The significantvariation in C σ , | ∆ C σ | ≃
3, can be used to identify theparameter-driven topological phase transition at finite U in the finite-size clusters. Moreover by adiabatic continu-ity to the non-interacting limit, we can also confidentlyidentity the two phases between the topological phasetransition.Next we present the finite-size analysis for the spinChern number in the GKM-Hubbard model, shown inFig.9. Since only three different sizes, 6 ×
6, 12 ×
12 and18 ×
18 are available, we are unable to fully capture thescaling behavior. However, the trends are sufficient toinfer the value of the thermodynamic spin Chern num-ber. On the other hand, the judgement can be also ar-rived at by the principle of adiabatic continuity to thenon-interacting limit of the GKM model. We tentatively C (a) t =0 t =0.25t t =0.33t t =0.34t t =0.45t t =0.6t (b) FIG. 9. (Color online) The tentative finite-size analysis of thespin-up Chern number C ↑ of the GKM-Hubbard model at (a) U = 3 t and (b) U = 4 t vs 1 /L . consider the Chern number scaling as 1 /L in Fig. 9 (or1 /L , not shown here); showing that as the value of t N isfar away from the critical point, the spin Chern numbersextrapolate well to whole integers, C ↑ = 1 or C ↑ = − C σ = 1 and C σ = − t N = 0 . t and0 . t ), the scaling analysis is less reliable and one needsbigger sizes to determine the behavior. However, it isstill helpful in determining the location of the topolog-ical phase transition. In Fig. 9 (a), we can recognizethat at t N = 0 . t the spin Chern number shows a dropwith increasing system size; thus it is a trivial state. Bycontrast, Fig. 9 (b) shows that the spin Chern numberat t N = 0 .
34 does not show a clear drop, suggestingthat it is still in the topological insulator regime. Thus,interactions stabilize the topological phase in the GKM-Hubbard model. Note that for t N values away from the critical value, the finite size scaling behavior is muchclearer in terms of how the thermodynamic limit is ap-proached.Although the values of the spin Chern number sufferfrom strong finite-size effects, the topological phase tran-sition boundary determined by the topological invariantin the GKM-Hubbard model has weak finite-size depen-dence. For U = 3 t , on the 6 ×
6, 12 ×
12, and 18 ×
18 clus-ters, t c N = 0 . t , 0 . t and 0 . t , respectively. For U = 4 t , t c N are 0 . t , 0 . t and 0 . t , respectively,suggesting the spin Chern number is a reliable meansto detect topological phase transitions in interacting sys-tems.These interaction effects that cause the critical bound-ary in phase space to shift must originate from the dy-namical quantum fluctuations, since the Hartree-Fockmean-field theory is unable to capture any phase bound-ary shift (for the U values we consider below the mag-netic phase transition). We were not able to developa perturbative argument for this shift, either. (ii) Dimerized Kane-Mele Hubbard model
We next turn to the DKM-Hubbard model: H = H DKM + U . Recall that at U = 0, the critical point oc-curs at t cd = 2 t and is independent of value of λ SO . Sim-ilar to the GKM-Hubbard model, correlation effects in-duce a shift of the phase boundary, but the critical valueof t d moves towards (into) the topological phase. In otherwords, correlation destabilizes the topological insulatorphase–a behavior opposite to the GKM-Hubbard model.With finite interactions at U = 2 t and λ SO = 0 . t , t cd isdetermined within 1 .
94 and 1 .
96 by observing the Green’sfunction behavior and the Z topological invariant. Here we also employ the QMC combined with the com-putation of the spin Chern number using Eq. (9) andthe Z index using Eq. (7). Likewise, the values of U weconsidered were below the magnetic transition. Figs. 10(a) and (b) show C ↑ and the Z index vs t d for U = 2 t and U = 4 t , respectively. In the 12 ×
12 cluster (redcircles), we can see that the spin Chern number jumpsat t d = 1 . t for U = 2 t and t d = 1 . t for U = 4 t ,and, simultaneously, the value of the Z index turnsfrom ( − ν = − ×
12 cluster, one sees that the topolog-ical phase transition occurs between the | C σ | = 1 stateto the | C σ | = 0 state and a variation | ∆ C σ | ≈ L = 12 cluster, t cd at U = 2 t is esti-mated to be 1 . t − . t , whereas at U = 4 t it lies within1 . t − . t . The critical point has shifted by roughly25%. The DKM-Hubbard model also has weaker finite-size effects on the topological phase boundaries. For the L = 6 cluster, t cd s are estimated around 1 . t − . t and1 . t − . t for U = 2 t and 4 t , respectively. The compar-0 (- ) C (a) U=2t (- ) t d /t C (b) U=4t FIG. 10. (Color online) The Z invariant ( − ν (upper panels,for L = 12 only) and C σ (lower panels) for the DKM-Hubbardmodel as a function of t d /t at (a) U = 2 t and (b) U = 4 t . λ SO = 0 . t and the systems sizes are chosen as 6 × ×
12 (red circles). At t d = t , the systemreduces to the standard KM model and t d = 2 t (vertical blueline) is the critical point for the non-interacting limit. For thesake of clarity, only the data for C ↑ is presented here. ison of the results for the two cluster sizes show similarlocations of the topological phase boundaries, thus sug-gesting weak finite-size effect on the critical points. (iii) t L - and t N -Dimerized Kane-Mele Hubbardmodels Lastly, we present QMC results for the on-site Hub-bard models of the t L -Kane-Mele model Eq. (4) andthe t N -dimerized KM model Eq. (5). These two mod-els represent polar opposites with regard to their non-interacting hopping Hamiltonians. The former like theGKM model preserves the full D point group, whilst thelater breaks it down almost completely to just the inver-sion subgroup Z ( i )2 . Thus, the t N -dimerized KM modelis even less symmetric than the DKM model. The moti-vation for considering these other variants is to demon-strate more examples of interacting TI phases and therole crystal symmetry or lack thereof might play and helpcontrast the different outcomes of explicitly breaking orpreserving the crystal symmetry of the underlying KMmodel. Given that for the hopping Hamiltonians thatwe have set out to study, the crystal symmetry alreadygreatly influences the low-energy character of the critical Z i n v a r i an t (- ) t L /t L=6 U=4 QMC L=12 U=4 QMC (a) -0.1 0.0 0.1 0.2-101 Z i n v a r i an t (- ) t /t L=6 U=4 QMC L=12 U=4 QMC (b)
FIG. 11. (Color online) The Z invariants for interacting andnoninteracting cases in the (a) t L -KM model and (b) t N -dimerized KM model at t d = 1 . t . For reference purposes, thenon-interacting Z invariant were computed and presented asthe blue lines using L × L = 1200 × Z invariant by the QMC for U = 4 t on6 × ×
12 (red circles). All calculationsare preformed with λ SO = 0 . t . theory – such as deciding the number of Dirac cones –between the topological insulator phase and the normalinsulator phase, it is reasonable to expect that crystalsymmetry will have a significant role to play in shiftingphase boundaries.The QMC results of the correlation effects on thesetwo models are displayed in Fig. 11 where we still used λ SO = 0 . t and U = 4 to compare with Fig. 8 and Fig.10. For simplicity, we only show the Z invariants asa function of the tight-binding parameters: t L and t N in Eq. (4) and Eq. (5), respectively. Note that, nearthe critical point, ( − ν shows a poor approximation toan integer value (not ≃ ≃ − stabilizes the topological insulator state in the t L -KM model. Forthe 6 × t cL = 0 . t − . t , and for the 12 × t cL = 0 . t − . t . Like the GKM model, the t cL has weak finite-size effect on the phase boundaries.In the next panel, Fig. 11 (b) exhibits the interact-ing Z invariant against the t N parameter for the t N -dimerized KM model at t d = 1 . t . The non-interactinglimit, t c N = 0 . t , is indicated by the blue line. Turning1on interaction, the phase boundary moves towards (into)the topological state regime; thus correlation destabilizes the topological insulator phase. We have numerically ex-amined that with finite bond dimerization, the topolog-ical critical points are always pushed to the topologicalinsulator regime under correlation.Our observations of the effects of Hubbard-type inter-actions on these KM model variants show a systematicpattern: The stability of the topological insulator phaseas measured by its occupied volume in the phase dia-gram is diminished when more of the symmetries of the D point group of the lattice are explicitly broken by theHamiltonian. A posteriori, we elevate our observationsto a speculative conjecture of a principle: in the absenceof any spontaneous symmetry breaking, the Hubbard in-teraction will displace the critical line of the interactingtopological quantum phase transition in favor of the nor-mal insulator phase if fewer crystal point group symme-tries – but which must include inversion – are present inthe non-interacting portion of the tight-binding Hamil-tonian. The condition regarding spontaneous symme-try breaking excludes competition with magnetic anddensity-waves phases. This is important to state sinceat very strong coupling either phase gives way to theN´eel ordered phase. It is worth reiterating that due tothe special form of these Hamiltonians at half-filling, theQMC methodology employed is free of sign-problems andis essentially an exact method for finite clusters up tostatistical noise; which can always be systematically im-proved with greater sampling. In a related study of theplaquette KM model, which is not too dissimilar from theones we have considered, qualitatively consistent resultsare obtained. V. DISCUSSIONS
We make a few remarks regarding low energy effectivetheories and present some related speculations. As waspreviously mentioned, mean-field calculations at the levelof the Hartree-Fock approximation in the KM model and also our own computations for models Eqs. (2)-(5) are unable to demonstrate a continuous shift in thetopological quantum phase transition boundaries at weakcoupling, although a transition to a magnetic state doesoccur at strong coupling. We rationalize this by not-ing that the Hubbard U in two dimensions for a lowenergy effective critical field theory of gapless linearlydispersing Dirac fermions is irrelevant under scaling. Infact, as is well known even in the case of long-rangeCoulomb interactions in graphene, which is a archetypefor this variety of field theory, a Renormalization Group(RG) analysis also produces the conclusion that the Diracnodes are perturbatively stable, albeit with anomalousscaling dimensions due to quantum fluctuations. Thusthe phase boundary shifts – significantly observable onlyat relatively large U ∼ t – that we have observed in ourQMC exact computations are effects at intermediately strong interactions, which is beyond the weak-couplinglow energy-effective theory description. The implicationsare that the standard field theoretic RG computations atone loop order would be unreliable in capturing the in-termediately strong coupling physics of interest.Nevertheless, we speculate that an explanation of thedichotomous behavior of the phase boundary shifts mustinvolve the fact that there are a different number of Diraccones present at the critical topological transition point( three in the GKM and t L -KM models but one in theDKM and t N -dimerized KM models) and that this is themain influence of point group symmetry to the low en-ergy physics. It is tempting to relate our observations to alarge- N study of graphene with Hubbard interactions,but we are cautious and reluctant to since N = 1 , N . In spite of this, Functional Renor-malization Group (fRG) computations or more recentdimensional regularization d = 3 − ǫ, ǫ = 1 studies applied to the Hubbard graphene system have been en-couraging in describing physics near strong coupling. Wewill leave these very interesting lines of investigations forfuture work, as these computations are by no means triv-ial undertakings. VI. SUMMARY AND CONCLUSIONS
In this work we have analyzed variants of the archetyp-ical model of a time-reversal symmetric topological insu-lator, the Kane-Mele model. These generalized modelsEqs. (2)-(5) exhibit various space group symmetries ofthe honeycomb lattice on which they are formulated. Inthe non-interacting limit, all of the models exhibit a topo-logical phase transitions between a Z topological insu-lator phase and a normal insulating phase. By means ofthe unbiased QMC method, we further study the inter-acting variants of these model systems by including on-site Hubbard interactions. The projective determinantQMC method that we employ is free of sign problems(at half-filling which is the only filling considered) andis essentially exact up to statistical sampling noise. Theregime that interests us most is the intermediately strong U regime before magnetic order sets in. We demonstratethat the topological phase of our numerically exact in-teracting ground states can be ascertained by computingeither the Z invariant or the spin Chern number C σ viathe zero-frequency single-particle Green’s function. Thusour work is a numerical implementation of the theoreticalproposal of Refs. [41 and 42] for finite-size clusters usingreliably accurate QMC. The spin Chern number had notbeen previously computed with QMC, and we argue ontechnical grounds that it is complementary to the inver-sion symmetry based expression for the Z invariant.Accompanied with finite-size scaling analyses and adi-abatic continuity to non-interacting limits, – which wenumerically observe – we argue that the spin Chern num-ber is a robust classification method of interacting TI’s.Moreover the spin Chern number may be utilized in cir-2cumstances where inversion and even time-reversal sym-metry are absent, and may be generalized to the casewhere S z conservation is absent. Our numerically ex-act QMC results suggest that quantum fluctuations fromintermediately strong interactions can act to either sta-bilize or destabilize the topological phase, depending onwhether the hopping terms preserve or break lattice sym-metries (when U is less than the value which inducesmagnetism). Although admittedly bold and speculative,we conjecture a principle that in the situations wherespontaneous symmetry breaking phases are excluded, on-site Hubbard interactions will destabilize the TI phase inhoneycomb models when the lattice point group D isexplicitly broken down to a subgroup containing Z ( i )2 in-version by the tight-binding Hamiltonian, and stabilizedwhen the full D M. W. Young, S.-S. Lee, and C. Kallin,Phys. Rev. B , 125316 (2008). D. Pesin and L. Balents, Nat. Phys. , 376 (2010). M. Kargarian and G. A. Fiete,Phys. Rev. B , 085106 (2010). R. Li, J. Wang, X.-L. Qi, and S.-C. Zhang, Nat. Phys. ,284 (2010). M. Hohenadler, T. C. Lang, and F. F. Assaad, Phys. Rev.Lett. , 100403 (2011). M. Hohenadler, Z. Y. Meng, T. C. Lang, S. Wessel, A. Mu-ramatsu, and F. F. Assaad, Phys. Rev. B , 115132(2012). D. Zheng, G.-M. Zhang, and C. Wu, Phys. Rev. B ,205121 (2011). B. Swingle, M. Barkeshli, J. McGreevy, and T. Senthil,Phys. Rev. B , 195139 (2011). J. Maciejko, X.-L. Qi, A. Karch, and S.-C. Zhang,Phys. Rev. Lett. , 246809 (2010). M. Kargarian and G. A. Fiete,Phys. Rev. Lett. , 156403 (2013). M. Levin and A. Stern, Phys. Rev. B , 115131 (2012). M. Levin and A. Stern,Phys. Rev. Lett. , 196803 (2009). J. E. Moore, Y. Ran, and X.-G. Wen,Phys. Rev. Lett. , 186805 (2008). T. Neupert, L. Santos, S. Ryu, C. Chamon, and C. Mudry,Phys. Rev. B , 165107 (2011). A. R¨uegg and G. A. Fiete,Phys. Rev. Lett. , 046401 (2012). X. Wan, A. M. Turner, A. Vishwanath, and S. Y.Savrasov, Phys. Rev. B , 205101 (2011). M. Hohenadler and F. F. Assaad,J. Phys. Cond. Matt. , 143201 (2013). G. A. Fiete, V. Chua, M. Kargarian, R. Lundgren,A. Ruegg, J. Wen, and V. Zyuzin, Physica E , 845(2012). J. Maciejko, V. Chua, and G. A. Fiete, Phys. Rev. Lett. , 016404 (2014). J. Maciejko and A. R¨uegg,Phys. Rev. B , 241101 (2013). J. E. Moore, Nature , 194 (2010). M. Z. Hasan and C. L. Kane,Rev. Mod. Phys. , 3045 (2010). X.-L. Qi and S.-C. Zhang,Rev. Mod. Phys. , 1057 (2011). M. K¨onig, S. Wiedmann, C. Brune, A. Roth, H. Buh-mann, L. Molenkamp, X.-L. Qi, and S.-C. Zhang,Science , 766 (2007). A. Roth, C. Br¨une, H. Buhmann, L. W. Molenkamp,J. Maciejko, X.-L. Qi, and S.-C. Zhang,Science , 294 (2009). D. Hsieh, D. Qian, L. Wray, Y. Xia, Y. Hor, R. J. Cava,and M. Z. Hasan, Nature , 970 (2008). S. Raghu, X.-L. Qi, C. Honerkamp, and S.-C. Zhang,Phys. Rev. Lett. , 156401 (2008). Y. Zhang, Y. Ran, and A. Vishwanath,Phys. Rev. B , 245331 (2009). J. Wen, A. R¨uegg, C.-C. J. Wang, and G. A. Fiete,Phys. Rev. B , 075125 (2010). Q. Liu, H. Yao, and T. Ma,Phys. Rev. B , 045102 (2010). A. R¨uegg and G. A. Fiete,Phys. Rev. B , 201103 (2011). A. Go, W. Witczak-Krempa, G. S. Jeon, K. Park, andY. B. Kim, Phys. Rev. Lett. , 066401 (2012). H.-H. Hung, L. Wang, Z.-C. Gu, and G. A. Fiete, Phys.Rev. B , 121113 (2013). J. C. Budich, R. Thomale, G. Li, M. Laubach, and S.-C.Zhang, Phys. Rev. B , 201407 (2012). J. C. Budich, B. Trauzettel, and G. Sangiovanni,Phys. Rev. B , 235104 (2013). L. Wang, X. Dai, and X. C. Xie,Euro. Phys. Lett.. , 57001 (2012). T. Yoshida, R. Peters, S. Fujimoto, and N. Kawakami,Phys. Rev. B , 085134 (2013). C. N. Varney, K. Sun, M. Rigol, and V. Galitski, Phys.Rev. B , 115125 (2010). C. N. Varney, K. Sun, M. Rigol, and V. Galitski,Phys. Rev. B , 241105 (2011). T. C. Lang, A. M. Essin, V. Gurarie, and S. Wessel, Phys. Rev. B , 205101 (2013). Z. Wang and S.-C. Zhang, Phys. Rev. X , 031008 (2012). Z. Wang, X.-L. Qi, and S.-C. Zhang,Phys. Rev. B , 165126 (2012). Q. Niu, D. J. Thouless, and Y.-S. Wu, Phys. Rev. B ,3372 (1985). T. Fukui and Y. Hatsugai,Phys. Rev. B , 121403 (2007). L. Fu and C. L. Kane, Phys. Rev. B , 045302 (2007). E. Prodan, Phys. Rev. B , 125327 (2009). C. L. Kane and E. J. Mele,Phys. Rev. Lett. , 226801 (2005). C. L. Kane and E. J. Mele,Phys. Rev. Lett. , 146802 (2005). B. A. Bernevig, T. L. Hughes, and S.-C. Zhang, Science , 1757 (2006). A. P. Schnyder, S. Ryu, A. Furusaki, and A. W. Ludwig,in
American Institute of Physics Conference Series , Vol.1134 (2009) pp. 10–21. A. Kitaev, in
ADVANCES IN THEORETICALPHYSICS: Landau Memorial Conference , Vol. 1134(AIP Publishing, 2009) pp. 22–30. X.-L. Qi, T. L. Hughes, and S.-C. Zhang, Phys. Rev. B , 195424 (2008). L. Fu, C. L. Kane, and E. J. Mele, Phys. Rev. Lett. ,106803 (2007). J. E. Moore and L. Balents, Phys. Rev. B , 121306(2007). We have also specialized to the case where there is also onlya single valence and conduction band per-spin species. D. N. Sheng, Z. Y. Weng, L. Sheng, and F. D. M. Haldane,Phys. Rev. Lett. , 036808 (2006). L. Sheng, D. N. Sheng, C. S. Ting, and F. D. M. Haldane,Phys. Rev. Lett. , 136602 (2005). M. A. N. Ara´ujo, E. V. Castro, and P. D. Sacramento,Phys. Rev. B , 085109 (2013). Y. Hatsugai, M. Kohmoto, and Y.-S. Wu,Phys. Rev. B , 4898 (1996). C. Xu and J. E. Moore, Phys. Rev. B , 045322 (2006). See the supplemental material in Ref. [33]. V. Gurarie, Phys. Rev. B , 085426 (2011). S. R. Manmana, A. M. Essin, R. M. Noack, and V. Gu-rarie, Physical Review B , 205119 (2012). Y. You, Z. Wang, J. Oon, and C. Xu, ArXiv e-prints(2014), arXiv:1403.4938 [cond-mat.str-el]. W. Brinkman and T. Rice, Physical Review B , 4302(1970). J. E. Avron, R. Seiler, and B. Simon,Phys. Rev. Lett. , 51 (1983). D. N. Sheng, L. Balents, and Z. Wang,Phys. Rev. Lett. , 116802 (2003). E. Prodan, New J. Phys. , 065003 (2010). H.-H. Lai and H.-H. Hung,Phys. Rev. B , 165135 (2014). A. Georges, G. Kotliar, W. Krauth, and M. J. Rozenberg,Rev. Mod. Phys. , 13 (1996). T. Fukui, Y. Hatsugai, and H. Suzuki, J. Phys. Soc. Jpn. , 1674 (2005). G. Sugiyama and S. E. Koonin, Ann. Phys. , 1 (1986). S. Sorella, S. Baroni, R. Car, and M. Parrinello, Europhys.Lett. , 663 (1989). S. R. White, D. J. Scalapino, R. L. Sugar, E. Y. Loh, J. E.Gubernatis, and R. T. Scalettar, Phys. Rev. B , 506(1989). F. F. Assaad,
Quantum Monte Carlo methods on lattices:The determinantal approach in Quantum Simulations ofComplex Many-Body Systems: ¿From Theory to Algo-rithms, Lecture Notes (NIC Series Vol. , 2002). F. F. Assaad, AIP Conf. Proc. , 117 (2003). S. Rachel and K. Le Hur, Phys. Rev. B , 075106 (2010). W. Wu, S. Rachel, W.-M. Liu, and K. Le Hur, Phys. Rev.B , 205102 (2012). J. Gonz´alez, F. Guinea, and M. Vozmediano, NuclearPhysics B , 595 (1994). I. F. Herbut, Physical review letters , 146401 (2006). L. Janssen and I. F. Herbut, ArXiv e-prints (2014),arXiv:1402.6277 [cond-mat.str-el]. I. F. Herbut, V. Juriˇci´c, and O. Vafek, Physical Review B , 075432 (2009). F. F. Assaad and I. F. Herbut, Physical Review X ,031010 (2013). Z. Y. Meng, T. C. Lang, S. Wessel, F. F. Assaad, andA. Muramatsu, Nature , 847 (2010). H.-H. Hung,
Exotic quantum magnetism and superfluidityin optical lattices (PhD thesis, University of California, SanDiego, 2011). N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H.Teller, and E. Teller, J. Chem. Phys. , 1087 (1953). J. E. Hirsch, Phys. Rev. B , 4403 (1985). F. F. Assaad and M. Imada, J. Phys. Soc. Jpn. , 189(1996). Appendix A: Quantum Monte Carlo
The projective quantum Monte Carlo method (QMC)is given by projecting an arbitrary trivial wave function | ψ T i (requiring h ψ T | ψ i 6 = 0) onto the ground state wavefunction | ψ i . The expectation value of an ob-servable A is obtained by h A i = lim Θ →∞ h ψ T | e − Θ2 H Ae − Θ2 H | ψ T ih ψ T | e − Θ H | ψ T i , (A1)where Θ is the projective parameter. To carry out theprocedures numerically, we need to discretize the projec-tion operator e − Θ H into tiny time propagators e − ∆ τH with Θ = ∆ τ M : e − Θ H = ( e − ∆ τH ) M where M is thenumber of time slices and ∆ τ is chosen as a small num-ber. The first-order Suzuki-Trotter decomposition canfurther decompose e − ∆ τH as e − ∆ τH ≃ e − ∆ τH e − ∆ τH U , (A2)where H is the tight-binding Hamiltonian, which couldbe equation (2) and equation (3), for the GKM modeland the DKM model, respectively; H U = U P i ( n i − is the repulsive Hubbard on-site interaction; n i = P σ c † i,σ c i,σ . To represent e − ∆ τH U in terms of the single-particle basis, we need to implement the SU (2)-invariantHubbard-Stratonovich transformation e − ∆ τ U ( n i − = 14 X l = ± , ± γ ( l ) e i √ ∆ τ U η ( l )( n i − + O (∆ τ ) , (A3)4where γ ( ±
1) = 1 + √ / γ ( ±
2) = 1 − √ / η ( ±
1) = ± q − √
6) and η ( ±
2) = ± q √
6) are4-component auxiliary fields. In the current literature,∆ τ t = 0 .
05 and Θ t = 40 are used through the content. Implementing equation (A2) and (A3), H turns out tobe τ -dependent since H U is associated with the auxiliaryfield configuration η ( l i,τ ); then e − Θ H = Q Mτ =1 e − ∆ τH τ .The denominator of equation (A1) (named the projectorpartition function), is evaluated as follows h ψ T | e − Θ H | ψ T i = h ψ T | M Y τ =1 e − ∆ τH τ | ψ T i ∼ = h ψ T | M Y τ =1 e − ∆ τH e − ∆ τH U,τ | ψ T i = ( 14 ) MN X { l i,τ } ((cid:16) Y i,τ γ ( l i,τ ) (cid:17) Y σ Tr (cid:16) M Y τ =1 e − ∆ τ P i,j c † i,σ [ H σ ] ij c j,σ e i √ ∆ τ U η ( l i,τ )( n i,σ − ) (cid:17)) = ( 14 ) MN X { l i,τ } ((cid:16) Y i,τ γ ( l i,τ ) (cid:17) det (cid:16) O ↑ [ η ( l i,τ )] (cid:17) det (cid:16) O ↓ [ η ( l i,τ )] (cid:17)) , (A4)where P l i,τ runs over possible auxiliary configurations η ( l i,τ ), where i = 1 ∼ N are site indices and τ = 1 ∼ M are imaginary time indices; H σ is the matrix kernel of H with spin σ . Each time propagator e − ∆ τH e − ∆ τH U is a N × N matrix and Tr( Q τ e ··· ) = det( O σ ) representsthe trace over fermion degrees of freedom.Given such a N -site and M -time slice system, thesummation in the above equation has a degree of 4 NM ,and generally, it is impossible to consider all configura-tions. The auxiliary field configuration, {· · · η ( l i,τ ) · · · } ,however, can be determined by Monte Carlo importancesamplings. For simplicity, we used the Metropo-lis algorithm in this paper. The physical meaningof (cid:16) Q i,τ γ ( l i,τ ) (cid:17) Q σ = ↑ , ↓ det (cid:16) O σ [ η ( l i,τ )] (cid:17) is the prob-ability weight at the given auxiliary field configura-tion { η ( l i,τ ) } . When this term is proven positive-definitive, QMC simulations are free-sign and the re-sults are numerically exact. This is always true in thehalf-filling Kane-Mele-type model without Rashba spin-orbital coupling.
To obtain zero-frequency Green’s function, we firstevaluate the time-displaced Green’s function G ( r , τ ).The unequal-time Green’s function is defined as G σ ( τ, r i , r j ) = h ψ | c σ ( τ, r i ) c † σ ( r j ) | ψ i = h ψ | e τH c σ ( r i ) e − τH c † σ ( r j ) | ψ i . (A5)Then we perform the Fourier transform from real spaceto momentum space r → k , and imaginary time to theMatsubara frequency τ → iω , G σ ( iω, k ) = 1 β Z β dτ e iωτ N X r i ,r j e i k · ( r i − r j ) G σ ( τ, r i , r j ) . The zero-frequency is given setting iω = 0. To calculatethe spin Chern number C σ , however, we need the single-particle Green’s functions for all momentum points andimplement equation (9). This procedure is slightly dif-ferent from the approach to evaluate the Z index, for which only time-reversal invariant momentum points arerequired. For sign-free QMC simulations, one can ac-curately calculate the zero-frequency Green’s functionsin system sizes which are larger than the small clustersin an exact diagonalization and then evaluate the spinChern numbers using the projection operators equation(9). This approach is useful to identify different topolog-ical phases in the interacting level without using twistedboundary conditions.Note that for more generic cases, the Green’s functionsare a 4 × , i.e., G = (cid:18) G ↑↑ G ↑↓ G ↓↑ G ↓↓ (cid:19) , (A6)where G ↑↑ = G ↑ ( G ↓↓ = G ↓ ) as defined in equation (A5),and G ↑↓ = h c ↑ ( τ ) c †↓ i . Without Rashba spin-orbital cou-pling, however, the G σσ ′ = 0 for σ = σ ′ . Therefore, forthe simplified Kane-Mele-type model, the Green’s func-tions reduce to 2 × × G ↑↑ ( iω = 0 , k i ) = α k i σ x , k i ∈ TRIM , (A7)where some coefficients multiply the σ x Pauli matrix, andin equation (7), ˜ η ( k i ) = ± U , ˜ η ( k i ) = ± k i = M , g ij = [ G ( iω = 0 , M )] ij , vs the number5of measurements ( m ) in Figs. 12. λ SO = 0 . t and U = 4 t are used. The test system size is 2 × . To recover equa-tion (A7), one should expect that Re[ g ] ≃ Re[ g ], andIm[ g ] = Im[ g ] = || g || ≃
0. It has been demon-strated that the values of Re[ g ] = α k i can be used toidentify the topological property . g ij
1/ m
Re[g ] Re[g ] Im[g ] Im[g ] ||g || ||g || (a) t =0.32t (b) t =0.37t
1/ m
FIG. 12. (Color online) The matrix elements of the zero-frequency Green’s functions G (0 , M ) vs the number of sam-plings m at (a) t N = 0 . t and (b) t N = 0 . t . λ SO = 0 . t and U = 4 t . Re[ g ij ] and Im[ g ij ] denote the real part andimaginary part of [ G (0 , M )] ij , respectively; || g ii || denotesthe diagonal component of G (0 , M ) in magnitudes. Fig. 12 (a) shows t N = 0 . t in the Z topologicalinsulator phase and (b) for t N = 0 . t in the trivial in-sulator. At small m , the real parts of g and g are notequal; furthermore, g and g have imaginary parts,and both of g and g are finite. However, one can seethat, upon sampling sufficient times, Re[ g ] ≃ Re[ g ],and meanwhile Im[ g ], Im[ g ], || g || go to zero.Thus, in the m → ∞ limit, equation (A7) is recov-ered, and then the value of h ˜ η ( k i ) i over QMC simulationsmonotonically approaches to ± k i and interacting case, equation (A7) doesnot hold. However, for the non-interacting case we foundthat the value of resulting spin Chern number C σ isnot sensitive to the number of samplings provided theyare large enough in number. Throughout our paper, wechoose the number of measurements large enough (mostlyover several thousands) to determine the 2 × Appendix B: Projection operator expression of theSpin-Chern number
In this section we provide a description of the projec-tion operator expression used to evaluate the spin Chernnumber for finite lattices and its practical numerical im-plementation. The expression for the Chern numbers us- ing the projection operators onto the occupied bands is, C σ = 12 πi Z B.Z.
Tr ( P σ dP σ ∧ dP σ )= 12 πi Z B.Z. Tr n P σ ( k ) h ∂ k x P σ ( k ) ∂ k y P σ ( k ) − ∂ k y P σ ( k ) ∂ k x P σ ( k ) io dk x dk y , (B1)where P σ ( k ) is the spectral projector operator con-structed using the Bloch eigenvectors (eigenspace) at k with energies below the Fermi energy ǫ F , i.e., E n ( k ) < ǫ F and for spin sector- σ . A merit of this formulation ofthe Chern number is the manifest independence of theU(1) phases of the Bloch states. The Bloch eigenstatesthemselves are obtained from diagonalizing the interact-ing zero-frequency Green’s functions G σ ( k , | µ i i = µ i | µ i i , (B2)and then P σ ( k ) = X µ i > | µ i ih µ i | , (B3)where choosing µ i > E n < ǫ F , i.e. R-zero of the G σ ( k , U (1) gaugeinvariant. The integral is over the Brillouin zone (BZ),and, in practical numerical application, the region of in-tegration over the BZ does not need to be a Wigner-Seitzunit cell in reciprocal lattice space, as long as the entirereciprocal lattice unit cell is covered.In a finite-size system, the set of k -vectors is dis-cretized, so we will need to replace the integral withthe summation over finite momentum points. For con-venience we can map the momentum points as a N = L x × L y square grid of spacing h and label each k withdiscrete coordinate indices { m, n } . Then we can approxi-mate the partial derivatives ∂ k x P σ ( k ) and ∂ k y P σ ( k ) usingthe symmetric finite difference as ∂ k x P σ ( k ) ≈ P σ,i +1 ,j − P σ,i − ,j h ,∂ k y P σ ( k ) ≈ P σ,i,j +1 − P σ,i,j − h . Thus, in equation (B1) we simplify P σ ( k ) h ∂ k x P σ ( k ) , ∂ k y P σ ( k ) i ≈ P σ,i,j h (cid:16) [ P σ,i +1 ,j , P σ,i,j +1 ] + [ P σ,i,j +1 , P σ,i − ,j ]+[ P σ,i − ,j , P σ,i,j − ] + [ P σ,i,j − , P σ,i +1 ,j ] (cid:17) . Note that due to periodic boundary conditions in BZ, P σ,L x +1 ,j ≡ P σ, ,j and P σ,i,L y +1 ≡ P σ,i, . Then the6Chern number is approximated as C σ = 12 πi Z B.Z.
Tr ( P σ dP σ ∧ dP σ ) ≈ πi N X i,j =1 P σ,i,j (cid:16) [ P σ,i +1 ,j , P σ,i,j +1 ] + [ P σ,i,j +1 , P σ,i − ,j ]+[ P σ,i − ,j , P σ,i,j − ] + [ P σ,i,j − , P σ,i +1 ,j ] (cid:17) ..