Entanglement Spectra of Interacting Fermions in Quantum Monte Carlo Simulations
EEntanglement Spectra of Interacting Fermions in Quantum Monte Carlo Simulations
Fakher F. Assaad, Thomas C. Lang, and Francesco Parisen Toldin Institut f¨ur Theoretische Physik und Astrophysik,Universit¨at W¨urzburg, Am Hubland, D-97074 W¨urzburg, Germany Department of Physics, Boston University, Boston, MA 02215, USA (Dated: September 20, 2018)In a recent article T. Grover [Phys. Rev. Lett. , 130402 (2013)] introduced a simple methodto compute Renyi entanglement entropies in the realm of the auxiliary field quantum Monte Carloalgorithm. Here, we develop this approach further and provide a stabilization scheme to computehigher order Renyi entropies and an extension to access the entanglement spectrum. The methodis tested on systems of correlated topological insulators.
PACS numbers: 02.70.Ss,03.67.-a,71.10.-w,73.43.-f
I. INTRODUCTION
The concept of entanglement spectra and entropies isemerging as a powerful tool to access universal quantitiessuch as the central charge of an associated conformal fieldtheory and characterize ground states and their degenera-cies beyond the Landau paradigm. Typically, the groundstate of a local Hamiltonian is entangled over short dis-tances. This short range entanglement leads to an arealaw (or more precisely, a perimeter law) of the Renyi andvon Neumann entanglement entropies [1]. Correctionsto the area law have implications for the nature of theground state. For instance, topologically ordered statesin two dimensions with anyonic elementary excitations [2]are characterized by a sub-leading correction to the arealaw [3, 4]. This signature of topological order has becomea numerical tool that has been implemented within theframeworks of variational quantum Monte Carlo (QMC)[5], stochastic series expansions [6, 7], and density ma-trix renormalization group (DMRG) [8–10]. Other cor-rections to the area law include sub-leading logarithmi-cally diverging corrections which signalize the presenceof Goldstone modes of spontaneously broken continuoussymmetries [11]. Entanglement entropies have also beenintroduced in the realm of the valence bond quantumMonte Carlo simulations for spin systems [12–14]. TheRenyi entropy itself does not necessarily provide an effi-cient tool to characterize states of matter. In fact, andas pointed out recently in [15], there is no practical wayof distinguishing a Chern insulator from a trivial bandinsulator by computing the entanglement entropy.More information on the nature of the ground state iscontained in the entanglement spectrum. The entangle-ment spectrum picks up the edge physics of topologicallyordered fractional quantum Hall states [16–19] as well asthat of topological insulators [20, 21]. In magnetically or-dered states with spontaneously broken SU(2) spin sym-metry, Anderson’s tower of states can be detected in theentanglement spectrum [22]. In a recent work Chandran et al. caution against the unconditional assumption thatthe low energy sector of the entanglement spectrum con-tains universal properties of the state and of quantum phase transition [23].The aim of this article is to extend a recent proposalto compute Renyi entropies in the realm of the auxiliaryfield QMC algorithm [24], in order to access the entangle-ment spectrum. In the next section, we will set the stageand briefly review the zero temperature formulation ofthe auxiliary field QMC method as well as the calcu-lation of Renyi entropies. Emphasis will be placed onstability issues we encountered when computing higherorder Renyi entropies on large subsystems. Our approachto compute the entanglement spectrum is based on spec-troscopy along the replica time axis. This notion will beintroduced in Section III. We test the method by comput-ing the single-particle entanglement spectral function forthe Kane-Mele Hubbard model supplemented by a single-particle bond dimerization term [25]. This dimerizationtriggers a topological phase transition to a topologicallytrivial state, which manifests plainly in the entanglementspectrum. In section IV we conclude and discuss theshortcomings and strong points of our approach.
II. RENYI ENTROPIES FROM THE ZEROTEMPERATURE AUXILIARY FIELD QMCMETHODA. The Projective auxiliary field QMC algorithm
Here we provide a short description of the zero tem-perature, or projective, auxiliary field QMC (PQMC) al-gorithm for the Kane-Mele Hubbard model:ˆ H = ˆ c † T ˆ c + U (cid:88) iii (ˆ n iii − ≡ ˆ H T + ˆ H U . (1)ˆ c † is a vector of fermionic creation operators with entriesˆ c † iiiσ . The operator ˆ c † iiiσ creates an electron at lattice site iii with a z -component of spin σ , and ˆ n iii = (cid:80) σ ˆ c † iiiσ ˆ c iiiσ is thecharge density at site iii . The kinetic energy is defined bythe N × N Hermitian matrix T , and N corresponds to thenumber of single-particle states, which we will henceforthlabel with super-index x ≡ ( iii, σ ). In this article, ˆ H T will a r X i v : . [ c ond - m a t . s t r- e l ] M a r FIG. 1. Honeycomb lattice and the regions A and B. In thefigure region A has a width of W A = 4. be taken to be the dimerized Kane-Mele model. Usingthe spinor notation ˆ c † i = (cid:0) ˆ c † iii ↑ , ˆ c † iii ↓ (cid:1) ,ˆ H T = (cid:88) iii,jjj ˆ c † iii [ t iiijjj + i λ iiijjj · σ ] ˆ c jjj . (2)The hopping matrix takes non-vanishing values be-tween nearest neighbors of the honeycomb lattice, iii − jjj = ± δδδ , ± δδδ , ± δδδ (see Fig. 1), and we have imple-mented the following dimerization: t iiijjj = − t if iii − jjj = ± δδδ , ± δδδ − t (cid:48) if iii − jjj = ± δδδ . (3)The intrinsic spin-orbit term is given by λλλ iiijjj = λ (cid:40) ( iii − rrr ) × ( rrr − jjj ) | ( iii − rrr ) × ( rrr − jjj ) | if iii, jjj are n.n.n.0 otherwise , (4)where rrr is the intermediate site involved in the next near-est neighbor (n.n.n.) hopping process from site iii to jjj .The PQMC algorithm is based on filtering out theground state from a Slater determinant trial wave func-tion: | Ψ T (cid:105) = N p (cid:89) n =1 (cid:32) N (cid:88) x =1 ˆ c † x P xn (cid:33) | (cid:105) . (5)The trial wave function is thus defined by the rectangu-lar N × N p matrix P with N p the number of particles.Given this trial wave function, and assuming that it isnon-orthogonal to the ground state, observables can beobtained from (cid:104) ˆ O (cid:105) = lim Θ →∞ (cid:104) Ψ T | e − Θ ˆ H ˆ O e − Θ ˆ H | Ψ T (cid:105)(cid:104) Ψ T | e −
2Θ ˆ H | Ψ T (cid:105) , (6)for large but finite values of the projection parameter Θ.To compute the imaginary time propagation one first dis-cretizes the imaginary time, L Θ ∆ τ = 2Θ, and then car-ries out a Trotter decomposition to isolate the Hubbard interaction term,e −
2Θ ˆ H = L Θ (cid:89) τ =1 e − ∆ τ ˆ H T / e − ∆ τ ˆ H U e − ∆ τ ˆ H T / + O (∆ τ ) . (7)Hereafter we neglect the systematic and controllableTrotter error [26]. The key point of the algorithm is touse a Hubbard-Stratonovitch (HS) transformation to re-formulate the many-body imaginary time propagator asa sum of one-body problems by introducing an auxiliaryfield. We have adopted the discrete decomposition [27],e − ∆ τ U (ˆ n iii − = γ (cid:88) s = ± e sα ( n iii ↑ − n iii ↓ ) , (8)with cosh( α ) = e ∆ τU/ and γ = e − ∆ τU/ . With thistransformation, the imaginary time propagation readse −
2Θ ˆ H = γ NL Θ (cid:88) sss ··· sss L Θ L Θ (cid:89) τ =1 e ˆ c † A ( sss τ )ˆ c e − ∆ τ ˆ c † T ˆ c , (9)with A ( sss τ ) xy = δ xy s iiiτ σ α . Recall that x = ( iiiσ ) and that σ takes the value 1 ( −
1) for the up (down) z -componentof the spin. For a given configuration of Ising variables, sss ≡ { sss · · · sss L Θ } , we now have to solve a free fermionproblem in an external space and time dependent field.Since under a single body propagator a Slater determi-nant remains a Slater determinant,e ˆ c † h ˆ c N p (cid:89) n =1 (cid:0) ˆ c † P (cid:1) n | (cid:105) = N p (cid:89) n =1 (cid:0) ˆ c † e h P (cid:1) n | (cid:105) , (10)and the overlap of two Slater determinants defined by P and P (cid:48) is a determinant, (cid:104) Ψ (cid:48) T | Ψ T (cid:105) = det (cid:0) P (cid:48)† P (cid:1) , (11)we can integrate out the fermions to obtain: (cid:104) ˆ O (cid:105) = (cid:88) sss P sss (cid:104) ˆ O (cid:105) sss . (12)Here, P sss = det ( U The n th Renyi entropy is defined as: S n = − n − H A [ˆ ρ n A ] . (18)To evaluate the above quantity, one introduces n replicas,(or n independent QMC simulations) such thate − ( n − S n = (cid:88) sss , ··· ,sss n P sss n · · · P sss Tr H A (cid:2) ˆ ρ A ( sss n ) · · · ˆ ρ A ( sss ) (cid:3) = (cid:88) sss , ··· ,sss n P sss n · · · P sss n (cid:89) m =1 det [ − G A ( sss m )] × det (cid:34) + n (cid:89) m =1 G A ( sss m ) − G A ( sss m ) (cid:35) . (19)Here we have made use of the identity:Tr (cid:104) e ˆ ccc † T ˆ ccc · · · e ˆ ccc † T n ˆ ccc (cid:105) = det (cid:2) + e T · · · e T n (cid:3) and thesuperscripts denote the replicas. At n = 2 one only S n Λ U / t = 2, L = 6, W A = 4 n = 8 n = 7 n = 6 n = 5 n = 4 n = 3 n = 2 S n nU / t = 2, L = 6, W A = 4, Λ = 5 10 -6 FIG. 2. Renyi entanglement entropies as a function of theregularization parameter Λ (a) and order n (b) for a 6 × U/t = 2, λ/t = 0 . t (cid:48) /t = 1and the region A defined in Fig. 1. Errorbars are smallerthan the symbol size. For these simulations, we have used aprojection parameter Θ t = 40. This guarantees convergencewithin our error-bars. For each replica, we have carried out500 × sweeps. Each sweep consists of sequentially updatingeach Ising spin in the space-imaginary time lattice. needs the knowledge of G A ( sss ) to compute the Renyientropy and one recovers Eq. (6) of Ref. [24]. For n > G A ( sss ) which will also be required tocompute the entanglement spectrum. However, G A ( sss )is in general a singular matrix. To show this, let us firstconsider the Green function in Eq. (14). As pointed outin [32] the Green function in the PQMC algorithm is aprojector operator, G ( sss ) = G ( sss ) , (20)such that the eigenvalues of G ( sss ) are given by 0 or1. Since the Tr [ G ( sss )] = N p , G ( sss ) contains N p non-vanishing and N − N p vanishing eigenvalues. Thereby G − ( sss ), and by the same token [ − G ( sss )] − , does notexist. In the absence of interactions, the above is merelystating that in the zero temperature limit the occupationnumber of single-particle eigenstates is given by 1 or 0.That G − ( sss ) does not exist does not necessarily meanthat G A ( sss ) is singular. In fact for the half-filled case,if the domain A (cf. Fig. 1) contains √ N s out of N s sites, G A ( sss ) turns out to be generically regular. On theother hand, if the domain A contains more sites than thenumber of particles N p , G A ( sss ) is indeed singular.At finite temperatures the Green function matrix isnever singular due to thermal broadening. This suggeststhe following regularization scheme. Let G ⊥ ( sss ) be a pro-jector on the vector space spanned by the elements of thekernel of G ( sss ). One can then consider the quantity˜ G ( sss ) = ( − Λ) G ( sss ) + Λ G ⊥ ( sss ) . (21)Essentially, Λ introduces a finite temperature since the oc-cupation number of occupied (unoccupied) single-particlestates is reduced (enhanced) from 1 (0) to 1 − Λ (Λ) [33].By definition ˜ G ( sss ) is regular. From ˜ G ( sss ) we can con-struct ˜ G A ( sss ) which equally proves to be regular. If Λis small, we do not expect this regularization scheme toalter the nature of the entanglement spectrum, nor theentanglement entropies. Furthermore, within the sameQMC run, it is possible to extrapolate to Λ → n = 8 on a 6 × U/t = 2, λ/t = 0 . t = t (cid:48) with region A being a ribbon ofwidth W A = 4 (see Fig. 1). For this specific problem,values of Λ = 5 × − suffice to guarantee stability andto render systematic errors negligible within error bars.From now on we will always use the presented regular-ization scheme and generically use the notation G A ( sss )instead of ˜ G A ( sss ).Instead of implementing the above regularizationscheme, one could use a finite temperature algorithm[28, 34]. There are however many cases in which an ad-equate choice of the trial wave function leads to consid-erable computational gain insofar as ground state prop-erties are concerned [35]. In this case, the above reg-ularization scheme will certainly be more useful than afinite-temperature simulation.Figure 2(b) plots the Renyi entropies as a function of n . It is interesting to note that extrapolation to n = 1,corresponding to the von Neumann entropy, is difficult.Detailed knowledge of the higher order Renyi entropieshas been used in Ref. [36] to access the entanglementspectra. In the next section we will present an alternativeapproach. III. ACCESSING THE ENTANGLEMENTSPECTRUM Generically, in quantum Monte Carlo simulations, wecan access the spectrum of the Hamiltonian by comput-ing imaginary time displaced correlation functions in var-ious channels. Similarly, for the entanglement Hamilto-nian, given by ˆ ρ A ≡ e − ˆ H E Tr H A (cid:104) e − ˆ H E (cid:105) , (22) -10-5 0 5 10 0 1 2 3 4 5 6 E n t a ng l e m e n t s p ec t r u m kaU / t = 0, t ’/ t = 1.5 L = W A = L = W A = -10-5 0 5 10 0 1 2 3 4 5 6 E n t a ng l e m e n t s p ec t r u m kaU / t = 0, L = 128, W A = 52 t ’/ t = 1.9 t ’/ t = 2.1 -10-5 0 5 10 0 1 2 3 4 5 6 E n t a ng l e m e n t s p ec t r u m kaU / t = 0, L = 72, W A = 48, t ’/ t = 2.1 FIG. 3. Entanglement spectrum of the dimerized Kane-Melemodel at λ/t = 0 . U/t = 0. In (a) and (b) the dimer-ization runs along the δδδ direction (See Fig. 1). In (c) thedimerization runs along the δδδ direction. The spectrum hasbeen obtained by inverting eq. (26). we can define the following replica time displaced corre-lation functions S EO ( τ E ) ≡ (cid:104) ˆ O † ( τ E ) ˆ O (cid:105) A ≡ Tr H A (cid:104) e − ( n − τ E ) ˆ H E ˆ O † e − τ E ˆ H E ˆ O (cid:105) Tr H A (cid:104) e − n ˆ H E (cid:105) , (23)for an observable O ∈ H A . Here τ E and n are integerswith τ E < n . In analogy to imaginary time, one can inserta complete set of eigenstates of ˆ H E , ˆ H E | Ψ i (cid:105) = E i | Ψ i (cid:105) toobtain S EO ( τ E ) = 1 π (cid:90) d ω e − τ E ω ± e − τ E ω A EO ( ω ) , (24a) with A EO ( ω ) = π (cid:80) i e − nE i (cid:88) i,j e − nE i |(cid:104) Ψ i | ˆ O | Ψ j (cid:105)| × δ ( ω + E i − E j ) (cid:0) ± e − βω (cid:1) . (24b)Here, the plus (minus) sign refers to single-particle(particle-hole) excitations. Hence, in principle entangle-ment spectral functions can be computed by using stan-dard maximum entropy methods [37–39]. Note however,that in contrast to imaginary time, replica time is quan-tized to integer values in our approach. Hence, with theuse of the maximum entropy method we will only be ableto resolve the low lying spectral features.The Monte Carlo evaluation of Eq. (23) is very similarto computing imaginary time displaced correlation func-tions in the realm of the finite temperature formulation ofthe auxiliary field determinant QMC algorithm [28]. Tohighlight the similarity we will introduce the notation:ˆ B A ( sss ) = e − ˆ aaa † ln [ G − ( sss ) − ] ˆ aaa , and B A ( sss ) = e − ln [ G − ( sss ) − ],such that (cid:104) ˆ O † ( τ E ) ˆ O (cid:105) A = (cid:80) sss , ··· ,sss n P sss n · · · P sss Z n (cid:0) sss n · · · sss (cid:1) (cid:104) ˆ O † ( τ E ) ˆ O (cid:105) A (cid:0) sss n · · · sss (cid:1)(cid:80) sss , ··· ,sss n P sss n · · · P sss Z n ( sss n · · · sss ) , (25a)with Z n (cid:0) sss n · · · sss (cid:1) = (cid:34) n (cid:89) m =1 det [ − G A ( sss m )] (cid:35) Tr H A (cid:104) ˆ B A ( sss n ) · · · ˆ B A ( sss ) (cid:105) , (25b)and (cid:104) ˆ O † ( τ E ) ˆ O (cid:105) A (cid:0) sss n · · · sss (cid:1) = Tr H A (cid:104) ˆ B A ( sss n ) · · · ˆ B A ( sss τ E +1 ) ˆ O † ˆ B A ( sss τ E ) · · · ˆ B A ( sss ) ˆ O (cid:105) Tr H A (cid:104) ˆ B A ( sss n ) · · · ˆ B A ( sss ) (cid:105) . (25c)The denominator in Eq.(25a) is merely e − ( n − S n . Tocompute the numerator one has to evaluate the replicatime displaced correlation function of Eq. (25c). Tech-nically speaking this is identical to computing imagi-nary time displaced correlation functions in the realmof the finite temperature auxiliary field QMC algorithm.For a given value of the HS fields Wick’s theorem holdssuch that it suffices to compute the single-particle replicatime displaced Green function. We refer the reader toRef. [28] for a detailed description. To implement theabove, we need to run n independent PQMC simulations(or replicas) [40]. At any one time, the configurationswill be distributed according to the probability distribu-tion P sss n · · · P sss . With this configuration of HS fields, wecompute the B A ( sss m ) and G A ( sss m ) matrices from whichcan evaluate the n-th Renyi entropy as well as the replicatime displaced correlation functions [41].We have tested the approach for the Kane-Mele Hub- bard model defined in Eq. (1). In the non-interactingcase, the entanglement Hamiltonian satisfies the relation − G A = 12 tanh ( H E / . (26)Using the spectral flattening trick [2, 20, 21] we can adi-abatically connect ( − G ) to the original Hamiltonian.Thereby both Hamiltonians share the same topologicalcharacter and restricting ( − G ) to a region A withedges will reveal the edge physics of the original Hamil-tonian. The eigenvalues of the matrix ( − G ) take values ± / 2. Thereby the eigenvalues of ( − G A ) are restrictedto the range ( − / 2; 1 / 2) and the spectrum of H E is notbounded. In the plots of Fig. 3 we consider an energywindow [ − t, t ]. Since the considered cut A preservestranslation symmetry along the a direction, ka = k · a remains a good quantum number. Alternative choices ofA which do not break any lattice symmetries have beenproposed [42].In the absence of interactions, the critical dimerizationwhere the topological phase transition between topologi-cal and trivial band insulators takes place reads t (cid:48) c /t = 2.As shown in Fig. 3(a), at t (cid:48) /t = 1 . < t (cid:48) c /t one observesthe typical edge signature of the quantum spin Hall(QSH) state: a single crossing at the time-reversal sym-metric point ka = π , which due to particle-hole symme-try, is pinned at zero energy. In the vicinity of the criticalpoint, the velocity of the edge state grows, and beyond t (cid:48) c /t a gap opens [cf. see Fig. 3(b)].It is interesting to note that the entanglement spec-trum depends on the direction in which the dimerizationis imposed. In particular, for a dimerization along the δδδ bond, the entanglement spectrum is given in Fig. 3(c).For t (cid:48) > t (cid:48) c this is still the spectrum of a trivial insula-tor: the surface state is disconnected from the bulk andan even number of crossings at ka = π and 0 are visible[43]. Hence, care has to be taken when analyzing the en-tanglement spectrum before drawing conclusions on thenature of the state.In the presence of interactions, we can compute thereplica time displaced single-particle Green function G E αβ ( k, τ E ) = (cid:88) σ (cid:104) ˆ a † kασ ( τ E ) ˆ a kβσ (cid:105) A , (27)where the indices α, β run over the sites across the widthof the ribbon defined by region A, and k is the conservedcrystal moment associated with translation invariancealong the a direction. As mentioned above, G E αβ ( k, τ E ),together with the replica time displaced correlation func-tions of other observables, provides the input for max-imum entropy methods which in principle, allow us todetermine the spectrum of ˆ H E . On the left side of Fig. 4we show G E αβ ( k, τ E ) at U/t = 2 and for various valuesof the dimerization. We consider the orbital α = 1 situ-ated on the edge of the ribbon as well as two momenta, ka = π and 0. At t (cid:48) /t = 1 . G E , ( π, τ E ) does not decayas a function of replica time τ E thus signaling a gaplesssingle-particle excitation at ka = π . In contrast, at k = 0a clear gap is visible. We hence conclude that we are ina topological phase. At t (cid:48) /t = 1 . 95 analysis of the samequantity shows large finite size effects, which one can in-terpret in terms of the vanishing of the single-particle gapat ka = π , and a robust gap at k = 0. Hence the data at t (cid:48) /t = 1 . 95 allows interpretation in terms of a topologicalphase. On the other hand, at t (cid:48) /t = 2 analysis of the dataas a function of system size shows that the gap remainsrobust at both considered k -points. Based on the entan-glement spectra of the single-particle Green function onecan conclude that the topological phase transition occursin the interval 1 . < t (cid:48) /t < 2. This is in accordance withindependent results of Ref. [25].Using the maximum entropy method and Eq. (24a)we can analytically continue the replica time single- particle spectral function to obtain the entanglement single-particle spectral function, A E1 , ( π, ω ). This quan-tity is plotted on the right side of Fig. 4. The analyt-ical continuation confirms the conclusions drawn fromthe replica time data. Deep in the topological phase at t (cid:48) /t = 1 . t (cid:48) /t = 2 [Fig. 4(c)]the edge mode acquires a gap. Close to the transition at t (cid:48) /t = 1 . 95 [Fig. 4(b)] we interpret the gap at k = π interms of a finite size effect. IV. CONCLUSIONS Based on the work of T. Grover [24] we have demon-strated how to compute in a numerically stable mannerRenyi entanglement entropies in the realm of the pro-jective auxiliary field quantum Monte Carlo algorithm.Furthermore we have introduced replica time displacedcorrelation functions from which the low-lying propertiesof the entanglement spectrum can be reconstructed. Thestrong point of such an approach is two-fold: One canstudy the entanglement spectrum in different channels,single-particle, particle-hole, or particle-particle, and onecan classify the spectrum in terms of symmetries allowedby the choice of the subsystem A. We have successfullytested this approach for the topological quantum phasetransition from the QSH insulator to a trivial band in-sulator in the weak coupling regime around U/t = 2.The weak point of our present implementation is thatthe replica time displaced correlation functions becomestatistically noisy at stronger couplings. This originatesfrom the fact that the stochastic evaluation of the Renyientropies becomes progressively noisy as a function ofcoupling strength, size of the subregion A and n. Sucha behavior can be understood by assuming an area law, S n ∝ α n ( U ) ∂A , where ∂A is the length of the bound-ary of subsystem A. The quantity we evaluate stochasti-cally is (cid:104) Z n (cid:105) ∝ e − ( n − α n ( U ) ∂A . Numerically we see that α n ( U ) is a growing function of U and weakly dependenton n . Hence (cid:104) Z n (cid:105) decreases exponentially with U , n and ∂A but the fluctuations do not decrease as quickly suchthat the signal to noise ratio increases strongly as a func-tion of the above mentioned parameters. As a conse-quence, and at this point, it is computationally very de-manding to pin down the Hubbard U driven transitionfrom the QSH state to the magnetically ordered phase[44–46]. Different sampling schemes, or improved esti-mators will be required to solve these issues.We would like to thank L. Bonnes, P. Br¨ocker,T. Grover, I. Herbut, A. L¨auchli, R. Melko and C. Mudryfor discussions. We thank the LRZ-M¨unich and theJ¨ulich Supercomputing center for generous allocationof CPU time. Financial support from the DFG grantAS120/9-1 as well as the DFG funded FOR1807 is ac-knowledged. G E , ( k , τ E ) τ E U / t = 2, t ’/ t = 1.5 L = ka = π , W A = L = ka = , W A = L = ka = π , W A = L = ka = , W A = ππ /2 -10 -5 0 5 10 k a ω / tA E1,1 ( k, ω ), t’ / t = 1.5 G E , ( k , τ E ) τ E U / t = 2, t ’/ t = 1.95 L = ka = π , W A = L = ka = , W A = L = ka = π , W A = L = ka = , W A = ππ /2 -10 -5 0 5 10 k a ω / tA E1,1 ( k, ω ), t’ / t = 1.95 G E , ( k , τ E ) τ E U / t = 2, t ’/ t = 2 L = ka = π , W A = L = ka = , W A = L = ka = π , W A = L = ka = , W A = ππ /2 -10 -5 0 5 10 k a ω / tA E1,1 ( k, ω ), t’ / t = 2 FIG. 4. Replica time displaced single-particle Green function (left) and the corresponding single-particle spectral function(right) for the Kane-Mele Hubbard model at U/t = 2, λ/t = 0 . t (cid:48) /t across the transitionfrom the QSH insulator to the topologically trivial band insulator. In computing the single-particle spectral function we haveconsidered the 12 × 12 lattice with W A = 8. For the largest system sizes, L = 12, we have used 40 × sweeps per replica andset the imaginary time propagation parameter to Θ t = 40. The spectral functions are obtained with the stochastic maximumentropy method [38, 39]. [1] J. Eisert, M. Cramer, and M. B. Plenio, Rev. Mod. Phys. , 277 (2010).[2] A. Kitaev, Annals of Physics , 2 (2006).[3] M. Levin and X.-G. Wen, Phys. Rev. Lett. , 110405(2006).[4] A. Kitaev and J. Preskill, Phys. Rev. Lett. , 110404(2006).[5] Y. Zhang, T. Grover, and A. Vishwanath, Phys. Rev. B , 075128 (2011).[6] M. B. Hastings, I. Gonz´alez, A. B. Kallin, and R. G.Melko, Phys. Rev. Lett. , 157201 (2010).[7] S. V. Isakov, M. B. Hastings, and R. G. Melko, Nat Phys , 772 (2011).[8] S. Yan, D. A. Huse, and S. R. White, Science , 1173(2011).[9] H.-C. Jiang, Z. Wang, and L. Balents, Nat Phys , 902(2012).[10] S. Depenbrock, I. P. McCulloch, and U. Schollw¨ock,Phys. Rev. Lett. , 067201 (2012).[11] M. Metlitski and T. Grover, arXiv:1112.5166 .[12] A. W. Sandvik, Phys. Rev. Lett. , 207203 ((2005)).[13] K. Beach and A. W. Sandvik, Nuclear Physics B ,142 (2006).[14] F. Alet, S. Capponi, N. Laflorencie, and M. Mambrini,Phys. Rev. Lett. , 117204 (2007).[15] J. Budich, J. Eisert, and E. J. Bergholtz, arXiv:1311.3309.[16] H. Li and F. D. M. Haldane, Phys. Rev. Lett. , 010504(2008).[17] A. M. L¨auchli, E. J. Bergholtz, J. Suorsa, and M. Haque,Phys. Rev. Lett. , 156404 (2010).[18] R. Thomale, A. Sterdyniak, N. Regnault, and B. A.Bernevig, Phys. Rev. Lett. , 180502 (2010).[19] X.-L. Qi, H. Katsura, and A. W. W. Ludwig, Phys. Rev.Lett. , 196402 (2012).[20] L. Fidkowski, Phys. Rev. Lett. , 130502 (2010).[21] A. M. Turner, Y. Zhang, and A. Vishwanath, Phys. Rev.B , 241102 (2010).[22] F. Kolley, S. Depenbrock, I. McCulloch, U. Schollw¨ock,and V. Alba, arXiv:1307.6592 .[23] A. Chandran, V. Khemani, and S. L. Sondhi,arXiv:1311.2946 .[24] T. Grover, Phys. Rev. Lett. , 130402 (2013).[25] T. C. Lang, A. M. Essin, V. Gurarie, and S. Wessel, Phys.Rev. B , 205101 (2013).[26] D. Rost, E. V. Gorelik, F. Assaad, and N. Bl¨umer, Phys.Rev. B , 155109 (2012).[27] J. E. Hirsch, Phys. Rev. B , 4059 (1983).[28] F. F. Assaad and H. G. Evertz, in Computational Many Particle Physics , Vol. 739 of Lecture Notes inPhysics , edited by H. Fehske, R. Schneider, and A. Weiße(Springer Verlag, Berlin, 2008), pp. 277–356.[29] I. Peschel, Journal of Physics A: Mathematical and Gen-eral , L205 (2003).[30] Note that ˆ ρ ( sss ) is not a density matrix but a general Gaus-sian operator. In particular, Tr (cid:2) ˆ ρ ( sss )ˆ c † x ˆ c x (cid:3) ≡ G x,x ( sss ) canbecome negative or even complex depending upon thechoice of the HS transformation.[31] With this definition of ˆ ρ A one will readily verify that forany observable ˆ A ∈ H A , Tr A (cid:104) ˆ ρ H A ˆ A (cid:105) = Tr (cid:104) ˆ ρ ˆ A (cid:105) .[32] M. Feldbacher and F. F. Assaad, Phys. Rev. B , 073105(2001).[33] Note that for the here considered half-filled case where N p = N/ 2, this regularization scheme does not alter theparticle number.[34] S. R. White, D. J. Scalapino, R. L. Sugar, E. Y. Loh,J. E. Gubernatis, and R. T. Scalettar, Phys. Rev. B ,506 (1989).[35] S. Capponi and F. F. Assaad, Phys. Rev. B , 155114(2001).[36] C.-M. Chung, L. Bonnes, P. Chen, and A. M. L¨auchli,arXiv:1305.6536 .[37] M. Jarrell and J. Gubernatis, Physics Reports , 133(1996).[38] A. Sandvik, Phys. Rev. B , 10287 (1998).[39] K. S. D. Beach, arXiv:0403055 (2004).[40] Alternatively, replica time displaced Green functions canbe obtained in a single run from measurements whichare sufficiently separated in imaginary time, such as toguarantee statistical independence.[41] For the Kane-Mele Hubbard model at half-band fill-ing, one can show that Z n (cid:0) sss n · · · sss (cid:1) is positive. Henceinstead of using a re-weighting scheme, one can inprinciple devise an algorithm which directly samples P sss n · · · P sss Z n (cid:0) sss n · · · sss (cid:1) .[42] M. Legner and T. Neupert, Phys. Rev. B , 115114(2013).[43] M. Z. Hasan and C. L. Kane, Rev. Mod. Phys. , 3045(2010).[44] M. Hohenadler, T. C. Lang, and F. F. Assaad, Phys. Rev.Lett. , 100403 (2011).[45] M. Hohenadler, Z. Y. Meng, T. C. Lang, S. Wessel, A.Muramatsu, and F. F. Assaad, Phys. Rev. B , 115132(2012).[46] F. F. Assaad, M. Bercx, and M. Hohenadler, Phys. Rev.X3