# PFT-SSR: Parallax Fusion Transformer for Stereo Image Super-Resolution

###### Abstract

Stereo image super-resolution aims to boost the performance of image super-
resolution by exploiting the supplementary information provided by binocular
systems. Although previous methods have achieved promising results, they did
not fully utilize the information of cross-view and intra-view. To further
unleash the potential of binocular images, in this letter, we propose a novel
Transformer-based parallax fusion module called Parallax Fusion Transformer
(PFT). PFT employs a Cross-view Fusion Transformer (CVFT) to utilize cross-
view information and an Intra-view Refinement Transformer (IVRT) for intra-
view feature refinement. Meanwhile, we adopted the Swin Transformer as the
backbone for feature extraction and SR reconstruction to form a pure
Transformer architecture called PFT-SSR. Extensive experiments and ablation
studies show that PFT-SSR achieves competitive results and outperforms most
SOTA methods. All code will be available.

Index Terms—  Stereo Image Super-Resolution, Parallax Fusion Transformer,
Stereo Cross Attention, SSR.

##  1 Introduction

00footnotetext: Corresponding author: Juncheng Li, Tieyong Zeng00footnotetext:
This work was supported in part by the National Key Research and Development
Program of China under Project no. 2021YFE0203700, 2018AAA0100102, and
2018AAA0100100, the National Natural Science Foundation of China under Grant
no. 61972212, and the YangFan Project of Shanghai under Grant no. 23YF1412800.

Binocular cameras have been widely employed to improve the perception
capabilities of vision systems in devices such as self-driving vehicles and
smartphones. With the rapid development of binocular cameras, stereo image
super-resolution (SSR) is becoming increasingly popular in academia and
industry. Specifically, SSR attempts to reconstruct a high-resolution (HR)
image from a pair of low-resolution (LR) images. With the help of additional
information from a pair of binocular images at the same physical location,
making full use of the information from both two images is crucial for stereo
image super-resolution (SSR).

The easiest way to implement stereo image SR is to perform single image SR
(SISR) methods [1, 2, 3, 4, 5, 6] on stereo image pairs, respectively. These
approaches, however, neglected the cross-view information between the pair of
images and are incapable of reconstructing high-quality images. To address
this problem, current strategies have focused on building novel cross-view
feature aggregation modules, loss functions, and so on, to improve the
efficiency with which image pair interaction features are used. For example,
[7] first combined depth estimation and image resolution tasks with multiple
image inputs. After that, StereoSR [8] took the lead in introducing CNN into
Stereo SR. iPASSR [9] suggested a symmetric bi-directional parallax attention
module (biPAM) and an inline occlusion handling scheme as its cross view
interaction module to exploit symmetry cues for stereo image SR. Recently,
several more advanced strategies for improving Stereo SR performance have been
introduced. For instance, NAFSSR [10] designed a new CNN-based backbone NAFNet
[11] and proposed a novel Stereo Cross Attention Module (SCAM) as parallax
fusion block. These network topologies typically included a CNN backbone for
obtaining intra-view information and a parallax fusion module for combining
cross-view attention. Since the existence of parallax, we discovered that it
is also highly crucial for cross-picture features and intra-picture features
to promote each other in the process of binocular feature fusion. However,
these two processes in existing works are often relatively independent, which
is not conducive to the full use of image features. Meanwhile, the quality of
the input features is vital for image fusion efficiency. However, existing
works never consider the degree of match between the backbone networks and
parallax fusion blocks. Therefore, the combination of these two pieces will be
sub-optimal.

Fig. 1: The complete architecture of the proposed Parallax Fusion Transformer
for Stereo Image Super-Resolution (PFT-SSR). This is a dual-stream network and
interacts through an interaction module. Due to page limit, please zoom in to
see details.

In this work, we address the aforementioned problems by introducing the
Transformer to stereo image SR. Recently, Transformer demonstrated strong
performance in various low-level tasks [12, 13, 14], which can learn global
information of images to further improve model performance. However, directly
merging current CNN-based parallax fusion modules (PFM) and Transformer will
not result in outstanding performance. This is because the CNN-based parallax
fusion modules and Transformer have different properties, leading to PFM that
cannot fully utilize the features from the Transformer backbone. To address
this issue, we designed a new parallax fusion module, named Parallax Fusion
Transformer (PFT). PFT contains a Stereo Cross Attention Module (SCAM) and a
Feature Refining Module (FRM). Among them, the SCAM gets the cross-view
attention and FRM will fuse the cross-view feature with the local window
features. The cross-view features and intra-view features (local window
features) will enhance each other to get a better representation for image
super-resolution task. With the help of PFT, the proposed model can well-adapt
the deep features with the parallax feature fusion blocks to fully utilize the
representational potential of the Transformer.

The contributions of this letter can be summarized as follows: 1) We propose a
novel Parallax Fusion Transformer (PFT) layer with a Cross-view Fusion
Transformer (CVFT) and an Intra-view Refinement Transformer (IVRT). 2) Based
on the proposed PFT, we design a pure Transformer network (named PFT-SSR) to
further improve the feature extraction ability of Transformer-based networks.
3) Extensive experiments have illustrated the effectiveness of PFT-SSR.

##  2 METHODOLOGY

In this paper, we propose a Parallax Fusion Transformer for Stereo Image
Super-Resolution, called PFT-SSR. As shown in Fig. 1, the proposed PFT-SSR
consists of three parts: stereo feature extraction, feature interaction, and
SR image reconstruction. For Stereo SR, the model takes two images
$x_{LR}^{L}$, $x_{LR}^{R}\in R^{B\times C_{in}\times H\times W}$ as inputs and
then outputs $x_{HR}^{L}$, $x_{HR}^{R}\in R^{B\times C_{out}\times S*H\times
S*W}$. Among them, $B$, $C_{in}$, $C_{out}$, $H$, and $W$ are the input batch
size, the number of channels, height, and weight, respectively. Meanwhile, $S$
is the upscaling factor, which is used to control the size of the output
images. Specifically, we first use two convolutional layers to extract shallow
features of the input images respectively. After that, we further extract the
deeper feature representations with SwinIR [13] backbone, which contains three
consequent Residual Swin Transformer Blocks (RSTBs)

Fig. 2: Visual results (x4) achieved by different methods on Flickr1024.

$\begin{split}I_{d}^{L}=f_{ex}(f_{s}(x_{LR}^{L})),\quad
I_{d}^{R}=f_{ex}(f_{s}(x_{LR}^{R})).\end{split}$ (1)

Then, the extracted features are fed into the proposed Parallax Fusion
Transformers (PFT) for cross-view interaction and intra-view refinement

$I_{f}^{L},I_{f}^{R}=f_{PFT}(I_{d}^{R},I_{d}^{L}).$ (2)

With fused features, we apply RSTBs again to obtain the refined features, with
a residual connection from the shallow image feature (ignored in formula for
simplicity).

$\begin{split}I_{r}^{L}=f_{cov}(f_{re}(I_{f}^{L})),\quad
I_{r}^{R}=f_{cov}(f_{re}(I_{f}^{R})).\end{split}$ (3)

Finally, a Reconstruction module that contains a single convolutional layer
and a PixelShuffle layer is used to reconstruct the final SR images.

Fig. 3: The architecture of stereo cross attention module.

###  2.1 Swin Transformer Backbone

In this work, we use Swin Transformer Blocks [15] to build the backbone of our
network. Specifically, a Swin Transformer Layer firstly reshapes the input
feature map $I_{in}$ to $\frac{HW}{M^{2}}\times M^{2}\times C$ and performs
standard self-attention locally on each window. For each of $\frac{HW}{M^{2}}$
feature maps, let input be $X\in{R}^{M^{2}\times C}$, then query, key, and
value should be

$Q=XP_{Q},\quad K=XP_{K},\quad V=XP_{V},$ (4)

where $P_{Q}$, $P_{K}$, and $P_{V}$ are linear projection matrices. Then, the
attention matrix is calculated within the local windows

$Attention(Q,K,V)=SoftMax(QK^{T}/\sqrt{d}+B)V,$ (5)

where $B$ is the positional encoding for Transformer. The model also apply an
MLP with two fully connected layers and GELU non-linearity on the attention
matrix for feature transformations. Meanwhile, the LayerNorm [16] layer is
added before both Attention Block and MLP with residual connection. Though
local attention can greatly reduce the amount of computation, there is no
connection across local windows. To solve this problem, Swin Transformer
proposed a shifted window mechanism to shifts the feature map by
$(\lfloor\frac{M}{2}\rfloor,\lfloor\frac{M}{2}\rfloor)$ pixels before
partitioning. The process can be expressed as

$\begin{split}X=MSA(LN(X))+X,X=MLP(LN(X))+X\end{split}$ (6)

where regular partitioning and shift partitioning are used alternately before
each MSA. With the help of this backbone, our model can extract sufficient
useful image features.

###  2.2 Parallax Fusion Transformer

In order to make full use of the features of the left and right images, we
propose a Parallax Fusion Transformer (PFT). As shown in Fig. 1, PFT contains
4 PFT blocks, and each PFT block consists of 6 PFT layers and a convolutional
layer. Meanwhile, each PFT layer has two different Transformer blocks, i.e.
Cross-view Fusion Transformer (CVFT) and Intra-view Refinement Transformer
(IVRT). Among them, CVFT adopts stereo cross-attention module (SCAM [10]) to
learn the features of another view and IVRT takes the local-window Transformer
to better merge features from the other view to its feature map. Specifically,
we first apply CVFT to achieve cross-view attention via SCAM. However, using
single-head SCAM to get the cross-view information cannot adapt to different
parallax. Therefore, we further use IVRT to make cross-view information from
the other branch better interact with intra-view features. With this
’Attention-Refine’ paradigm, our PFT-SSR shows a compelling effect on cross-
view attention.

Table 1: Quantitative comparison on different datasets. PSNR$/$SSIM values
achieved on both the left images (i.e., Left) and a pair of stereo images
(i.e., $\left(\textit{Left}+\textit{Right}\right)/2$) are reported. Among
them, the best results are highlighted.

Method | Scale | Left | $\left(\textit{Left}+\textit{Right}\right)/2$
---|---|---|---
KITTI 2012 | KITTI 2015 | Middlebury | KITTI 2012 | KITTI 2015 | Middlebury | Flickr1024
EDSR [3] |  $\times$2 |  30.83$/$0.9199 |  29.94$/$0.9231 |  34.84$/$0.9489 |  30.96$/$0.9228 |  30.73$/$0.9335 |  34.95$/$0.9492 |  28.66$/$0.9087
RCAN [17] |  $\times$2 |  30.88$/$0.9202 |  29.97$/$0.9231 |  34.80$/$0.9482 |  31.02$/$0.9232 |  30.77$/$0.9336 |  34.90$/$0.9486 |  28.63$/$0.9082
StereoSR [18] |  $\times$2 |  29.42$/$0.9040 |  28.53$/$0.9038 |  33.15$/$0.9343 |  29.51$/$0.9073 |  29.33$/$0.9168 |  33.23$/$0.9348 |  25.96$/$0.8599
PASSRnet [19] |  $\times$2 |  30.68$/$0.9159 |  29.81$/$0.9191 |  34.13$/$0.9421 |  30.81$/$0.9190 |  30.60$/$0.9300 |  34.23$/$0.9422 |  28.38$/$0.9038
iPASSR [9] |  $\times$2 |  30.97$/$0.9210 |  30.01$/$0.9234 |  34.41$/$0.9454 |  31.11$/$0.9240 |  30.81$/$0.9340 |  34.51$/$0.9454 |  28.60$/$0.9097
SSRDE-FNet [20] |  $\times$2 |  31.08$/$0.9224 |  30.10$/$0.9245 |  35.02$/$0.9508 |  31.23$/$0.9254 |  30.90$/$0.9352 |  35.09$/$0.9511 |  28.85$/$0.9132
PFT-SSR (Ous) |  $\times$2 |  31.15$/$0.9166 |  30.16$/$0.9187 | 35.08$/$0.9516 |  31.29$/$0.9195 |  30.96$/$0.9306 | 35.21$/$0.9520 |  29.05$/$0.9049
EDSR [3] |  $\times$4 |  26.26$/$0.7954 |  25.38$/$0.7811 |  29.15$/$0.8383 |  26.35$/$0.8015 |  26.04$/$0.8039 |  29.23$/$0.8397 |  23.46$/$0.7285
RCAN [17] |  $\times$4 |  26.36$/$0.7968 |  25.53$/$0.7836 |  29.20$/$0.8381 |  26.44$/$0.8029 |  26.22$/$0.8068 |  29.30$/$0.8397 |  23.48$/$0.7286
StereoSR [18] |  $\times$4 |  24.49$/$0.7502 |  23.67$/$0.7273 |  27.70$/$0.8036 |  24.53$/$0.7555 |  24.21$/$0.7511 |  27.64$/$0.8022 |  21.70$/$0.6460
PASSRnet [19] |  $\times$4 |  26.26$/$0.7919 |  25.41$/$0.7772 |  28.61$/$0.8232 |  26.34$/$0.7981 |  26.08$/$0.8002 |  28.72$/$0.8236 |  23.31$/$0.7195
SRRes+SAM |  $\times$4 |  26.35$/$0.7957 |  25.55$/$0.7825 |  28.76$/$0.8287 |  26.44$/$0.8018 |  26.22$/$0.8054 |  28.83$/$0.8290 |  23.27$/$0.7233
iPASSR [9] |  $\times$4 |  26.47$/$0.7993 |  25.61$/$0.7850 |  29.07$/$0.8363 |  26.56$/$0.8053 |  26.32$/$0.8084 |  29.16$/$0.8367 |  23.44$/$0.7287
SSRDE-FNet [20] |  $\times$4 |  26.61$/$0.8028 |  25.74$/$0.7884 |  29.29$/$0.8407 |  26.70$/$0.8082 |  26.45$/$0.8118 |  29.38$/$0.8411 |  23.59$/$0.7352
PFT-SSR (Ours) |  $\times$4 |  26.64$/$0.7913 |  25.76$/$0.7775 | 29.58$/$0.8418 |  26.77$/$0.7998 |  26.54$/$0.8083 | 29.74$/$0.8426 |  23.89$/$0.7277

