Introduction

Twisted bilayer graphene (TBG) burst on the scene as a tunable two carbon-atom layers thick system realizing a remarkable multitude of interaction-driven macroscopic quantum phenomena1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20. Although significant progress has been achieved in understanding the nontrivial topology of the narrow bands, as well as the correlated electron states in the magic-angle TBG21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42, many important questions remain open. One of the most fascinating question is the origin and the mechanism of the quantum anomalous Hall (QAH) state with Chern number C = ±16,7 at three-quarters filling of the system, aligned with the hexagonal boron nitride (hBN), and the insulating state which replaces the QAH in devices without the hBN alignment.

Currently, the prevailing opinion is that the QAH can be obtained from narrow band models with large Coulomb interactions28,43,44,45,46, but that the nontrivial topology of the narrow bands prevents a faithful construction of local “Hubbard-like” tight-binding models that locally respect all the symmetries23. Although there exists no a priori Wannier obstruction, as the narrow bands’ total Chern number vanishes, there is yet no clear understanding of how the QAH could arise within such correlated lattice model, even in principle, in the limit where the Coulomb interactions dominate the kinetic energy.

Precisely such a state was sought by Raghu et al. in an entirely different context47, coining the term topological Mott insulator (TMI), which we define to be a QAH in a strong coupling limit of a local lattice model with a vanishing ratio of the bandwidth to Coulomb interaction. However, the original proposal47 was subsequently shown not to host a QAH, and therefore not TMI either48,49. More recent works have found the interaction-induced QAH state in a different model, but it is stabilized by the kinetic energy and necessitates sizable bandwidth50,51,52. Because it gives way to more conventional Mott insulators in the strong coupling regime52, these models do not host a TMI.

Here we show that the TMI is realized in a simple lattice model introduced by two of the authors as a local description of the correlations within the TBG narrow bands26,53,54. The key new ingredients are the off-site terms appearing alongside the usual on-site terms in the projected density operator. Physically, such terms originate in the extended multi-peak nature of the maximally localized Wannier states22,24 arising from the nontrivial topology25,25,35,55,56,57,58 of the narrow bands, and, importantly, remain finite even when the bandwidth vanishes.

Results

Honeycomb moiré lattice model

In the strong coupling limit, the aforementioned model (as illustrated in the upper panels of Fig. 1) is

$$H = U_{0}{\mathop{\sum}\limits_{\hexagon}}(Q_{\hexagon}+\alpha T_{\hexagon}-1)^{2},$$
(1)

where U0 constitutes the overall energy scale in the problem (≈40 meV in TBG and set to unity henceforth). \(Q_{\hexagon}\equiv \,\frac{1}{3}\mathop{\sum }\nolimits_{l = 1}^{6}{c}_{{{{{{{{\bf{R}}}}}}}}+{\delta }_{l}}^{{{{\dagger}}} } {c}_{{{{{{{{\bf{R}}}}}}}}+{\delta }_{l}}^{\,}\) represents the cluster charge term22,25,53,54,59,60 (c.f. Fig. 1c), and \(T_{\hexagon}\equiv\, \mathop{\sum }\nolimits_{l = 1}^{6}[{(-1)}^{l}{c}_{{{{{{{{\bf{R}}}}}}}}+{\delta }_{l+1}}^{{{{\dagger}}} }{c}_{{{{{{{{\bf{R}}}}}}}}+{\delta }_{l}}^{\,}+h.c.]\) represents the Coulomb induced hopping with alternating sign (c.f. Fig. 1d). Fermion annihilation and creation operators \({c}_{{{{{{{{\bf{R}}}}}}}}+{\delta }_{l}}^{\ }\) and \({c}_{{{{{{{{\bf{R}}}}}}}}+{\delta }_{l}}^{{{{\dagger}}} }\) are defined at the sites of the honeycomb lattice R + δl, where R = m1L1 + m2L2 with integer m1,2 spans the triangular Bravais lattice. The hexagon centers, over which we sum in Eq. (1), are connected to the six nearest honeycomb lattice sites l = 1, 2, ⋯ 6 through δl (c.f. Fig. 1e). As we focus on the three-quarters filling of the TBG, where the spin and orbital degrees of freedom are assumed to be polarized, Eq. (1) thus constitutes a simplification to the full Hamiltonian of ref. 26. The parameter α controls the relative strength of charging and assisted-hopping of the projected Coulomb interaction. It originates from the overlap of two neighboring Wannier states in the continuum model and thus depends on the lattice relaxation. Due to the background charge from the remote bands, which is approximated to be uniform in Eq. (1), the projected Coulomb interaction is in the form of density-density repulsion43,61,62, instead of being normal ordered. Although the projected interaction contains other terms such as next-nearest neighbor interaction, the more detailed calculations at the chiral limit have shown that the interaction-induced dispersion of the charged excitation at the charge neutrality point is dominated by α, the nearest neighbor assisted hopping63.

