Realization of topological Mott insulator in a twisted bilayer graphene lattice model

Magic-angle twisted bilayer graphene has recently become a thriving material platform realizing correlated electron phenomena taking place within its topological flat bands. Several numerical and analytical methods have been applied to understand the correlated phases therein, revealing some similarity with the quantum Hall physics. In this work, we provide a Mott-Hubbard perspective for the TBG system. Employing the large-scale density matrix renormalization group on the lattice model containing the projected Coulomb interactions only, we identify a first-order quantum phase transition between the insulating stripe phase and the quantum anomalous Hall state with the Chern number of ±1. Our results not only shed light on the mechanism of the quantum anomalous Hall state discovered at three-quarters filling, but also provide an example of the topological Mott insulator, i.e., the quantum anomalous Hall state in the strong coupling limit.

In the present paper, the authors applied the state of art DMRG method to study the phases in twisted bi-layer graphene (TBG) based on a tight binding model constructed by the Wannier functions derived from a single valley. Over all speaking, I agree this is a timely research on a very hot topic in condensed matter by a reliable numerical method, which has been developed by some of the authors for years. Overall speaking, I agree that their paper has captured some features of the TBG system and deserves to be published somewhere. However, I don't think the present paper meets the scientific standard of nature communication due to an obvious weak point, the questionable relevance to the real materials. It is no doubt that TBG is an itinerant system and its Wannier obitals for the flat bands actually extend to quite some range away from the center. In the theoretical model adopted by the present paper, only the interaction terms involving particles in the neighbouring super cells are considered in their DMRG calculation. Similar Wannier functions have also been discussed in detail in PHYSICAL REVIEW X 8, 031087 (2018). In table 1 of the above reference, the authors list the strength of the projected Coulomb interaction at different range, from which we can see that from on-site to the 5th neighbour the interaction strength only drops to about 1/3 of the on-site value. The above results suggest that in order to provide a faithful description for the correlation effects in TBG we definitely need to include further range interaction terms in the Wannier representation. Another problem of the strong coupling model discussed here is that the finite energy dispersion of the flat bands has been completely neglected, which I don't think is a good approximation. Comparing to the interaction strength U0 (about 40meV) considered in the present paper, the energy dispersion of the flat bands (which is about 8meV at magic angle) is small but definitely not irrelevant. Based on the above two points, I don't agree to accept this paper for publication in Nature communication.
Reviewer #3 (Remarks to the Author): In this work, the authors numerically and analytically study a local 'Hubbard-like' model on the hexagonal lattice which is derived from Wannier localizing the flat bands of twisted bilayer graphene. The Hamiltonian contains a single tuning parameter alpha, which is shown to drive a first order phase transition from a unidirectional charge density wave phase to a uniform quantum anomalous Hall state which spontaneously breaks time reversal symmetry. Within the field of moire materials the results presented in this work are of great interest because they show that even despite the fragile topology of the flat bands in the Bistritzer-MacDonald Hamiltonian, the electron interactions in the Wannier basis can be truncated to produce tractable real-space models which correctly capture the strong-coupling physics. The appearance of a quantum anomalous Hall phase in this model is also of general interest even after stripping away the connections to twisted bilayer graphene, because it is the first instance of a 'topological Mott insulator'. In his seminal paper, Haldane has shown that one does not need a magnetic field to have electrons occupy topologically non-trivial extended states which respond to Laughlin's adiabatic flux insertion by producing a quantized Hall current. The results in this work show that one does not even need a kinetic energy term in the Hamiltonian for this to occur.
The numerical results presented by the authors are detailed and convincing, and can in certain regimes also be matched to perturbative and mean-field calculations. This leaves little room for doubt about the correctness of the results. For this reason, I support publication of the manuscript.
There is perhaps only one statement in the paper that I would like to ask the authors to clarify. In particular, in the right column of page 3 it is stated that 'the NNN hopping amplitude is purely imaginary'. How is this hopping amplitude defined? And how do I see that it is indeed purely imaginary?