Cross-view Fusion Transformer (CVFT): The core component of CVFT is SCAM, and
the whole process of SCAM is shown in Fig. 3. Given input image features
$X_{L},X_{R}\in{R}^{H\times W\times C}$, we first perform layer normalization
to get scaled features. Due to the nature of stereo images, we use the same
$Q$ and $K$ for representing intra-view features. Then, we get cross-view
attention both from right to left and from left to right by

$\begin{split}F_{R\rightarrow
L}=Attention(T_{1}^{L}\overline{X_{L}},T_{1}^{R}\overline{X_{R}},T_{2}^{R}\overline{X_{R}}),\\\
F_{L\rightarrow
R}=Attention(T_{1}^{R}\overline{X_{R}},T_{1}^{L}\overline{X_{L}},T_{2}^{L}\overline{X_{L}}),\end{split}$
(7)

where $Attention$ is defined same as Eq. (5). Besides, $T_{1}^{L}$,
$T_{1}^{R}$, $T_{2}^{L}$, and $T_{2}^{R}$ are linear projection matrices.
After getting the cross-view attention feature, we use a weighted residual
connection to merge it to the corresponding image feature, which are
formulated as

$\begin{split}Y_{L}=\alpha_{L}F_{R\rightarrow
L}+X_{L},Y_{R}=\alpha_{R}F_{L\rightarrow R}+X_{R},\end{split}$ (8)

where $\alpha_{L}$ and $\alpha_{R}$ are learnable scalars. After observing the
corrected features, we apply MLP and LayerNorm to get the final outputs and
the whole process can be expressed as

$\begin{split}X=SCAM(X)+X,X=MLP(LN(X))+X.\end{split}$ (9)

Intra-view Refinement Transformer (IVRT): One key difficulty of SSR is the
different parallax brought by various stereo systems. Although SCAM shows
great cross-view attention ability, it cannot adapt various parallax. After
observing this, we used a Transformer with local-window attention for feature
refinement. Regular partitioning is adopted before the MSA so that the
features after the interaction of the two views can be further fused and
enhanced, which is helpful for the final SR image reconstruction.

##  3 Experiment

###  3.1 Experimental Settings

800 images from Flickr1024[21] and 60 images from Middlebury[22] are chosen
for training. To make the Middlebury dataset matches the spatial resolution of
the Flickr1024 dataset, we perform bicubic downsampling by a factor of 2 on
each image. And then, we use bicubic downsampling to these GT images by the
factors of 2 and 4 to get the input images. We follow previous works [9, 10,
20] on this setting to make comparison fair. During training, we use the L1
loss function for supervision, PSNR and SSIM as quantitative metrics to make
easy comparison with previous methods. These metrics are calculated on RGB
color space with a pair of stereo images. To evaluate SR results, we use KITTI
2012 [23], KITTI 2015 [24], Middlebury [22], and Flickr1024 [21] for test.

###  3.2 Comparison to state-of-the-art methods

We compare our proposed PFT-SSR with several state-of-the-art methods,
including SISR methods (e.g., EDSR [3], RCAN [17]) and stereo image SR methods
(e.g., StereoSR [18], PASSRnet [19], iPASSR [9], and SSRDE-FNett [20]).
According to TABLE 1, we can clearly observe that our PFT-SSR achieves
outstanding results and outperforms most other SOTA methods, especially on
Flickr102. Meanwhile, we also show the qualitative comparisons in Figs. 2.
Obviously, our PFT-SSR can reconstruct more accurate SR images with more
accurate edges and texture details. This fully demonstrates the effectiveness
of the proposed PFT-SSR.

###  3.3 Ablation Study

Cross-view interaction is the key part in Stereo SR. In this part, we do
ablation on the choice of this technology to show the strong stereo image
fusion ability of the proposed PFT. We use Swin Transformer [15] Blocks as
backbones and take the same number of Swin Transformer, biPAM [9], and our
proposed PFT as the cross-view interaction module in this part. According to
TABLE 2, it is obviously that the proposed PFT can improve the model
performance more effectively, which fully illustrates the effectiveness of
PFT.

Table 2: Ablation study on PFT under Flickr1024. Backbone | Module | PSNR (x4) | SSIM (x4)
---|---|---|---
Swin Transformer | None | 23.54 | 0.7120
Swin Transformer | RSTB (SwinIR) | 23.65 | 0.7164
Swin Transformer | BiPAM | 23.42 | 0.7068
Swin Transformer | PFT (Ours) | 23.83 | 0.7268

##  4 Conclusion

In this paper, we proposed a PFT-SSR for stereo image super-resolution, which
contains a well-designed Parallax Fusion Transformer (PFT). PFT consists of a
Cross-view Fusion Transformer (CVFT) and an Intra-view Refinement Transformer
(IVRT), specially designed for cross-view interaction. It is worth mentioning
that PFT can better merge different parallaxes to utilize the features of the
left and right images fully. Meanwhile, PFT can also better adapt to the
current popular Transformer-based backbone. Extensive experiments show that
PFT-SSR outperforms most current models and achieves promising outcomes.

## References

  * [1] Chao Dong, Chen Change Loy, Kaiming He, and Xiaoou Tang.  Learning a deep convolutional network for image super-resolution.  In ECCV, pages 184–199, 2014.
  * [2] Jiwon Kim, Jung Kwon Lee, and Kyoung Mu Lee.  Accurate image super-resolution using very deep convolutional networks.  In CVPR, pages 1646–1654, 2016.
  * [3] Bee Lim, Sanghyun Son, Heewon Kim, Seungjun Nah, and Kyoung Mu Lee.  Enhanced deep residual networks for single image super-resolution.  In CVPR Workshops, pages 136–144, 2017.
  * [4] Tao Dai, Jianrui Cai, Yongbing Zhang, Shu-Tao Xia, and Lei Zhang.  Second-order attention network for single image super-resolution.  In CVPR, pages 11065–11074, 2019.
  * [5] Yukai Shi, Haoyu Zhong, Zhijing Yang, Xiaojun Yang, and Liang Lin.  Ddet: Dual-path dynamic enhancement network for real-world image super-resolution.  IEEE Signal Processing Letters, 27:481–485, 2020.
  * [6] Yongsong Huang, Zetao Jiang, Rushi Lan, Shaoqin Zhang, and Kui Pi.  Infrared image super-resolution via transfer learning and psrgan.  IEEE Signal Processing Letters, 28:982–986, 2021.
  * [7] Arnav V Bhavsar and AN Rajagopalan.  Resolution enhancement in multi-image stereo.  IEEE Transactions on Pattern Analysis and Machine Intelligence, 32(9):1721–1728, 2010.
  * [8] Longguang Wang, Yulan Guo, Yingqian Wang, Juncheng Li, Shuhang Gu, Radu Timofte, Liangyu Chen, Xiaojie Chu, Wenqing Yu, Kai Jin, et al.  NTIRE 2022 challenge on stereo image super-resolution: Methods and results.  In CVPR Workshop, pages 906–919, 2022.
  * [9] Yingqian Wang, Xinyi Ying, Longguang Wang, Jungang Yang, Wei An, and Yulan Guo.  Symmetric parallax attention for stereo image super-resolution.  In CVPR Workshop, pages 766–775, 2021.
  * [10] Xiaojie Chu, Liangyu Chen, and Wenqing Yu.  Nafssr: Stereo image super-resolution using nafnet.  In CVPR, pages 1239–1248, 2022.
  * [11] Liangyu Chen, Xiaojie Chu, Xiangyu Zhang, and Jian Sun.  Simple baselines for image restoration.  arXiv preprint arXiv:2204.04676, 2022.
  * [12] Alexey Dosovitskiy, Lucas Beyer, Alexander Kolesnikov, Dirk Weissenborn, Xiaohua Zhai, Thomas Unterthiner, Mostafa Dehghani, Matthias Minderer, Georg Heigold, Sylvain Gelly, et al.  An image is worth 16x16 words: Transformers for image recognition at scale.  arXiv preprint arXiv:2010.11929, 2020.
  * [13] Jingyun Liang, Jiezhang Cao, Guolei Sun, Kai Zhang, Luc Van Gool, and Radu Timofte.  Swinir: Image restoration using swin transformer.  In ICCV, pages 1833–1844, 2021.
  * [14] Zhisheng Lu, Juncheng Li, Hong Liu, Chaoyan Huang, Linlin Zhang, and Tieyong Zeng.  Transformer for single image super-resolution.  In CVPR, pages 457–466, 2022.
  * [15] Ze Liu, Yutong Lin, Yue Cao, Han Hu, Yixuan Wei, Zheng Zhang, Stephen Lin, and Baining Guo.  Swin transformer: Hierarchical vision transformer using shifted windows.  In ICCV, pages 10012–10022, 2021.
  * [16] Jimmy Lei Ba, Jamie Ryan Kiros, and Geoffrey E Hinton.  Layer normalization.  arXiv preprint arXiv:1607.06450, 2016.
  * [17] Yulun Zhang, Kunpeng Li, Kai Li, Lichen Wang, Bineng Zhong, and Yun Fu.  Image super-resolution using very deep residual channel attention networks.  In ECCV, 2018.
  * [18] Daniel S Jeon, Seung-Hwan Baek, Inchang Choi, and Min H Kim.  Enhancing the spatial resolution of stereo images using a parallax prior.  In CVPR, pages 1721–1730, 2018.
  * [19] Longguang Wang, Yingqian Wang, Zhengfa Liang, Zaiping Lin, Jungang Yang, Wei An, and Yulan Guo.  Learning parallax attention for stereo image super-resolution.  In CVPR, pages 12250–12259, 2019.
  * [20] Qinyan Dai, Juncheng Li, Qiaosi Yi, Faming Fang, and Guixu Zhang.  Feedback network for mutually boosted stereo image super-resolution and disparity estimation.  In ACMMM, pages 1985–1993, 2021.
  * [21] Yingqian Wang, Longguang Wang, Jungang Yang, Wei An, and Yulan Guo.  Flickr1024: A large-scale dataset for stereo image super-resolution.  In ICCV Workshops, 2019.
  * [22] Daniel Scharstein, Heiko Hirschmüller, York Kitajima, Greg Krathwohl, Nera Nešić, Xi Wang, and Porter Westling.  High-resolution stereo datasets with subpixel-accurate ground truth.  In GCPR, pages 31–42, 2014.
  * [23] Andreas Geiger, Philip Lenz, and Raquel Urtasun.  Are we ready for autonomous driving? the kitti vision benchmark suite.  In CVPR, pages 3354–3361, 2012.
  * [24] Moritz Menze and Andreas Geiger.  Object scene flow for autonomous vehicles.  In CVPR, pages 3061–3070, 2015.


# Correlation effects in magic-angle twisted bilayer graphene:
An auxiliary-field quantum Monte Carlo study

Zhi-Yu Xiao  Department of Physics, College of William & Mary, Williamsburg,
Virginia 23187, USA  Shiwei Zhang  Center for Computational Quantum Physics,
Flatiron Institute, New York, NY 10010, USA

###### Abstract

Magic angle twisted bilayer graphene (MATBG) presents a fascinating platform
for investigating the effects of electron interactions in topological flat
bands. The Bistritzer-MacDonald (BM) model provides a simplified quantitative
description of the flat bands. Introducing long-range Coulomb interactions
leads to an interacting BM (IBM) Hamiltonian, a momentum-space continuum
description which offers a very natural starting point for many-body studies
of MATBG. Accurate and reliable many-body computations in the IBM model are
challenging, however, and have been limited mostly to special fillings, or
smaller lattice sizes. We employ state-of-the-art auxiliary-field quantum
Monte Carlo (AFQMC) method to study the IBM model, which constrains the sign
problem to enable accurate treatment of large system sizes. We determine
ground-state properties and quantify errors compared to mean-field theory
calculations. Our calculations identify correlated metal states and their
competition with the insulating Kramers inter-valley coherent state at both
half-filling and charge neutrality. Additionally, we investigate one- and
three-quarter fillings, and examine the effect of many-body corrections beyond
single Slater determinant solutions. We discuss the effect that details of the
IBM Hamiltonian have on the results, including different forms of double-
counting corrections, and the need to establish and precisely specify many-
body Hamiltonians to allow more direct and quantitative comparisons with
experiments in MATBG.

###### pacs:

Valid PACS appear here

##  I INTRODUCTION

The interaction-driven insulating phase and superconducting phase discovered
in magic-angle twisted bilayer graphene (MATBG) [1, 2] are inspiring intense
efforts, in both the theoretical and experimental domains, to understand the
mechanism behind them [3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17,
18, 19]. However, due to the interplay between electron-electron interactions
and flat bands with non-trivial topology, accurate, ab initio modeling of
TABLG presents many challenges [13, 14, 20, 21, 22]. The exploration of
effective models for MATBG falls loosely into several categories: Hubbard-like
models [5, 23] based in real-space; the so-called interacting Bistritzer and
MacDonald (BM) model based in momentum space, by including long-range
electron-electron interactions in the low-energy continuum BM model [3]; and
more recently, a heavy fermion description [24, 25].

