Topological invariants in interacting Quantum Spin Hall: a Cluster Perturbation Theory approach

Using Cluster Perturbation Theory we calculate Green's functions, quasi-particle energies and topological invariants for interacting electrons on a 2-D honeycomb lattice, with intrinsic spin-orbit coupling and on-site e-e interaction. This allows to define the parameter range (Hubbard U vs spin-orbit coupling) where the 2D system behaves as a trivial insulator or Quantum Spin Hall insulator. This behavior is confirmed by the existence of gapless quasi-particle states in honeycomb ribbons. We have discussed the importance of the cluster symmetry and the effects of the lack of full translation symmetry typical of CPT and of most Quantum Cluster approaches. Comments on the limits of applicability of the method are also provided.

Topological invariants are by now widely recognized as a powerful tool to characterize different phases of matter; in particular they turn out to be useful in the classification of topological insulators. 1,2 In the topological insulator phase, solids are characterized by gapped bulk bands but present gapless edge states that allow charge or spin conductivity on the boundaries. The presence of such gapless edge states is linked to the emergence of non-vanishing topological invariants via a bulk-boundary correspondence. 3 This topological feature ensures the robustness of the edge states against disorder. 4,5 A two-dimensional honeycomb lattice with spin-orbit coupling has been identified as a remarkable and paradigmatic example of topological insulator. This system is the prototype of the so-called Quantum Spin Hall (QSH) system presenting a quantized spin-Hall conductance at the boundaries. The topological nature of QSH insulators is identified by a time-reversal (T ) -topological invariant Z 2 6-8 . In the same way as the Thouless-Kohmoto-Nightingale-den Nijs 9 (TKNN) topological invariant was defined for the integer quantum Hall effect, the above Z 2 invariant was defined for the topological insulator in terms of band eigenvectors and, as such, only applies to noninteracting systems. On the other hand, in the presence of electron-electron interaction, there are generalizations of the TKNN invariant based on twisted boundary conditions 10 and on many-body Green's functions, [11][12][13] and unlike the former, the latter construction can be straightforwardly extended to the Z 2 invariant to classify interacting QSH systems.
The field of interacting topological insulators is attracting growing interest (see Refs. 14 and 15 for recent reviews on two-dimensional systems) and the definition of theoretical and computational tools to evaluate topological invariants in the presence of e-e interaction is extremely timely. The approach that seems most promis-ing is the one developed in Ref. 12 where it has been demonstrated that topological invariants are determined by the behaviour of the one-particle propagator at zero frequency only; more precisely it has been shown that the eigenvectors of the single particle Hamiltonian that yields topological invariants for non-interacting systems, should be replaced in the interacting case by the eigenvectors of the operator G −1 (k, ω) at ω = 0 and momentum k. The extension of topological invariants to interacting systems is in this sense straightforward, the only demanding task remaining the determination of the dressed Green's function. This concept has been recently applied to identify the topological character of heavy fermion mixed valence compounds [16][17][18][19] and of the half-filled honeycomb lattice with an additional bond dimerization 20 .
In recent years a new class of many-body approaches has been developed to calculate the one-particle Green's function of extended systems solving the many body problem in a subsystem of finite size and embedding it within an infinite medium. These methods gather under the name of Quantum Cluster theories 21 and include Cluster Perturbation Theory 22,23 (CPT), Dynamical Cluster Approach 24 (DCA), Variational Cluser Approximation 25 (VCA), Cellular Dynamical Mean Field Theory 26 (CDMFT). They have found an unified language within the variational scheme based on the the Self Energy Functional approach 27 . These methods, with different degrees of accuracy, give access to non trivial many body effects and have been applied both to model systems and to realistic solids 28,29 .
In this paper we consider the Kane-Mele-Hubbard model 30-33 describing a 2-dimensional honeycomb lattice with both local e-e interaction and spin-orbit coupling and we adopt an approach based on CPT to determine the one-particle propagator, the topological hamiltonian G −1 (k, ω = 0) and its eigenvectors. This allows us to arXiv:1404.1287v3 [cond-mat.str-el] 2 Feb 2015 identify a general procedure that can be extended to any Quantum Cluster approach and to investigate how Green's function-based topological invariants can be effectively calculated.
The paper is organized as follows: in section I we recall how topological invariants can be obtained in terms of G −1 (k, ω = 0). Section II describes how topological invariants are obtained by CPT and section III reports the results in terms of topological invariants and spectral functions for the Kane-Mele-Hubbard model of 2D and 1D honeycomb lattices.

