Topological insulator on honeycomb lattices and ribbons without inversion symmetry

We study the Kane-Mele-Hubbard model with an additional inversion-symmetry-breaking term. Using the topological Hamiltonian approach, we calculate the $\mathbb{Z}_2$ invariant of the system as function of spin-orbit coupling, Hubbard interaction $U$, and inversion-symmetry-breaking on-site potential. The phase diagram calculated in that way shows that, on the one hand, a large term of the latter kind destroys the topological non-trivial state. On the other hand, however, this inversion-symmetry-breaking field can enhance the topological state, since for moderate values the transition from the non-trivial topological to the trivial Mott insulator is pushed to larger values of interaction $U$. This feature of an enhanced topological state is also found on honeycomb ribbons. With inversion symmetry, the edge of the zigzag ribbon is magnetic for any value of $U$. This magnetic moment destroys the gapless edge mode. Lifting inversion symmetry allows for a finite region in interaction strength $U$ below which gapless edge modes exist.

However, the influence of interactions onto the topological classification is still not fully understood. Just recently, new phase transitions in strongly correlated topological insulators have been reported [22,23]. The most used quantity to characterize topological order, namely the Z 2 invariant introduced by Fu, Kane, and Mele [1,[24][25][26], relies on defined Bloch bands and is thus not directly applicable for interacting systems. A generalization is possible using the so-called topological Hamiltonian [27][28][29], an artificially noninteracting system determined by the Green's function.
Determining the topological phase becomes more difficult if an inversion-symmetry-breaking term such as a staggered on-site potential [1,43], a Rashba coupling [1,39], or site-dependent hoppings [36,42] are included. A possibility to analyze topological phases is to * Robert Triebl: robert.triebl@tugraz.at calculate the spin Chern number C S [36,[42][43][44][45][46]. This approach requires spin to be a good quantum number and has the drawback that due to numerical artifacts a good quantization of C S is not given close to phase transitions. Another approach is to look directly for gapless edge states and use bulk-boundary correspondence [37,[39][40][41].
In this paper, we calculate the Z 2 invariant of the KMH model with an inversion-symmetry-breaking on-site potential by combining the topological Hamiltonian with a method introduced by Soluyanov and Vanderbilt [47,48] that is based on maximally localized Wannier charge centers [49]. This enables a precise calculation of invariants without restricting the systems to certain symmetries. Furthermore, we investigate bulk-boundary correspondence by calculating the spectral functions of a zigzag ribbon. We show that bulk-boundary correspondence has to be treated with care in strongly interacting systems since time-reversal symmetry might be lifted locally at the edges due to spontaneous symmetry breaking. The Green's functions in our approach are obtained by a twosite dynamical impurity approximation [50][51][52][53][54].

II.1. Kane-Mele-Hubbard model
The Kane-Mele-Hubbard Hamiltonian is used exemplary since it is a toy model for strongly correlated topological insulators. The noninteracting part as proposed by Kane and Mele [1,2] is given by on a honeycomb lattice, where c † i is the creation operator of a spinor c † i↑ , c † i↓ , · denotes nearest neighbors, and arXiv:1605.07855v2 [cond-mat.str-el] 13 Dec 2016 · next-nearest neighbors. The first term is a tightbinding nearest-neighbor hopping term, which is commonly used to model the Dirac cones of graphene up to first order. The second is the intrinsic spin-orbit coupling, leading to the quantum spin Hall topological insulating state as it opens a gap [1,2,45]. The third term is a staggered on-site potential, where ξ i = 1 if site i belongs to sublattice A of the honeycomb lattice, and −1 if it belongs to sublattice B. This distinction between the two sublattices breaks inversion symmetry and has a crucial influence on topology: the KM model with a sublattice potential is a topological insulator for any finite λ SO , as long as |λ ν | < 3 √ 3λ SO . The gap closes for |λ ν | = 3 √ 3λ SO and reopens for |λ ν | > 3 √ 3λ SO , but the topology becomes trivial in that case.
Interaction effects can be introduced by a Hubbard interaction U n ↑ n ↓ on each site, leading to the Kane-Mele-Hubbard Hamiltonian [30] Throughout the paper, the energy scale is defined by t ≡ 1, and the length scale by the lattice parameter a ≡ 1.

II.2. Calculation of topological invariants
Topological systems are classified by their dimension and symmetries, as summarized in the periodic table of topological matter [3][4][5]. The important symmetry in case of the KM model is time reversal, leading to the topological class AII, specified by a Z 2 invariant ν. A possibility to define this invariant is via time-reversal polarizations of a one-dimensional system that depends on an additional pumping parameter, as introduced by Fu and Kane [24]. In case of a noninteracting two-dimensional system, this definition is applicable if k y is considered as the pumping parameter. For the actual calculation of ν, inversion-symmetric and non-inversion-symmetric systems are treated differently, as discussed in the following.
If inversion symmetry is present, the four TRIMs Γ i contain the whole topological information. The Z 2 invariant ν can be obtained by where ξ n (Γ i ) is the eigenvalue of the parity operator at momentum k = Γ i of Kramer's pair n [25]. If, on the other hand, inversion symmetry is broken, one needs information on how the Bloch states evolve continuously between the TRIMs. Soluyanov and Vanderbilt suggested [47,48] to use hybrid Wannier functions which are maximally localized [49]. The topology is determined by tracking the maximally localized Wannier charge centers (WCCs) assigned to the occupied bands along the pumping parameter k y , which are given bȳ x n (k y ) = 0k y n| x |0k y n . This function is defined modulo a lattice constant that is chosen to be 1, sox(k y ) has a periodicity of 2π in k y , and a period of 1 alonḡ x. The KM model has only two occupied bands which form a Kramers pair because of time-reversal invariance, so Kramer's degeneracy enforces the two WCCs to be equal at k y = 0 and π. Tracking the WCCs continuously from k y = 0 to 2π, the system is trivial if the very same WCCs intersect at both points, and nontrivial if there is a shift which is a multiple of the lattice constant. Examples are given in Fig. 1. If the spin in z direction S z is conserved, the continuous tracking is straight forward since each WCC can be assigned to a certain spin. If no conserved quantity helps identifying the respective WCCs, the two cannot be sorted and some more advanced method has to be applied, as for example tracking the difference of the WCCs [47].
We now turn to the determination of topological states for a system with electron-electron interactions. Here, topological invariants cannot be defined as described above since one-electron Bloch functions are not eigenstates. A more general definition of the first Chern num-ber uses Green's functions [55][56][57], with k 0 = iω, which gives the integer coefficient of the quantum Hall effect of a two dimensional system. If spin is a good quantum number, the Chern invariant can be evaluated separately for each spin. C ↑ is then evaluated from the spin up block of the Green's function, C ↓ from the spin down block. This leads to a quantized spin Chern number C S = (C ↑ − C ↓ )/2, which is integer for time-reversal invariant Hamiltonians. Modulo 2, this quantity can be used as a Z 2 invariant. In the general case, a Z 2 invariant ν is obtained from a dimensional reduction of the second Chern number which describes the response of a four dimensional insulator [58,59]. Starting from definitions (5) and (6), Wang et al. showed that the topological information is already captured in the Green's function at zero frequency [27][28][29]. They conclude that minus the inverse Green's function at zero frequency can be considered as the Bloch Hamiltonian of an artificial noninteracting system which has the same Chern invariant and the same Z 2 invariant as the interacting one, as long as they are continuously connected. Thus, this Bloch Hamiltonian is called topological Hamiltonian [29] H t (k) = −G −1 (ω = 0, k) of the interacting system. A main consequence is that methods devised for noninteracting Hamiltonians are sufficient to calculate topological numbers related to the more complicated integrals (5) and (6), as for example C S and ν. A direct evaluation of (5) or (6) is therefore not necessary.
If the system obeys inversion symmetry, G −1 (ω, k) commutes at the TRIMs k = Γ i with the parity transformation matrix P , and as a consequence, ther are simultaneous eigenstates |α(ω, Γ i ) of G −1 and P : The topological invariant ν can be calculated from these eigenvalues through [27] (−1) ν = R zeros Here, the convention (−1) 1/2 = +i is used. In the noninteracting case, this equation reduces to the Fu-Kane formula (3) [27]. The direct evaluation of topological invariants through Eq. (9) became already a standard procedure in case of interacting systems with inversion symmetry [22,23,35,38,39]. In this work, we are interested in the topological invariants of an interacting system without inversion symmetry, where Eq. (9) cannot be applied. For this case we propose to use a combination of the topological Hamiltonian with the Soluyanov-Vanderbilt method of WCCs as described in the beginning of this section. In practice, we first calculate the Green's function at zero frequency using a dynamical impurity approximation as explained in the next section. The obtained topological Hamiltonian can then be used just like a Bloch Hamiltonian to determine the Z 2 invariant. This in turn is done with Wannier charge centers as proposed by Soluyanov and Vanderbilt [47], just using the eigenstates of the topological Hamiltonian |α(ω = 0, k) instead of the Bloch functions |ψ nk of the noninteracting case.

II.3. Variational cluster approach
As described in the previous section, the one-electron Green's function is needed to determine the topological Hamiltonian. Since an exact solution of the full many-body problem is not possible, an approximative method has to be chosen. Here we apply the Variational cluster approach (VCA) [50,52], because the Kane-Mele-Hubbard model is known to have an antiferromagnetic moment [30-33, 35, 38, 39] which can efficiently be treated by the VCA with symmetry-breaking Weiss fields [53,54].
The VCA is based on the self-energy functional approach, which uses the fact that the grand potential of an arbitrary interacting system H = H 0 (t) + H 1 (U ) has to be a stationary point of the self-energy functional where F [Σ] denotes the Legendre transform of the Luttinger-Ward functional Φ[G] [50,60]. The approximation of this method is to restrict the space of selfenergies Σ. This subset S of self-energies is spanned by all Σ(t ) that are the exact self-energies of a so-called reference system H = H 0 (t ) + H 1 (U ). The interaction parameters U are the same as in the original system, but H and H can differ in the one-particle parameters. The one-particle parameters t of the reference system H are chosen such that the self-energy of the reference system can be calculated exactly. To obtain the approximative physical self-energy Σ ∈ S, a stationary point of Ω t [Σ(t )] has to be found as t is varied. The parametrized functional can be reduced to and can thus be calculated if the Green's function of the reference system is known. Quite generally, reference systems in the VCA are clusters of finite size, which can be treated by exact diagonalization techniques [50,52,53]. In case of the KMH model, several cluster sizes have already been analyzed [37][38][39]. However, the tiling of the lattice into clusters of finite sizes breaks artificially some symmetries, which can change the topological phase diagram [61]. That is why we choose as a reference system for VCA single-site clusters, which are coupled to one additional bath site by a hopping V . This rather simple approach, called two-site dynamical impurity approximation (DIA) [51], has two advantages. First, despite its simplicity, it gives accurate results for the transition towards an antiferromagnetic insulator for two-dimensional Hubbard models [51]. Second, which is even more important, the lattice symmetries are trivially satisfied. A drawback of this method is the locality of the self-energy. We will show below, however, that for known cases we get very good agreement with existing results obtained by numerically much more expensive methods.
Since the honeycomb lattice has two distinct sites, the unit cell is tiled by two clusters, which are coupled by the noninteracting part of the Hamiltonian, as shown in Fig. 2. On-site energies on both impurity and bath site, as well as the connecting hopping between them, give in total three variational parameters per cluster. However, in the inversion-symmetric case (λ ν = 0), the on-site energies are fixed by particle-hole symmetry and only one parameter remains.
In order to capture symmetry breaking necessary for the emerging antiferromagnetic moment, a Weiss field has to be added [54]. Without any symmetry considerations, these fields on both A and B sites give in total 6 variational parameters. Due to the inversion-symmetry breaking on-site potential λ ν , a second Weiss field is used to enable unequal electron densities on the two sublattices. As in Eq. (1), ξ i = ±1, depending on the sublattice. This Weiss field is basically a renormalisation of λ ν in the reference system, which is caused by the interplay of the sublattice potential and Hubbard interaction.
The method described so far considers bulk properties. Introducing an edge destroys translational symmetry and influences therefore local magnetization. As known from field theoretical investigations, mean-field approximation gives a finite magnetization on the zigzag edge for every finite interaction strength [31]. This could lead to a breakdown of the bulk-boundary correspondence and may cause problems for calculating topological invariants using the existence of gapless edge states as a proof for nontrivial topology, which has so far been used in some cases of interacting systems without inversion symmetry [37,39,41]. Vice versa, a nontrivial topological invariant in the bulk may not result in gapless edge states due to locally broken time-reversal symmetry caused by spontaneous symmetry breaking. Therefore, we additionally implemented the DIA on the zigzag ribbon in order to compare the topological invariants defined by the bulk Green's function to the existence of gapless edge states. The ribbon is translationally invariant in the x direction, whereas the sites along the width of the ribbon are distinct. If a unit cell contains N pairs of A and B sites, 2N clusters containing each a bath and an impurity site have to be solved and effectively coupled by the noninteracting part of the Hamiltonian (see Fig. 3). In order to keep the number of parameters manageable, the on-site energies and hybridisations are chosen to be constant along the ribbon. To allow for edge magnetization, the antiferromagnetic Weiss fields for each pair of sites A and B is varied independently, only assuming a mirror symmetry y → −y.

III.1. Bulk
As mentioned in the methods section, the hopping to the bath sites, the magnetic Weiss fields, and the sublattice potential Weiss field have to be determined in the VCA. For all stationary points, the ferromagnetic part of the Weiss field vanishes, hence only an antiferromagnetic ordering h A = −h B is possible. Without spin-orbit coupling, the system has full SU(2) symmetry, so only the absolute value of the Weiss field has to be determined. When spin-orbit coupling is included, only the xy-plane is still degenerate, but the degeneracy of the z direction is lifted. This means that we have to deal with two antiferromagnetic Weiss fields, h z and h x . To analyze the direction of the antiferromagnetic moment, we calculate a two-dimensional surface of the self-energy functional Ω(h z , h x ), where all other variational parameters are optimized for each set of variables (h z , h x ). The stationary points, i.e. extrema and saddle points, of this twodimensional surfaces are physical solutions, where the stable solution is the one with lowest potential Ω. Fig. 4 shows the value of the self-energy functional as a function of both in-plane and out-of-plane AF symmetry-breaking field. Depending on the KMH model parameters, up to three different stationary points exist: A saddle point of Ω if h points in z direction; a minimum if it is in the xy plane; the nonmagnetic solution, which can be both maximum or minimum, depending on the parameters. This is consistent with the results of other cluster geometries [37,39]. The local minimum h ẑ is never the physically realized solution with the lowest grand potential Ω for all sets of parameters considered here. Hence, only one variational quantity is needed for the AF Weiss field, namely the in-plane antiferromagnetic component. As mentioned above, the on-site energy levels of both impurity and bath are fixed by particle hole symmetry and the given chemical potential. Therefore, in total three cluster parameters have to be optimized: The hopping V between impurity and bath, the in-plane antiferromagnetic Weiss field h x , and the potential difference between the two sublattices ∆.
Directly from the two-site DIA one can distinguish two phases, the antiferromagnetic insulator for large U and the nonmagnetic insulator for small U . The system reduces to the ordinary Hubbard model on the honeycomb lattice if λ SO = 0 and λ ν = 0. In this case, the magnetization direction is not important since SU (2) symmetry is not broken. The mean-field critical interaction is U c = 2.23 [30,62], which is lower as compared to more accurate methods. Quantum Monte Carlo simulations show that it is actually slightly above 4 [32][33][34][35]62]. The two-site DIA considered in this work is expected to give similar results as other variational methods. VCA gives critical interactions between 2.4 and 4, depending on the cluster geometries [37][38][39], which coincides with our DIA results of U c = 3.7, where we observe a second order phase transition. With increasing λ SO , all methods show that U c increases as well. Mean-field [30], however, overestimates here the slope in comparison with the more elaborate methods [32][33][34][35][37][38][39]. The reason for that is analyzed in the Appendix A. Our results show a similar behaviour as VCA with different cluster geometries [39].
To sum up, in the inversion-symmetric case the two-site DIA is in good agreement with other methods. We can therefore expect that the method is suitable to explore the model when inversion symmetry is broken.
Using the topological Hamiltonian defined in Eq. (7) in combination with the Soluyanov-Vanderbilt method, information on the topological properties can be obtained in addition to the magnetic ordering. In the noninteracting case, a topological phase transition occurs at λ ν = 3 √ 3λ SO , as known from the original work by Kane and Mele [1,2]. Including a Hubbard interaction U , the topological Hamiltonian has the same structure as the noninteracting Hamiltonian, as long as the antiferromagnetic moment vanishes. However, both self-energy and staggered on-site Weiss field renormalize the energy scales. The interplay of interaction and on-site energy can be seen as follows: Without interaction, the sublattice with the lower on-site energy has a higher double occupancy. A finite Hubbard U punishes double occupancies, and reduces as a result the double occupancy on the sublattice with lower on-site energy. Hence, the sublattice potential λ ν is effectively lowered in case of a finite U , stabilizing the topological phase, and shifting the critical λ ν to higher values. The resulting phase diagram is shown in Fig. 5. This stabilization effect is also captured in mean-field, although with quantitative differences [43]. We want to note that we cross-checked the validity of our WCC approach by calculating the spin Chern number C S directly from Eq. (5) for the selected value of U = 1. We found perfect quantitative agreement. This reasoning for the stabilisation of the topological phase is only valid in case of weak interactions where the antiferromagnetic Weiss field is zero. In the strongly interacting regime, the non-vanishing Weiss field causes a time-reversal symmetry breaking term proportional to σ x ⊗ τ z (σ acts in spin space, τ in sublattice space) in the topological Hamiltonian. As a consequence, the topological invariant in the sense of Fu and Kane [25] is not defined. This can also be seen in the WCC, where the lifted Kramer's degeneracy does not enforce the two WCCs to be identical at half the period of the pumping parameter. Examples of the WCCs are shown in Fig. 6. In this regime, not just quantitative, but also qualitative differences compared to a standard Hartree-Fock mean-field arise, as discussed in Appendix A. To sum up, three phases exist for a given spin-orbit coupling: (i) a topological insulator continuously connected to the quantum spin Hall phases of the non-interacting KM-model if both λ ν and U are small enough; (ii) a trivial band insulator if λ ν is large; (iii) an antiferromagnetic insulator with in-plane magnetization for large U . The phase boundaries are shown in Fig. 5. Interestingly, similar results of an enhanced topological phase have been reported for the Kane-Mele model including long-ranged Coulomb interactions [63]. There, the Coulomb interaction induces charge-density-wave fluctuations, while our model shows static charge ordering through staggered potentials.

III.2. Ribbon
In order to analyze the robustness of the topological phases presented in the last section and to investigate the bulk-boundary correspondence, we calculate directly the edge properties on a zigzag ribbon of finite width.
We first consider the inversion-symmetric case, λ ν = 0.  Mean-field results have shown different magnetizations at the edge than in the middle of the ribbon [31]. This agrees with our results, and an example of the structure of the Weiss fields across the ribbon profile is shown in the inset of Fig. 7. The larger field at the edges decays quickly to the bulk value. The optimized values of both edge and midpoint antiferromagnetic fields as a function of U are shown in Fig. 7 for λ SO = 0.1. At the edges, any finite U results in a finite antiferromagnetic field. Sites that are not at the edges have a Weiss field comparable to the bulk values. Just below the bulk magnetic transition at U ≈ 3.8 they become finite, though small, which is a finite-size effect caused by the increasing correlation length as the magnetic transition is approached. The main consequence of the non-vanishing Weiss field is that the finite magnetization at the edges breaks time-reversal symmetry for any U and gaps therefore the edge states.
As the interaction is below the critical value for the bulk magnetic transition, topological analysis of the bulk suggests a topological insulator with gapless edge states, but a local symmetry breaking at the edges causes the edge states to gap. This local effect, namely that local timereversal symmetry breaking by a magnetic field causes states to gap, cannot be captured within a topological invariant of the two-dimensional (2D) system. However, at what point in the phase diagram this local symmetry breaking occurs, depends both on the specific model and also on the edge geometry. For example, for the armchair ribbon there is a region at small U with vanishing edge magnetization and therefore gapless edge states, even in the inversion-symmetric case λ ν = 0. In the last paragraph it is demonstrated that gapless edge states are impossible on a zigzag ribbon for any finite U , as long as λ ν = 0. This picture changes if inversion symmetry is broken. From the bulk calculations we know that λ ν tends to suppress magnetic ordering, where it increases the critical value of interaction U c for the magnetic transition (Fig. 5). The same principle is observed looking at the edge magnetization as a function of λ ν . For given U and λ SO , the Weiss field at the edges changes only marginally as λ ν is increased, and the edge is magnetic. However, at a critical value λ c ν , the magnetic moment drops to 0 in a first-order phase transition. This critical value λ c ν strongly depends on U . For λ SO = 0.1, for example, we get λ c ν = 0.006 as U = 1, and it raises by an order of magnitude to λ c ν = 0.07 for U = 2 and to λ c ν = 0.35 for U = 3. This argument can of course be turned around. Fixing the sublattice potential λ ν and varying the interaction strength U , one finds a critical value U c for the magnetic transition with finite magnetic moment only for U > U c . This critical value U c raises continuously with increasing sublattice potential λ ν , starting from U c = 0 at λ ν = 0.
Exemplary spectral functions are shown in Fig. 8, where we use spin-orbit coupling strength λ SO = 0.1 and interaction strength U = 2.5. If the sublattice potential λ ν is below the critical value, as in the top panel of Fig. 8, the edge is magnetic and the edge states are gapped. For λ ν > λ c ν there is no magnetization at the edge, and gapless states occur. We want to stress again that gapless edge states do not occur at any finite U in the inversion-symmetric case. To sum up, an inversionsymmetry-breaking term can stabilize the gapless edge state.

IV. CONCLUSION AND DISCUSSION
We have investigated the topological properties of the Kane-Mele-Hubbard model, comparing cases with and without inversion symmetry. For the calculation of the topological invariants we apply a combination of the topological Hamiltonian approach and the Wannier charge center method. This approach allowed to calculate the phase diagram of the KMH model in the U -λ ν plane. The inversion-symmetry-breaking term λ ν has a two-fold effect. First, for large values the topological order is destroyed and a trivial insulator obtained. Second, in combination with interactions the topological order is enhanced, pushing the phase boundaries towards the antiferromagnetic insulator to larger critical values of U .
This effect can also be seen in the surface properties of the honeycomb lattice. In agreement with previous studies, our calculations on the zigzag ribbon geometry have shown that with inversion symmetry any finite value of U results in a finite edge magnetization, which in turn produces a finite gap in the edge states. Introducing an inversion-symmetry-breaking field, this critical value U c is shifted to finite values, below which the whole ribbon including the edge is nonmagnetic, and a gapless surface state exists. As a result, one can find gapless edge states on the zigzag ribbon only when inversion symmetry is lifted and the interaction strength U is small enough, such that no ordered magnetic moments can form.
Our study is based on the Kane-Mele Hamiltonian, which was introduced as the low-energy Hamiltonian for graphene. Since the bulk gap in graphene is minute, the effects that we propose here are difficult to see in this material. However, there is increasing interest in artificial honeycomb systems using heavy atoms, such as bismuthene on SiC substrate [64]. Since these systems are grown artificially, it might be possible to modify their structure such that inversion symmetry is broken and the influence of this symmetry breaking on the topological properties can be studied. ACKNOWLEDGMENTS R.T. thanks Georg W. Winkler for helpful discussion. We acknowledge financial support from the Austrian Science Fund FWF, START program Y746.
Appendix A: Discussion -Comparison to mean-field As mentioned in Sec. III.1, the basic structures of the topological Hamiltonian could also be found in a mean-field approximation since the self-energy is diagonal. Usually, the z axis is chosen as the axis of mean-field decomposition [30]. The resulting matrix is then qualitatively different from the topological Hamiltonian of the DIA, since the mean-field magnetic moment points in the z direction. In order to respect that the easy axis is inplane, we did a mean-field decoupling in the x direction n i↑ n i↓ ≈ ( n i← n i→ + n i→ n i← − n i← n i→ ) , (A1) where | → ← = 1/ √ 2 (|↑ ± |↓ ). Within this framework, the same phases as in the DIA appear, where the meanfield one-electron Bloch Hamiltonian corresponds to the topological Hamiltonian. The phase boundaries, however, will shift since a bare mean-field approach does not capture quantum dynamics as the DIA.
In case of the Hubbard model on a honeycomb lattice λ SO = λ ν = 0, the magnetization direction is not important since SU(2) symmetry is not broken. The mean-field critical interaction for any quantization axis is U c = 2.23 [30,62]. If λ SO = 0, the difference between the two mean-field methods is important. Since the in-plane magnetic moment is always favorable, a restriction of the magnetization direction to be out-of-plane requires stronger interactions for the stability of the antiferromagnetic solution. This is the case in a conventional mean-field theory [30,43], hence, U c is overestimated in comparison with an in-plane mean-field approach (A1). Consequently, the slope of the U c -λ SO phase boundary is higher if z is used as a quantization axis.
In addition to the magnetic transition considered so far, using Wannier charge centers as an analytical tool allows again to extract topological information. The DIA results are described in the previous sections, showing the phase diagram of three different phases in figure 5. As mentioned above, the mean-field decoupling in the x direction gives qualitatively the same phases since the MF Bloch Hamiltonian has the same structure as the DIA topological Hamiltonian, but underestimates U c . New phases appear, however, in the standard Hartree-Fock approach where the z axis is the quantization direction. The Hamiltonian splits into spin up and spin down parts, which are decoupled if neither Rashba coupling nor inplane magnetization are present. Hence, even though time-reversal symmetry is broken in the presence of an antiferromagnetic moment, a Z 2 invariant can be defined using the spin Chern number ν S = C S mod 2, C S = (C ↑ −C ↓ )/2 as introduced by Sheng et al. [45]. The Chern numbers of the two spin categories are determined with the Wannier charge centers: Because of the conservation of S z , the two WCC can be labeled by their spin. The Chern number C S is then given by the difference of the WCCsx ↑ andx ↓ as they evolve continuously from 0 to 2π.
In the inversion-symmetric case, the only mean-field parameter that has to be determined self-consistently is the antiferromagnetic moment M AF = n A↑ − n B↑ = n B↓ − n A↓ . A change of both Chern numbers C ↑ and C ↓ occurs when the gap closes at a critical moment M c AF = 12 √ 3/U , which follows from diagonalizing the mean-field Bloch Hamiltonian. Since M AF rises continuously from 0 as U is increased, magnetic and topological transition do not coincide, leading to an antiferromagnetic quantum spin Hall phase between the two transitions.
If additionally inversion symmetry is broken, both onsite energy and occupation of A and B sites are different. Together with the magnetic order, this leads to different M c AF for spin up and spin down electrons. If C ↑ = 0 and C ↓ = 1 or vice versa, the total Chern number C = C ↑ +C ↓ is nontrivial. Hence, for a certain parameter range, an antiferromagnetic Chern insulator is realized (see Fig.  9). Both Chern insulator and antiferromagnetic quantum spin Hall insulator have also been found recently for cases where the symmetry breaking is not due to an on-site potential, but due to a spin-dependent hopping [44]. These phases are stable since for certain parameter regions the out-of-plane magnetization is energetically favorable.
The topological properties of the Chern insulator are not bound to time-reversal symmetry but related to the spin structure only. The number of edge states is directly determined by the Chern numbers of spin up and spin down electrons. As an example, the bands of a zigzag ribbon in the Chern insulator phase with only one edge state are shown in Fig. 9. Hence, bulk boundary correspondence is fully satisfied if the antiferromagnetic moment is in the z direction, but not if it is in-plane.