In this work, we take the second approach, and perform a systematic many-body
study using the interacting BM (IBM) model. Previous studies of the IBM model
have sought to explain the interaction-driven insulating phase with
perturbation theory [13, 15] and Hartree Fock (HF) approaches [13, 6, 26, 27,
28]. Many-body studies, including exact diagonalization (ED) [15, 29],
determinantal quantum Monte Carlo (DQMC) [30, 31, 32, 33], density matrix
renormalization group (DMRG) [34, 35, 32, 36, 37, 37] and quantum chemistry
methods [38] have shown that the proper treatment of electron-electron
interactions and their interplay with dispersion are indispensable for
determining ground state properties. Much progress has been made in
characterizing the properties of this model. Important obstacles remain,
however, for establishing comparisons and connections with experiments [1, 39,
40, 41, 42, 43, 19].

The obstacles are several fold. The relevant model parameters and electron
fillings lie within the most difficult regime, and require very high
resolution to capture accurately the interplay between band dispersion and
electron-electron interactions. This is at the heart of strong correlation
physics, and presents challenges of various forms to every computational
method. Here the presence of topological bands and long-range interactions
necessitates large system sizes in order to properly reach the thermodynamic
limit. It has also become clear that the answers can be sensitive to model
details including double counting corrections, the effect of remote bands, and
accounting for strain, etc. The lack of accurate solutions to the model makes
it more challenging to pinpoint and understand discrepancies with experimental
observations.

Here we tackle this problem by employing the phaseless auxiliary-field quantum
Monte Carlo (AFQMC) method [44] to investigate the ground state of the IBM
model at integer fillings. Our approach is based on the momentum space
description and captures the non-trivial topology by retaining eight “flat”
bands. Long-range (screened) electron-electron interaction is treated
explicitly and rigorously. The phaseless AFQMC method mitigates the sign/phase
problem by a constrained path approximation [45, 46]; it enables simulations
of large lattice sizes (here calculations are performed up to $12\times 12$
cells), from which properties at the thermodynamic limit can be inferred. The
accuracy of the contrained path approximation (more precisely the phaseless
approximation [44] because of the long-range interaction) depends on the form
of Hubbard-Stratonovich transformation. The particular form to decompose a
long-range interaction written in momentum space has been applied in many
AFQMC studies of solids [47]. A more general form to decompose a 4-index
Coloumb interaction written in an arbitrary basis is also available [48, 49],
and has been applied for two decades in quantum chemistry [50, 51]. The
extensive benchmarks in these show that the phaseless AFQMC approach allows
accurate treatment of correlated materials with long-range interactions.

Using this approach we perform a systematic and comprehensive exploration of
the ground-state phase diagram of the IBM model. We focus on four integer
fillings that are specified by the number of electrons per
${\mathbf{k}}$-point $n_{f}$ (charge neutrality $n_{f}=4$, three quarter
$n_{f}=3$, half-filling $n_{f}=2$, and quarter-filling $n_{f}=1$). Our study
scans a wide range of tunneling ratio $\kappa$ that captures the lattice
relaxation between the intra-sublattice (AA) and inter-sublattice (AB$\&$BA)
regions.

At charge neutrality and half-filling, a transition from correlated insulating
states (i.e., Kramers inter-valley coherent, K-IVC) to metallic states is
observed as the tunneling ratio $\kappa$ is increased. This transition is
characterized by the correlation energy, charge gap, and related order
parameters. Our results indicate a loss of K-IVC long-range order which
contributes to the insulator-metal transition. At quarter-filling, a
transition from the correlated insulating state to a semi-metallic state is
observed as a function of tunneling ratio $\kappa$. This transition is
characterized by charge gap, and C2T symmetry. (Here, C2T indicates 180-degree
in-plane rotation symmetry C2 with time reversal T.) The AFQMC results confirm
that the ground states in both phases can be well described by a single
determinant (i.e., solution of HF), in agreement with previous works using
quantum chemistry methods (treating no valley degree of freedom) [38]. We also
find a degenerate state with correlations between two valleys in the semi-
metallic regime.

These results, in comparison with experimental observations, will also shed
light on the modeling of MATBG. Although much needs to be done in this model,
including fractional fillings, the systematic and accurate results at integer
filling provide an important step in determining what additional ingredients,
if any, will need to be considered to quantitatively describe MATBG.

The rest of the paper is organized as follows. In Sec. VII.2, we briefly
introduce the IBM model. In Sec. III, we introduce our AFQMC algorithm for
studying Coulomb interaction with the topological band and its typical
application in the IBM model. In Sec. IV, we describe the zero-temperature
phase diagram of the IBM model at integer fillings. Finally, in Sec. V, we
discuss connections with experiments comment on possible future work.

##  II IBM model

Our model consists of the continuum BM model $\hat{H}_{BM}$, which describes
the moire bands at charge neutrality, and long-range screened Coulomb
interactions $\hat{V}$. We briefly introduce IBM model here while elaborating
more details in App. VII.

The continuum BM model $\hat{H}_{BM}$ is a tight-binding model that describes
electrons hopping between $2p_{z}$ orbitals of carbon atoms in twisted bilayer
graphene. By taking Bloch vector $\bm{p}$, sub-lattice $\sigma$, and layer
$\mu$, spin $s$ to identify the Bloch state (see Eq. 11), the continuum BM
model can be written as:

$\displaystyle\hat{H}_{BM}$
$\displaystyle=\sum_{\bm{p}}T^{\textup{intra}}_{\bm{p},\bm{p}}c^{\dagger}_{\bm{p}}c_{\bm{p}}+T^{\textup{inter}}_{\bm{p}^{\prime},\bm{p}}c^{\dagger}_{\bm{p}^{\prime}}c_{\bm{p}}$
(1) $\displaystyle T^{\textup{intra}}_{\bm{p},\bm{p}}$
$\displaystyle=\hbar\upsilon_{F}(p_{x}\sigma_{x}\tau_{z}+p_{y}\sigma_{y})$
$\displaystyle T^{\textup{inter}}_{\bm{p}^{\prime},\bm{p}}$
$\displaystyle=\frac{1}{2}\sum_{l=0,1,2}\delta_{\bm{p}^{\prime},\bm{p}-\tau_{z}\bm{q}_{l}}(\mu_{x}-i\mu_{y})T_{l}+\textup{h.c.}$
$\displaystyle T_{l}$
$\displaystyle=\omega_{0}e^{-i\theta\sigma_{z}\mu_{z}\tau_{z}/2}+\omega_{1}e^{2\pi
il\sigma_{z}\tau_{z}/3}\sigma_{x}e^{-2\pi il\sigma_{z}\tau_{z}/3}$

where $\theta=1.05^{\circ}$ is the twist angle,
$\nu_{F}=2700*\frac{3a_{0}}{2\hbar}$ meV is the Dirac velocity, and
$\omega_{0}$, $\omega_{1}$ describes lattice relaxation between A-A and A-B
stacks, respectively. Transition momentum $\bm{q}_{l}=R(2\pi l/3)\bm{q}_{0}$,
where $R(\psi)=e^{i\psi\sigma_{y}}$ is a rotation matrix.
$\bm{q}_{0}=\frac{8\pi\textup{sin}(\theta/2)}{3\sqrt{3}a_{0}}(0,-1)^{T}$ and
$a_{0}=1.42\AA$ is the intra-layer distance between carbon atoms.

As two graphene layers stack with a twisted angle $\theta$, the Brillouin
zones (BZ) of mono-layer graphene rotate with angle $\theta$ related to the
other layer, and this rotation leads to periodic moire pattern and folds BZ to
mini-BZ (MBZ) [3]. That is we use $\varepsilon_{\bm{k},i}$ and
$|\varphi_{\bm{k},i}\rangle$ to specify each band $n=0,...,N_{G}-1$ in MBZ:

$\hat{H}_{BM}|\varphi_{\bm{k},i}\rangle=\varepsilon_{\bm{k},i}|\varphi_{\bm{k},i}\rangle$
(2)

with $i=(n,\tau,s)$ to label bands $n$, valley $\tau$, and spin $s$. $N_{G}$
is truncated for numerical implementation.

We mention here that a correction term is included in addition to the BM model
Hamiltonian $\hat{H}_{BM}$. The additional correction term,
$[\hat{V}]_{\phi_{BM}}$, is a mean-field approximation of the Coulomb
interactions $\hat{V}$ with respect to the ground state $\phi_{BM}$ of the BM
model at charge neutrality. This correction is adopted to remove or counter
the double counting of the interaction effect [13, 28, 30], as parameters used
in the BM model are based on first-principle calculations [3] which already
account for the effect of Coulomb interactions at the density-functional
theory level. There are somewhat different ways to treat the correction [38,
31, 32, 33], and the results turn out to be quite sensitive to the details of
the correction. As we illustrate in Fig. 1, the band-structure of the
dispersion part of the Hamiltonian,
$\hat{H}_{0}\equiv\hat{H}_{BM}-[\hat{V}]_{\phi_{BM}}$, can be significant
different from that of $\hat{H}_{BM}$.

Figure 1: (Top) Comparison of bandwidths between $\hat{H}_{BM}$ and
$\hat{H}_{0}$ versus $\kappa$ in a $6\times 6$ system. $\Gamma$ and $M$ label
the ${\mathbf{k}}$-point where the largest gap is observed. (Bottom) Band
structures of $\hat{H}_{0}$, $\hat{H}_{BM}$ and $\hat{H}_{HF}$ for
$\kappa=0.0$ and $\kappa=0.8$. $\hat{H}_{HF}$ indicates the IBM model band
structure obtained from Hartree-Fock (HF). Solid and dashed lines specify the
two valleys respectively (degenerate at $\kappa=0$).

It is also important to mention that the model is parameterized by the ratio
$\kappa=\omega_{0}/\omega_{1}$, which describes the lattice relaxation between
AA regions (characterized by $\omega_{0}$) and AB$\&$BA regions (characterized
by $\omega_{1}$). This ratio $\kappa$ is responsible for the emergence of flat
bands [10] and the ground state symmetry [14, 13]. To explore the effects of
lattice relaxation in TBG, we fix $\omega_{1}=105$meV and study
$\kappa\in(0,1)$, which covers the entire range from the chiral limit
($\kappa=0$) through the more realistic case ($\kappa\sim 0.8$).

The IBM model is built based on eight “flat” bands (i.e., two bands with two
valleys and two spins that are closest to the Fermi surface ) from the BM
model. Effect from the remote bands are neglected. The Hamiltonian of our IBM
model is

$\displaystyle\hat{H}_{IBM}=\hat{H}_{0}(\kappa)+\hat{V}$ (3)
$\displaystyle=\sum_{\bm{k},i}\varepsilon_{\bm{k},i}c^{\dagger}_{\bm{k},i}c_{\bm{k},i}-\sum_{\bm{k},i,j}[V]_{\bm{k},i;\bm{k},j}c^{\dagger}_{\bm{k},i}c_{\bm{k},j}+\frac{1}{2}\sum_{\bm{q}}\bar{\rho}_{\bm{q}}V_{\bm{q}}\bar{\rho}_{-\bm{q}}.$

where $\bm{k}$ is on a hexagonal lattice with size $L$ to discretize the first
MBZ (see Fig. 9), $\bar{\rho}_{\bm{q}}$ is the density operator at momentum
$\bm{q}=\bm{k}-\bm{k}^{\prime}$ with a background shift [49], and $V_{\bm{q}}$
is the Fourier transformed Coulomb potential. The transfer momentum $\bm{q}$
remains unrestricted, but is truncated for numerical implementation. The
Coulomb potential is two-gate screened:

$V_{\bm{q}}=\frac{e^{2}}{2\epsilon\epsilon_{0}}\frac{1}{A}\frac{\textup{tanh}(|\bm{q}|d)}{|\bm{q}|}.$
(4)

which is the Fourier-transformed Coulomb potential with the distance $d=20$ nm
between two layers and relative permittivity $\epsilon=10$ for vacuum
permittivity $\epsilon_{0}$. We comment on the choice of relative permittivity
in App. XI. $A=N_{M}A_{0}$ is the total area of the system with
$A_{0}=\sqrt{2}a^{2}_{M}/2$,
$a_{M}=\sqrt{3}a_{0}/2\textup{sin}\frac{\theta}{2}=13.4$ nm and $N_{M}$ is the
number of moire unite cells of this system.

##  III AFQMC with phaseless constraint

The IBM model is a much simplified model, which still contains both a non-
trivial band-structure, including corrections, and long-range Coulomb
interaction. The phaseless AFQMC [44] has the capacity to study such a general
Hamiltonian while controlling the sign problem to restore low-polynomial
computational scaling with lattice sizes. To give a brief self-contained
description, we adopt the notation in Ref. [49] and start from the Hamiltonian
we treat, written in “Monte Carlo form”:

$\hat{H}=\hat{T}+\frac{1}{2}\sum^{\Gamma}_{\gamma}\hat{L}_{\gamma}^{2}$ (5)

where $\hat{T}$ and $\hat{L}_{\gamma}$ denote sets of one-body operators whose
matrix elements are explicitly specified in the chosen basis. To study the IBM
model, we take
$\hat{L}_{\bm{q}}=\sqrt{\frac{V_{\bm{q}}}{2}}(\bar{\rho}_{\bm{q}}+\bar{\rho}_{\bm{-q}})$
and
$\hat{L}^{\prime}_{\bm{q}}=i\sqrt{\frac{V_{\bm{q}}}{2}}(\bar{\rho}_{\bm{q}}-\bar{\rho}_{\bm{-q}})$
[47].

