Topological minibands and interaction driven quantum anomalous Hall state in topological insulator based moiré heterostructures

The presence of topological flat minibands in moiré materials provides an opportunity to explore the interplay between topology and correlation. In this work, we study moiré minibands in topological insulator films with two hybridized surface states under a moiré superlattice potential created by two-dimensional insulating materials. We show the lowest conduction (highest valence) Kramers’ pair of minibands can be \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathbb{Z}}}_{2}$$\end{document}Z2 non-trivial when the minima (maxima) of moiré potential approximately form a hexagonal lattice with six-fold rotation symmetry. Coulomb interaction can drive the non-trivial Kramers’ minibands into the quantum anomalous Hall state when they are half-filled, which is further stabilized by applying external gate voltages to break inversion. We propose the monolayer Sb2 on top of Sb2Te3 films as a candidate based on first principles calculations. Our work demonstrates the topological insulator based moiré heterostructure as a potential platform for studying interacting topological phases.

In this section, the topology of inversion symmetric moiré system in Fig. 2(a) of the main text with ϕ = 0, α = 1, V 0 /E 0 = 0 is studied under the perturbation of moiré potential strength ∆ 1 .
By the Fu-Kane parity criterion [1], the Z 2 invariant ν can be determined by the parities λ i at time reversal (T ) invariant momenta one Γ and three M of moiré Brillouin Zone (MBZ) for one of the degenerate states by (−1) ν = i λ i At ϕ = 0, α = 1, V 0 /E 0 = 0, the crystal symmetry of this system is described by the point group D 6h with six-fold rotation C 6z about the z-axis, the inversion I, the y-directional mirror M y and the z-directional mirror M z .We label the original basis of our model Hamiltonian Eq. 1 in main text by |k, J z , α⟩, where k is the momentum, J z = ±1/2 labels two spin states of surface states and α = t, b labels the top and bottom surfaces.The z-directional mirror M z transforms the top surface to the bottom surface and thus it relates the basis wave-functions on two surfaces by M z |k, J z , t⟩ = e −iπJz |k, J z , b⟩ M z |k, J z , b⟩ = e −iπJz |k, J z , t⟩. ( We may transform the basis wave-functions to the bonding and anti-bonding states of two surface states as with I = ± labels the transformation property under the inversion parity and the eigen-values of the M z operator M z |k, J z , ±⟩ = ±e −iπJz |k, J z , ±⟩.
On these bonding and anti-bonding basis the Hamiltonian H 0 in the main text Eq.1 can be written in a block diagonal form, with where k± = −i(∂ x ±i∂ y ), m = m 0 +m 2 (−∂ 2 x −∂ 2 y ), I 4×4 is a 4×4 identity matrix, and m z = ±i labels the eigen-values of the mirror operator M z .H M is the moiré potential with ϕ = 0, α = 1, V 0 /E 0 = 0.
We next determine the parities of lower-energy minibands at T invariant momenta, including one Γ and three M in the moiré BZ, of the Hamiltonian H 0 in the limit |∆ 1 | ≪ |vb M  1 |, |m| via perturbation theory.As the H 0 is block diagonal in the m z = ±i subspace, we may perform the perturbation calculation for the m z = −i block while the miniband parity of the m z = +i block can be related by TR symmetry.For the m z = −i block, the unperturbed Hamiltonian is H TI mz=−i while H M is treated as the perturbation.We choose the eigen-wavefunctions of H TI mz=−i to possess a well-defined gauge at Γ in the moiré BZ, which can be written as for m > 0 with the eigen-energies for m < 0 with the eigen-energies , where cos θ k = m/ √ m 2 + v 2 k 2 and ke iϕ k = k x + ik y .
I|ψ TI I,−i (k)⟩ = I|ψ TI I,−i (−k)⟩, (10) and the expression for the eigen-energy can be unified as so the inversion parity I also labels different eigen-energies of our model Hamiltonian.The second lower-index −i in the eigen-state labels the M z eigen-values.This definition of the eigen-states |ψ TI I,mz ⟩ is also used in the main text.From the expression of the eigen-energies (11), the higher energy state with E TI CB (k) = E TI I=+sgn(m) (k) that corresponds to the conduction bands should be given by |ψ TI  CB,−i (k)⟩ = |ψ TI +sgn(m),−i (k)⟩ (12) and the lower energy state with E TI VB (k) = E TI I=−sgn(m) (k) for the valence bands should be For m z = +i subspace, we use TR symmetry operator, given by in the basis Eq. ( 5) with K for complex conjugate, to define |ψ TI I,+i (k)⟩ = −iT |ψ TI I,−i (−k)⟩, (15) and the commutation relation [T , I] = 0 leads to the same inversion parity for two degenerate states |ψ TI I,mz=±i (Γ i )⟩ at any TR-invariant momentum Γ i .
The band gap and the inversion parity of minibands at Γ of the moiré BZ are determined by the hybridization term m in the limit |∆ 1 | ≪ |m|, for which the moiré potential does not play a role.Thus, we only need to consider the unperturbed Hamiltonian H TI mz=−i in Eq. (7), which is diagonal, and the eigen-state |ψ TI +,−i (Γ)⟩ = (1, 0) T has the eigen-energy m and |ψ TI −,−i (Γ)⟩ = (0, 1) T has the eigen-energy −m.The lower index I directly gives the parity of the eigen-state at Γ, namely I|ψ TI I,−i (Γ)⟩ = I|ψ TI I,−i (Γ)⟩.The parity of the CB1 for the eigen-state |ψ TI CB,−i (Γ)⟩ is λ Γ = +sgn(m) and that of the VB1 for the eigen-state |ψ TI VB,−i (Γ)⟩ is λ Γ = −sgn(m), depending on the sign of m.Therefore, the parities of CB1 and VB1 at Γ are opposite, Different from the Γ point, the moiré potential is essential in determining the parities of the minibands at M in the moiré BZ.Without moiré potential, the eigen-states so the spectrum is gapless at M, even with a finite m.The moiré potential will couple these two states at M and −M as both belong to the same momentum in moiré BZ.By projecting the full Hamiltonian H 0 (r) into the subspace spanned by these two states |ψ TI CB,−i (±M)⟩, we find the effective Hamiltonian H CB ef f (M) for CB1 and CB2 is given through the degenerate perturbation by where with the parity λ CB1 M = (+sgn(m))(−sgn(∆ 1 )).For VB1 and VB2, the effective Hamiltonian at M is given by where with the parity λ VB1 M = (−sgn(m))(+sgn(∆ 1 )).Thus, The parity at M for CB1 and VB1 are the same.Because the Z 2 invariant is (−1) ν = λ Γ (λ M ) 3 , (−1) νCB1 = −sgn(∆ 1 ) and (−1) νVB1 = +sgn(∆ 1 ), so that ν CB1 and ν VB1 are differed by 1.Thus, we conclude ν CB1 + ν V B1 = 1 mod 2, namely one of CB1 and VB1 has nonzero Z 2 invariant and the other has trivial Z 2 invariant.The above conclusion of topology of CB1 and VB1 can also be understood from the chiral symmetry operator C of H TI , defined by C = τ z s z , when the chemical potential is at the charge neutrality point, where τ acts on the top/bottom surface degrees of freedom and s acts on spin.The emergence of the chiral symmetry requires dropping higher-order k terms, e.g.k 2 terms, in H TI , which are not important at the moiré energy scale.This operator has the commutation relations On the basis of Eq. ( 5), the form of chiral symmetry operator is transformed into which mixes the eigen-states with opposite M z eigen-values, namely This implies At Γ, the opposite parities between |ψ CB,−i (Γ)⟩ and |ψ VB,+i (Γ)⟩ (+sgn(m) for |ψ CB,−i (Γ)⟩ and −sgn(m) for |ψ VB,+i (Γ)⟩ directly come from the anti-commutation relation {C, I} = 0.
At M , the CB1 (VB1) and CB2 (VB2) are degenerate for H TI , so we need to consider the first order perturbation from H M .For the convenience of the discussion, we introduce the inversion adapted basis functions for CB1, CB2, VB1, and VB2 as with the parity They are related by chiral symmetry As I, H M = 0, the first order perturbation correction from H M is diagonal.For CB1 and CB2 |Ψ CB,I,mz (M)⟩, we find the perturbation Hamiltonian is with cos θ kM = m/ m 2 + v 2 k 2 M , while for VB1 and VB2 |Ψ VB,I,mz (M)⟩, the perturbation Hamiltonian is The eigen-energy of the system at M after taking into first order perturbation is which is equivalently So, E CB,+ (M) > E  21) for all cases.In the above analysis, the key is that H M commutes with C and leads to the same parity at M, different from the case at Γ where H TI anti-commutes with C and results in opposite parities.This leads to one of CB1 and VB1 to be topologically non-trivial while the other to be trivial.
Additionally, from the perturbation perspective, the bandwidth of CB1 can be reduced by increasing m, ∆ 1 , and moiré lengths as shown in Fig. 2. Increasing m and ∆ 1 opens larger gaps at Γ and M and thus flattens CB1.A larger moiré unit cell can reduce b M 1 , as well as E 0 , and thus increases both m/E 0 and ∆ 1 /E 0 , leading to the narrowing of the bandwidth of CB1.

B. Topological phase transition when varying ϕ
In this section, we study the topological phase transition of our system when varying ϕ of the moiré potentials.For general ϕ, there is no inversion symmetry.As shown in Fig. 3(a), ν CB1 changes from 1 to 0 when ϕ varies from 0 to 1/6.Between the two phases, there is a gap closing around ϕ ≈ 1/48 at K and K ′ .The gap closing at K, K ′ happens between two states |u Jz=−1/2 (K)⟩ and |u Jz=3/2 (K)⟩, belonging to K6 and K4 irreducible representations as summarized in Tab.I [2,3], respectively, with different angular momenta J z = −1/2 and J z = 3/2 under threefold rotation C 3 , where |u Jz (K)⟩ are eigen-states of H 0 (K) as Eq. ( 44).The effective Hamiltonian on the basis From the phase diagram in Fig. 2(a)(b) of the main text, we notice that the Z 2 topological property of the system shows a periodicity when ϕ varies by 1/3.Indeed, one can show that the moiré potential ∆(r) with ϕ and ∆ ′ (r) with ϕ + 1/3 (with the same ∆ 1 parameter) are related by a constant shift as As a constant shift of potential term cannot change the band topology of the system, ν must keep the same for ϕ and ϕ + 1/3 while keeping other parameters.For NI phase, the Wyckoff position of Wannier orbitals should also shift accordingly by a M 1 /3 + 2a M 2 /3, as shown in Fig. 8. Similar topological phase transitions happen for α = 0 and V 0 /E 0 = 1.2 by a Dirac-type gap closing at K and K ′ between two states with different angular momenta when ϕ varies from 0 to 1/6 to 1/3, as shown by Fig. 3(b).

C. Atomic limits at m2 → ∞
In this section, we will provide theoretical understanding of the non-trivial morié minibands from the atomic limits of the CB1 and CB2 with a large m 2 term (the quadratic term of the inter-surface coupling m), and discuss how the realistic models with a small m 2 are connected to this atomic limit.For α = 1, V 0 /E 0 = 0, ϕ = 0 with the inversion symmetry in Fig. 2(d) of the main text, the energy spectra for increasing m 2 are shown in Fig. 5(a)-(c).We focus on CB1 and CB2 as a whole for atomic limits because they together have ν CB1 + ν CB2 = 0 and are topologically trivial.When increasing m 2 , we do not find any gap closing between CB1, CB2 and other valence bands or higher conduction bands.Thus, the topological properties of CB1 and CB2 remain the same, and the CB1 and CB2 are adiabatically connected to those corresponding bands in the large m 2 limit.When the m 2 term dominates in H 0 , for ϕ = 0, we may consider the Hamiltonian in the m z = ±i basis, Eq. ( 6), and drop the linear term ±iv k± in the off-diagonal component first.Then, the remaining part of the Hamiltonian just describes the 2D electron gas (2DEG) with a simple parabolic dispersion on a hexagonal potential,

C6
My T Supplementary Table II.(a)(b) symmetry operators in the irreducible representation at high symmetry momenta Γ and K for the double space group 183 P 6mm corresponding to the point group C6v.
orbitals localized at the moiré hexagonal potential minima of the Wyckoff positions 1b and 1c for the point group D 6h and give rise to a Dirac cone at K and K', similar to the case of graphene.The off-diagonal linear term in Eq. ( 6) represents the strong spin-orbit-coupling (SOC) of TI thin films, which gives rise to a small gap opening for the dispersion in Fig. 5(c) and can be treated perturbatively.We perform a k • p type of perturbation expansion of the full Hamiltonian H 0 (k) around K. The basis wave functions are chosen to be the eigen-states of H 0 (K) in Eq. 1 of the main eigenstates without SOC (v f = 0) for CB1 and CB2 with the irreps K6 for |u 1/2,+i (K)⟩, |u −1/2,−i (K)⟩ and K4 , K5 for |u 3/2,+i (K)⟩, |u 3/2,−i (K)⟩ (Fig. 5(c)), the detailed forms of which can be numerically evaluated.The relevant symmetry operators are with σ acts on the different m z , τ acts on different J z in one m z , and K is the complex conjugate.The SOC couples |ũ Jz,mz (K)⟩ and valence bands and contributes a k-independent term from the first order Löwdin perturbation [4] by The effective Hamiltonian H ef f around K to the first order in k with where C 0 , ∆ KM , v f are material dependent parameters and can be obtained numerically from the perturbation expansion.The above effective Hamiltonian H ef f (k) resembles the Kane-Mele model [5] with the SOC term ∆ KM σ z τ z , which provides another understanding of the non-trivial Z 2 topology of the CB1 and CB2 in our moiré system.For α = 0, V 0 /E 0 = 1.2, ϕ = 0, similar procedure can be applied to find the atomic limits of CB1 and CB2 at a large m 2 .The point group in this case is C 6v group.For m 2 = 1.5E 0 /|k M | 2 in Fig. 6(e), the effective Hamiltonian on the same basis as Eq. ( 37) is given by Besides Kane-Mele SOC term ∆ KM , there is another Rashba type of SOC term ∆ R (σ x τ y − σ y τ x ) as the inversion symmetry is broken for α = 0 [5].The Rashba term couples two basis functions |ũ 3/2,±i (K)⟩ ( K4 and K5 irreps) and opens the gap between these two states, as schematically shown in Fig. 6 6(e).In this limit, the topology of the CB1 and CB2 is ν CB1 + ν CB2 = 0, as the CB1 and CB2 together form an atomic limit.
With decreasing m 2 to m 2 = 0.1E 0 /|k M | 2 , we notice the nodes at K between CB1 and CB2 remains, but there is another band crossing between CB2 and higher conduction bands at Γ in Fig. 6(b).This band crossing at Γ changes the overall Z 2 topology of CB1 and CB2 to ν CB1 + ν CB2 = 1 for a smaller m 2 .In Fig. 7, we also show the band dispersion and the irreducible representations at high symmetry momenta for other higher-energy minibands (labelled by CB3, CB4, CB5 and CB6).We find the minibands of CB3, CB4 and CB5 are always touching each other and their total Another transition between CB5 and CB6 occurs at m 2 ≈ 0.035E 0 /|k M | 2 (See Fig. 7e), and after this transition, ν CB3 + ν CB4 + ν CB5 becomes zero while the other non-trivial Z 2 number is moved to even higher energy minibands.For m 2 < 0.1E 0 /|k M | 2 , these additional transitions only occur for higher-energy minibands, while the Z 2 topology of CB1 and CB2 remains the same (ν CB1 + ν CB2 = 1).For m 2 < 0.04E 0 /|k M | 2 , we find a gap between CB1 and CB2 opens at K due to the interchange between the K6 and K5 minibands.Thus, the isolated CB1 with ν CB1 = 1 and CB2 with ν CB2 = 0 states can be found in Fig. 6(a) for m 2 = 0.

D. Normal insulator phases of atomic limits
We construct the maximally localized Wannier functions [6] for the topologically trivial region for the CB1 as shown in Fig. 8.The locations of Wannier functions show the NI phase of CB1 has localized orbitals at Wyckoff positions 1b for ϕ = 1/6, 1a for ϕ = 1/2, 1c for ϕ = 5/6, as indicated in the phase diagram Fig. 2(a)(b) of the main text.Comparing the Wannier functions with the moiré potentials, they are located at minima of moiré potentials and correspond to the lowest conduction bands as expected.Since the minima of potentials change from one to another when tuning ϕ, the localized orbitals shift from one location to the other.The phase transition between two NI phases with orbitals at different Wyckoff positions has gap closing [7], shown as the semi-metal phase in Fig. 2(a)(b) of the main text, as they belong to different atomic limits.

Supplementary Note 2. HARTREE FOCK METHODS FOR COULOMB INTERACTION A. Eigenbasis projection
In this section, we project the Coulomb interaction into the eigenbasis of the non interacting Hamiltonian H 0 (k) [8,9].The non-interacting moiré Hamiltonian in the second quantization form is where u n G,α (k) satisfies the eigen equation for H 0 (k) with energies E n 0 (k).By replacing G with G + G 0 and G ′ with G ′ + G 0 in Eq. ( 44), we obtain which can be viewed as the eigen equations for u n G,α (k + G 0 ) by replacing k with k + G 0 in Eq. (44), Thus, we can fix the periodic gauge for the eigen-state as As u n G,α (k) is a set of orthonormal basis, we can take the inverse of the above expansion as and To improve the efficiency of the numerical calculations, we need to further fix the gauge freedom of eigenstates.An important step is to choose the real gauge for the Hamiltonian and eigenbasis due to the space-time inversion symmetry C 2z T in 2D for moiré potential with ϕ = 0. Take C 2z T = U CT K with K as complex conjugate.U CT is unitary and satisfies the U CT * U CT = 1 from (C 2z T ) 2 = 1.Under the basis transformation and the corresponding Hamiltonian and eigenbasis can be chosen to be real.There is still a SO(2) gauge freedom left for eigenstates for Fig. 2(d) in the main text with inversion and ± gauge freedom for Fig. 2(e) in the main text without inversion.
In the eigenbasis, the non-interacting Hamiltonian is The dual-gated Coulomb interaction potential is [8, 10] where S is the area, d is the dual-gate distance, ϵ 0 ϵ r are permittivity, e is electron charge.The Coulomb interaction Hamiltonian in second quantization form is k1,k2,q,G G1,α1,G2,α2 V (q + G) k1,k2,q,G m1,n1,m2,n2 with the form factor The form factor satisfies In the real eigenbasis, the form factors are all real.

B. Self-consistent Hartree-Fock mean field Theory
In this section, we treat the Coulomb interaction under the Hartree-Fock (HF) approximation [8].The basic idea is the decoupling of four-fermion operators by The expectation value of the two-fermion operator is the density matrix determined by with n F as the Fermi distribution function and ψ HF j,m (k),E HF j (k) as the j-th eigenstates and eigenenergies of Hartree-Fock Hamiltonian where H HF nm [ρ](k) is defined in Eq.Eq. ( 61) below.We always choose E HF j=1 (k) < E HF j=2 (k) < ...., so the mean field ground state is given by the eigen-state ψ HF j=1 (k).Here, we do not consider non-uniform order parameters in real space with the form c † m (k)c n (k + q) for q ̸ = 0.The Coulomb interaction under Hartree-Fock approximation is with the Hartree term with n as the number of bands projected.Since the H 0 (k) comes from DFT with Hartree-Fock interaction, the non-interacting states ψ HF j,m (k) = δ j,m or ρ 0 (k) would be a solution to the Hartree-Fock mean-field Hamiltonian.To achieve this, the H I [ρ 0 ] is subtracted from H I [ρ] [8,11].We define the Hartree-Fock Hamiltonian to be We solve H HF [ρ](k) self-consistently in the following standard procedures.We first choose an initial guess of the density matrix, denoted as ρ ini (k), as the order parameter for the filling of one band (half filling in two-band model and one quarter filling for four-band model).Based on ρ ini (k), we can construct H HF [ρ ini ](k) from Eq. ( 61) and calculate the corresponding new eigenstates that allow us to construct the new density matrix, denoted as ρ new (k).We reset ρ ini (k) = ρ new (k) and continue the iterative process until the convergence is achieved.The criterion for the convergence is taken as the spectra ẼHF j (k) of H HF [ρ ini ] and with max taken for all bands j in H HF and k on the high symmetry lines Γ−K −M as shown in Fig. 9(a).E 0 = v|b M 1 |.The final self-consistent solution for the density matrix is denoted as ρ HF (k) which is determined by the eigen wavefunctions |ψ HF j (k)⟩ by Eq. ( 57) and (58).The energy per particles for each self-consistent solution to the meanfield Hamiltonian is with N as the number of electrons for the filling.
C. Two-band model of CB1 In this section, we discuss the self-consistent solutions of H HF (k) at half filling of two-band model for CB1.Below we will discuss both the inversion-symmetric and asymmetric cases.
We first describe our gauge choice of the non-interacting eigen-states for the case α = 1, V 0 /E 0 = 0, ϕ = 0 with inversion symmetry, which is important to simplify the numerical calculations.The non-interacting states are |u CB1 ±i (k)⟩ with the m z eigenvalues ±i.The mirror Chern number [12] C = ±1 can be defined for |u CB1 ±i (k)⟩.The relative phase between two M z states are fixed by C 2z T through It turns out that the Hartree-Fock calculations can be simplified by taking the real gauge due to the C 2z T symmetry and thus we transform the basis wavefunctions into the real-gauge form where φ k is the remaining relative U (1) phase between eigen-states opposite m z (spin U (1) symmetry).The real eigen-states with different φ k can be related by a SO(2) transformation which shifts φ k to φ k + φk .The other symmetry operators can be taken as and which transforms Pauli matrices σ as We performed the self-consistent calculations on the basis |u R,CB1 ± (k)⟩ and generally consider the following two types of order parameters: Supplementary Table III.A summary of symmetries preserved (✓) or broken (×) by the mirror polarized states with ρ R,HF y (k) and the mirror coherent states with ρ R,HF zx (k).ρ R,HF (k) are the self-consistent solutions from the mean-field Hamiltonian H HF (k).The symmetry operators are written in two basis.τ, s are Pauli matrices for the surface and spin basis as Eq.1 in the main text.σ are the Pauli matrices for the real basis |u R,CB1 Different symmetry properties of ρ R y (k) and ρ R zx (k) under the C 2z T and SO(2) symmetry guarantee that they will not mix with each other.We may start from the initial density matrix ρ R ini (k) = ρ R y (k) with certain forms of f 0y and f y , which preserves the R( φk ) symmetry, [ρ R ini (k), R( φk )] = 0.As the Hartree-Fock Hamiltonian ] = 0 for any φ k .From Eq. (66) of R( φk ), the Hamiltonian has to take the form where h 0 (k), h y (k) are some functions of k which can be determined numerically.From the above form of the Hamiltonian, the new density matrix can be evaluated as is the Fermi distribution function.Thus, the R( φk ) symmetry is preserved in the self-consistent calculation process and thus the Pauli matrices σ x and σ z cannot be generated in the final ρ R,HF y (k).Similar argument can be applied to the initial density matrix As a result, the Hamiltonian form has to be and the new density matrix is which can be obtained from the states |u CB1 −i (k)⟩.Although the initial density matrix ρ R ini is independent of k, the H HF (k) in Eq. ( 61) depends on k and the self-consistent density matrix should in principle depend on k.The self-consistent solutions are shown in Fig. 9, in which we evaluate the overlap m (k)⟩ in Fig. 9(a) for the filled bands at half-filling.Furthermore, the Chern number for the band j can be evaluated by where the Berry curvature is calculated by [13] Ω with δk x , δk y as the momenta connecting neighboring momentum grid points in the x, y direction.Our calculation shows C = +1 for the filled band |ψ HF j=1 (k)⟩ with the Berry curvature distribution shown in Fig. 9(d).For ρ R zx (k), the initial density matrices are taken as for a certain uniform value of φ, which corresponds to states cos φ|u R,CB1 + (k)⟩ − sin φ|u R,CB1 − (k)⟩.The HF energy spectrum from this initial ρ R ini in Fig. 10 shows nodes at K, K ′ .These nodes can be understood from nonzero Euler number, denoted as µ, a topological invariant defined for a two-band model with the C 2z T symmetry [14][15][16].The non-interacting eigen-state of CB1 has non-trivial Z 2 number ν CB1 = 1, and the Euler number can be related to the Z 2 number by ν CB1 = µ mod 2 [16].Thus, when ν CB1 = 1, µ has to be an odd number, which gives rise to 2µ of gapless Dirac nodes in the spectrum.Because C 2z T is preserved for the initial density matrix ρ R ini , this symmetry remains throughout the whole self-consistent calculation process, so Euler class is still well-defined for the final self-consistent Hatree-Fock ground state.We evaluate the Wannier center flow for the final Hatree-Fock ground state, which is shown in Fig. 10(d which are denoted as mirror coherent states. The true ground state of the system is obtained by comparing the energies E I [ρ R ] of two self-consistent density matrices in Fig. 10(e).Above the critical interaction value around 0.05E 0 ≈ 2 meV, our calculation shows that the mirror polarized state with ρ R y (k) has lower energies than the non-interacting ground state and the mirror coherent states with ρ R zx (k).This is because non-interacting ground state and mirror coherent state have gapless excitations in their spectrum, while the mirror polarized states are fully gapped.Thus, we conclude that the true ground state is a mirror polarized Chern insulator.
For the case with α = 0, V 0 /E 0 = 1.2, ϕ = 0 without inversion, the mirror symmetry M z is broken so we cannot characterize the non-interacting eigen-state with mirror eigen-values and mirror Chern number.However, the C 2z T symmetry remains, so we can still choose the real gauge for non-interacting states as |u R,CB1 1 (k)⟩, |u R,CB1 2 (k)⟩, which satisfies where n = 1, 2 labels two spin-split bands for the Kramers' pair of CB1.Consequently, two types of order parameters, ρ R y (k) that breaks C 2z T and ρ R zx (k) that breaks spin U (1) symmetry, do not mix with each other.The self-consistent solutions with two types of order parameters are shown in Fig. 11.ρ R y (k) breaks T symmetry and one band of CB1 with nonzero Chern number is gapped from the other band.ρ R zx (k) has Dirac nodes in spectra at K, K ′ with energies E I [ρ R zx (k)] higher than the other case.The ground state is an interaction-driven Chern insulator, same as the inversion symmetric case.

D. Four-band model with CB1 and CB2
In the main text, we have discussed the important role of the band mixing between CB1 and CB2 induced by the Coulomb interaction, which can result in the interacting ground state varying from the QAH state to a trivial insulator state for the realistic Coulomb interaction strength for the inversion symmetric case (V 0 = 0), while the QAH state remains for the realistic Coulomb interaction when a large asymmetric potential V 0 is applied.The difference between the inversion symmetric and asymmetric cases is that both CB1 and CB2 carry non-trivial Z 2 number, ν CB1 = ν CB2 = 1, for inversion symmetric case, while a strong asymmetric potential V 0 gives a trivial insulator phase for CB2, ν CB1 = 1 and ν CB2 = 0, for inversion asymmetric case.This effect can only be taken into account when considering both CB1 and CB2, and thus it is important to go beyond the two-band model discussed above and consider a four-band model with both CB1 and CB2.In this section, we will provide more details of our numerical self-consistent calculations of the interacting ground state within the HF approximations for the four-band model.Below we always assume the 1/4 filling of four bands, which corresponds to the 1/2 filling of CB1.

F. Spin density wave
In this section, we discuss the nonuniform magnetic order parameters within the Hartree-Fock mean field approximation following Ref.[19].
To construct the nonuniform magnetic order parameters, we first write down the Coulomb interaction in Eq. (53) using spin operators under the Hartree Fock approximation as where τ and s are Pauli matrices for the top/bottom surfaces and spin degrees of freedom as in Eq. ( 1) of the main text, respectively.i = x, y, z.E C is the condensation energy.Eq. ( 81) can be restored to Eq. ( 59) by projecting the electron operators into the eigenstate basis (Eq.48), imposing q = 0 and making use of the equality The nonuniform magnetic order parameter or spin density wave has the form of M i (Q) = ⟨S i (Q)⟩ in Ref. [19] with where we have used the eigen-state projection Eq.48 in the second line.As in Ref. [19], Q takes the values of 16(a) that connects K and K ′ in the moiré BZ because of the large density of states around K and K ′ Fig. 16(c).Λ si n,n ′ (k + Q, k) is the form factor of the spin operator defined as and satisfies A real order parameter M(r) = Q e iQ•r M(Q) in the real space requires The order parameters have the C 3z symmetry by where R 3z is the rotation matrix about z axis.C 3z symmetry adapted basis unit vectors for M(Q) is shown in Fig. 16(b) satisfying where µ = 1, 2, 3 and ẑ is the unit vector in the z direction.We decompose M(Q) into basis unit vectors as Here we only consider the non-uniform order parameter with C 3z symmetry, so M µ is independent of Q.The interaction Hamiltonian in the M(Q) order parameter channel is with The following approximation is also adopted The last equation holds because V (Q) only depends on |Q| in Eq. ( 52) and are same for Because M(Q) breaks the moiré translation symmetry, it folds moiré BZ as shown in Fig. 16(a) and couples states at κ, κ + K 1 , κ + K ′ 1 where κ is the momentum in the folded moiré BZ (FMBZ).The explicit form of the Hartree-Fock mean field Hamiltonian in the basis C As H HF is hermitian.In the last equation, we have used Eq. ( 85) and the fact that ) is a moiré reciprocal lattice vector vector.To perform the self-consistent calculations, we choose the initial condition of the magnetic order parameters as corresponding to the magnetic meron lattice discussed in Ref. [19].With H HF [M] and the initial condition, the self-consistent calculations are done until with M(K 1 ), M(K 1 ) as the order parameters for ath and a + 1th iterations, respectively.We next discuss the Hartree-Fock mean-field solutions at half-filling of CB1 for the inversion symmetric case.In the single-particle spectrum, there are Von-Hove singularities near K and K ′ at filling 0.8 of CB1, which is away from the half filling for our self-consistent calculations, as shown in Fig. 16(c).The Von-Hove singularities correpsond to a peak in density of states (DOS) calculated by with broadening γ = 0.005E 0 .Here we choose the Lorentz form instead of δ function for the numerical calculations.The critical Coulomb interaction strength for nonzero magnetic meron texture is around 0.18E 0 , as shown in Fig. 16(d), which is much larger than the critical interaction 0.05E 0 for the QAH phase discussed in the main text.In the moiré TI system, the Coulomb interaction is estimated to be 0.13E 0 , which is smaller than that for magnetic meron texture but larger than that of QAH phase.Thus, we conclude that the QAH state, instead of magnetic meron texture, is the ground state at half filling of the lowest conduction band.We also expect that the magnetic meron texture in Ref Ref. [19] can be stabilized for the filling around 0.8 of the CB1.