Fig. 1: The honeycomb moire lattice model and phase diagram.
figure 1

a YC and b XC geometries with PBC along vertical (L1 − L2 for XC and 2L1 − L2 for YC) and OBC along horizontal direction. The number of sites on the cylinders is N = W × L × 2, with length L (the number of vertical armchair/zigzag chains, c.f. the gray-shaded lines) and W is the number of 2-site unit cells (c.f. the red-shaded rectangles) along those chains. c Shows the cluster charge operator Q⎔, which counts the electron number in a hexagon and d demonstrates the assisted hopping term T with alternating-sign structure. e The labeling of six sites within hexagon R. f The phase diagram contains two distinct insulating phases, i.e., the stripe phase for α < αc, and the QAH state for α > αc ≃ 0.12. g The schematic plot of the emergent current through a mean-field tight-binding analysis of the QAH state.

Full size image

The original bandwidth W ~ 8 meV24 is much smaller than U0, suggesting the system is in the strong coupling regime. Furthermore, after the states on the remote bands are integrated out, the superexchange interaction (≲5 × 10−3e2/(ϵLm)) is found to be negligible compared with the projected Coulomb interaction61; this justifies neglecting additional fermion bilinear (kinetic) terms in Eq. (1). The kinetic term, as well as the further-range assisted hopping terms, may shift the critical value αc of the phase transition but do not qualitatively change the phase diagram in Fig. 1f. In addition, we do not include the additional symmetry breaking term produced by the possible hBN alignment that favors the QAH phase64, but focus on the topological phase transitions purely driven by interactions.

It is worth emphasizing that Eq. (1) corresponds to the leading order terms when the distance to the gates lg is about the same as the moiré lattice constant ∣L1∣, and thus the electron-electron repulsion decays exponentially when the inter-electron separation is larger than ∣L1∣26. With larger lg, the longer range aspect of the Coulomb repulsion will have to be included, but because currently there is no experimental indication that there are significant changes in the nature the insulating states for different lg1,3,65, it is reasonable to neglect the longer range terms in Eq. (1). We should note that terms in Eq. (1) are purely real, and because the two QAH states with opposite Chern numbers transform into each other under complex conjugation, the QAH state is not a priori favored by this model. In what follows, we will demonstrate that, for a range of α, Eq. (1) naturally leads to the TMI ground state via spontaneous symmetry breaking without including any other interactions or kinetic terms.

Phase diagram