AFQMC performs ground-state calculations through the imaginary-time
projection:

$\displaystyle|\Psi_{0}\rangle$
$\displaystyle\propto\underset{\beta\to\infty}{\lim}e^{-\beta\hat{H}}|\phi_{I}\rangle$
(6) $\displaystyle\approx(e^{-\tau\hat{H}})^{n}|\phi_{I}\rangle$
$\displaystyle\approx\int
d\\{\textbf{x}\\}\prod^{n}_{i=1}p(\textbf{x}_{i})\hat{B}(\textbf{x}_{i})|\phi_{I}\rangle\,,$

which propagates the given initial state $|\phi_{I}\rangle$ to the ground
state $|\Psi_{0}\rangle$ of the Hamiltonian $\hat{H}$, if the overlap
$\langle\Psi_{0}|\phi_{I}\rangle$ is non-zero. The second line of Eq. (6) is
obtained by using the Suzuki-Trotter decomposition [52, 53] to break the
imaginary-time evolution operator into time slices with $n=\beta/\tau$. The
third line is then obtained by applying the Hubbard-Stratonovich
transformation [54, 55] to decouple the short-time propagator in each time
slice with probability

$p(\textbf{x})=\prod_{\gamma}\frac{1}{\sqrt{2\pi}}e^{-x_{\gamma}^{2}/2},$ (7)

and the one-body propagator

$\hat{B}(\textbf{x})=e^{-\tau\hat{T}/2}e^{\sum_{\gamma}x_{\gamma}\sqrt{-\tau}\hat{L}_{\gamma}}e^{-\tau\hat{T}/2}\,,$
(8)

where x indicates a series of auxiliary-fields (AFs):
$\textbf{x}=\\{x_{\gamma}\\}$ ($x_{\gamma}\in\mathbb{R}$). Without loss of
generality, we can take $|\phi_{I}\rangle$ as a single determinant. Since a
one-body propagator propagates a single determinant to another single
determinant: $|\phi^{\prime}\rangle=\hat{B}(\textbf{x})|\phi\rangle$, the
integration in Eq. (6) can then be performed by Monte Carlo, mapping the
sampling of $\\{x_{\gamma}\\}$ into random walks in the manifold of Slater
determinants along the direction of imaginary time evolution [46].

Because of fermion antisymmetry, the Monte Carlo sampling of Eq. (6) in
general suffers from a phase problem [44]. The phase problem in AFQMC can be
understood by considering [46] a single determinant, through random chance in
the time evolution, becoming perpendicular to the ground state, after which
propagation of this determinant effectively contributes nothing but
statistical fluctuations to the ground state estimation. Since the number of
such determinants in general grows exponentially with projection time, the
signal of Monte Carlo sampling will eventually be overwhelmed by statistic
noise. Eliminating these determinants that are perpendicular to the ground
state is an exact condition which removes the sign or phase problem [45, 44].
To implement this condition, AFQMC uses a trial wave function
$|\Psi_{T}\rangle$ to constrain the random walk paths. For the general form of
interaction, the orbitals in the Slater determinant become complex,
effectively introducing an overall gauge degree of freedom that is redundant
in the projection. The constraint in this case is a phaseless approximation
that projects the weight of each Slater determinant sample to the real axis in
the complex plane [49]. If the trial wave function is the exact ground state,
then this constraint is unbiased. Many community benchmark studies have shown
that AFQMC with constraint is highly accurate, both in model systems and in
real materials [56, 57, 58].

To further validate the accuracy and robustness of the phaseless AFQMC
approach in this study, we have taken several additional measures. In the IBM
model different low-lying states of different symmetry characters can be
(nearly) degenerate. The AFQMC method we use is in general much less prone to
ergodicity problems in the sampling, because it takes the form of many (almost
independent) streams of branching random walks. Still, we perform multiple
AFQMC calculations using different low-lying HF states of different symmetry
as trial/initial wave function. The final energy is compared to gauge the
preferred symmetry of the many-body ground state. (Although the ground-state
energy computed by AFQMC is not strictly variational [59], the energy is very
accurate, and we only draw conclusions about the ordering of the states when
the energy separation is much more significant than the expected resolution of
the phaseless approximation.) For the study of $n_{f}=2,4$, C2T-symmetry-
breaking states and K-IVC states are barely distinguishable in energy around
the chiral limit. In these cases, we apply yet another step, using a
constraint release method, AFQMC with Metropolis release constraint (AFQMC-
MRC) [60], to minimize the constraint bias and identify the actual symmetry of
the ground state. Further details are given in the App. IX.

To summarize, in this work, we employ AFQMC to systematically investigate the
IBM model. Compared with previous quantum Monte Carlo work [30, 31, 33], this
method is not limited to special fillings, high temperatures, or small system
sizes. It allows us to access the whole range of $\kappa$ to determine the
complete phase diagram (with the largest lattice up to $12\times 12$), while
treating both valley and spin degrees of freedom.

##  IV Results

In this section, we describe our numerical results for the ground state of the
IBM model from AFQMC, divided into several subsections by filling fractions:
charge neutrality ($n_{f}=4$) and half-filling ($n_{f}=2$), quarter-filling
($n_{f}=1$), and three-quarter-filling ($n_{f}=3$). The tunneling ratio
$\kappa$ for realistic cases varies sample by sample and can affect the ground
state behaviors. To present a systematic study and facilitate comparison with
experiments, this work scans $\kappa$ from $0$ to $1$. At the flat band limit,
analytical studies [14, 13] indicate that the candidate ground states are
correlated insulating states: quantum anomalous Hall (QAH) state, valley Hall
(VH) state and K-IVC state for fillings $n_{f}=4$ and $2$, and QAH for filling
$n_{f}=1$. HF results (on the IBM model but with no double counting
correction) [26] support the above analysis and suggest the ground state at
$n_{f}=3$ filling is a mixture of the states at $n_{f}=4,2$ and $n_{f}=1$
fillings. As discussed in the last paragraph of Sec. III, we can use each of
the candidate symmetry states above as the trial wave function in AFQMC. In
this sense, our calculations can be viewed as using an unbiased (or at least
less biased than other possibilities) method, namely AFQMC, to judge the
validity of the different candidates as the many-body ground state.

Our results show that the HF solution is in some cases quite accurate, while
in other situations qualitatively incorrect. It is worth emphasizing that,
while a weakly interacting system is often well described by a single
determinant solution from a mean-field calculation, the reverse is not always
true. That is, a system which can be described well by a single Slater
determinant is not necessarily weakly interacting. An example of this is the
Wigner crystal state in the electron gas at the low-density limit, an
essentially classical state whose energy can be well approximated by that of a
single-determinant solution in the spirit of unrestricted HF, but which is a
strongly interacting state, with the interaction energy dominating the kinetic
(zero-point) energy. In the IBM, the system is often more analogous to this
situation when the HF solution is accurate. We quantify the errors from mean-
field theory calculations throughout the parameter regimes.

###  IV.1 Breaking of K-IVC insulating state at $n_{f}=4$ and $n_{f}=2$

At $n_{f}=4$ and $n_{f}=2$ fillings, we observe a transition from a metallic
state to a K-IVC insulating state. At large $\kappa$ ($\kappa>\sim 0.8$), we
find the ground states are K-IVC states characterized by a spontaneous
hybridization between the two valleys that preserves an effective time-
reversal symmetry. The K-IVC state occupies both spins at charge neutrality
and is polarized at one spin at half-filling. This is also the picture from
the HF solution, but quantitative discrepancies are seen between HF and AFQMC.
In contrast, at small $\kappa$ (i.e., $\kappa<\sim 0.2$) we observe
qualitative differences from the HF solution. The system is characterized by
much larger correlation energy, a loss of long range K-IVC order, and
vanishing charge gap. Below we first show our observation of this insulator-
metal transition through the ground-state energy, including correlation
energy, followed by results on the charge gap. Then we study the K-IVC order
parameter and band occupancy to complete this discussion.

Figure 2:  Ground-state energies and correlation energies. The top panel shows
the HF energy per electron, $E_{\rm HF}$, as a function of $\kappa$, for both
$n_{f}=4$ and $n_{f}=2$. The bottom panel shows the corresponding correlation
energies, $E_{c}$. Most results are obtained from a $6\times 6$ lattice in
momentum space, but at three values of $\kappa$ results are also shown for the
thermodynamic limit (TDL), which is obtained via extrapolation with $6\times
6$, $9\times 9$, and $12\times 12$ lattices. The inset in the bottom panel
shows the correlation energy relative to the total, $|E_{c}|/E_{\rm HF}$.

In Fig. 2, we show the HF energy per electron, $E_{\rm HF}$, and the AFQMC
correction over HF (i.e., the correlation energy) versus $\kappa$. The
correlation energy is an important metric widely adopted in electronic
structure calculations and quantum chemistry, defined as the difference
between the true many-body ground-state energy and that from HF, which in this
case is $E_{c}\equiv E_{\rm AFQMC}-E_{\rm HF}$, where $E_{\rm AFQMC}$ is the
ground-state energy per electron from AFQMC (calculated by the mixed estimator
[44, 49], with the HF state as the trial wave function). For $n_{f}=2$, HF
yields a spin-polarized K-IVC state; for $n_{f}=4$, it finds a superposition
of two spin-polarized K-IVC states with opposite spins and no correlation
between them. The K-IVC states at the two fillings have exactly the same HF
energy per electrons, since the mean-field background of Coulomb interaction
is removed in Eq. (3) by shifting the charge-density operators according to
the chemical potential at the corresponding filling $n_{f}$. In the chiral
limit ($\kappa=0$), HF yields an emergent C2T-symmetry-breaking (QAH, VH)
state that preserves $U_{V}(1)$ symmetry at both fillings and is degenerate
with the K-IVC state. Both the K-IVC and the C2T-symmetry-breaking states are
correlated insulating states. Details of the above analysis are elaborated in
App. VIII. The above observation agrees with previous analytical studies [14,
13] and the numerical HF results [26]. To simplify our studies, we focus on
the K-IVC state here, but briefly discuss the close competition between it and
the C2T-symmetry-breaking state in App. IX.

The amount of correlation energy $|E_{c}|$ grows as $\kappa$ decreases, as
seen in the bottom panel of Fig. 2. Beyond $\kappa=0.9$ virtually no
correlation energy is present, making the HF description essentially exact. At
lower $\kappa$ many-body correlations make significant corrections to the HF
results, which could render perturbative analysis more problematic.
Interestingly, although the correlation energy magnitude decreases
monotonically with $\kappa$, the percentage of the correlation energy relative
to the total ground-state energy (here measured with respect to the HF energy)
is not monotonic, reaching a maximum at intermediate $\kappa$ for $n_{f}=4$,
as seen in the inset of the bottom panel. We also observe a strong
manifestation of correlation effects in the separation of the $n_{f}=4$ and
$n_{f}=2$ curves, in contrast with previous analytical predictions [14, 13]
around the “flat” band limit. Correlation between different spins leads to a
dramatic lowering of the $n_{f}=4$ correlation energy from the spin-polarized
state.

Figure 3:  Charge gaps versus $\kappa$ at $n_{f}=4$ (top panel) and $n_{f}=2$
(bottom panel). Similar to Fig. 2, a $6\times 6$ lattice is used in the
calculations, but spot checks are performed at $\kappa=0.0$,$0.5$, and $0.9$
by extrapolating through $6\times 6$, $9\times 9$, and $12\times 12$ lattices
for the TDL. Shaded pink and green backgrounds separating the semi-metallic
and insulating phases are a rough guide to the eye.

We next examine the ground state properties in some detail. We compute the
charge gap of the system via the add/remove energy:

$\Delta_{c}=\frac{E(N+1)+E(N-1)-2\,E(N)\,,}{2}$ (9)

where $N$ is the total number of electrons in the simulation cell at the
targeted filling fraction, and $E(N)$ is the total ground-state energy. The
results are shown in Fig. 3 as a function of $\kappa$ for both $n_{f}=4$ and
$n_{f}=2$. At $\kappa\approx 0.9$, AFQMC finds the same finite charge gaps as
HF at both fillings, which is consistent with the picture from Fig. 2 and
confirms the K-IVC ground states in this regime. As $\kappa$ decreases, we
observe a quantitative difference between the two filling fractions, with the
charge gap increasing for $n_{f}=4$ and remaining essentially constant for
$n_{f}=2$, until $\kappa$ reaches an intermediate value around $0.4$, when a
sharp drop in $\Delta_{c}$ is seen at both fillings. As $\kappa$ further
decreases, the AFQMC charge gap is much reduced compared to the HF gap. A more
careful finite-size scaling study, as illustrated for $\kappa=0$, shows that
$\Delta_{c}$ from AFQMC becomes vanishingly small at the TDL, while the gap
from from HF remains very large. Thus, contrary to the HF prediction, the
many-body state is metallic at small $\kappa$ at both fillings, with a metal-
insulator transition taking place roughly around $\kappa\sim 0.3$.

Figure 4:  K-IVC order parameter (top) and band occupancy (bottom) versus
$\kappa$ at $n_{f}=2$ filling. Results are obtained with a $6\times 6$
lattice. Shaded pink and green backgrounds separate the metallic/semi-metallic
and insulating phases as in Fig. 3.

