Topological invariants in interacting quantum spin Hall: a cluster perturbation theory approach

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 (2D) honeycomb lattice with spin–orbit coupling has been identified as a remarkable and paradigmatic example of a 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 ($\mathcal{T}$)—topological invariant ${{\mathbb{Z}}_{2}}$ [6–8]. In the same way as the Thouless–Kohmoto–Nightingale–den Nijs (TKNN) [9] topological invariant was defined for the integer quantum Hall effect, the above ${{\mathbb{Z}}_{2}}$ invariant was defined for the topological insulator in terms of band eigenvectors and, as such, only applies to non-interacting 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–13], and unlike the former, the latter construction can be straightforwardly extended to the ${{\mathbb{Z}}_{2}}$ invariant to classify interacting QSH systems.

The field of interacting topological insulators is attracting growing interest (see [14, 15] for recent reviews on 2D 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 promising is the one developed in [12] where it has been demonstrated that topological invariants are determined by the behavior 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,\omega )$ 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–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 (CPT) [22, 23], dynamical cluster approach (DCA) [24], variational cluser approximation (VCA) [25], cellular dynamical mean field theory (CDMFT) [26]. 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 2D 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,\omega =0)$ and its eigenvectors. This allows us to 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 2 we recall how topological invariants can be obtained in terms of ${{G}^{-1}}(k,\omega =0)$. Section 3 describes how topological invariants are obtained by CPT and section 4 reports the results in terms of topological invariants and spectral functions for the Kane–Mele–Hubbard model of 2D and 1D honeycomb lattices.

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–13, 34]. As shown in [11, 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 non-interacting topological Hamiltonian [35]

$$\begin{eqnarray}&&{{h}_{{\rm topo}}}(k)\equiv -{{G}^{-1}}(k,0)\end{eqnarray} \tag{ 1 }$$

and its eigenvectors

$$\begin{eqnarray}&&{{h}_{{\rm topo}}}(k)\mid k,n,s\rangle ={{\epsilon }_{ns}}(k)\mid k,n,s\rangle \end{eqnarray} \tag{ 2 }$$

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=\uparrow \downarrow$). 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

$$\begin{eqnarray*}&&\Theta =\mathbb{I}\;\otimes \;{\rm i}{{\sigma }_{y}}\;K,\end{eqnarray*}$$

where ${{\sigma }_{y}}$ acts on the spin indices, K denotes complex conjugation and $\mathbb{I}$ is the identity for the sublattice indices. The matrix

$$\begin{eqnarray}&&{{w}_{ns,n^{\prime} s^{\prime} }}(k)\equiv \langle -k,n,s\mid \Theta \mid k,n^{\prime} ,s^{\prime} \rangle \end{eqnarray} \tag{ 3 }$$

is thus a block-diagonal matrix, and is antisymmetric at time-reversal invariant momenta (TRIM) ${{\Gamma }_{i}}$ defined by the condition that $-{{\Gamma }_{i}}$ = ${{\Gamma }_{i}}+\mathcal{G}$ with $\mathcal{G}$ a reciprocal lattice vector. The generalized ${{\mathbb{Z}}_{2}}$ topological invariant can thus be defined [7, 12] as the exponent Δ in the expression

$$\begin{eqnarray}&&{{(-1)}^{\Delta }}\equiv \mathop{\prod }\limits_{{\rm TRIM}}\frac{\sqrt{{\rm det} [w({{\Gamma }_{i}})]}}{{\rm P}{{{\rm f}}^{{}}}[w({{\Gamma }_{i}})]}\end{eqnarray} \tag{ 4 }$$

and used to classify trivial insulators ($\Delta =0$, mod 2) from topological QSH insulators ($\Delta =1$, mod 2). In the presence of inversion symmetry this definition is even simpler, involving just the parity eigenvalues ${{\eta }_{n}}({{\Gamma }_{i}})=\pm 1$ of the occupied bands at ${{\Gamma }_{i}}$ for any of the two spin sectors

$$\begin{eqnarray}&&{{(-1)}^{\Delta }}=\mathop{\prod }\limits_{{\rm TRIM}}{{\mathop{\prod }\limits_{n=1}^{N}}^{{}}}{{\eta }_{n}}({{\Gamma }_{i}}).\end{eqnarray} \tag{ 5 }$$