We solve the TBG lattice model in Eq. (1) using DMRG on long cylinders of XC (zigzag, Fig. 1a) and YC (armchair, Fig. 1b) geometries, with widths W up to 6 and lengths L up to 24. The details of DMRG implementation and finite-size analysis are given in the Methods and Supplementary Note 1. The obtained ground state phase diagram, as a function of α, is shown in Fig. 1f. We identify two gapped insulating phases: a stripe phase with charge density wave (CDW) for small α, and a TMI phase for α > αc ≈ 0.12. These two ground states are separated by a first-order quantum phase transition (QPT). In Fig. 2, we show results for various quantities, including the ground state energy eg, entanglement entropy SE, charge structure factor Cn, and the imaginary part of the equal time correlation \(\langle J\rangle \equiv \frac{{i}}{2}\langle ({c}_{l}^{{{{\dagger}}} }{c}_{l^{\prime} }^{\ }-{c}_{l^{\prime} }^{{{{\dagger}}} }{c}_{l}^{\ })\rangle\). As shown in Fig. 2a, the eg curve exhibits a discontinuity in the slope (a kink) at αc, indicating the first-order QPT. In Fig. 2b, we calculate the entanglement entropy \({S}_{E}(x)\equiv -{{{{{{{\rm{Tr}}}}}}}}[{\rho }_{{{{{{{{\mathcal{A}}}}}}}}}(x){{{{{{\mathrm{ln}}}}}}}\,({\rho }_{{{{{{{{\mathcal{A}}}}}}}}}(x))]\), with \({\rho }_{{{{{{{{\mathcal{A}}}}}}}}}(x)\) the reduced density matrix of the subsystem \({{{{{{{\mathcal{A}}}}}}}}\) consisting of the first x columns (c.f. Fig. 1a, b). By setting x = L/2 (for even L), i.e., cutting at the very center of the system, we compute SE(L/2) and show it vs. α in Fig. 2b, where an evident “jump” takes place right at the QPT. In addition, for α < αc, the negligibly small SE(L/2) indicates the existence of a nearly direct product state with virtually no charge fluctuations in the CDW pattern. On the other hand, the sizable SE(L/2) for α > αc indicates a finite amount of quantum entanglement in the ground state. In the insets of Fig. 2b, SE(x) vs. subsystem length x shows a flat plateau in the bulk of the system, indicating that both phases in Fig. 1f are gapped, consistent with the exponentially decaying single-particle Green’s functions also obtained by our DMRG (see the Supplementary Note 1).

Fig. 2: Identification of two insulating phases.
figure 2

a The ground-state energy per site \({e}_{g}\equiv \frac{1}{N}\langle {\psi }_{g}|\hat{H}|{\psi }_{g}\rangle\), shown as a function of α, with total number of sites N = 2WL and \(|{\psi }_{g}\rangle\) the DMRG ground state. b Entanglement entropy SE, c stripe order parameter Cn(M), d both correlations \(\langle J\rangle_{{{\mathrm{NN}}}}\) and \(\langle J\rangle_{{{\mathrm{NNN}}}}\), are shown versus α, all showing abrupt changes of behavior at αc ≃ 0.12. The mean-field energies for both phases are as well shown in a. The detailed entanglement profile SE vs. subsystem x is shown in the inset of (b), and Cn(k) vs. k in the first Brillouin zone (BZ) shown in the inset of (c).

Full size image

Stripe and QAH insulators

The emergence of the stripe phase at small α can be understood from a perturbative analysis26. Up to second-order corrections (c.f. Supplementary Note 2), we find the ground-state energy eg/U0 ≃ α2, and plot it together with the DMRG results in Fig. 2a, where the high accuracy of this analytical calculation can be clearly seen. The CDW order can be characterized by the structure factor, \({C}_{n}({{{{{{{\bf{k}}}}}}}})\equiv \frac{1}{N}\mathop{\sum }\nolimits_{\lambda = 1}^{2}{\sum }_{{{{{{{{\bf{R}}}}}}}}}{e}^{-i{{{{{{{\bf{k}}}}}}}}\cdot ({{{{{{{\bf{R}}}}}}}}+{\delta }_{\lambda })}{\tilde{n}}_{{{{{{{{\bf{R,\lambda }}}}}}}}}\), where the quantity \({\tilde{n}}_{{{{{{{{\bf{R,\lambda }}}}}}}}}=\langle {c}_{{{{{{{{\bf{R}}}}}}}}+{\delta }_{\lambda }}^{{{{\dagger}}} }{c}_{{{{{{{{\bf{R}}}}}}}}+{\delta }_{\lambda }}^{\ }\rangle -1/2\) counts the number of electrons (with respect to the half filling) on the honeycomb site R + δλ. In Fig. 2c, we find that Cn(k) peaks at \({{{{{{{\bf{M}}}}}}}}=(0,\frac{2\pi }{\sqrt{3}| {{{{{{{{\bf{L}}}}}}}}}_{1}| })\) for α < αc, and drops abruptly to 0 for α > αc, confirming that the small-α regime has a CDW order, while for α > αc the insulating phase has no charge order. Remarkably, this α > αc regime turns out to be a topological phase with spontaneous time-reversal symmetry (TRS) breaking and a quantized Hall conductance, i.e., a QAH phase.