Nick Bultinck Title: Realization of Topological Mott Insulator in a Twisted Bilayer Graphene Lattice Model
Dear Reviewers, We thank you very much for the very positive assessment and insightful comments that have helped us to further improve our manuscript. We are more than happy to see that both Reviewers #1 and #3 have rather positive assessment of our work. They think our results constitute "an important addition to the growing literature" (Reviewer #1), and are "of general interest even after stripping away the connections to twisted bilayer graphene, because it is the first instance of a topological Mott insulator" (Reviewer #3). Besides, Reviewer #2 has raised some concern on the relevance to real materials. In the following, we give a pointby-point response to the comments, with the reviewers' text being cited in blue and our subsequent response in normal format. Note that all bibliographic citations below make direct reference to the bibliography in our manuscript.
Overall, by taking the reviewers' valuable comments and suggestions in to account, we have made careful revisions accordingly and believe our manuscript is now further improved.

Response to Reviewer #1
Reviewer #1: The author employ DMRG to study the lattice model describing twisted bilayer graphene introduced by some of the authors (Kang and Vafek). The main finding of the work is a phase diagram for the model in the single valley spineless limit at half-filling which is relevant to twisted bilayer graphene at filling -3 or +3. At the latter filling, experiments have already observed quantum anomalous Hall states as well as Chern zero states which is compatible with the two phases observed by the authors. The quantum anomalous Hall state has already been identified in several previous theoretical studies based on momentum space continuum model approaches including works by some of the authors (arxiv:2002.10360). The main novelty here is to be able to obtain this phase from a real space model. The paper is well written and represent an important addition to the growing literature on the topic. I recommend for publication provided the authors address a few comments that I detail below: Reply: We thank Reviewer #1 for his/her concise summary of our work and its relevance towards the recent developments of the field, we also thank Reviewer #1 for the support and the constructive comments.
Reviewer #1: 1. The authors make the claim that this is the first study to find a quantum anomalous Hall state. While this is likely true, a recent study by the same authors (Ref. 35) have already identified a strong coupling topological valley Hall state, which is a very close relative to the state identified here, using quantum Monte Carlo simulation of the same model at charge neutrality. It would be useful for the authors to clarify this point.
Reply: We thank the reviewer for the valuable comment. Indeed, the quantum Monte Carlo study by some of us revealed the quantum valley Hall (QVH) state in the real space model at charge neutrality ⌫ = 0 for a specific choice of kinetic energy hopping terms. There are several differences between the QVH and the quantum anomalous Hall (QAH) state discovered in the present work. First, the two states are realized in different physical regimes. The QVH phase (Ref. [52] now in the updated reference list) is found in the intermediate coupling regime, and turns into the intervalley-coherence state (IVC) at the strong coupling limit of the model, so the QVH is more of theoretical interest than the realistic TBG experiment. In contrast, the QAH state discovered here is in the strong coupling limit at the filling of ±3, with both the spin and valley polarized. It should be closer to the real TBG system than that of the QVH. Second, these states are obtained in the strong coupling regime at very different fillings. As arXiv:2009.13530 illustrates, the QAH occurs at the odd fillings, and the ground state at the even fillings is usually the IVC state with zero Chern number.
We have incorporated these clarification into the revised manuscript as Reviewer #1 suggested.
Reviewer #1: 2. The paper will benefit from relating the model parameters to physical quantities, even on a qualitative level. For instance, most samples which observe the quantum anomalous Hall state are aligned with the hexagonal boron nitride (hBN) substrate. This induces a sublattice potential which breaks C2 symmetry. The lattice model studied by the authors was derived for the pristine case without hBN alignment. What would change if this effect was included? is the alpha parameter going to increase or decrease? Similarly, it is unclear what physical quantities control the parameter alpha. Is it sensitive to lattice relaxation? to interaction screening?
Reply: We thank the reviewer for the insightful comment. Since the hBN alignment does not close the band gap between the flat bands and the remote bands, the Hilbert space of the flat bands are almost intact. Therefore, the Wannier basis and the associated projected Coulomb interaction are not affected by the alignment. However, the alignment breaks the C 2 symmetry and thus introduces the additional symmetry breaking terms in the Hamiltonian. Although such terms are small compared with the interaction and thus the system is still in the strong coupling regime, they favor the QAH state with a particular combination of the Chern number and the valley.
The parameter ↵ controls the relative strength between the cluster charging and assisted hopping interactions. As mentioned above, the projected interaction including the value of ↵ is not changed by the hBN alignment. As explained in Ref. [26], ↵ is originated from the overlaps of two neighboring Wannier states, and therefore, it can be calculated by the Bistritzer-MacDonald model whose parameters vary with the lattice relaxation. More detailed discussions on the real space interaction model and the associated ↵ will be presented in a separate work.
The form of the interaction in Eq. (1) is obtained when the screening distance is comparable with the moiré lattice constant, and thus contains only the cluster charging and the assisted nearest neighbor hopping terms. With the screening distance larger than moiré lattice constant, the interactions would include more terms for the sizable longer-range repulsion, but experimentally, barely change the qualitative properties of the insulating phases. Motivated by this observation, we argue our Hamiltonian in Eq.
(1) provides a minimal model for the explanation of the interaction dominated physics and the emergence of TMI.
We have incorporate these suggestions of the reviewer into the revised manuscript.
Reviewer #1: 3. The QAH spontaneously breaks two-fold rotation symmetry which is an emergent symmetry of twisted bilayer graphene. This symmetry is however implemented non-locally in the lattice model studied which means that any finite system will already introduce some explicit symmetry breaking that may slightly alter the energy competition by introducing a bias to C2-breaking states. While I do not expect any qualitative change to the phase diagram, I think the authors should comment on the issue and present some quantitative estimates for such explicit C2 breaking due to finite size effects.
Reply: While the C 2 symmetry is broken by the truncation over the interaction range, the Hamiltonian in Eq. (1) is invariant under the combination of the chiral particle-hole, the unitary particle-hole, and C 2 T symmetries. This combined symmetry can be locally implemented and belong to the U (4)⇥U (4) symmetry group in the chiral limit. Since this symmetry transforms the two QAH states with opposite Chern numbers into each other, it has to be broken for the system to develop a QAH state. Therefore, our model does not favor a priori the QAH state even if the C 2 symmetry is broken by truncating the interaction range. More detailed discussion on the real space model and the associated symmetries will be presented elsewhere in a separate work.