To better characterize the ground-state phases, we also compute the K-IVC
order parameter and band occupancy. Expectation values of operators that
commute with the Hamiltonian, such as the ground-state energies and gaps
discussed above, can be computed by the so-called mixed estimator [44] in
AFQMC. Other expectations require additional technical steps to compute, for
example by back-propagation [45, 61] or via linear response [62, 63]. We take
the latter approach here, which requires separate calculations for different
expectations but tends to be more accurate. More specifically, we calculate
the expectation value of an operator $\hat{O}$ with respect to the ground
state from AFQMC, $\langle\hat{O}\rangle$, via Hellman-Feynman theorem

$\langle\hat{O}\rangle=\underset{\lambda\to
0}{\lim}\frac{\mathrm{d}\langle\hat{H}_{IBM}+\lambda\hat{O}\rangle}{\mathrm{d}\lambda},$
(10)

which is then evaluated by finite difference with several ground-state energy
calculations at small $\lambda$. Further details are given in App. X.

In Fig. 4, we show the computed K-IVC order parameter
$\langle\tau_{x}\,n_{x}\rangle$ and band occupancy $\langle-n_{z}\rangle$ for
different $\kappa$ at $n_{f}=2$ filling. The subscripts of $n_{*}$ and
$\tau_{*}$ denote the related Pauli operator on bands and valleys (i.e.,
$\tau_{x}n_{x}=\frac{1}{L^{2}}\sum_{\textbf{k},s}c^{\dagger}_{\textbf{k},i^{\prime}=(-n,-\tau,s)}c_{\textbf{k},i=(n,\tau,s)}$,
and $n_{z}=-n\frac{1}{L^{2}}c^{\dagger}_{\textbf{k},i}c_{\textbf{k},i}$ with
$\tau=1$ for $+$ valley and $\tau=-1$ for $-$ valley, $n=1$ for top band and
$n=-1$ for bottom band.) The K-IVC order parameter is used to identify the
strength of K-IVC order (see App. VIII) and the band occupancy shows the
difference in electron distributions between the top and bottom bands. At
$\kappa=0.9$, a strong signal of K-IVC order is seen, and electrons are almost
evenly distributed over the top and bottom bands. These observations again
confirm a K-IVC ground state. As $\kappa$ decreases, we find a decay of the
K-IVC order in both the HF and AFQMC results. However, the K-IVC order in
AFQMC decays much faster, and vanishes at $\kappa\sim 0.2$, while in HF it
stays at a significant finite value. In the bottom panel the band occupancy
increases as $\kappa$ decreases. A large $\langle-n_{z}\rangle$ implies a
higher occupancy of the bottom bands. The computed value from AFQMC is always
larger than from HF, indicating that the electrons tend to stay in the bottom
bands more. Since electrons forming K-IVC pairs tend to be evenly distributed
over the top/bottom band, the results on $\langle-n_{z}\rangle$ suggest a
decay of K-IVC order with decreasing $\kappa$, consistent with
$\langle\tau_{x}n_{x}\rangle$. By symmetry [13], the above analysis and
observations can be directly applied to $n_{f}=4$ filling.

We also examined the stability of a fully spin-polarized ground state for
$n_{f}=2$ and non-spin-polarized ground state for $n_{f}=4$ by computing the
energy cost for flipping spins. These results are included in App. IX.

###  IV.2 Breaking of C2T symmetry at $n_{f}=1$

Figure 5:  Ground-state energies and correlation energies for spin valley
polarized (SVP) state and inter-valley coherent (IVC) state, at $n_{f}=1$. The
top panel shows the HF energy per electron, $E_{\rm HF}$, as a function of
$\kappa$, for the two types of states considered. The bottom panel shows the
corresponding correlation energies, $E_{c}$. Most results are obtained from a
$6\times 6$ lattice in momentum space, but at three values of $\kappa$ results
are shown for the SVP state at the TDL, using the same extrapolation procedure
as in Fig. 2. Beyond $\kappa=0.8$ (indicated by the vertical black dashed
lines) the SVP and IVC states are degenerate within HF, while small energy
differences are seen in AFQMC. However when extrapolated to the TDL at
$\kappa=0.85$, the energy difference in AFQMC is zero within statistical error
bar.

At $n_{f}=1$, we find that the ground state can be reasonably approximated by
HF through a comparison between HF and AFQMC results in both correlation
energy and charge gap. This is consistent with results from a previous quantum
chemistry study [38]. Based on HF, we observed a transition along $\kappa$
from the C2T-symmetry-breaking correlated insulating state (i.e., Chern
polarized QAH state [13]) to a C2T-symmetry-restored semi-metallic state and
an emerged spin-polarized inter-valley-correlated semi-metallic state.

In Fig. 5, we show HF energy in the top panel and the AFQMC correction over HF
in the bottom panel for different $\kappa$ at $n_{f}=1$. Our HF study found a
spin-valley-polarized (SVP) state for all $\kappa$, and a spin-polarized
inter-valley coherent (IVC) state at $\kappa\geq 0.8$. In the bottom panel,
the same correlation energy behaviors are seen as observed for $n_{f}=4,2$,
where many-body effects are found to be strong at the chiral limit (i.e.,
$\kappa=0$) and become vanishingly small at large $\kappa$. However, the
correlation energy for $n_{f}=1$ is much smaller than in $n_{f}=4,2$, and
AFQMC only leads to a quantitative correction over HF.

Figure 6: Charge gap for $n_{f}=1$ filling versus $\kappa$. Charge gaps are
calculated by measuring the energy required to add or remove an electron from
the system. AFQMC proceeds with the corresponding HF state as the trial wave
function, whose HF charge gap is illustrated in the figure. The thermodynamic
limits (TDL) are extrapolated through 6×6, 9×9, and 12×12 lattice
calculations. Shaded pink and green backgrounds separate the metallic/semi-
metallic and insulating phases.

In Fig. 6, we show the charge gap for different $\kappa$. A finite charge gap
is seen for $\kappa<0.8$ at the thermodynamic limit, which suggests an
insulating ground state. A sudden drop of the charge gap is observed around
$\kappa=0.8$, indicating a transition to semi-metallic states. This behavior,
however, is quite different from what is observed at $n_{f}=4,2$, where mean-
field calculations do not capture the metallic states that emerge.

Figure 7: K-IVC order parameter $\tau_{x}n_{x}$ and C2T-symmetry-breaking
order parameters $C(gK)$ at $n_{f}=1$ for HF ground states with spin-valley-
polarized (SVP) or Inter-valley coherent (IVC) symmetry breaking. In the
degenerate region (i.e., between $\kappa=0.8$ and $\kappa=0.95$), HF states
are determined by imposing the spin-polarized/inter-valley-coherent initial
density, respectively. Shaded green and pink backgrounds separate the
insulating and metallic/semi-metallic phases as in Fig. 6.

In Fig. 7, we show that the above insulating-semi-metallic transition is
accompanied by a restoration of C2T-symmetry and an emerge of IVC state. Since
our energy and charge gap results show only minor differences between AFQMC
and HF here, we consider the HF state as a qualitatively good approximation
for the ground state, and we use HF to study the behavior of IVC and
C2T-symmetry by focusing on the K-IVC order parameter
$\langle\tau_{x}n_{x}\rangle$ and the C2T-symmetry-breaking order parameter
$C(gK)$ proposed in [38], following identical notation ($g$ is unitary
symmetry and $K$ is complex conjugation). $C(gK)$ vanishes when C2T symmetry
is completely restored. At small $\kappa$, the ground state is well described
by a spin-valley polarized QAH state where C2T-symmetry is significantly
broken, and no K-IVC order is observed. As $\kappa$ increases, a degeneracy of
the SVP state and spin-polarized IVC state is seen. For the SVP state, a
sudden drop of $C(gK)$ indicates a restoration of C2T-symmetry. For the spin-
polarized IVC state, we find a significant K-IVC order parameter with
significant C2T symmetry breaking. Remarkably, $\tau_{x}n_{x}\sim 1$ suggests
all electrons at $n_{f}=1$ form K-IVC pairs, but this is not a K-IVC state
found in $n_{f}=4,2$ as there are not enough K-IVC pairs to fill entire
valence bands.

###  IV.3  Mixture at $n_{f}=3$

At $n_{f}=3$ filling, we observed a “phase diagram” as a function of $\kappa$
that shows features from both $n_{f}=2$ and $n_{f}=1$. In Fig. 8, we show the
charge gap for different $\kappa$. For the HF study, we observed a K-IVC state
polarized at one spin flavor and a C2T-symmetry-breaking state polarized at
the other spin and one valley, which agrees with previous results [26]. For
the AFQMC study, we also observed a combination of two states in the many-body
sense. We observed a metallic state for small $\kappa$, and a transition to an
insulating state as $\kappa$ is increased. As $\kappa\sim 0.8$ is approached,
the drop in the charge gap is observed, which suggests a semi-metallic state.
This behavior can be understood as a result of the direct combination between
the state at $n_{f}=2$ and $n_{f}=1$.

Figure 8: Charge gap for $n_{f}=3$ filling versus $\kappa$. AFQMC calculations
use the corresponding HF state as the trial wave function, whose charge gap is
also shown in the figure. Shaded pink and green backgrounds separate the
metallic/semi-metallic and insulating phases.

##  V  Summary and discussion

In this work, we investigated a model for twisted bilayer graphene that
combines the BM model with realistic electron-electron Coulomb interactions.
To minimize the difference between real materials and the model, we consider a
mean-field correction to remove the double counting of Coulomb interactions
and tune the ratio $\kappa$ to include the effects of lattice relaxation. We
employ HF and AFQMC to study this model at three typical integer fillings, in
order to probe and quantify the effect of electron correlation.

We implemented state-of-the-art phaseless AFQMC to enable an accurate and
systematic study of the ground state of the IBM model for a range of
parameters and fillings. By controlling the sign problem with the phaseless
constraint, we perform calculations in the ground state and in lattices as
large as $12\times 12$ to reach the thermodynamic limit. Numerical Hartree
Fork is first carried out to produce mean-field candidates of the ground
states. Based on symmetry breaking in perturbation theories, these mean-filled
states are correlated insulating states with either K-IVC or C2T symmetry
breaking. We then use the HF states as trial wave functions for AFQMC and
determine the actual ground states. The work presented here paves the way for
systematic studies of realistic models of MATBG. The same AFQMC technology we
have set up can be directly applied to reach beyond the integer fillings
studied here.

The outcome of these calculations is summarized as ground state phase diagrams
for the IBM model for different integer fillings as a function of the
tunneling ratio $\kappa$. Our results at charge neutrality are consistent with
those from previous sign-problem-free QMC studies which focused on limited
parameters [30, 31, 33]. At charge neutrality and half-filling, we observe an
insulator-metal transition with a loss of spin gap and K-IVC long-range order.
The results suggest that this transition is the result of many-body effects
driven by competition between topological band and Coulomb interaction. At
quarter-filling, we observe an insulator-metal transition that comes with a
restoration of C2T symmetry and an emerged inter-valley-correlated state. Our
results indicate that the transition can be qualitatively described by HF. We
observe the restoration of C2T symmetry confirming conclusions of a previous
study which did not include valley degree of freedom [38]. An extra spin-
polarized IVC state is seen in our study. We also investigate the ground state
phase diagram at three-quarter-filling. As a function of $\kappa$ it shows
both features from both charge neutrality/half-filling and quarter-filling.

We find that the ground state of the IBM model is very sensitive to the
corrections of the bare band structure. This has also been seen in other many-
body studies [38, 30, 36]. Though conclusions from many-body studies,
including ours, appear to yield qualitatively consistent results, the detailed
model in each work involves different mean-field corrections, which lead to
quantitative discrepancies. We believe it would be very valuable to have a
single model to study. Agreement between different many-body calculations
would then help “converge” a set of reliable solutions to the model.
Discrepancies between different methods would help identify areas of
improvement needed and spur method development and synergistic studies [58,
57].

For example, recent experiments [19] suggest the ground state of twisted
bilayer graphene is in a time-reversal symmetric inter-valley coherent (T-IVC)
phase, while a small strain can lead to incommensurate Kekulé spiral (IKS)
order. Are these states captured by the IBM model but not by the detailed
Hamiltonian we (and others) have studied? Or are there limitations in our
computations? The aforementioned unique model and converged solution would
allow more direct comparisons with experimental observations, which can then
identify model limitations, and improve our understanding of essential
physics.

##  VI  Acknowledgments

We thank L. Lin, A. H. MacDonald, O. Vafek for helpful conversations. Z.X.
acknowledges partial support from the U.S. Department of Energy (DOE) under
grant DE-SC0001303. Z.X. is also grateful for the support and hospitality of
the Center for Computational Quantum Physics (CCQ), where this work was
performed. Computing was carried out at the computational facilities at
Flatiron Institute. The Flatiron Institute is a division of the Simons
Foundation.