Supplementary Figure 3 . 1 K5 e −iπ/ 3 K6 e iπ/ 3 (
(a) Spectra for different ϕ in Fig. 2(a) of the main text.(b) Spectra for different ϕ in Fig. 2(b) of the main text.Insets are real space moiré potentials for different ϕ.States are labelled with the irreducible representations by the little groups C3v at Γ and C3 at K. K5, K6, K4 represent the angular momentum states under C3 with Jz = 1/2, −1/2, 3/2 respectively.ν is the Z2 invariant for the lowest conduction bands CB1.b) K Supplementary Table I. (a)(b) symmetry operators in the irreducible representation at high symmetry momenta Γ and K for the double space group 156 P 3m1 corresponding to the point group C3v.Supplementary Figure 4. (a)(b)(c) The Wannier center flows for CB1 with ϕ = 0, 1/6, 1/3, corresponding to Fig. 3(a)(c)(d), respectively.(d)(e)(f) The Wannier center flows for CB1 with ϕ = 0, 1/6, 1/3, corresponding to Fig. 3(e)(g)(h), respectively.|u Jz=−1/2 (K)⟩ and |u Jz=3/2 (K)⟩ has the Dirac fermion form H ef f (k) = v K (k x σ x + k y σ y ) + m K σ z , up to the linear order, with σ x,y,z are Pauli matrices for the two band basis.The gap closing can be captured by one parameter, namely the Dirac mass m K that is controlled by ϕ, corresponding to the co-dimension 1 case.The T symmetry guarantees the gap closing also occurring at K ′ , and the gap closings at K and K ′ lead to the change of Z 2 number ν by 1.The normal insulator (NI) states are localized at moiré potential minima of the Wyckoff position 1b shown by the insets of spectrum for ϕ = 1/6 in Fig. 3(a) (See SM Sec.I.D).From ϕ = 1/6 to ϕ = 1/3 in Fig. 3(a), another Dirac-type gap closing should happen at K and K ′ , and we find the system with ϕ = 1/3 has ν CB1 = 1.