The definition of ${{\mathbb{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.

We are interested in the Kane–Mele–Hubbard model for a 2D honeycomb lattice

$$\begin{eqnarray}&&\hat{H}=\mathop{\sum }\limits_{il,i^{\prime} l^{\prime} s}{{t}_{il,i^{\prime} l^{\prime} }}(s)\hat{c}_{ils}^{\dagger }{{\hat{c}}_{i^{\prime} l^{\prime} s}}+U\mathop{\sum }\limits_{il}\hat{c}_{il\uparrow }^{\dagger }{{\hat{c}}_{il\uparrow }}\hat{c}_{il\downarrow }^{\dagger }{{\hat{c}}_{il\downarrow }}.\end{eqnarray} \tag{ 6 }$$

The hopping term ${{t}_{il,i^{\prime} l^{\prime} }}(s)$ includes both the first-neighbor spin-independent hopping and the Haldan–Kane–Mele second-neighbor spin–orbit coupling [6, 36] given by $\imath {{t}_{{\rm KM}}}{{s}_{z}}{{({{d}_{1}}\times {{d}_{2}})}_{z}}$, where d₁ and d₂ are unit vectors along the two bonds that connect site $il$ with site $i^{\prime} l^{\prime}$. Here $i,i^{\prime}$ run over the M atomic positions within the unit cell (cluster) and $l,l^{\prime}$ 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

$$\begin{eqnarray}&&{{G}_{ij}}(k,\omega )=\langle {{\Psi }_{0}}\mid \hat{c}_{ki}^{\dagger }\hat{G}{{\hat{c}}_{kj}}\mid {{\Psi }_{0}}\rangle +\langle {{\Psi }_{0}}\mid {{\hat{c}}_{ki}}\hat{G}\hat{c}_{kj}^{\dagger }\mid {{\Psi }_{0}}\rangle ,\end{eqnarray} \tag{ 7 }$$

where $\hat{G}$ = $\frac{1}{\omega -\hat{H}}$ and

$$\begin{eqnarray*}&&\hat{c}_{ki}^{\dagger }=\frac{1}{\sqrt{L}}\mathop{\sum }\limits_{l}^{L}{{{\rm e}}^{-{\rm i}k\cdot ({{R}_{l}}+{{r}_{i}})}}\hat{c}_{li}^{\dagger };\ \ {{\hat{c}}_{ki}}=\frac{1}{\sqrt{L}}\mathop{\sum }\limits_{l}^{L}{{{\rm e}}^{{\rm i}k\cdot ({{R}_{l}}+{{r}_{i}})}}{{\hat{c}}_{li}}\end{eqnarray*}$$

with R_l the lattice vectors (L $\to \;\infty$) 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-energy—in 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

$$\begin{eqnarray}&&{{G}_{ij}}(k,\omega )=G_{ij}^{c}(\omega )+\mathop{\sum }\limits_{{{i}^{\prime }}}^{M}{{B}_{ii^{\prime} }}(k,\omega ){{G}_{i^{\prime} j}}(k,\omega ).\end{eqnarray} \tag{ 8 }$$

Here $G_{ij}^{c}$ 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^{\prime} }}(k,\omega )$ is given by

$$\begin{eqnarray*}&&{{B}_{ii^{\prime} }}(k,\omega )=\mathop{\sum }\limits_{l}^{L}{{{\rm e}}^{{\rm i}k\cdot {{R}_{l}}}}\mathop{\sum }\limits_{{{i}^{\prime\prime }}}^{M}G_{i{{i}^{\prime\prime }}}^{c}(\omega ){{t}_{{{i}^{\prime\prime }}0,i^{\prime} l}}(s),\end{eqnarray*}$$

where ${{t}_{{{i}^{\prime\prime }}0,i^{\prime} l}}$ is the hopping term between site i' and i'' belonging to different clusters.

Equation (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 equation (1). The diagonalization of the topological Hamiltonian is then required to obtain the eigenvectors to be used for the calculation of ${{\mathbb{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 quasi-particle spectral functions

$$\begin{eqnarray}&&A(k,\omega )=\frac{1}{\pi }\mathop{\sum }\limits_{n}{\rm Im}\;G(k,n,\omega ),\end{eqnarray} \tag{ 9 }$$

where

$$\begin{eqnarray*}&&G(k,n,\omega )=\frac{1}{M}\mathop{\sum }\limits_{ii^{\prime} }{{{\rm e}}^{-{\rm i}k\cdot ({{r}_{i}}-{{r}_{{{i}^{\prime }}}})}}\alpha _{i}^{n*}(k)\alpha _{{{i}^{\prime }}}^{n}(k){{G}_{ii^{\prime} }}(k,\omega )\end{eqnarray*}$$

with n the band index and $\alpha _{i}^{n}(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.

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 six-site cluster (figure 1 (a)) commonly used in quantum cluster calculations [37–40] has been adopted. In order to check the role of cluster size and geometry we have also considered the eight-site cluster of figure 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' two-site unit cell or any larger unit cell (four, six, eight 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 eight-site tiling gives rise to a wrong k-dispersion. This appears quite clearly by comparing spectral functions (cf equation (9)) obtained with six-site and eight-site tilings at k-points along the border of the 2D Brillouin zone. In particular, for the six-site tiling the quasi-particle energies display the correct symmetry, and energies at any k-point and its rotated counterpart—$\vec{k}$ and $R\vec{k}$, with R a point group rotation—coincide (see figures 2 (a), (c) ). This well-known basic rule is violated for the eight-site tiling and the quasi-particle energies at K, $K$ do not coincide with the values at K' and K'' (see figures 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 eight-site tiling has a preferred direction so that the dispersions along $K-K^{\prime}$ and $K^{\prime} -K^{\prime\prime}$ are different.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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 behavior for ${{t}_{{\rm KM}}}=0$ and $U/t\leqslant 3.5$ that we find for the eight-site tiling in agreement with [38, 42] should be considered an artifact due to the wrong cluster symmetry and not to be used to infer any real behavior 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 six-site cluster⁵ . These considerations are quite general and quasi-particle states, topological invariants and phase diagrams obtained by quantum cluster approaches using tilings with the wrong symmetry (two-, four-, eight-site clusters) [38, 41, 42] are for this reason not reliable.

We focus now on the Green's function at $\omega =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 basis (equations (7) and (8)) and the corresponding topological Hamiltonian ${{({{h}_{{\rm topo}}})}_{ij}}=-G_{ij}^{-1}(k,\omega =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 six atomic sites, giving rise to six topological bands. The product over the first three occupied bands at the TRIM points corresponding to the six-site cluster provides, according to equation (5), the ${{\mathbb{Z}}_{2}}$ invariant. Figure 3 reports the resulting $U-{{t}_{{\rm KM}}}$ phase diagram showing the parameter range where the system behaves as either a topologically trivial insulator (TTI, $\Delta =0$) or a (QSH, $\Delta =1$)⁶ .

Zoom In Zoom Out Reset image size

Few comments are in order: since we are monitoring the topological phase transition by the ${{\mathbb{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 ${{\mathbb{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\to 0$ is enough to destroy the semi-metallic behavior at ${{t}_{{\rm KM}}}=0$; see figure 4. This is at variance with quantum Monte Carlo (QMC) results [43, 44] that at ${{t}_{{\rm KM}}}=0$ identify a semimetallic behavior up to $U/t\sim 3.5$ ⁷ . It has been recently shown [45] that the existence at ${{t}_{{\rm KM}}}=0$ of an excitation gap down to $U\to 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 behavior that is found for the Kane–Mele–Hubbard model by quantum cluster approaches such as CDMFT [38] and VCA [42] 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 behavior 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 [45]. 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.

Zoom In Zoom Out Reset image size

By calculating spectral functions in the same parameter range we observe that at the transition points the single particle excitation gap ${{\Delta }_{{\rm 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 [46, 47] and QMC calculations [48]; 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. Figure 4 shows the behavior of ${{\Delta }_{{\rm sp}}}$ at different values of ${{t}_{{\rm KM}}}/t$ as a function of $U/t$. Figures 5 (a)–(c) shows as an example the quasi-particle band structure obtained for ${{t}_{{\rm 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).

Zoom In Zoom Out Reset image size

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 obtained as the sum of the spectral functions over a large sample of k-points (figure 6).

Zoom In Zoom Out Reset image size

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 figures 5 (d)–(f) where a zoom of the quasi-particle band structure around the point K is compared with the eigenvalues ${{\epsilon }_{nk\uparrow }}$ of h_topo.

We may then conclude, in agreement with QMC calculations [20, 49], that a change in the ${{\mathbb{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 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 figure 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 figure 7 for ${{t}_{{\rm 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.

Zoom In Zoom Out Reset image size

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, naturally leads to the topological Hamiltonian matrix that enters in the computation of the ${{\mathbb{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 ${{\mathbb{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.