I. TOPOLOGICAL HAMILTONIAN AND TOPOLOGICAL INVARIANTS
In the search for an extension of topological invariants from the non-interacting to the interacting case, the Green's function has proved to be the fundamental tool [11][12][13]34 . As shown in Refs. 11 and 12 the dressed one-particle Green's function at zero frequency contains all the topological information that is required to calculate topological invariants: the inverse of the Green's function at zero frequency defines a fictitious noninteracting topological hamiltonian 35 and its eigenvectors h topo (k)|k, n, s = ns (k)|k, n, s are the quantities to be used to compute the topological invariants for the interacting system. Here n, s are band and spin indices respectively (s =↑↓). The latter is a good quantum number if -as in the model we study below-the spin orbit interaction only involves the z component of the spin. Hence, we can take the time-reversal operator to be where σ y acts on the spin indices, K denotes complex conjugation and I is the identity for the sublattice indices. The matrix w ns,n s (k) ≡ −k, n, s|Θ|k, n , s is thus a block-diagonal matrix, and is antisymmetric at time-reversal invariant momenta (TRIM) Γ i defined by the condition that −Γ i = Γ i +G with G a reciprocal lattice vector. The generalized Z 2 topological invariant can thus be defined 7,12 as the exponent ∆ in the expression and used to classify trivial insulators (∆ = 0, mod 2) from topological QSH insulators (∆ = 1, mod 2). In the presence of inversion symmetry this definition is even simpler, involving just the parity eigenvalues η n (Γ i ) = ±1 of the occupied bands at Γ i for any of the two spin sectors The definition of Z 2 for an interacting system is thus formally identical to the non-interacting case, involving in both cases the eigenstates of a single particle hamiltonian; in the presence of e-e interaction the difficult task remains the calculation of the topological hamiltonian in terms of the interacting Green's function. In the next section we will describe how this can be done within the CPT paradigm.

II. KANE-MELE-HUBBARD MODEL AND CPT
We are interested in the Kane-Mele-Hubbard model for a 2D honeycomb latticê The hopping term t il,i l (s) includes both the firstneighbor spin-independent hopping and the Haldane-Kane-Mele second-neighbor spin-orbit coupling 6,36 given by ıt KM s z (d 1 × d 2 ) z , where d 1 and d 2 are unit vectors along the two bonds that connect site il with site i l .
Here i, i run over the M atomic positions within the unit cell (cluster) and l, l refer to lattice vectors identifying the unit cells of the lattice. The on-site e-e repulsion is described by the U -Hubbard term.
In order to solve the eigenvalue problem (2), in strict analogy with what is done in any standard Tight-Binding scheme for non-interacting hamiltonians, a Bloch basis expression of the topological hamiltonian, namely of the dressed Green's function and of its inverse, is required with R l the lattice vectors (L → ∞) and r i the atomic positions inside the unit cell. (These relations hold in any spin sector and we have therefore intentionally omitted the spin index). In the following we will adopt a many body technique to calculate the one-particle dressed Green's function based on the CPT 22,23 . This method shares with other Quantum Cluster formalisms the basic idea of approximating the effects of correlations in the infinite lattice with those on a finite-size cluster. Different Quantum Cluster approaches differ for the strategy adopted to embed the cluster in the continuum and to express the lattice Green's function -or the corresponding self-energyin terms of the cluster one. The common starting point is the choice of the M -site cluster used to tile the extended lattice.
In CPT the Green's function (7) for the extended lattice is calculated by solving the equation Here G c ij is the cluster Green's function in the local basis obtained by exact diagonalization of the interacting hamiltonian for the finite cluster; we separately solve the problem for N, N-1 and N+1 electrons and express the cluster Green's function in the Lehmann representation at real frequencies. The matrix B ii (k, ω) is given by where t i 0,i l is the hopping term between site i and i belonging to different clusters. Eq. (8) is solved by a M × M matrix inversion at each k and ω. A second M × M matrix inversion is needed to obtain the topological hamiltonian according to eq. (1). The diagonalization of the topological hamiltonian is then required to obtain the eigenvectors to be used for the calculation of Z 2 according to (5). It is worth recalling that the eigenvalues of h topo in principle have nothing to do with the quasi-particle excitation energies: only the topological information is encoded in G ij (k, 0), but the full Green's function is needed to calculate quasiparticle spectral functions A(k, ω) = 1 π n Im G(k, n, ω) where with n the band index and α n i (k) the eigenstate coefficients obtained by the single-particle band calculation. 29 In the next section, analyzing in the detail all the information that can be deduced from the explicit calculation of the interacting Green's function, we will also be able to investigate more closely the relations between the eigenstates of the topological hamiltonian and the quasi-particle energies.