## References

  * Cao  _et al._ [2018a] Y. Cao, V. Fatemi, A. Demir, S. Fang, S. L. Tomarken, J. Y. Luo, J. D. Sanchez-Yamagishi, K. Watanabe, T. Taniguchi, E. Kaxiras, R. C. Ashoori, and P. Jarillo-Herrero, Correlated insulator behaviour at half-filling in magic-angle graphene superlattices, Nature 556, 80 (2018a).
  * Cao  _et al._ [2018b] Y. Cao, V. Fatemi, S. Fang, K. Watanabe, T. Taniguchi, E. Kaxiras, and P. Jarillo-Herrero, Unconventional superconductivity in magic-angle graphene superlattices, Nature 556, 43 (2018b).
  * Bistritzer and MacDonald [2011] R. Bistritzer and A. H. MacDonald, Moiré bands in twisted double-layer graphene, Proceedings of the National Academy of Sciences 108, 12233 (2011), https://www.pnas.org/doi/pdf/10.1073/pnas.1108174108 .
  * Xu and Balents [2018] C. Xu and L. Balents, Topological superconductivity in twisted multilayer graphene, Phys. Rev. Lett. 121, 087001 (2018).
  * Kang and Vafek [2018] J. Kang and O. Vafek, Symmetry, maximally localized wannier states, and a low-energy model for twisted bilayer graphene narrow bands, Phys. Rev. X 8, 031088 (2018).
  * Xie and MacDonald [2020] M. Xie and A. H. MacDonald, Nature of the correlated insulator states in twisted bilayer graphene, Phys. Rev. Lett. 124, 097601 (2020).
  * Dodaro  _et al._ [2018] J. F. Dodaro, S. A. Kivelson, Y. Schattner, X. Q. Sun, and C. Wang, Phases of a phenomenological model of twisted bilayer graphene, Phys. Rev. B 98, 075154 (2018).
  * Zou  _et al._ [2018] L. Zou, H. C. Po, A. Vishwanath, and T. Senthil, Band structure of twisted bilayer graphene: Emergent symmetries, commensurate approximants, and wannier obstructions, Phys. Rev. B 98, 085435 (2018).
  * Bi  _et al._ [2019] Z. Bi, N. F. Q. Yuan, and L. Fu, Designing flat bands by strain, Phys. Rev. B 100, 035448 (2019).
  * Tarnopolsky  _et al._ [2019] G. Tarnopolsky, A. J. Kruchkov, and A. Vishwanath, Origin of magic angles in twisted bilayer graphene, Phys. Rev. Lett. 122, 106405 (2019).
  * Kang and Vafek [2019] J. Kang and O. Vafek, Strong coupling phases of partially filled twisted bilayer graphene narrow bands, Phys. Rev. Lett. 122, 246401 (2019).
  * Seo  _et al._ [2019] K. Seo, V. N. Kotov, and B. Uchoa, Ferromagnetic mott state in twisted graphene bilayers at the magic angle, Phys. Rev. Lett. 122, 246402 (2019).
  * Bultinck  _et al._ [2020] N. Bultinck, E. Khalaf, S. Liu, S. Chatterjee, A. Vishwanath, and M. P. Zaletel, Ground state and hidden symmetry of magic-angle graphene at even integer filling, Phys. Rev. X 10, 031034 (2020).
  * Khalaf  _et al._ [2021] E. Khalaf, S. Chatterjee, N. Bultinck, M. P. Zaletel, and A. Vishwanath, Charged skyrmions and topological origin of superconductivity in magic-angle graphene, Science Advances 7, eabf5299 (2021), https://www.science.org/doi/pdf/10.1126/sciadv.abf5299 .
  * Lian  _et al._ [2021] B. Lian, Z.-D. Song, N. Regnault, D. K. Efetov, A. Yazdani, and B. A. Bernevig, Twisted bilayer graphene. iv. exact insulator ground states and phase diagram, Phys. Rev. B 103, 205414 (2021).
  * Kwan  _et al._ [2021] Y. H. Kwan, G. Wagner, T. Soejima, M. P. Zaletel, S. H. Simon, S. A. Parameswaran, and N. Bultinck, Kekulé spiral order at all nonzero integer fillings in twisted bilayer graphene, Phys. Rev. X 11, 041063 (2021).
  * Utama  _et al._ [2021] M. I. B. Utama, R. J. Koch, K. Lee, N. Leconte, H. Li, S. Zhao, L. Jiang, J. Zhu, K. Watanabe, T. Taniguchi, P. D. Ashby, A. Weber-Bargioni, A. Zettl, C. Jozwiak, J. Jung, E. Rotenberg, A. Bostwick, and F. Wang, Visualization of the flat electronic band in twisted bilayer graphene near the magic angle twist, Nature Physics 17, 184 (2021).
  * Lisi  _et al._ [2021] S. Lisi, X. Lu, T. Benschop, T. A. de Jong, P. Stepanov, J. R. Duran, F. Margot, I. Cucchi, E. Cappelli, A. Hunter, A. Tamai, V. Kandyba, A. Giampietri, A. Barinov, J. Jobst, V. Stalman, M. Leeuwenhoek, K. Watanabe, T. Taniguchi, L. Rademaker, S. J. van der Molen, M. P. Allan, D. K. Efetov, and F. Baumberger, Observation of flat bands in twisted bilayer graphene, Nature Physics 17, 189 (2021).
  * Nuckolls  _et al._ [2023] K. P. Nuckolls, R. L. Lee, M. Oh, D. Wong, T. Soejima, J. P. Hong, D. Călugăru, J. Herzog-Arbeitman, B. A. Bernevig, K. Watanabe, T. Taniguchi, N. Regnault, M. P. Zaletel, and A. Yazdani, Quantum textures of the many-body wavefunctions in magic-angle graphene (2023), arXiv:2303.00024 [cond-mat.mes-hall] .
  * Ahn  _et al._ [2019] J. Ahn, S. Park, and B.-J. Yang, Failure of nielsen-ninomiya theorem and fragile topology in two-dimensional systems with space-time inversion symmetry: Application to twisted bilayer graphene at magic angle, Phys. Rev. X 9, 021013 (2019).
  * Po  _et al._ [2019] H. C. Po, L. Zou, T. Senthil, and A. Vishwanath, Faithful tight-binding models and fragile topology of magic-angle bilayer graphene, Phys. Rev. B 99, 195455 (2019).
  * Song  _et al._ [2021] Z.-D. Song, B. Lian, N. Regnault, and B. A. Bernevig, Twisted bilayer graphene. ii. stable symmetry anomaly, Phys. Rev. B 103, 205412 (2021).
  * Koshino  _et al._ [2018] M. Koshino, N. F. Q. Yuan, T. Koretsune, M. Ochi, K. Kuroki, and L. Fu, Maximally localized wannier orbitals and the extended hubbard model for twisted bilayer graphene, Phys. Rev. X 8, 031087 (2018).
  * Song and Bernevig [2022] Z.-D. Song and B. A. Bernevig, Magic-angle twisted bilayer graphene as a topological heavy fermion problem, Phys. Rev. Lett. 129, 047601 (2022).
  * Rai  _et al._ [2023] G. Rai, L. Crippa, D. Călugăru, H. Hu, L. de’ Medici, A. Georges, B. A. Bernevig, R. Valentí, G. Sangiovanni, and T. Wehling, Dynamical correlations and order in magic-angle twisted bilayer graphene (2023), arXiv:2309.08529 [cond-mat.str-el] .
  * Zhang  _et al._ [2020] Y. Zhang, K. Jiang, Z. Wang, and F. Zhang, Correlated insulating phases of twisted bilayer graphene at commensurate filling fractions: A hartree-fock study, Phys. Rev. B 102, 035136 (2020).
  * Zhang  _et al._ [2022] S. Zhang, X. Lu, and J. Liu, Correlated insulators, density wave states, and their nonlinear optical response in magic-angle twisted bilayer graphene, Phys. Rev. Lett. 128, 247402 (2022).
  * Liu  _et al._ [2021] S. Liu, E. Khalaf, J. Y. Lee, and A. Vishwanath, Nematic topological semimetal and insulator in magic-angle bilayer graphene at charge neutrality, Phys. Rev. Res. 3, 013033 (2021).
  * Potasz  _et al._ [2021] P. Potasz, M. Xie, and A. H. MacDonald, Exact diagonalization for magic-angle twisted bilayer graphene, Phys. Rev. Lett. 127, 147203 (2021).
  * Hofmann  _et al._ [2022] J. S. Hofmann, E. Khalaf, A. Vishwanath, E. Berg, and J. Y. Lee, Fermionic monte carlo study of a realistic model of twisted bilayer graphene, Phys. Rev. X 12, 011061 (2022).
  * Zhang  _et al._ [2021] X. Zhang, G. Pan, Y. Zhang, J. Kang, and Z. Y. Meng, Momentum space quantum monte carlo on twisted bilayer graphene, Chinese Physics Letters 38, 077305 (2021).
  * Da Liao  _et al._ [2021] Y. Da Liao, J. Kang, C. N. Breiø, X. Y. Xu, H.-Q. Wu, B. M. Andersen, R. M. Fernandes, and Z. Y. Meng, Correlation-induced insulating topological phases at charge neutrality in twisted bilayer graphene, Phys. Rev. X 11, 011014 (2021).
  * Huang  _et al._ [2024] C. Huang, X. Zhang, G. Pan, H. Li, K. Sun, X. Dai, and Z. Y. Meng, Evolution from quantum anomalous hall insulator to heavy-fermion semimetal in magic-angle twisted bilayer graphene, Phys. Rev. B 109, 125404 (2024).
  * Soejima  _et al._ [2020] T. Soejima, D. E. Parker, N. Bultinck, J. Hauschild, and M. P. Zaletel, Efficient simulation of moiré materials using the density matrix renormalization group, Phys. Rev. B 102, 205111 (2020).
  * Chatterjee  _et al._ [2022] S. Chatterjee, M. Ippoliti, and M. P. Zaletel, Skyrmion superconductivity: Dmrg evidence for a topological route to superconductivity, Phys. Rev. B 106, 035421 (2022).
  * Parker  _et al._ [2021] D. E. Parker, T. Soejima, J. Hauschild, M. P. Zaletel, and N. Bultinck, Strain-induced quantum phase transitions in magic-angle graphene, Phys. Rev. Lett. 127, 027601 (2021).
  * Wang  _et al._ [2022] T. Wang, D. E. Parker, T. Soejima, J. Hauschild, S. Anand, N. Bultinck, and M. P. Zaletel, Kekulé spiral order in magic-angle graphene: a density matrix renormalization group study (2022), arXiv:2211.02693 [cond-mat.str-el] .
  * Faulstich  _et al._ [2022] F. M. Faulstich, K. D. Stubbs, Q. Zhu, T. Soejima, R. Dilip, H. Zhai, R. Kim, M. P. Zaletel, G. K.-L. Chan, and L. Lin, Interacting models for twisted bilayer graphene: a quantum chemistry approach (2022), arXiv:2211.09243 [cond-mat.str-el] .
  * Polshyn  _et al._ [2019] H. Polshyn, M. Yankowitz, S. Chen, Y. Zhang, K. Watanabe, T. Taniguchi, C. R. Dean, and A. F. Young, Large linear-in-temperature resistivity in twisted bilayer graphene, Nature Physics 15, 1011 (2019).
  * Lu  _et al._ [2019] X. Lu, P. Stepanov, W. Yang, M. Xie, M. A. Aamir, I. Das, C. Urgell, K. Watanabe, T. Taniguchi, G. Zhang, A. Bachtold, A. H. MacDonald, and D. K. Efetov, Superconductors, orbital magnets and correlated states in magic-angle bilayer graphene, Nature 574, 653 (2019).
  * Saito  _et al._ [2020] Y. Saito, J. Ge, K. Watanabe, T. Taniguchi, and A. F. Young, Independent superconductors and correlated insulators in twisted bilayer graphene, Nature Physics 16, 926 (2020).
  * Stepanov  _et al._ [2020] P. Stepanov, I. Das, X. Lu, A. Fahimniya, K. Watanabe, T. Taniguchi, F. H. L. Koppens, J. Lischner, L. Levitov, and D. K. Efetov, Untying the insulating and superconducting orders in magic-angle graphene, Nature 583, 375 (2020).
  * Saito  _et al._ [2021] Y. Saito, J. Ge, L. Rademaker, K. Watanabe, T. Taniguchi, D. A. Abanin, and A. F. Young, Hofstadter subband ferromagnetism and symmetry-broken chern insulators in twisted bilayer graphene, Nature Physics 17, 478 (2021).
  * Zhang and Krakauer [2003] S. Zhang and H. Krakauer, Quantum monte carlo method using phase-free random walks with slater determinants, Phys. Rev. Lett. 90, 136401 (2003).
  * Zhang  _et al._ [1997] S. Zhang, J. Carlson, and J. E. Gubernatis, Constrained path monte carlo method for fermion ground states, Phys. Rev. B 55, 7464 (1997).
  * Zhang [2019] S. Zhang, Auxiliary-field quantum monte carloat zero- and finite-temperature, in  _Many-Body Methods for Real Materials_ , edited by E. Erik Koch and Shiwei Zhang (eds.) (Verlag des Forschungszentrum Jülich, 2019).
  * Suewattana  _et al._ [2007] M. Suewattana, W. Purwanto, S. Zhang, H. Krakauer, and E. J. Walter, Phaseless auxiliary-field quantum monte carlo calculations with plane waves and pseudopotentials: Applications to atoms and molecules, Phys. Rev. B 75, 245123 (2007).
  * Purwanto  _et al._ [2011] W. Purwanto, H. Krakauer, Y. Virgus, and S. Zhang, Assessing weak hydrogen binding on Ca+ centers: An accurate many-body study with large basis sets, The Journal of Chemical Physics 135, 164105 (2011), https://pubs.aip.org/aip/jcp/article-pdf/doi/10.1063/1.3654002/15444035/164105_1_online.pdf .
  * Shi and Zhang [2021] H. Shi and S. Zhang, Some recent developments in auxiliary-field quantum monte carlo for real materials, The Journal of Chemical Physics 154, 024107 (2021), https://doi.org/10.1063/5.0031024 .
  * Lee  _et al._ [2022] J. Lee, H. Q. Pham, and D. R. Reichman, Twenty years of auxiliary-field quantum monte carlo in quantum chemistry: An overview and assessment on main group chemistry and bond-breaking, Journal of Chemical Theory and Computation 18, 7024 (2022), pMID: 36255074, https://doi.org/10.1021/acs.jctc.2c00802 .
  * Shee  _et al._ [2023] J. Shee, J. L. Weber, D. R. Reichman, R. A. Friesner, and S. Zhang, On the potentially transformative role of auxiliary-field quantum Monte Carlo in quantum chemistry: A highly accurate method for transition metals and beyond, The Journal of Chemical Physics 158, 10.1063/5.0134009 (2023), 140901, https://pubs.aip.org/aip/jcp/article-pdf/doi/10.1063/5.0134009/16824716/140901_1_5.0134009.pdf .
  * Suzuki [1976] M. Suzuki, Relationship between d-Dimensional Quantal Spin Systems and (d+1)-Dimensional Ising Systems: Equivalence, Critical Exponents and Systematic Approximants of the Partition Function and Spin Correlations, Progress of Theoretical Physics 56, 1454 (1976), https://academic.oup.com/ptp/article-pdf/56/5/1454/5264429/56-5-1454.pdf .
  * Trotter [1959] H. F. Trotter, On the product of semi-groups of operators, Proceedings of the American Mathematical Society 10, 545 (1959).
  * Hubbard [1959] J. Hubbard, Calculation of partition functions, Phys. Rev. Lett. 3, 77 (1959).
  * Stratonovich [1957] R. L. Stratonovich, On a Method of Calculating Quantum Distribution Functions, Soviet Physics Doklady 2, 416 (1957).
  * Williams  _et al._ [2020] K. T. Williams, Y. Yao, J. Li, L. Chen, H. Shi, M. Motta, C. Niu, U. Ray, S. Guo, R. J. Anderson, J. Li, L. N. Tran, C.-N. Yeh, B. Mussard, S. Sharma, F. Bruneval, M. van Schilfgaarde, G. H. Booth, G. K.-L. Chan, S. Zhang, E. Gull, D. Zgid, A. Millis, C. J. Umrigar, and L. K. Wagner (Simons Collaboration on the Many-Electron Problem), Direct comparison of many-body methods for realistic electronic hamiltonians, Phys. Rev. X 10, 011041 (2020).
  * Motta  _et al._ [2019] M. Motta, C. Genovese, F. Ma, Z.-H. Cui, R. Sawaya, G. K. Chan, N. Chepiga, P. Helms, C. Jimenez-Hoyos, A. J. Millis, _et al._ , Ground-state properties of the hydrogen chain: insulator-to-metal transition, dimerization, and magnetic phases, arXiv preprint arXiv:1911.01618 (2019).
  * LeBlanc  _et al._ [2015] J. P. F. LeBlanc, A. E. Antipov, F. Becca, I. W. Bulik, G. K.-L. Chan, C.-M. Chung, Y. Deng, M. Ferrero, T. M. Henderson, C. A. Jiménez-Hoyos, E. Kozik, X.-W. Liu, A. J. Millis, N. V. Prokof’ev, M. Qin, G. E. Scuseria, H. Shi, B. V. Svistunov, L. F. Tocchio, I. S. Tupitsyn, S. R. White, S. Zhang, B.-X. Zheng, Z. Zhu, and E. Gull (Simons Collaboration on the Many-Electron Problem), Solutions of the two-dimensional hubbard model: Benchmarks and results from a wide range of numerical algorithms, Phys. Rev. X 5, 041041 (2015).
  * Carlson  _et al._ [1999] J. Carlson, J. E. Gubernatis, G. Ortiz, and S. Zhang, Issues and observations on applications of the constrained-path monte carlo method to many-fermion systems, Phys. Rev. B 59, 12788 (1999).
  * Xiao  _et al._ [2023] Z.-Y. Xiao, H. Shi, and S. Zhang, Interfacing branching random walks with metropolis sampling: Constraint release in auxiliary-field quantum monte carlo, Journal of Chemical Theory and Computation 19, 6782 (2023), pMID: 37661928, https://doi.org/10.1021/acs.jctc.3c00521 .
  * Motta and Zhang [2017] M. Motta and S. Zhang, Computation of ground-state properties in molecular systems: Back-propagation with auxiliary-field quantum monte carlo, Journal of Chemical Theory and Computation 13, 5367 (2017), pMID: 29053270, https://doi.org/10.1021/acs.jctc.7b00730 .
  * Xu  _et al._ [2022] H. Xu, H. Shi, E. Vitali, M. Qin, and S. Zhang, Stripes and spin-density waves in the doped two-dimensional hubbard model: Ground state phase diagram, Phys. Rev. Res. 4, 013239 (2022).
  * Qin  _et al._ [2020] M. Qin, C.-M. Chung, H. Shi, E. Vitali, C. Hubig, U. Schollwöck, S. R. White, and S. Zhang (Simons Collaboration on the Many-Electron Problem), Absence of superconductivity in the pure two-dimensional hubbard model, Phys. Rev. X 10, 031016 (2020).