Response to Reviewer #2
Reviewer #2: In the present paper, the authors applied the state of art DMRG method to study the phases in twisted bi-layer graphene (TBG) based on a tight binding model constructed by the Wannier functions derived from a single valley. Over all speaking, I agree this is a timely research on a very hot topic in condensed matter by a reliable numerical method, which has been developed by some of the authors for years.
Reply: We thank the reviewer for the nice evaluation of our work, especially for his/her acknowledgement of our efforts in developing reliable lattice model of TBG and unbiased large-scale numerical methodologies in addressing them, over the years.
Reviewer #2: However, I don't think the present paper meets the scientific standard of nature communication due to an obvious weak point, the questionable relevance to the real materials. It is no doubt that TBG is an itinerant system and its Wannier obitals for the flat bands actually extend to quite some range away from the center. In the theoretical model adopted by the present paper, only the interaction terms involving particles in the neighbouring super cells are considered in their DMRG calculation. Similar Wannier functions have also been discussed in detail in PHYSICAL REVIEW X 8, 031087 (2018). In table 1 of the above reference, the authors list the strength of the projected Coulomb interaction at different range, from which we can see that from on-site to the 5th neighbour the interaction strength only drops to about 1/3 of the on-site value. The above results suggest that in order to provide a faithful description for the correlation effects in TBG we definitely need to include further range interaction terms in the Wannier representation.
Reply: We thank the reviewer for the insightful comment and would need to go several different angles to address it.
First, we would like to point out that the reference PRX 8, 031087 (2018) (the reference [24] in our revised manuscript) has assumed that the Coulomb interaction is / 1/r and projected it to the basis of the constructed Wannier states. Experimentally, however, the interaction is always suppressed at the long distance by the metallic gates. As a consequence, the long-range density-density interaction exponentially decays at the distance r exceeding the distance to the gates. Numerically, we have already found in the reference PRL 122, 246401 (2019) (Ref. [26] in our revised manuscript) that the projected interaction can be well approximated by Eq. (1) when the distance between double gates is about the order of the moiré lattice constant and this is case for the model in our present study.
Secondly, on the other hand, we completely agree with Reviewer #2 that in the TBG system the Wannier orbitals for the flat bands extend to quite some distance, and that is why over the past several years, we have been working very actively in developing reliable lattice model beyond the conventional on-site interaction and more importantly, developing unbiased numerical methodologies to gradually include more extended interactions into the computations. As the reviewer is well aware of, what have been presented in this work, with the third nearest neighbor interaction at the same strength of the on-site interaction, is already at the best of the present computational capability. While we are still developing better numerical methods that could take even longer range interaction into account in an unbiased manner, we would like to quote the words from Reviewer #3 to emphasize such an understanding, that, "Within the field of moire materials the results presented in this work are of great interest because they show that even despite the fragile topology of the flat bands in the Bistritzer-MacDonald Hamiltonian, the electron interactions in the Wannier basis can be truncated to produce tractable real-space models which correctly capture the strong-coupling physics." Thirdly, this work is primarily theoretical in that, although we tried to be qualitatively consistent with the TBG systems at ⌫ = ±3 fillings, our scope is actually beyond the TBG. We find our real-space model with interaction-only can give rise to the long-pursued topological Mott insulator state (i.e., the QAH state in the present case). This is certainly an important theoretical discovery in its own right, as such state has been investigated by the communities for more than a decade without success. We would again quote the text of Reviewer #3 to stress the point "The appearance of a quantum anomalous Hall phase in this model is also of general interest even after stripping away the connections to twisted bilayer graphene, because it is the first instance of a 'topological Mott insulator'. In his seminal paper, Haldane has shown that one does not need a magnetic field to have electrons occupy topologically non-trivial extended states which respond to Laughlin's adiabatic flux insertion by producing a quantized Hall current. The results in this work show that one does not even need a kinetic energy term in the Hamiltonian for this to occur." Reviewer #2: Another problem of the strong coupling model discussed here is that the finite energy dispersion of the flat bands has been completely neglected, which I don't think is a good approximation. Comparing to the interaction strength U0 (about 40meV) considered in the present paper, the energy dispersion of the flat bands (which is about 8meV at magic angle) is small but definitely not irrelevant.
Reply: We thank the reviewer for the comment. We believe the system is in the strong coupling regime for the following two reasons. First, the original bandwidth W ⇠ 8 meV is much smaller than the interaction U 0 , and thus suggest the system is in the strong coupling regime. Second, in a recent work (Ref. [61]) by some of the authors, the wavefunction and the associated dispersion of the flat bands have been carefully calculated by integrating out the states on the remote bands. The corresponding superexchange interaction is found to be . 10 2 e 2 /(✏L m ), that is tiny compared with U 0 and thus provides strong argument to neglect the kinetic terms. We have added such discussion in the revised manuscript.
Reviewer #2: Based on the above two points, I don't agree to accept this paper for publication in Nature communication.
Reply: With the detailed responses above and corresponding revisions in the manuscript, we humbly believe that our work both provides a microscopic mechanism of the observed QAH and, even more intriguingly, also goes beyond TBG by realizing the topological Mott insulator. We therefore sincerely hope Reviewer #2 can agree with us that it meets the standard of Nature Communications.