Supplementary Figure 5 .
(a)(b)(c) Spectra with increasing m2 for α = 1, V0/E0 = 0, ϕ = 0 of Fig. 2(d) in the main text.Green (Orange) dots denote even (odd) parities at Γ and M .(d) Spectrum of 2DEG on the moiré potential with ϕ = 0 shown in Fig. 1(c) of the main text.Spectrum in (c) is labelled with irreps by the little group C6v at Γ and C3v at K.

)
Supplementary Figure 6.(a)-(e) Spectra with increasing m2 for α = 0, V0/E0 = 1.2, ϕ = 0 of Fig. 2(e) in the main text.Different colorful dots Γ represents different irreps of the little group C6v.Spectrum in (c) is labelled with irreps by the little group C6v at Γ and C3v at K. (f) spectrum around K before and after breaking the inversion symmetry.Jz is the angular momentum of the state at K under C3.
with ∆(r) the hexagonal potential, as shown in Fig.1(c) in the main text.The corresponding conduction band dispersion with m z = −i is shown Fig. 5(d), while the m z = +i conduction bands are degenerate with m z = −i bands.The lowest two conduction bands of the Hamiltonian H 2DEG can be viewed as coming from two s-wave atomic Supplementary Figure 7. (a)-(f) Spectra with reducing m2 for α = 0, V0/E0 = 1.2, ϕ = 0 of Fig. 2(e) in the main text.Different colorful dots Γ represents different irreps of the little group C6v shown in Tab.II.Cyan bands are CB1 and CB2.Orange ones are Cb3-5.Black ones are CB6 and higher energy bands.Insets in (d)(f) are enlargement of spectra around Γ in the dashed boxes.