##  VII IBM model

###  VII.1 BM model

We recap key features of the BM model [3, 30]. To establish a tight-binding
description of TBG, we first consider electrons on two graphene layers with no
inter-layer interactions (i.e., the distance between two layers
$a_{0}\longrightarrow\infty$). We use unit cell index $l$, sub-lattice
$\sigma$, and layer $\mu$ to identify the position of a selected carbon atom.
We use $c^{\dagger}_{l,m,s}$, where $m=(\sigma,\mu)$, to specify the operator
that creates an electron with spin $s$ in the $2p_{z}$ orbital
$\phi(\bm{r}-\bm{R}_{l,m})$ around the carbon atom at $\bm{R}_{l,m}$.

The Bloch state

$\displaystyle\Phi_{\bm{p},m}(\bm{r})=e^{i\bm{p}\cdot\bm{r}}u_{\bm{p},m}(\bm{r})$
(11)

where $u_{\bm{p},m}(\bm{r})=u_{\bm{p},m}(\bm{r}+\bm{R}_{l,m})$ is a periodic
function, under tight-binding approximation, can be written as:

$\Phi_{\bm{p},m}(\bm{r})=\sum_{l}e^{i\bm{p}\cdot\bm{R}_{l,m}}\phi(\bm{r}-\bm{R}_{l,m}).$

In this appendix, $c^{\dagger}_{\bm{p},m,s}$ indicates creating a spin $s$
electron with Bloch state $\Phi_{\bm{p},m}(\bm{r})$.

As demonstrated in Fig. 9, if two layers stack with a twisted angle $\theta$,
the Brillouin zones (BZ) of one mono-layer of graphene rotate with angle
$\theta$ relative to the other layer, and this rotation leads to a periodic
moire pattern and a mini-BZ (MBZ). With the above notations, the tight binding
Hamiltonian of bilayer graphene with no inter-layer hopping (for $\bm{p}$
measured around and from $K^{\mu}_{\pm}$ points) can be written as below:

$\displaystyle\hat{H}$
$\displaystyle=\sum_{\bm{p}}c^{\dagger}_{\bm{p}}T^{\textup{inter}}_{\bm{p},\bm{p}}c_{\bm{p}}\,$
(12) $\displaystyle T^{\textup{intra}}_{\bm{p},\bm{p}}$
$\displaystyle=\hbar\upsilon_{F}(p_{x}\sigma_{x}\tau_{z}+p_{y}\sigma_{y})e^{-i\theta\sigma_{z}\mu_{z}\tau_{z}/2}\,,$

where $\nu_{F}$ is Dirac velocity,
$c^{\dagger}_{\bm{p}}=(...,c^{\dagger}_{\bm{p},m,s},...)$ is a row vector of
electron creation operators. $\bm{p}=(p_{x},p_{y})$ run over the 1st BZ of
graphene at layer $\mu$ and the choice of the $K^{\mu}_{\pm}$ point is denoted
by valley $\tau$. $\tau_{\star}$, $\sigma_{\star}$, and $\mu_{\star}$ are the
Pauli matrices for the related flavors.

Figure 9: Two twisted graphene Brillouin zones (BZ) with twisted angle
$\theta$ and the corresponding mini Brillouin zones (MBZ) for TBG moire
lattice generated from the twist. The momentum space hexagonal lattice (i.e.,
yellow point) discretizes MBZ. $G_{1}$ and $G_{2}$ are primitive vectors of
the MBZ.

By incorporating inter-layer hopping into the tight-binding Hamiltonian and a
$U_{g}=e^{-i\theta\sigma_{z}\mu_{z}\tau_{z}/4}$ basis transformation (i.e.,
$c_{\bm{p}}\rightarrow U_{g}c_{\bm{p}}$ ), we have:

$\displaystyle\hat{H}_{BM}$
$\displaystyle=\sum_{\bm{p}}T^{\textup{intra}}_{\bm{p},\bm{p}}c^{\dagger}_{\bm{p}}c_{\bm{p}}+T^{\textup{inter}}_{\bm{p}^{\prime},\bm{p}}c^{\dagger}_{\bm{p}^{\prime}}c_{\bm{p}}$
(13) $\displaystyle T^{\textup{intra}}_{\bm{p},\bm{p}}$
$\displaystyle=\hbar\upsilon_{F}(p_{x}\sigma_{x}\tau_{z}+p_{y}\sigma_{y})$
$\displaystyle T^{\textup{inter}}_{\bm{p}^{\prime},\bm{p}}$
$\displaystyle=\frac{1}{2}\sum_{l=0,1,2}\delta_{\bm{p}^{\prime},\bm{p}-\tau_{z}\bm{q}_{l}}(\mu_{x}-i\mu_{y})T_{l}+\textup{h.c.}$
$\displaystyle T_{l}$
$\displaystyle=\omega_{0}e^{-i\theta\sigma_{z}\mu_{z}\tau_{z}/2}+\omega_{1}e^{2\pi
il\sigma_{z}\tau_{z}/3}\sigma_{x}e^{-2\pi il\sigma_{z}\tau_{z}/3}\,,$

where $\theta=1.05^{\circ}$ is the twist angle,
$\nu_{F}=2700*\frac{3a_{0}}{2\hbar}$ meV is the Dirac velocity, and
$\omega_{0}$, $\omega_{1}$ describes lattice relaxation between A-A and A-B
stacks, respectively. In this work, we take $\omega_{1}=105$meV and tune
$\kappa=\frac{\omega_{0}}{\omega_{1}}$ to scan $\omega_{0}$. Transition
momentum $\bm{q}_{l}=R(2\pi l/3)\bm{q}_{0}$, where
$R(\psi)=e^{i\psi\sigma_{y}}$ is a rotation matrix.
$\bm{q}_{0}=\frac{8\pi\textup{sin}(\theta/2)}{3\sqrt{3}a_{0}}(0,-1)^{T}$ and
$a_{0}=1.42\AA$ is the intra-layer distance between carbon atoms.

Since the model $\hat{H}_{BM}$ [3] above are described in “microscopic basis”
(i.e., unit cell index $l$, sub-lattice $\sigma$ and layer $\mu$), we
transform it to “band basis” that simplifies the implementation of AFQMC. As
demonstrated in Fig. 9, we discretize MBZ at valley $\tau$ as $L\times L$
lattice, and each lattice point is labeled by
$\bm{k}+\bm{G}=(k_{1}+m_{1}*G_{1},k_{2}+m_{2}*G_{2})$ with $\bm{k}$ running
over the 1st MBZ and $G_{1},G_{2}$ are the primitive vectors of MBZ. $\Gamma$
point of the hexagon is set to be the origin $\bm{k}=(0,0)$. We shift $\bm{p}$
with $\bm{p}+K^{\mu}_{\pm}=\bm{k}+\bm{G}+\Gamma$, then $\hat{H}_{BM}$ is
block-diagonalized for $\bm{k}$, spin $s$, valley $\tau$ where only terms that
connect $\bm{k}$ and $\bm{k}+\bm{G}$ are non-zero. In the numerical study, we
truncate $\bm{G}$ with the closest $N_{G}$ vectors ($-4\leq m_{1},m_{2}\leq 4$
in this work). Then we can use $\varepsilon_{\bm{k},i}$ and
$|\varphi_{\bm{k},i}\rangle$ to specify eigenvalues and eigenstates of each
band $n=0,...,N_{G}-1$ in this model:

$\hat{H}_{BM}|\varphi_{\bm{k},i}\rangle=\varepsilon_{\bm{k},i}|\varphi_{\bm{k},i}\rangle$
(14)

which is in band basis labeled by $i=(s,\tau,n)$ and Bloch wave functions
$|\varphi_{\bm{k},i}\rangle=\sum_{\bm{G},m}\alpha_{\bm{k}+\bm{G},i,m}|\Phi_{\bm{k}+\bm{G},m}\rangle$.

We limit this model to two flat bands closest to the Fermi level to allow a
practical numerical study. We only consider one-body and two-body terms
involving two flat bands $n=\pm$ where $\pm$ specify the top and bottom flat
bands. When $|\varphi_{\bm{k},i}\rangle$ are degenerated, we fix the gauge
with C2T symmetry as used in previous DQMC work [30].

###  VII.2 Interacting BM model

The two-body interaction is already included in the BM model at an
independent-electron level through the original ab initio density-functional
computations. To include the screened Coulomb interaction, the IBM model in
Eq. (3) also needs to include a background subtraction to correct this double-
counting effect. The Coulomb interaction $\hat{V}$ is limited in two flat
bands, and the charge-density operators are shifted according to the chemical
potential at $n_{f}$ filling:

