{'fraction_non_alphanumeric': 0.039986177617613666, 'fraction_numerical': 0.011749025028385249, 'mean_word_length': 4.252268602540835, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Quantum Spin Hall Insulators with Interactions and Lattice Anisotropy

We investigate the interplay between spin-orbit coupling and electron-electron interactions on the honeycomb lattice combining the cellular dynamical mean-field theory and its real space extension with analytical approaches. We provide a thorough analysis of the phase diagram and temperature effects at weak spin-orbit coupling. We systematically discuss the stability of the quantum spin Hall phase toward interactions and lattice anisotropy resulting in the plaquette-honeycomb model. We also show the evolution of the helical edge states characteristic of quantum spin Hall insulators as a function of Hubbard interaction and anisotropy. At very weak spin-orbit coupling and intermediate electron-electron interactions, we substantiate the existence of a quantum spin liquid phase.


I. INTRODUCTION
Time-reversal invariant topological insulators 1 -bulk insulators with metallic surfaces -are characterized by a Z 2 invariant 2,3 and cannot be adiabatically connected to trivial band insulator phases unless the single particle gap closes. While Z 2 topological insulators (TIs) are robust against disorder 4 , rigorous and general results about the fate of TIs in the presence of prominent electron-electron interactions are limited 5 . Strongly correlated TIs as well as exotic time-reversal invariant Mott insulator phases have been predicted [6][7][8][9][10][11][12] apart from more conventional magnetically ordered phases.
By analogy to the quantum Hall effect twodimensional TIs are also named quantum spin Hall (QSH) insulators. They were originally proposed to be realized in graphene 2 and later also in HgTe/CdTe quantum wells 13 where subsequent experiments 14 measuring a quantized conductance established the field of TIs. They possess an odd number of pairs of time-reversal conjugate counter-propagating edge states (the helical edge states) 2, [13][14][15] . There are other promising proposals to stabilize the QSH effect in real materials: graphene endowed with heavy adatoms like indium and thallium 16 and synthesized silicene 17 was shown to exhibit a stable QSH phase. Particularly interesting are monolayers or thin films of Iridium-based materials X 2 IrO 3 (X=Na or Li) which have been debated to possibly host QSH phase 18,19 since both spin orbit coupling (SOC) and electron-electron interactions are quite strong in such materials. All these systems have in common the underlying honeycomb lattice where recently a gapped quantum spin liquid for intermediate interactions was found using quantum Monte Carlo (QMC) [20][21][22] .
Very recent progress within ultracold atoms in tunable optical lattices establishes a second class of systems (beside Iridium based materials) where topological interacting phases might be realized. This progress is two-fold: (i) tunable hexagonal lattices have been realized 23 and (ii) different types of spin-orbit interactions are feasible now 24 . Additional electron-electron onsite interactions are considered as a standard tool in optical lattices 25 . All these achievements in such a rapidly evolving field indicate the demonstration of topological interating phases within cold atoms in the very near future.
In this paper, we investigate the interplay between SOC and electron-electron interactions on the honeycomb lattice and combine two very paradigmatic models: to capture the non-trivial band topology we consider the Kane-Mele (KM) model 2 without Rashba term and to describe interaction effects the Hubbard model, merging to the Kane-Mele-Hubbard model (KMH) 8,21,22,[26][27][28] . Our goal is to combine the cellular dynamical mean-field theory (CDMFT) 29   filling including temperature effects, and thoroughly address the fate of helical edge states as a function of interactions and SOC coupling. Additionally, we introduce the plaquette-honeycomb model which illustrates that the QSH phase is also stable toward lattice anisotropy.

II. KANE-MELE-HUBBARD MODEL
The KMH Hamiltonian on the honeycomb lattice reads (1) where i, j, label the sites on the honeycomb lattice, c iσ is the electron annhilation operator, n iσ = c † iσ c iσ , t is the hopping integral (t is our reference energy scale, and hence we set t ≡ 1), λ the SOC, U the onsite interaction, and ν ij = sgn[(d i ×d j ) z ] whered i/j are the two vectors along the links from j to i (ν ij = ±1) 2 . While the Hubbard model respects SU(2) spin and C 6 lattice symmetry, the SOC breaks the spin symmetry down to U(1) and the lattice to C 3 while it leaves time reversal symmetry invariant. Details of the used CDMFT method and its real space extension are provided in Appendix A.
For large interactions, the phase is an easy-plane antiferromagnet 8 . First, it is useful to start with a mean-field (MF) consideration of the TBI-SDW transition which can also be described within slaverotor theory 8 where the condensation of magnetic monopoles results in an XY -instability. At λ = 0, Eq. (1) is spin rotationally invariant and a staggered magnetization pointing in any direction will be a good order parameter. For finite λ, one can learn from the where the operators a kσ and b kσ are associated with the two sublattices. Eventually one obtains the MF spectrum ε MF ± = ± |g(k)| 2 + γ(k) 2 + U 2 |m| 2 + const. and from there we find the MF equation. Note that g(k) is the nearest neighbor hopping and γ(k) the second neighbor (spin-orbit) Haldane term 36 . The transition line differs slightly from the analog calculation with the order parameter S z i . Mainly, the form of ε MF ± reveals that the mean field approximation causes another mass term which does not compete with the SOC mass γ(k). Consequently, when passing from the TBI to the SDW phase we do not expect closing of the single particle gap. Indeed, within CDMFT no closing of the single particle gap is observed (while the gap has a local minimum at the transition). We also computed S x i which becomes finite at the TBI-SDW transition (see also App. A); note that this transition is of 3D XY universality at T = 0.
For very weak SOC we identify a small phase at intermediate U which is non-magnetic, exhibits a finite spin gap, and is separated from the TBI phase by closing of the single-particle gap (see Fig. 3). This phase is reminiscent of the recently found quantum spin liquid (SL) phase 20,21 , in particular, position and shape of this nonmagnetic insulator phase essentially coincides with the SL phase found within QMC 21 . Here, we shall mention that since the CDMFT method approximates the selfenergy by restricting it to the chosen cluster (see inset of to observe the SL phase, see Fig. 2. We have plotted single particle gap ∆ sp and magnetization m for λ = 0 and various cluster sizes C. For C = 2 and C = 4 single particle gap and magnetization become finite simultaneously while for C ≥ 6 the onset of both quantities occurs for different U indicating the SL phase. The existence of such possible spin-gapped phases on the honeycomb lattice has also been addressed theoretically 37 . In the phase diagram of Fig. 1 there is an additional semi-metal (SM) phase where the Dirac cones of graphene are preserved. Using QMC as impurity solver, we systematically discuss temperature effects at weak SOC. We observe that the SL phase appears at very low temperatures (see Fig. 2) and that the SM domain prominently develops with temperature as the TBI gap is also very small when λ → 0. For T → 0, the SM phase reduces to a line pointing along the λ = 0 axis which survives until relatively large interactions. From Renormalization Group point of view, this SM line at λ = 0 is known to be stable for weak U 38 . Hence, we show the SM region for T = 0.025, 0.0125, and 0.005 in Fig. 1. Remember that both the SL and TBI phases at λ ≤ 0.03 possess a small single-particle gap. In the inset of Fig. 3 we present the single particle gap and magnetization for λ = 0.02 and T = 0.025. The evolution of the phase diagram with temperature is shown in Fig. 3 for λ = 0.02.
We are using continuous time quantum Monte Carlo (CTQMC) as an impurity solver which enables us to access finite temperatures. This is clearly an advantage of the method since experiments are performed at finite temperatures but a careful interpretation of results is still needed. In Fig. 3, for instance, the SDW phase at finite T violating Mermin-Wagner theorem is clearly an artifact of the CDMFT method. In contrast, the stability of the TBI phase with increasing temperature is reliable and important for experimentalists. From an experimental perspective, it may be interesting to know how T c can be raised and for what parameter settings T c is conveniently large. Since the spin orbit gap for λ = 0.02 is tiny, T c < 0.05. Increasing λ (and the bulk gap), however, leads to an (approximately linear) increase of T c up to T c > 0. 16  line is a linear fit. We corroborate that the critical temperature T c follows the gap at T = 0. For finite values of U a similar behavior of T c is obtained. Note that T c does not necessarily indicate a phase transition from TBI to metal; to clarify the precise nature of this boundary would require further work.
Within CDMFT, we also access the double occupancy D occ = n i↑ n i↓ = ∂F /∂U where F is the free energy. In Fig. 5 we show that D occ which measures the itinerant properties of electrons is an insightful quantity to detect the magnetically ordered phase 31,32 . The double occupancy is bounded as 1/4 ≥ D occ ≥ 0. For U = 0, the ground state is a TBI with occupied lower and empty upper bands, |TBI = k l † k↑ l † k↓ | 0 , where (a kσ , b kσ ) T = T σ (l kσ , u kσ ) T and the matrices T σ are given by with Here g(k) is the nearest neighbor hopping term, γ(k) is the Haldane term (second neighbor spin orbit hopping), and ε(k) = |g(k)| 2 + γ(k) 2 is the single-particle Kane-Mele energy spectrum. Essentially, the matrices T σ contain the eigenvectors of the single-particle eigenstates. The Fermi level lies in the gap, ε F = 0. Choosing a site i that belongs to the sublattice "a" we obtain, TBI| n i↑ n i↓ |TBI = where we used in the second last line l † k1↑ l k2↑ l † k3↓ l k4↓ = δ k1k2 δ k3k4 as well as u † kσ u kσ = 0, u † kσ l kσ = 0 and l † kσ l kσ = 1. N refers to the number of unit cells. Using the relations we confirm that independent of λ in agreement with CDMFT, see Fig. 5. Above we considered without loss of generality n iσ = a † iσ a iσ . Identical calculations performed with the b iσ operators which belong to the other sublattice yields of course the same result. In the opposite limit, U → ∞, we clearly find D occ = 0 as a fingerprint of Mott physics since we impose half filling.

III. HELICAL EDGE STATES
To describe the edge states associated with QSH insulators, first we apply the concept of the helical Luttinger liquid (HLL) 15 .
Hence, we linearize the spectrum around the Fermi points and switch from lattice operators to field operators ψ R↑ (x) and ψ L↓ (x) which are right-and left-moving fields, respectively; we obtain for the non-interacting part Note that a standard single-particle (disorder) backscattering term ψ † R↑ ψ L↓ + h.c., which opens up a mass gap in the spinless Luttinger liquid, is not allowed since the model is odd under time-reversal symmetry, T 2 = −1. Only two time-reversal invariant interactions are allowed: the forward scattering as well as the umklapp scattering ∼ ψ † R↑ ψ † R↑ ψ L↓ ψ L↓ which is not intrinsically present in the Hubbard model. Instead we shall include the forward interaction H I = U dx ψ † R↑ ψ R↑ ψ † L↓ ψ L↓ . As long as there is no magnetic order in the bulk, (H 0 + H I ) can be solved exactly resorting to bosonization which results in power-law decaying spin and charge correlations in the HLL. This result is also obtained through a spin wave analysis at the edges for weak interactions 39 . It is worth mentioning that, in contrast, spin-spin correlation functions in the bulk decay very rapidly 8  ing the bosonization procedure and introducing the Luttinger parameter K as usual, the Hamiltonian becomes: (7) Here, v is the renormalized plasmon velocity and a the lattice spacing, the field φ contains both spin and charge degrees of freedom and ∂ x θ and φ are conjugate variables 15 . The Sine-Gordon term is a relevant perturbation for repulsive interactions (since K 2) and hence, through the pinning of the φ-field, the edge modes acquire a charge gap at the SDW transition. We can also check that the spin sector at the edges is described by an Ising order characterized by S x i = 0 and S y i = 0. This shows that the SDW transition affects the charge sector of the HLL in a non-trivial way resulting in a metalinsulator transition of Kosterlitz-Thouless type. Note that the disordering-ordering transition in the spin sector also influences the charge properties of the edges.
In addition, we performed CDMFT simulations on honeycomb cylinders with length x = (3a/4)L with both armchair and zigzag edges (L = number of sites in xdirection). By computing the spectral function A(k y , ω) we extract the edge state spectrum in the presence of interactions. In contrast to previous work 21 we treat the full microscopic Hamiltonian. We observe the following: (i) The plasmon velocity v associated with the edge modes slightly decreases for increasing U as expected from the HLL. (ii) The intensity of the spectral func- We further performed computations for the spectral function on zigzag-ribbons. Now, in contrast to the armchair-case, the periodic boundary conditions are imposed in x-direction yielding k x as the good quantum number. For zigzag edges one extracts the more familiar spectrum 2 of the KM model at small λ (see Fig. 7, top): the edge states connect the upper band of a gapped Dirac cone K with the lower band of the other gapped Dirac cone K and vice versa. In addition, we have also shown the situation for stronger SOC λ = 0.2 (see Fig. 7, bottom) where the bulk spectrum is relatively flat. Both spectral functions are computed for U = 2.0. The main findings (i) -(iii) obtained for armchair-ribbons also apply for zigzag-ribbons.
We also computed the spectral function on a cylinder for parameters U and λ which belong to the spin liquid phase (see Fig. 6). While the single particle gap is very small we did not find edge states inside the gap. While it is not known whether the true SL exhibits edge states, the non-magnetic insulator phase found within CDMFT is topologically trivial and absence of edge states expected.

IV. PLAQUETTE HONEYCOMB MODEL
Next, we show that the QSH phase is not only stable towards relatively large interactions but also to lattice anisotropy resulting in the plaquette-honeycomb model, as sketched in the inset of Fig. 8. Hopping amplitudes which connect yellow plaquettes are t = αt while hopping amplitudes inside a cluster are t. Spin-orbit couplings λ and λ = αλ are analogously defined. In the absence of spin-orbit couplings, note that this type of anisotropy immediately opens a gap at the Dirac points. Mostly when α → 0, we check that the spectrum exhibits flat bands associated with localized states on a given hexagon, hence resulting in the plaquette insulator (PI) in Fig. 8. Also for large α, the system is an insulator (I) at the Fermi level. For finite SOC λ, we find that the QSH phase is stable over a large region of the phase diagram (containing the KM model at α = 1). The gap closes at the phase transition between the PI/I and the QSH phase yielding a SM phase boundary.
These SM lines are stable in the presence of finite electron-electron interactions U up to a critical U c (α, λ) what is of comparable size as U c (1, λ) (i.e., for the isotropic KMH model). Preliminary results for some parameters indicate that for U > U c (α = 1, λ) there is a non-magnetic insulator phase like in the α = 1 case. The detailed investigation of the interacting phase diagram of this plaquette honeycomb model is left for future studies.

V. CONCLUSION
To summarize, through a combination of CDMFT and complementary approaches, we have substantiated that the QSH phase is robust toward prominent electronelectron interactions and lattice anisotropies. We have provided a quantitative analysis of the edge state properties which shows that the ordering transition for the spins is also accompanied by a metal-insulator transition at the edges. For very weak SOC, we have also confirmed the existence of a spin-gapped insulating phase (SL). Finally, we have introduced an anisotropic version of the KM model with the potential to host a rich phase diagram when including the Hubbard interaction.

Acknowledgments
We acknowledge discussions with F. F. Assaad, G. Fiete, M. Hohenadler, T. Lang, and G. C. Liu In order to investigate edge states within CDMFT we use the real-space extension of the homogeneous CDMFT method 29,30 . We consider (nano-) ribbons (i.e., cylinder geometry) with translation symmetry in direction along the edges (which we aligned for armchairribbons along the y-axis).
Within CDMFT we map the original honeycomb lattice onto a N c -site effective cluster embedded in a selfconsistent medium. The effective cluster model is obtained via an iterative procedure which can be started with an initial guess of the cluster self-energy Σ(iω). The effective medium represented by the dynamical meanfield, which is also known as Weiss function g c (iω), is determined by the cluster self-energy Σ(iω) via the coarsegrained Dyson equation 31 . Since the translation symmetry in the x-direction is broken, we shall use the realspace CDMFT coarse-grained Dyson equation for the inhomogeneous system: (A1) Both g c and Σ are block diagonal matrices with block size N c and dimension ((N c ×N x )×(N c ×N x )). I, J are block indices, N c is the cluster size, and N x is the number of the clusters after performing the dimensional reduction, i.e., N x corresponds to the height of the ribbon.
As far as g c (iω) is determined, the impurity solver 40 can be used to compute the effective many-body impurity problem to obtain the full Green's function G c .
Eventually, by using Dyson equation of the cluster system Σ(iω) = g −1 c (iω) − G −1 c (iω), we recalculate the selfenergy Σ(iω) to finish the iterative loop. The CDMFT self-consistent iterative loop should be repeatedly carried out, until the numerical convergence has been reached.
Note that we use the spinor notation, Ψ † = {c † i↑ , c † i↓ } (i ∈ cluster) and Ψ = (Ψ † ) † since a spontaneous transverse magnetization mixes the ↑ and ↓ spin-parts of the system. Therefore the Green's function in the CDMFT framework takes the following form, From the above procedure we can easily obtain the quantites D occ 31,32 and S x . The double occupancy D occ = n i↑ n i↓ = ∂F/∂U , as the first derivative of the free energy, is an important quantity in studies of strongly correlated systems. As D occ is directly related to the potential energy this is an indicator for the transition order. In the following, we suppress the site index i of the electron operators. Using the notation of Eq. (A2) we find where 0 < τ 1 < τ 2 < τ 3 , and we used Wick's theorem. The transverse magnetic order parameter S x i reads S x i = 1 2 c † i↑ c i↓ + h.c. which can be expressed through the above defined Green's functions, In Fig. 9 we show both S x = S x i and S z = S z i vs. U for λ = 0.1. Since the SOC breaks the spin-rotational invariance the magnetization is inside the XY -plane; consequently we find a finite value of S x at the transition into the SDW phase while S z remains zero.
In the case of armchair-ribbons and PBC calculations we used a six-site cluster, for zigzag-ribbons an eightsite cluster, see inset of Fig. 1. In addition, we performed some PBC calculations with the eight-site cluster as well as with a 24-site cluster (see inset of Fig. 8). We found good quantitative agreement between the results performed with the different clusters which underlines the reliability of the CDMFT method. But we would like to emphasize that a six-site cluster is needed to observe the spin liquid phase. For smaller clusters sizes (2 and 4) the spin liquid phase cannot be found as shown in Fig. 2.