To reveal the TRS breaking in the large-α QAH phase, in Fig. 2d we show the correlation \(\langle J\rangle\) on both the nearest-neighbor (NN) and next-nearest-neighbor (NNN) \((l,l^{\prime} )\) pairs. We find a finite value of \(\langle J\rangle_{{{\mathrm{NN}}}}\) ~ 0.22 and \(\langle J\rangle_{{{\mathrm{NNN}}}}\) ~ 0.1 in the bulk of the cylinder for large-α phase, while they vanish in the stripe phase. In the QAH phase, the real part of \(\langle {c}_{l}^{{{{\dagger}}} }{c}_{l^{\prime} }^{}\rangle\) is negligibly [O(10−7~−8)] smaller compared to its imaginary part, and thus \(\langle {c}_{l}^{{{{\dagger}}} }{c}_{l^{\prime} }^{}\rangle\) emerging from interactions is virtually purely imaginary. The corresponding hopping process thus acquires a π/2 phase (labeled as i in Fig. 1g), rendering a 3π/2 flux for a circulating triangular loop current, which resembles the Haldane model66. The difference is that the TRS breaking NNN hopping term is introduced explicitly in the Haldane model, while here it emerges spontaneously due to electron interactions, a typical feature of TMIs. We also note that in a recent quantum Monte Carlo simulation applied at charge neutrality53 (i.e. even integer filling), a quantum valley Hall state is found at intermediate coupling for a specific choice of kinetic energy terms. Such a state is different from the QAH found at odd integer filling here as it preserves the TRS with helical valley edge modes and undergoes a first-order phase transition into intervalley coherent insulator at strong coupling, consistent with the exact results obtained in ref. 26.

Quantized Hall conductance

To reveal the topological properties in the large-α phase, we perform a flux insertion experiment on the cylindrical geometry (c.f. the inset of Fig. 3a) and compute the Hall conductance. We thread a ϕ-flux along the cylinder by modifying the boundary condition \({c}_{{{{{{{{\bf{R}}}}}}}}+W({{{{{{{{\bf{L}}}}}}}}}_{1}-{{{{{{{{\bf{L}}}}}}}}}_{2})+{\delta }_{\lambda }}\equiv {c}_{{{{{{{{\bf{R}}}}}}}}+{\delta }_{\lambda }}\) to \({c}_{{{{{{{{\bf{R}}}}}}}}+W({{{{{{{{\bf{L}}}}}}}}}_{1}-{{{{{{{{\bf{L}}}}}}}}}_{2})+{\delta }_{\lambda }}\equiv {e}^{-i\phi }{c}_{{{{{{{{\bf{R}}}}}}}}+{\delta }_{\lambda }}\) for XC geometry and \({c}_{{{{{{{{\bf{R}}}}}}}}+W({{{{{{{{\bf{L}}}}}}}}}_{1}-{{{{{{{{\bf{L}}}}}}}}}_{2}/2)+{\delta }_{\lambda }}\equiv {c}_{{{{{{{{\bf{R}}}}}}}}+{\delta }_{\lambda }}\) to \({c}_{{{{{{{{\bf{R}}}}}}}}+W({{{{{{{{\bf{L}}}}}}}}}_{1}-{{{{{{{{\bf{L}}}}}}}}}_{2}/2)+{\delta }_{\lambda }}\equiv {e}^{-i\phi }{c}_{{{{{{{{\bf{R}}}}}}}}+{\delta }_{\lambda }}\) for YC geometry. During the process of the flux insertion, ϕ is adiabatically increased from 0 to 2π in the DMRG calculations. One thereafter obtains the Hall conductance \({\sigma }_{H}=\frac{{e}^{2}}{h}{{\Delta }}Q\) by measuring the net charge pumping ΔQ from one edge of the cylinder to the other. In DMRG, we calculate the net charge transfer as \({{\Delta }}Q=\mathop{\sum }\nolimits_{x = L-l+1}^{L}[{\tilde{n}}_{x}^{{{{{{{{\rm{col}}}}}}}}}(\phi )-{\tilde{n}}_{x}^{{{{{{{{\rm{col}}}}}}}}}(0)]\), i.e. the pumped charge to the rightmost l columns (chosen as l = 3–4 in practice) where \({\tilde{n}}_{x}^{{{{{{{{\rm{col}}}}}}}}}(\phi )\) is the deviation of the charge number of the x-th column measured in the ϕ-flux inserted ground state \(|{\psi }_{\phi }\rangle\) from the half filling. For instance, we have \({\tilde{n}}_{x}^{{{{{{{{\rm{col}}}}}}}}}(\phi )=\mathop{\sum }\nolimits_{y = 1}^{W}\mathop{\sum }\nolimits_{\lambda = 1}^{2}\langle {\psi }_{\phi }|{\hat{n}}_{(x-1){{{{{{{{\bf{L}}}}}}}}}_{{{{{{{{\bf{1}}}}}}}}}+y({{{{{{{{\bf{L}}}}}}}}}_{{{{{{{{\bf{1}}}}}}}}}-{{{{{{{{\bf{L}}}}}}}}}_{2})+{\delta }_{\lambda }}-\frac{1}{2}|{\psi }_{\phi }\rangle\) for the XC geometry, and similar expressions for YC.