Supplementary Figure 8 .
(a)(b) The real-space maximally localized Wannier functions w R (r) for the lowest conduction bands with ϕ = 1/6 corresponding to Fig. 2(a) of the main text.(c) The real space moiré potentials with ϕ = 1/6.(d)(e)(f) Those for ϕ = 1/2 and (g)(h)(i) Those for ϕ = 5/6.where α = 1, ..., 4 labels both spin and layer index, f † α (k + G) is a fermion creation operator, k is within the first moiré BZ and G is Moiré reciprocal lattice vectors.The creation operators for eigenstates of H 0 (k) are are defined as

2 Supplementary
Figure 10.(a) The spectrum of H HF (k) for the mirror coherent states with ρ R zx (k) (black lines) .(b)(c) The overlap between mirror coherent states |ψ HF j=1 (k)⟩ for the filled band in (a) and the non interacting states |u CB1 ±i ⟩(k).(d) The Wannier center flow for both eigen-states of H HF (k).(e) The energy per particle EI[ρ R ] with non-interacting energy EI[ρ0] subtracted for ρ R y (k) and ρ R zx (k).Here the calculation is for the two-band model with the parameters ϕ = 0, α = 1, V0/E0 = 0. which has the C 2z T symmetry, [ρ R new (k), C 2z T ] = 0.So the Pauli matrix σ y cannot be mixed into the density matrix ρ R,HF zx (k) in the above procedure.Based on this symmetry argument, we can discuss the self-consistent solutions for the density matrix form ρ R y (k) and ρ R zx (k), separately, below.For ρ R y (k), we choose the initial density matrix as ).The nonzero Euler class with µ = 1 from the Wannier center flow guarantees the existence of 2 Supplementary Figure 11.(a)(b) The spectra (orange) of H HF (k) with ρ R y (k) for (a) and ρ R zx (k) for (b).The cyan lines are non-interacting spectrum.(c) Berry curvature for the lower band of H HF (k) in (a).(c) The energy per particle EI[ρ R ] with non-interacting energy EI[ρ0] subtracted for ρ R y (k) and ρ R zx (k).Here the calculation is for the two-band model with the parameters ϕ = 0, α = 0, V0/E0 = 1.2.Dirac nodes in the Hartree-Fock spectrum.Fig. 10(b)(c) shows that the Hartree-Fock solutions |ψ HF j=1 (k)⟩ shown in Fig. 10(a) are superposition of two m z states with the same probability