$\displaystyle\bar{\rho}_{\bm{q}}=\sum_{\bm{k},i,j}d^{\dagger}_{\bm{k}+\bm{q},i}\Lambda(\bm{k}+\bm{q},\bm{k})_{i,j}d_{\bm{k},j}-\frac{n_{f}}{8}\\!\\!\sum_{\bm{k}}\delta_{\bm{q},\bm{G}}\textup{tr}\\{\Lambda(\bm{k}+\bm{q},\bm{k})\\}$
(15)
$\displaystyle\Lambda(\bm{k}+\bm{q},\bm{k})_{i,j}\equiv\sum_{m}\alpha^{\star}_{\bm{k}+\bm{q},i,m}\alpha_{\bm{k},j,m}$

where $d^{\dagger}_{\bm{k},i}$ create state
$|\varphi_{\bm{k},i}\rangle=\sum_{\bm{G},m}\alpha_{\bm{k}+\bm{G},i,m}|\Phi_{\bm{k}+\bm{G},m}\rangle$
and momentum $\bm{k}$ that runs over MBZ. The transfer momentum $\bm{q}$
remains unrestricted but is truncated for numerical implementation which is
$(-3*G_{1},-3*G_{2})\leq\bm{q}\leq(3*G_{1},3*G_{2})$ in this work. The
potential

$V_{\bm{q}}=\frac{e^{2}}{2\epsilon\epsilon_{0}}\frac{1}{A}\frac{\textup{tanh}(|\bm{q}|d)}{|\bm{q}|}$
(16)

is the Fourier-transformed two-gate screened Coulomb interaction.
$A=N_{M}A_{0}$ is the total area of the system with
$A_{0}=\sqrt{2}a^{2}_{M}/2$,
$a_{M}=\sqrt{3}a_{0}/2\textup{sin}\frac{\theta}{2}=13.4$ nm and $N_{M}$ is the
number of moire unite cells of this system. We chose a distance $d=20$ nm
between two layers and relative permittivity $\epsilon=10$ ($\epsilon_{0}$ is
vacuum permittivity). We study the choice of relative permittivity in App. XI.

The inclusion of the Coulomb interaction in $\hat{V}$ should be accompanied by
a removal of the double-counting of Coulomb interacting in $\hat{H}_{BM}$, as
discussed earlier. The “bare” dispersion in TBG without Coulomb interaction
involved is approximated with

$\displaystyle H_{0}$ $\displaystyle=\hat{H}_{BM}-[\hat{V}]_{\phi_{BM}}$ (17)
$\displaystyle=\sum_{\bm{k},i}\varepsilon_{\bm{k},i}c^{\dagger}_{\bm{k},i}c_{\bm{k},i}-\sum_{\bm{k},i,j}[V]_{\bm{k},i;\bm{k},j}c^{\dagger}_{\bm{k},i}c_{\bm{k},j}.$

with a mean-field approximation $[V]$ of the Coulomb interaction. The choice
of mean-field state to estimate the double counting correction
$[V]_{\bm{k},i;\bm{k},j}$ can be varied. This paper uses the ground state of
$H_{BM}$ at $n_{f}=4$ as our mean-field reference state $\phi_{BM}$. The
discrepancies between bare dispersion $H_{0}$ and BM model $H_{BM}$ along
$\kappa$ are presented in Fig. 1.

In practice, we first write $\hat{H}_{BM}$ (Eq. 13) in “microscopic basis”
where basis are labeled with spin $s$, valley $\tau$, layer $\mu$ and
sublattice $\sigma$. We then diagonalize $\hat{H}_{BM}$ for flat bands closest
to the fermi surface and identify two bands through energy. For two flat bands
that are degenerated (i.e., $\epsilon_{\bm{k}}=0$), we break the degeneracy
with C2T and “small angle approximation” [30] by taking their eigenvectors as
our specified bands. We take the set of bands and related spin, valley, as our
band basis (labeled by $n$, $s$, $\tau$) where $H_{BM}$ is diagonalized.

##  VIII Symmetry in IBM model

According to a previous study [13, 14], at the spinless limit, the ground
states of IBM (without double counting correction) at chiral and flat band
limit are single determinants:

  * •

( Chern number $\pm 2$) fill two bands in the same Chern sector
$\longrightarrow$ quantum anomalous Hall states (QAH $\sigma_{z}\tau_{z}$).

  * •

( Chern number $0$) fill as linear combinations of the two bands in both Chern
sectors $\longrightarrow$ Valley Polarized states(VP $\tau_{z}$), Valley Hall
states(VH $\sigma_{z}$), T-symmetric Inter-valley-coherent states (T-IVC
$\sigma_{x}\tau_{x,y}$)and Kramers inter-valley-coherent states(K-IVC
$\sigma_{y}\tau_{x,y}$).

The limit away from flat band can be considered perturbatively as an emergence
of dispersion term:

$h(\bm{k})=h_{0}(\bm{k})\tau_{z}+h_{x}(\bm{k})\sigma_{x}+h_{y}(\bm{k})\sigma_{y}\tau_{z}$

that tunnel between opposite Chern (sublattice) bands belonging to the same
valley, and we have

  * •

At $\kappa=0$, QAH, VH, and K-IVC are exact ground states if there is no
dispersion. If dispersion is non-zero, the perturbation of the dispersion
indicates the candidates of ground state prefer QAH, VH, and K-IVC.

  * •

At $\kappa>0$, the combination of interaction terms (an extra sublattice
exchange term should be considered when $\kappa>0$) and the perturbation of
the dispersion prefer K-IVC.

HF [13, 26, 28] can produce the above correlated insulating states
respectively by imposing corresponding symmetry breaking in the initial
density. If the dispersion term is small, the breaking of $U_{V}(1)$ valley
charge conservation indicates K-IVC state and the unbroken of $U_{V}(1)$
valley charge conservation (no inter-valley terms) combined with C2T symmetry
helps us to distinguish VH and QAH (C2T symmetry is broken) with VP and SP
state (C2T symmetry is preserved). However, for a non-zero dispersion term,
especially when perturbation fails, the solution of HF has no need to be
eigenstates of QAH, VH, K-IVC correlation (i.e., exactly VH, QAH, and K-IVC
state at flat band limit). In this paper, the term “K-IVC state” specifies HF
ground state that have non-zero IVC, and “C2T-symmetry-breaking state”
summarize all HF ground states that have no IVC and breaks C2T symmetry (i.e.,
QAH and VH) as our study does not distinguish VH and QAH.

##  IX Low-lying states around chiral limit for $n_{f}=4,2$ fillings

###  IX.1 Spin gaps

We take the chiral limit for $n_{f}=4,2$ fillings. To determine the spin
polarization, we calculate the spin gaps

$\Delta_{s}=E(N_{\uparrow}-\delta N,N_{\downarrow}+\delta
N)-E(N_{\uparrow},N_{\downarrow})\,,$ (18)

where $\delta N=1,2,\cdots$. In Fig. 10, we show the spin gaps as a function
of spin flips ($\delta N$) at $\kappa=0.0$ from both HF and AFQMC. The
reference states, with zero spin-flip, are the states with the lowest HF
energy, namely, the non-spin polarized ground state
($N_{\uparrow}=N_{\downarrow}=N/2$) at $n_{f}=4$ and fully spin-polarized
state ($N_{\uparrow}=N$, $N_{\downarrow}=0$) for $n_{f}=2$ fillings. The AFQMC
calculations used the corresponding HF state as a trial wave function. We show
results from several lattice sizes, up to $12\times 12$, to help gauge finite-
size effects. Our data support a vanishing spin gap around zero spin-flip.

###  IX.2 Constraint robustness and further correction with constraint release

Many studies [58, 56, 57] have shown that AFQMC can recover from a
constraining trial wave function which is incorrect. However there has been no
previous application of AFQMC to the MATBG system. Here we illustrate how in
the most ambiguous cases we can reduce the dependence on which low-lying
states as the trial wave function by constraint release.

Figure 10:  Spin gap at the chiral limit ($\kappa=0.0$) for fillings $n_{f}=4$
(left panel) and $n_{f}=2$ (right panel). $\Delta_{s}$ is shown as a function
of $\delta N$ for a small number of spin flips. States with zero spin-flip
corresponds to the lowest HF energy. Dashed lines connect results from the
largest lattice sizes $12\times 12$.

To distinguish low-lying states that are closely competing (i.e., K-IVC state
and C2T-symmetry-breaking state, which are degenerated at the chiral limit in
the HF study), we detect ground state preference between these two states by
comparing the energy difference between AFQMC with corresponding trial wave
functions. That is the AFQMC run with lower energy gives a better description
of the ground state. To resolve such a small difference, AFQMC with Metropolis
release constraint [60] is performed to detect – and eliminate – any residual
systematic bias from the constrained path approximation. AFQMC with Metropolis
release can exponentially reduce the systematic bias by increasing released
imaginary time $\beta$. In Fig. 12, we show the detailed convergence in energy
for $\kappa=0$. In Fig. 11, we summarise the comparison of converged energy in
AFQMC with Metropolis release for different trial wave functions and different
$\kappa$. We show that AFQMC with K-IVC trial wave function has lower energy
for both $n_{f}=4,2$ fillings. Note that the competition between the ground
states in K-IVC and C2T-symmetry-breaking manifolds is only relevant when
$\kappa$ is small. The energy of ground states in K-IVC manifolds is much
lower than the others for $\kappa\gg 0$.

Figure 11:  Energy difference per electron between the ground states in K-IVC
and C2T-symmetry-breaking manifolds at different $\kappa$ and $6\times 6$
lattice for $n_{f}=4$ fillings (left) and $n_{f}=2$ fillings (right). AFQMC
proceeds with the corresponding HF state as trial/initial wave functions.
AFQMC with Metropolis release is performed to eliminate minor constraint bias.

Figure 12:  Quantifying and reducing the constrained path error by AFQMC with
Metropolis release constraint [60]. Energy per electron vs. released imaginary
time $\beta$ at $\kappa=0$ for charge neutrality $n_{f}=4$ (left) and half-
filling $n_{f}=2$ (right). Energy from AFQMC with K-IVC state or C2T-symmetry-
breaking state as the trial wave function is identified with orange and blue
lines, respectively. Dashed lines indicate the converged energy for the
Metropolis release procedure.

##  X order parameter

In AFQMC we can use the Hellman-Feynman theorem to calculate an order
parameter $\langle\hat{O}\rangle$ by taking the derivative of energy over
$\lambda$, [63]:

$\langle\hat{O}\rangle=\underset{\lambda\to 0}{\lim}\frac{\mathrm{d}\langle
H_{IBM}+\lambda\hat{O}\rangle}{\mathrm{d}\lambda}.$ (19)

For a small range of $\lambda$ around zero, we then use the numeric derivative
of the energy to estimate the average of observables in the ground state.

In Fig. 13, we show more details about our calculations for K-IVC order
parameter $\tau_{x}n_{x}$ and band occupancy $n_{z}$. To detect the K-IVC
order, there is a two-fold degeneracy of the K-IVC pairs that favor
$\tau_{x}n_{x}$ and $-\tau_{x}n_{x}$ (i.e., a direct estimation of
$\tau_{x}n_{x}$ in the ground state without breaking the symmetry will always
be zero). As presented in the top figures of Fig. 13, adding a positive
$\lambda$ will favor one of them to get lower energy and disfavor the other,
which gives higher energy. We can produce a HF solution for one vs the other
by imposing a corresponding initial density matrix. For $\kappa=0$ (i.e., in
the metallic phase), AFQMC with different HF trial wave functions gives
similar energy with zero derivatives around $\lambda=0$, supporting a loss of
K-IVC order. As $\kappa$ increases to the insulating phase, the linear
behavior and the increase of derivatives in AFQMC indicate the appearance of
K-IVC order in the ground state. There is no such degeneracy in the
calculation of band occupancy $n_{z}$, as the top and bottom bands have a non-
zero gap. We observed a consistent linear behavior across $\lambda=0$, and the
derivatives decrease as $\kappa$ increases, which indicates a reduction of the
gap between the top and bottom band.

Figure 13: Energy per electron versus $\lambda$ for K-IVC order parameter
$\tau_{x}n_{x}$ (top) and band occupancy $n_{z}$ (bottom) . HF results are
presented in blue, while AFQMC results are highlighted in orange. Dashed
lines/curves are linear/quadratic fitting to guide eyes. In the top figure,
blue star and circle marks identify the energy of HF states that prefer
$\tau_{x}n_{x}$ or $-\tau_{x}n_{x}$. The AFQMC results with corresponding HF
states as trial wave functions are denoted by orange stars and circle marks.
The derivative along $\lambda$ indicates the preference of the corresponding
order.

##  XI  The choice of the permittivity $\epsilon$

We explore the effects of the relative permittivity $\epsilon$. In our study,
we chose $\epsilon\approx 10$. To briefly verify the effects of different
$\epsilon$ in determining the ground state, we present a direct comparison of
two $\epsilon$ in terms of charge gap at $n_{f}=2$ filling in Fig. 14. Our
results show that by choosing a smaller $\epsilon=5$, the charge gap is almost
scaled by two, but the metal-insulator transition is still maintained.

These tests are obviously limited in scope and cannot rule out the possibility
that some choice of $\epsilon$ might qualitatively change the ground state.

Figure 14: Charge gap for $n_{f}=2$ fillings versus $\kappa$ for different
choice of $\epsilon$. Charge gaps are calculated by measuring the energy
required to add or remove an electron from the system. AFQMC proceeds with the
corresponding HF state as the trial wave function, the HF charge gap of which
is illustrated in the figure. Shaded pink and green backgrounds separate the
metallic/semi-metallic and insulating phases.