Fig. 3: Quantized hall conductance and QAH state.
figure 3

In the systems of both width W = 4, 6, a flux ϕ ∈ [0, 2π] is threading through the cylinder, a One electron is pumped from one edge to the other for α = 0.15 (QAH phase), while no charge response is observed for α = 0.1 (stripe phase). b In the real-space charge distribution of YC4 cylinder, no accumulations are observed in the bulk, i.e., only the charge near the left edge is pumped. c Entanglement spectrum computed at the central bond (between two columns) shows a two-fold degeneracy. For a typical QAH state with α = 0.25, we show in d the charge density nλ(k), with λ labeling the two eigenvalues of the 2 × 2\(\tilde{G}({{{{{{{\bf{k}}}}}}}})\) matrix associated with two sublattices. In e the von Neumann Entropy \({\bar{S}}_{{{{{{{{\rm{vN}}}}}}}}}\) averaged over all the k points, is shown versus α, where the SvN distribution in BZ is shown in the inset (with also α = 0.25).

Full size image

As shown in Fig. 3a, for both XC and YC systems (with widths W = 4 and 6) in the QAH phase (e.g., α = 0.15), we find a net charge transfer ∣ΔQ∣ = 1 through a 2π flux insertion, showing that the Chern number C = ±1. In addition, Fig. 3b shows the column charge distribution \({\tilde{n}}_{x}^{{{{{{{{\rm{col}}}}}}}}}\), where a half-charge \(\pm\kern-2pt \frac{1}{2}\) appears in two edges in \(|{\psi }_{\phi = 0}\rangle\). As ϕ gradually increases, the left/right-end charge smoothly reduces/increases from \(\pm \kern-2pt\frac{1}{2}\)  to \(\mp \frac{1}{2}\), which corresponds to an end-to-end pumping of a unit charge ΔQ = 1, without “disturbing” the charge distribution in the bulk. We note that there is two-fold degenerate QAH ground state (apart from the additional degeneracy due to half-charge zero edge modes, see discussion below), and the charge pumping could be ΔQ = ±1, corresponding to the spontaneous TRS breaking states with C = ±1.

Understanding the TMI phase