Supplementary Figure 16
CB1 ±i (k)⟩ with C = ±1 and |u CB2 ±i (k)⟩ with C = ∓1, where C denotes the Chern number of the minibands in the m z = −i .(a) The moiré BZ (MBZ) and the folded moiré BZ (FMBZ).(b) The C3z symmetry adapted basis unit vectors for M(Q).(c) The single-particle spectrum and density of states (DOS) for CB1 (blue lines).The black dotted line is the Fermi energy with Von-Hove singularities while the red dotted line is the Fermi energy for half filling of CB1.(d) The magnetic texture order parameters |M(K1)| versus different Coulomb interaction strength U (a M 1 )/E0.as U (a M 1 ) = 10 meV.The spectrum is shown in Fig. 15(a).The parameters for the spectrum are ∆ 1 = 14meV, ϕ = 1/3, and m 0 = 30meV.The energy scale is E ′ 0 = v|b M 1 | = 77meV.The bandwidth of CB1 is 7.4meV and the direct gap between CB1 and CB2 is 4meV.The ratio between U (a M 1 ) and bandwidth is smaller for the smaller |a M 1 |.When the Coulomb interaction is considered for CB1 with the density matrix ρ R y (k) as shown in Fig. 15(b), the mirror polarized states with C = +1 can be induced.