Response to Reviewer #3
Reviewer #3: In this work, the authors numerically and analytically study a local 'Hubbard-like' model on the hexagonal lattice which is derived from Wannier localizing the flat bands of twisted bilayer graphene. The Hamiltonian contains a single tuning parameter alpha, which is shown to drive a first order phase transition from a unidirectional charge density wave phase to a uniform quantum anomalous Hall state which spontaneously breaks time reversal symmetry. Within the field of moire materials the results presented in this work are of great interest because they show that even despite the fragile topology of the flat bands in the Bistritzer-MacDonald Hamiltonian, the electron interactions in the Wannier basis can be truncated to produce tractable real-space models which correctly capture the strong-coupling physics.
Reply: We thank the respected reviewer for his concise summary of our work and its relevance with the field.
Reviewer #3: The appearance of a quantum anomalous Hall phase in this model is also of general interest even after stripping away the connections to twisted bilayer graphene, because it is the first instance of a 'topological Mott insulator'. In his seminal paper, Haldane has shown that one does not need a magnetic field to have electrons occupy topologically non-trivial extended states which respond to Laughlin's adiabatic flux insertion by producing a quantized Hall current. The results in this work show that one does not even need a kinetic energy term in the Hamiltonian for this to occur.
Reply: Thanks for this insightful comment and deep understanding about the value of our work and the field of topological quantum matter.
Reviewer #3: The numerical results presented by the authors are detailed and convincing, and can in certain regimes also be matched to perturbative and mean-field calculations. This leaves little room for doubt about the correctness of the results. For this reason, I support publication of the manuscript.
Reply: We thank Reviewer #3 for the firm support for the publication of our work.
Reviewer #3: There is perhaps only one statement in the paper that I would like to ask the authors to clarify. In particular, in the right column of page 3 it is stated that 'the NNN hopping amplitude is purely imaginary'. How is this hopping amplitude defined? And how do I see that it is indeed purely imaginary?
Reply: Thanks for the careful reading, and indeed we did not give the definition clearly enough. What we have calculated with the DMRG is not the NNN hopping constant but the equal time NNN correlation. On the other hand, the hoppings of the approximated tight binding model are obtained based on the variational analysis, as shown in the SM. We have added corresponding description in the revised manuscript, with the precise definition of NNN hopping now provided.