With DMRG calculations, we can also calculate the single-particle Green’s function \({G}_{\lambda ,\lambda ^{\prime} }({{{{{{{\bf{R}}}}}}}}-{{{{{{{\bf{R}}}}}}}}^{\prime} )=\langle {c}_{{{{{{{{\bf{R}}}}}}}}+{\delta }_{\lambda }}^{{{{\dagger}}} }{c}_{{{{{{{{\bf{R}}}}}}}}^{\prime} +{\delta }_{\lambda ^{\prime} }}^{\ }\rangle\), from which we can find the electron occupation \({n}_{\lambda ,\lambda ^{\prime} }({{{{{{{\bf{k}}}}}}}})\) in the momentum space. Due to the two-sublattice structure, \({G}_{\lambda ,\lambda ^{\prime} }({{{{{{{\bf{R}}}}}}}}-{{{{{{{\bf{R}}}}}}}}^{\prime} )\) and its Fourier transformation \({\tilde{G}}_{\lambda ,\lambda ^{\prime} }({{{{{{{\bf{k}}}}}}}})\) are both 2 × 2 matrices (cf., Supplementary Note 1). The two eigenvalues {n1(k), n2(k)} of \(\tilde{G}({{{{{{{\bf{k}}}}}}}})\) are shown in Fig. 3d. We find for all allowed k points, the larger eigenvalue n2(k) ≃ 1 and the smaller value n1(k) ≃ 0, representing the “two-orbit” electronic structure with one orbit filled while the other left empty. Albeit small, charge fluctuations between the two orbits are still present. We compute the von Neumann entropy \({S}_{{{{{{{{\rm{vN}}}}}}}}}({{{{{{{\bf{k}}}}}}}})\equiv -\mathop{\sum }\nolimits_{\lambda = 1}^{2}{n}_{\lambda }({{{{{{{\bf{k}}}}}}}}){{{{{{\mathrm{ln}}}}}}}\,{n}_{\lambda }({{{{{{{\bf{k}}}}}}}})\) that measures the deviation of the DMRG ground state from a Slater determinant of Bloch states. In Fig. 3e, we show the calculated SvN averaged over the first BZ, which decreases as α increases, and becomes very small for large α cases. For example, we show the detailed k-dependent profile for the α = 0.25 case, in the inset of Fig. 3e. The relatively small SvN values suggest the QAH state, emerging in the interacting TBG model as revealed by DMRG calculations, actually very much resembles the Slater determinant ground state of the Haldane model and thus can be captured by a mean-field description.

To be specific, for small α, a second-order perturbation shows the charging term \({\mathop{\sum}\limits_{\hexagon}}(Q_{\hexagon}-1)^{2}\) favors the insulating phases in which each hexagon of the honeycomb lattice contains exactly one electron, i.e. \(Q_{\hexagon}=1\) for every hexagon. Among all the states satisfying this requirement, the first- and second-order corrections from the cross terms \(T_{\hexagon}(Q_{\hexagon}-1)\) vanish. The stripe phase is selected from such states because it minimizes the contribution of \(\langle{\mathop{\sum}\limits_{\hexagon}} T_{\hexagon}^{2}\rangle\), with the energy \(\langle H\rangle_{{{\mathrm{stripe}}}}\) ≈ α2U0 (c.f. Supplementary Note 2).

For large α, motivated by the resemblance of the DMRG ground state to the Slater determinant, we perform a variational mean-field calculation that approximates the true ground state with the ground state of a tight-binding model containing various hoppings (see Methods and Supplementary Note 3). In particular the Fig. 1g demonstrates the emergence of NNN currents which constitute a loop in each hexagon, spontaneously choosing either the left- or right-chiral direction (here the right chirality). We find that the cross terms, i.e. \(\langle T_{\hexagon} (Q_{\hexagon} -1)\rangle_{{{{{\mathrm{QAH}}}}}}\) become negative and thus favor the QAH phase. Including both the charging terms and \({\mathop{\sum}\limits_{\hexagon}} T_{\hexagon}^{2}\), the variational mean-field analysis results in \(\langle H\rangle_{{{\mathrm{QAH}}}}\) ≈ U0(0.037 − 0.27α + 0.71α2). Therefore, as α continuously increases from 0, the mean-field theory also finds the first-order phase transition from the stripe phase to the QAH, in agreement with the DMRG result mentioned earlier. The mean-field energy is shown in Fig. 2a as indicated by the blue and red dashed line for the stripe and QAH phases respectively. Both lines provide a good approximation to the DMRG energy curve, and the intersection of two mean-field energies also provides a very good estimate of the QPT value \({\alpha }_{c}^{{{{{{{{\rm{MF}}}}}}}}}\simeq 0.125\). Interestingly, the energy difference between the mean-field approximation and the DMRG calculation decreases as α moves away from the QPT, reflecting the suppression of the quantum fluctuations for large ∣α − αc∣, also illustrated by the SvN in Fig. 3e.