Supplementary Figure 19 . 2 . 1 , β 1 |
(a) The side view of the Sb2/Sb2Te3 heterostructure with AA stacking.The solid black lines mark the unitcell used in DFT calculations.Sb atoms in Sb2 monolayer are marked as gray, Sb atoms and Te atoms in Sb2Te3 films are marked as red and blue.The atoms in the region with the green background are frozen when we relax the lattice structures, the x and y coordinates of Sb atoms in the region with the yellow background are fixed.The value of the van der Waals gap inside 2QL Sb2Te3 films is 2.078 Å. (b-d) The top view for heterostructures with the stacking of AAmAA-X, AAmAB-X, and AAmBA-X (X = I, II, III).Corresponding dR is shown by the green arrows.The black arrows show the lattice vector of Sb2 monolayer.(e) Calculated band structures (cyan lines) of the heterostructure with the stacking of AAmAA, AAmAB, and AAmBA, respectively.The Fermi levels are set as zero.The orange dashed lines are fitted spectra by the Hamiltonian Eq. (100) for a uniform shift.The effective Hamiltonian ĤDFT (d R ) is then given by⟨k 1 , β 1 | ĤDFT (d R )|k 2 , β 2 ⟩ = δ(k 1 − k 2 )( δ(k 1 − k 2 )H DFT (k 1 , d R ),(100) where H DFT (k, d R ) is just Eq.(4) in the main text, h t/b D (k) are the top/bottom Dirac surface states same as Eq.(1) of the main text, s 0 are the identical matrix in spin space, m is the tunnelling between two surfaces, and α captures the difference in the potentials on two surfaces created by the Sb 2 layer.|k, β⟩ is the atomic Bloch states for the Sb 2 Te 3 and Sb 2 lattice with a constant shift.β 1,2 = 1, ..., 4 represents both the spin and layer degrees of freedom.The spectra of this model are given byE DFT η,ξ (k, d) = 1 + α 2 ∆(d) + η m 2 + 1 − α 2 ∆(d) + ξv 2 k 2 ) (101) with η = ±, ξ = ±.By fitting E DFT η,ξ (k, d) to the spectrum calculated from DFT in Fig.4(b) of the main text and Fig. 19(f), the parameters α, ∆(d), v and m in the model Hamiltonian can be obtained and documented in Tab.V for |k, β⟩ are atomic Bloch states for each layer underlying the moiré superlattice.As the moiré Hamiltonian H 0 does not preserve atomic lattice translation, the crystal momentum k is not a good quantum number and H 0 can mix different k states.The atomic Bloch wave function |k, β⟩ is related to atomic Wannier function |R, β⟩ by |k, β⟩ = R e ik•R |R, β⟩,(104)so the moiré Hamiltonian is transformed into the form on the atomic Wannier function basis as⟨k 1 , β 1 | Ĥ0 |k 2 , β 2 ⟩ = R1R2 e −ik1•R1 ⟨R 1 , β 1 | Ĥ0 |R 2 , β 2 ⟩e ik2•R2 .(105)Here ⟨R 1 , β 1 | Ĥ0 |R 2 , β 2 ⟩ describes the Hamiltonian matrix element between atomic Wannier functions located at R 1 and R 2 in the superlattice shown in Fig.20(a).As the overlap between atomic Wannier functions decays quickly as the distance increases, we only consider the local hopping within the length scale|R 2 − R 1 | ∼ O(|ã 1 |), where ã1 is the atomic primitive lattice vector for the Sb 2 Te 3 layer.In this atomic length scale, the Hamiltonian matrix element between two Wannier orbitals near R on the superlattice structure with twist angle θ in Fig.20(a) can be approximated locally by the Hamiltonian matrix element for two atomic layers with a constant shift d R in Fig.20(b)[33], whered R = R(θ)R − R(106)and R(θ) as the rotation operator for the Sb 2 layer with the rotating angle θ.This approximation is valid for a small twist angle θ because the local shift vector d R is almost uniform at the atomic length scale,d R2 ≈ d R1(107)for|R 2 −R 1 | ∼ O(|ã 1 |).The Hamiltonian matrix element between two atomic Wannier orbitals for the commensurate lattice is captured by ⟨R 1 , β 1 | ĤDFT (d R2 )|R 2 , β 2 ⟩, so we make the approximation⟨R 1 , β 1 | Ĥ0 |R 2 , β 2 ⟩ ≈ ⟨R 1 , β 1 | ĤDFT (d R2 )|R 2 , β 2 ⟩(108)and⟨k Ĥ0 |k 2 , β 2 ⟩ = R1R2 e −ik1•R1 ⟨R 1 , β 1 | ĤDFT (d R2 )|R 2 , β 2 ⟩e ik2•R2 .(109)To extract R 2 in ĤDFT (d R2 ) for the summation, we transform ĤDFT (d R ) to the momentum-space byĤDFT (d R ) = G e −i G•d R ĤDFT ( G)(110)as ĤDFT (d R + xã 1 + yã 2 ) = ĤDFT (d R ) is periodic for atomic lattice vectors (x, y are integers here), as shown in Fig.20(b).We also denote the atomic reciprocal lattice vector G = Gwz = w b1 + z b2 with integers w, z, so the summation over G is equivalent to the summation over w, z. b1,2 are atomic reciprocal lattice vectors satisfying bi • ãj = δ ij for i, j = 1, 2 and shown in Fig.20(c).SinceGwz • d R = Gwz • (R(θ)R − R) = ( Gwz − R(θ) Gwz ) • (R(θ)R) = G wz • (R + d R ) ≈ G wz • R(111) with the moiré reciprocal lattice vectors G wz given by G wz = Gwz − R(θ) Gwz = (w b1 + z b2 ) − R(θ)(w b1 + z b2 ) = w( b1 − R(θ) b1 ) + z( b2 − R(θ) b2 ) = w( b1 − b′ 1 ) + z( b2 − b′ 2 ) = wb M 1 + zb M 2 , (112) we have e −i G•d R ≈ e −iG•R , ĤDFT (d R ) ≈ w,z e −iGwz•R ĤDFT ( Gwz ).(113) (f).The other two states |ũ 1/2,i (K)⟩, |ũ 1/2,i (K)⟩ ( K6 irrep) remain degenerate and form a 2D irrep under the little group C 3v at K. When this energy splitting ∆ R is larger than the Kane-Mele SOC gap ∆ KM , the degenerate states with the 2D irrep K6 lies between the K4 and K5 state, leading to the band touching between CB1 and CB2 bands at K for m 2 = 1.5E 0 /|k M | 2 in Fig.