Summary of changes
All major changes are marked in blue color in the resubmitted text. Amongst them, main revisions are as follows.
• To respond the suggestions from Reviewers #1 and #2, on Page 2 we clarified the physical origin of the parameters in our model and discussed the irrelevance of the energy dispersion of the flat bands towards the strong coupling results revealed in this work, as well as the spontaneous nature of our QAH state without the need of including complex hopping terms to take the effect of hBN alignment into consideration.
• To respond the suggestion of Reviewer #3, on Page 3 (right column) we further explained the equaltime correlator hc † l c l 0 i that are found to be purely imaginary in the QAH state, and thus defined its imaginary part hJi = i 2 h(c † l c l 0 c † l 0 c l )i as the order parameter of QAH phase.
• To respond the suggestion of Reviewer #1, on Page 4 (left column) we explained the difference between the QAH state discovered in this work and the QVH state discovered by some of us in a recent work.
The authors have satisfactorily addressed all my concerns. I have only a minor issue regarding a statement in the revised manuscript "In addition, we do not include the additional small symmetry-breaking term produced by the possible hBN alignment that favors the QAH phase with a particular Chern number". The estimate for this term is between 15-20 meV (1901.08110) which is around half the value of U0 used in this work and can be larger if a different value of the screening epsilon is used, so it is not clear it can be dismissed as a small contribution. Maybe the authors should qualitatively discuss what is the likely effect of this term on the phase diagram they obtain.