Moreover, as shown in Fig. 3b, there exist half-charge zero modes on both edges of the cylinder with even W, which also coincide with the Haldane model wrapped on the cylinder (for more details, see the Supplementary Note 4). We also compute the entanglement spectrum (ES), defined as \({\xi }_{i}\equiv -{{{{{{\mathrm{ln}}}}}}}\,({\rho }_{i})\) with ρi the eigenvalues of the reduced density matrix. As shown in Fig. 3c, when we cut at the center of the system, a two-fold degeneracy in the ES is evident, which accounts for the half-charge zero modes in the edge (c.f. Fig. 3b), through the bulk-edge correspondence.

Discussion

As we mentioned, the QAH can be obtained from narrow band models of TBG with large Coulomb interactions, but these models are built in the basis of extended states27,28,43,44,45 making the interaction potential rather unwieldy. The results indeed show that several phases: QAH, strongly correlated topological semimetal, and insulating stripe phases, are energetically competitive for the ground states at odd integer fillings28,35,44,46,67.

The common belief, however, is that the nontrivial symmetry-protected topology of the narrow bands prevents a faithful construction of models within exponentially localized basis even when the bands’ total Chern number vanishes23. On the other hand, as first shown in the context of the Z2 topological insulators68, the obstruction is not as severe as in the case of a nonzero Chern band (or band composite). If the total Chern number vanishes, the exponentially localized Wannier states can be constructed69, but some of the protecting symmetries do not have a simple on-site implementation68,70,71. Because the transformation from the Bloch to Wannier basis is unitary and no information is lost in the process, it is therefore expected that the lattice tight-binding description should also result in the same ground state as found in unobstructed, extended states, basis. However, any practical implementation of this program needs to truncate the expansion of the interaction to on-site and few nearest neigbour sites. What is not obvious, therefore, is whether all the terms need to be included in the expansion or whether it can be truncated to recover the ground state.

The results presented here show that the truncation at just the nearest neighbor, parameterized by α in Eq. (1), is sufficient to recover the insulating and the topologically nontrivial phases. In addition, the main features of the single-particle excitation dispersion of the strong coupling correlated ground states at the charge neutrality point53 from the model in Eq. (1) match those computed exactly in the extended basis61,67. This demonstrates the practicality of Wannier description even for such symmetry-obstructed bands. Our real-space interaction-only model therefore establishes the microscopic mechanism of the evolution between the insulating stripe and QAH phases. Our effective model and its unbiased numerical solution therefore revealed the essence of the physics in this particular regime, and is also consistent with other theoretical calculations27,28,44,46.

As for relevance of our model towards the real system, it is understood that other than the \(Q_{\hexagon}\) and \(T_{\hexagon}\) terms, we do not include all the other projected interactions nor the small kinetic terms, i.e., the detailed feature of the TBG material, which will surely modify the specific value of αc. Apart from that they should not qualitatively alter the two phases and thus also the main conclusion of the present work. In addition to the ground states given above, the dispersion of the charged excitations produced by Eq. (1) is also found to be qualitatively consistent with more detailed calculation by two of the authors in refs. 61,63. Reference63 has also explicitly shown that the dispersion at the charge neutrality point is dominated by the α term in the chiral limit. For systems away from the chiral limit, it is expected that the inclusion of other terms may only quantitatively change the dispersion.

Methods

Density matrix renormalization group

We employ the DMRG method, realized in the matrix product state form and with U(1) charge symmetry implemented, to accurately find the ground state of the TBG model. Following standard 2D DMRG calculations, we map the cylindrical geometries through a snake-like path, i.e., a quasi-1D structure, where highly controllable and efficient simulations can be performed. In practice, we retain up to D = 512(1024) for W = 4(6) cylinders, with truncation errors ϵ < 5 × 10−5, for an accurate large-scale calculations. The detailed convergence check of the TBG model calculations can be seen in the Supplementary Note 1.

Mean-field analysis

We also applied the mean-field theory to approximate the interactions by a tight-binding model with variational hopping constants. The hopping amplitudes are obtained by minimizing the expectation value of the interactions in Eq. (1) for the state produced by the tight-binding model. In practice, the tight-binding model includes hopping amplitudes up to the 5th nearest neighbor. The details are presented in the Supplementary Note 3.