III. RESULTS
We have used the CPT formalism to calculate the dressed Green's function of the Kane-Mele-Hubbard model spanning a whole set of spin-orbit couplings t KM and U parameters. For the 2D honeycomb lattice the 6-site cluster ( Fig. 1 (a)) commonly used in Quantum Cluster calculations [37][38][39][40] has been adopted. In order to check the role of cluster size and geometry we have also considered the 8-site cluster of Fig. 1 (b). Both clusters represent a tiling for the honeycomb lattice but with very different cluster symmetries. Obviously this has no influence in the non-interacting case where either the "natural" 2-site unit cell or any larger unit cell (4, 6, 8 sites etc.) produce the same band structure. This is no more so if the e-e interaction is switched on: in any Quantum Cluster approach where the extended system is described as a periodic repetition of correlated units, the translation periodicity is only partially restored (it is preserved only at the superlattice level). This inevitably affects the quasi-particle band structure and for the 8-site tiling gives rise to a wrong k-dispersion. This appears quite clearly by comparing spectral functions (cfr. eq. (9)) obtained with 6-site and 8-site tilings at k-points along the border of the 2D Brillouine zone. In particular, for the 6-site tiling the quasi-particle energies display the correct symmetry, and energies at any k-point and its rotated counterpart-− → k and R − → k , with R a point group rotation-coincide (see fig. 2 (a), (c) ). This well-known basic rule is violated for the 8-site tiling and the quasiparticle energies at K, K do not coincide with the values at K and K (see fig. 2 (b), (d) ). Indeed the gap closes down around K and K but not at K and K . This is due to the fact that the 8-site tiling has a preferred direction so that the dispersions along K − K and K − K are different.
The dependence on the cluster's size and symmetry of the quasiparticle band dispersion in the Kane-Mele-Hubbard model is shown here for the first time and appears to be crucial in order to identify the accuracy and appropriateness of the results: any band structure of a non-interacting system violating the point group symmetries should be disregarded as wrong and unphysical; the same should be done for interacting systems since e-e repulsion does not affect the lattice point group symmetry. For this reason the semimetal behaviour for t KM = 0 and U/t ≤ 3.5 that we find for the 8-site tiling in agreement with Refs. 38 and 43 should be considered an artifact due to the wrong cluster symmetry and not to be used to infer any real behaviour of the model system. Only clusters that preserve the point group symmetries of the lattice should be used 21 and this criterion restricts the choice for the 2D honeycomb lattice to the 6-site cluster. 41 . These considerations are quite general and quasiparticle states, topological invariants and phase diagrams obtained by Quantum Cluster approaches using tilings with the wrong symmetry (2-, 4-, 8-site clusters) 38,42,43 are for this reason not reliable.
We focus now on the Green's function at ω = 0 and on the topological properties that can be deduced from it. As discussed in the previous section, the key quantity is the dressed Green's function expressed in a Bloch ba-  (Color online) Comparison between the quasiparticle band structures obtained for 6-site and 8-site tiling assuming U/t = 2. Panel (a) and (b) ( (c) and (d) ) show the results obtained with tKM /t = 0 ( tKM /t = 0.1) for 6-site and 8-site tiling respectively. Notice that the gapless band structure obtained for tKM /t = 0 with the 8-site tiling is a consequence of a band dispersion that violates the rotational symmetry of the lattice. sis (eqs. (7) and (8)) and the corresponding topological hamiltonian (h topo ) ij = −G −1 ij (k, ω = 0). The eigenvalue problem associated to h topo -a 6 × 6 matrix diagonalization at each k-point-is equivalent to a standard single particle Tight-Binding calculation for a unit cell containing 6 atomic sites, giving rise to 6 topological bands. The product over the first 3 occupied bands at the TRIM points corresponding to the 6-site cluster provides, according to eq. (5), the Z 2 invariant. Fig. 3 reports the resulting U − t KM phase diagram showing the parameter range where the system behaves as either a topologically trivial insulator (TTI, ∆ = 0) or a Quantum Spin Hall insulator (QSH, ∆ = 1) 44 .
Few comments are in order: since we are monitoring the topological phase transition by the Z 2 invariant, we are implicitly assuming both time reversal and parity invariance and an adiabatic connection between QSH and TTI phase to persist in all regimes. Anti-ferromagnetism that breaks both time-reversal and sublattice inversion symmetry is therefore excluded and we are assuming the system to remain non-magnetic at any U . In this sense the value of Z 2 can be considered an indicator of the topological properties of the system of interest only for a parameter range that excludes antiferromagnetism.
The behavior for t KM close to zero is worth noticing. According to our CPT calculation, at low t KM the QSH regime does not survive the switching on of e-e interaction: a value U → 0 is enough to destroy the semimetallic behavior at t KM = 0; see Fig. 4. This is at variance with Quantum Monte Carlo (QMC) results 45,46 that at t KM = 0 identify a semimetallic behavior up to U/t ∼ 3.5 47 . It has been recently shown 48 that the existence at t KM = 0 of an excitation gap down to U → 0 is characteristic of all Quantum Cluster schemes with the only exception of DCA. This is due to the aforementioned violation of translational symmetry in Quantum Cluster methods such as CPT, VCA and CDMFT, regardless of the scheme being variational or not, and independent on the details of the specific implementations (different impurity solvers, different temperatures). We stress here again that the semimetal behaviour that is found for the Kane-Mele-Hubbard model by Quantum Cluster approaches such as CDMFT 38 and VCA 43 is actually an artifact due to the choice of clusters with wrong symmetry. The only Quantum Cluster approach that is able to reproduce a semimetal behaviour at finite U is DCA. DCA preserves translation symmetry and has been shown to describe better the small U regime; it becomes however less accurate at large U values where it overemphasizes the semimetallic behavior of the honeycomb lattice. 48 In this sense DCA and the other Quantum Cluster approaches can be considered as complementary and it would be interesting to compare their results also in terms of parity invariants.
By calculating spectral functions in the same parameter range we observe that at the transition points the single particle excitation gap ∆ sp , namely the minimum energy separation between hole and particle excitations, closes down. The possibility for e-e interaction to induce a metallic behavior in a band insulator has been recently analysed by DMFT 49,50 and QMC calculations 51 ; here we observe that, in agreement with previous results 37, 38 , the same effect occurs in the honeycomb lattice made semiconducting by intrinsic spin-orbit interaction. Fig.  4 shows the behavior of ∆ sp at different values of t KM /t as a function of U/t. Fig. 5 (a-c) shows as an example the quasi-particle band structure obtained for t KM /t = 0.1 and for U/t = 2 (QSH regime), U/t = 3.5 (transition point from QSH to TTI, the gap closes down) and U/t = 4 (TTI regime).
Other effects are due to the e-e correlation, namely a band width reduction and the appearance of satellite structures below (above) valence (conduction) band. These effects are more clearly seen by looking at the density of states (DOS) obtained as the sum of the spectral functions over a large sample of k-points (Fig. 6).
Even if the eigenvalues of h topo cannot be identified with excitation energies, they exhibit a behavior similar to the quasi-particle energies. In particular the same gap closure appears in h topo eigenvalues at the transition  points. This is shown in Fig. 5 (d-f) where a zoom of the quasi-particle band structure around the point K is compared with the eigenvalues nk↑ of h topo .
We may then conclude, in agreement with QMC calculations 20,52 , that a change in the Z 2 invariant is associated to the closure of both the single-particle excitation gap and of the energy separation between filled and empty states of h topo : in strict analogy with the non-interacting case, a change of topological regime of the interacting systems is associated to a gap closure followed by a gap inversion in the fictitious band structure associated to h topo .
According to the bulk-boundary correspondence, 1D non-interacting systems should exhibit gapless edge states once the 2D system enters the QSH regime. In  (Color online) Density of states of the half filled Kane-Mele-Hubbard model for tKM /t = 0.1 at U/t = 2 (red), 3.5 (blue), 4 (green) compared with the non interacting results (dashed line). Satellite structures appear below and above the valence and conduction band respectively and the band width is reduced, an effect that is more evident for larger U .
the presence of e-e interaction this may not be true and a gap may open in edge states before the time-reversal Z2 invariant switches off 31 . We have calculated within CPT the spectral functions for a honeycomb ribbon with zigzag termination using the tiling shown in Fig. 1 (c). For any given value of t KM we have systematically found that gapless edge states persist up to a critical value of U that coincides with the one previously identified as the transition point from QSH to TTI regime in the 2D system. This is shown in Fig. 7 for t KM /t = 0.1. Here we notice that at critical value U/t = 3.5 a tiny gap exists between filled and empty states. We have checked that, increasing the ribbon width, this gap becomes smaller and smaller, and we may then attribute it to the finite width of the ribbon. In conclusion we have studied the topological properties of the Kane-Mele-Hubbard model by explicitly calculating the Green's function-based topological invariant. The CPT scheme, using Bloch sums as a basis set, nat-urally leads to the topological hamiltonian matrix that enters in the computation of the Z 2 topological invariant. The approach gives direct access to the dressed Green's function at real frequencies, avoiding the problem of analytic continuation, and does not require the extraction of a self-energy. We have shown that the interplay between the Hubbard interaction and the intrinsic spin-orbit coupling is coherently described by the values of Z 2 invariant and by 2D and 1D quasi-particle energies (gap closing at the transitions, edge states in the QSH phase). We have discussed the importance of the cluster symmetry and the effects of the lack of full translation symmetry typical of CPT and of most Quantum Cluster approaches. Comments on the limits of applicability of the method are also provided.