11institutetext: Faculty of Mathematics and Information Science, Warsaw
University of Technology, Koszykowa 75, 00-662 Warszawa, Poland
11email<EMAIL_ADDRESS>22institutetext: Systems Research Institute Polish Academy of Sciences
22email<EMAIL_ADDRESS>
# Using adversarial images to improve outcomes of federated learning for non-
IID data

Anastasiya Danilenka 11 0000-0002-3080-0303 Maria Ganzha 11
0000-0001-7714-4844 Marcin Paprzycki 22 0000-0002-8069-2152 Jacek Mańdziuk 11
0000-0003-0947-028X

###### Abstract

One of the important problems in federated learning is how to deal with
unbalanced data. This contribution introduces a novel technique designed to
deal with label skewed non-IID data, using adversarial inputs, created by the
I-FGSM method. Adversarial inputs guide the training process and allow the
Weighted Federated Averaging to give more importance to clients with
“selected” local label distributions. Experimental results, gathered from
image classification tasks, for MNIST and CIFAR-10 datasets, are reported and
analyzed.

###### Keywords:

Federated learning Adversarial attack Non-IID data.

##  1 Introduction

The concept of federated learning (FL) was introduced in [10]. It aims at
creating a shared/global model, combining information from multiple computing
nodes, without sharing individual source datasets. In FL, training proceeds in
rounds. In each round, clients receive the current version of the global
model, perform training using local data, and return the updated model
parameters. Thus, the clients in FL pipelines act as both data sources and
computational nodes, while the learning process is controlled, and
orchestrated, by the server/manager (for needed details see [10, 3]).

The process of combining updates, received from clients, into the new version
of the shared model, is a subject of intensive research. However, its most
straightforward realization averages all model parameters, returned by all
participating clients. This algorithm (FedAvg) has been summarized as
Algorithm 1 (see [7] for additional details).

Algorithm 1 Federated Averaging (FedAvg) algorithm. The number of server-side
epochs is $N$, the number of clients participating in the training round is
$K$, $m$ are parameters of the model, $s$ is a set of parameters needed for a
client to train the model (learning rate, number of local epochs, batch-size,
model-dependent parameters)

global model initialization

clients ready for training

for $i\neq N$ do

for $client_{k}$ in $clients[1..K]$ in parallel do

$m_{i}^{client_{k}}\leftarrow RequestModelUpdate(m_{i},s_{i})$

end for

$m_{i}\leftarrow{\frac{1}{K}}\times\sum_{k=0}^{K}m_{i}^{client_{k}}$

end for

A natural extension to FedAvg is weighting, i.e. assigning individual
importance to each client, based on some additional measure/knowledge. For
example, clients with more local data can be treated as delivering more
valuable input to the training process, and thus receive larger weights for
updates they “suggest” [10].

Federated learning pipeline utilizes client devices’ local datasets, while the
complete dataset is never gathered “in one place”. Therefore, statistical
properties of both, the local datasets and of the global dataset, are unknown.
As a result, data may turn out to be not identically independently distributed
(non-IID). It has been shown that non-IID datasets can (very) negatively
affect the quality of a global model produced by the FL [18, 19, 17].

From “statistical point of view”, non-IID data can be divided into categories,
based on the source of heterogeneity [19], i.e.: (1) data quantity skew (local
datasets differ in size), (2) label distribution skew (different devices have
different subsets of labels inside local dataset or have data relying on
preference), (3) attribute skew (local data has unique characteristic
features, noise, perturbations, etc.), (4) temporal skew (local data
distributions differ over time or data was collected in different time
periods). Obviously, combinations of skews can materialize in local datasets.

This work concentrates on non-IID datasets with label distribution skew. This
type of skew naturally occurs in (image) classification tasks. For instance,
individual devices may or may not have the full set of labels in their local
datasets. Depending on how many labels are represented inside the local
dataset, the skew can be extreme with only $1$ label per local dataset, or a
($C-1$)-label skew, where $C$ is the number of labels in classification task,
thus, assuming that each local dataset lacks data for one of the labels.

In this context, a novel method, to overcome the problems caused by label
distribution skew in non-IID data, called Adversarial Federated Learning
(AdFL), is proposed. It is conceptually based on adversarial attacks and is
applied to neural network based training with image-based data. Thus, we
proceed as follows. In Section 2, we introduce key concepts of adversarial
attacks. We follow, in Section 3, with the description of the AdFL method.
Section 4 presents the experimental setup and results of experiments,
conducted for MNIST [5] and CIFAR-10 [4] datasets. Conclusions and directions
for future research are covered in the last section.

It should be noted that, due to the space limitation, and the fact that a
number of fresh and comprehensive literature reviews already exist for
federated learning (see, for instance, [19, 3, 6]), we omit the standard
state-of-the-art discussion. Only references directly needed, in the context
of this contribution, are cited.

##  2 Adversarial inputs

In 2014 it was shown that a deep neural network can be misled and inclined to
make mistakes when classifying slightly modified images that still remain very
easy to categorize by humans. Moreover, these images have been initially
correctly classified by the model itself [13]. Specifically, it was shown that
by altering the original input image with small non-random perturbations, one
can create an “adversarial image” that is indistinguishable from the original
one by humans, but forces a classifier to make a mistake.

Overall, adversarial attacks can be divided into two categories: targeted (an
initial image is modified in order to be misclassified as a certain class),
and untargeted (modifications are made to mislead the model and misclassify
the image as “any” incorrect class) [8].

The source image can be modified using various techniques. For example,
methods belonging to the family of gradient-based techniques, like Fast
Gradient Sign Method (FGSM), or Projected Gradient Method (PGM) [1, 9] can be
used. These methods rely on having an access to model’s gradients. FGSM can be
described using the following formula (cf. [1] for more details):

$\eta=\epsilon sign(\nabla_{x}L(\theta,x,y))$

Here, $x$ stands for an initial input, $y$ – for the input’s label (for tasks
that assume labels), $\theta$ – model’s parameters, $L(\theta,x,y)$ – a loss
function (for example, cross entropy for multi-class classification task),
$\nabla_{x}$ – a gradient computed with back propagation, $\epsilon$ – a
constant chosen to control the amount of perturbations added to the input.

Moreover, the property of transferability of adversarial inputs was discovered
[13]. Specifically, it has been established that adversarial inputs generated
by one model can still mislead another similar model, even if their training
subsets are non-overlapping. Conducted studies showed that adversarial inputs
are most likely derived from the adversarial subspaces [14]. A number of
contributions has been devoted to discovering ways of improving
transferability, understanding it’s nature and boundaries [16, 15, 2]. The
transferability property of adversarial images also provides a possibility for
decision-based attacks. Here, a local model is trained on a data subset and
tries to learn how to mimic the decision boundary of the model that is a
target of the attack [12]. Overall, the success rate for adversarial attacks,
depending on various conditions, is in the range from 63% to 100%, including
cross-model attacks scenarios (where the source model and the target model
have different architectures) [11].

In what follows, features of adversarial attacks, i.e.: (1) transferability,
and (2) connectivity with decision boundary of the model, will be used in the
context of FL pipeline to deal with label distribution skew non-IID data.

##  3 Description of the proposed AdFL algorithm

Overall, AdFL utilizes transferability of adversarial inputs, and uses them to
improve quality of FL. Here, adversarial inputs are created not as threats to
misguide models into making mistakes, but are instantiated to measure the
“cohesion” among federated clients. This allows to modify the training process
without sharing local data.

As stated in Section 2, for gradient based attacks, having only an access to
the model parameters would be sufficient to successfully perform a targeted
adversarial attack. Interestingly, an access to the model parameters is
exactly what the server in the federated learning pipeline is granted, while
knowing nothing about the local (client) data. In particular, it is possible
for the server to use updates of individual model parameters, returned by
clients, to produce an adversarial image for every label that is defined in
given classification task. Adversarial images, in this case, serve as a
“communication tool” for updated models to check if the produced adversarial
samples are able to fit inside decision boundaries of the updated models.

Without access to the source data, images consisting of noise, or fully
black/white images, are used as the starting point for the adversarial input
generation. Examples of images produced from a totally black starting image
for MNIST and CIFAR-10 datasets are shown in Figure 1.

Figure 1: Adversarial images produced for MNIST and CIFAR-10 datasets starting
from a totally black image.

Note that, the adversarial image generation step occurs on the server side,
i.e. for the client devices no changes to the default FL pipeline are needed
(no additional information or communication is required). The clients operate
“as always”, i.e. they participate in the training of the global model, using
their local data, and send the resulting update proposals back to the server.
The enhanced algorithm only changes the way the server behaves, after model
updates are gathered from the client devices. Specifically, the server
behaviour can be summarized as follows:

Figure 2: Server-side training epoch routine in AdFL. Green block indicates
starting point of each epoch. Validation and testing routines are omitted.

  * •

Based on obtained parameter updates, the server recreates local models for
each client.

  * •

Obtained models are used to generate images for each of the labels, according
to the iterative version of the Fast Gradient Signed Method (FGSM). Starting
from a random noise image, it iteratively modifies the image, using gradients
retrieved from the client’s model, until the specified number of iterations is
reached.

  * •

Individual models make predictions of all adversarial inputs generated by all
models, including the ones generated “for itself”. Therefore, each client’s
model makes $K$ * $C$ predictions, where $K$ – is the number of clients
participating in a given training round and $C$ – is the number of labels
involved in given classification task.

  * •

Based on the results from the previous step, each client receives a weight
that serves as a measure of its ability to predict adversarial inputs produced
by other models, as well as ability to produce such images. The evaluation
process, used to calculate weigth, is described in Section 3.1.

  * •

A weighted version of the FedAvg algorithm (WFedAvg) is used to generate a
global model from the local models, where weighting takes place according to
the values calculated during the previous step.

The server-side algorithm is summarized in Algorithm 2. For clarity, the
detailed descriptions of adversarial images generation, and cross-prediction,
are omitted, leaving only the high-level representation of steps that are
introduced by the proposed algorithm. The diagram featured in Figure 2 shows
how each server-side epoch is designed in AdFL.

Algorithm 2 AdFL algorithm. The number of server-side epochs is $N$, the
number of clients participating in the training round is $K$, $m$ are
parameters of the model, $s$ is a set of parameters needed for a client to
train the model (learning rate, number of local epochs, batch-size, model-
dependent parameters), $w$ are weights of the updated versions of global
model, $C$ is the number of labels in the considered task

global model initialization

clients ready for training

for $i\neq N$ do

for $client_{k}$ in $clients[1..K]$ in parallel do

$m_{i}^{client_{k}}\leftarrow RequestModelUpdate(m_{i},s_{i})$

end for$\triangleright$ Produce adversarial images and use them for prediction

for $client_{k}$ in $clients[1..K]$ in parallel do

for $label_{c}$ in $labels[1..C]$ in parallel do

$images_{label_{c}}^{client_{k}}\leftarrow
CreateAdversarialImage(m_{i}^{client_{k}},label_{c})$

for $client_{j}$ in $clients[1..K]$ in parallel do

$predictions_{label_{c}}^{client_{j}}\leftarrow
PredictAdversarialImage(m_{i}^{client_{j}},images_{label_{c}}^{client_{k}})$

end for

end for

end for$\triangleright$ Measure each model performance and use it for WFedAvg

$w\leftarrow DefineModelsWeights(predictions)$

$m_{i}\leftarrow\sum_{k=0}^{K}m_{i}^{client_{k}}\times w_{k}$

end for

###  3.1 Weights calculation

As stated above, WFedAvg is to be used on the basis of predictions of
adversarial images for all classes and for all clients. However, it should be
noted that in this process two aspects of image recognition materialize: (a)
How well a given client recognizes adversarial images generated by other
models (clients)? and (b) How well adversarial images that are produced by a
given model (client) are recognized by the other models (clients)? Therefore,
the weight assigned to updates originating from a given client will consist of
two parts, corresponding to these two aspects.

####  3.1.1 The ability of the updated model to predict adversarial images
produced by other models

is measured according to Equation (1), where each multiplication consists of a
binary flag indicating whether or not the prediction for label $c$ generated
by model $k$ was correct and the probability returned by the model. The
results obtained for the model predicting its own adversarial inputs are
omitted.

$\begin{split}\mbox{predicted
others}=\sum_{k=1}^{K}{\sum_{c=0}^{C}\mbox{prediction
correct${}_{k,c}$}*{\mbox{returned probability${}_{k,c}$}}}\end{split}$ (1)

####  3.1.2 The ability of the updated model to produce adversarial images
that are recognized by the other models

is defined by Equation (2). Again, results for the model “predicting itself”
are omitted. The formula consists of the binary flag, indicating whether the
current model input for label $c$ was correctly predicted by model $k$, and
the probability assigned by model $k$ while predicting label $c$.

$\begin{split}\mbox{was
predicted}=\sum_{k=1}^{K}{\sum_{c=0}^{C}\mbox{prediction
correct${}_{k,c}$}*{\mbox{returned probability${}_{k,c}$}}}\end{split}$ (2)

####  3.1.3 The total weight

is a result of a simple summation of the scores from Equation (1) and Equation
(2). So, the absolute weight is calculated as:

$\begin{split}\mbox{weight total}=\mbox{predicted others}+\mbox{was
predicted}\end{split}$ (3)

After the total weights for all clients that participated in the training
epoch, are calculated they are normalized. The final weights for the current
epoch are used by the WFedAvg algorithm.

##  4 Experimental setup and results

###  4.1 Experimental setup

Initial set of experiments has been conducted using two standard benchmark
datasets – MNIST [5] and CIFAR-10 [4]. The main interest lied in exploring the
performance of the proposed AdFL method, compared to the default FL (with no
modifications to the server-side routine). The experiments structure was
designed as follows:

  * •

In each server-side epoch 5 clients were picked from the training pool.

  * •

Clients performed training for 5 epochs and returned results back to the
server.

  * •

A hold-out IID test dataset was stored on the server, and used to track IID
accuracy and per-class accuracy for the global model. The test set was
composed of samples that are not presented in any clients’ local datasets.

  * •

Validation clients were used to modify the learning rate during training.
Tracking average validation clients loss was done to adjust the learning rate
during the training process, by reducing the learning rate in inverse square
root of the number of epochs passed manner.

Note that the aim of experiments was not to find the set of best
“hyperparameters” and the most effective “experimental setup”. Instead, the
focus was on proving the potential of AdFL and the premises for its further
exploration. This was also the reason for using classic benchmark datasets,
which are “well understood” and simple enough to focus on the effects of the
proposed approach, rather than having them “hidden” due to the complexty of
the dataset itself.

Table 1: Label probability during experiment Label |  0 |  1 |  2 |  3 |  4 |  5 |  6 |  7 |  8 |  9
---|---|---|---|---|---|---|---|---|---|---
Probability % |  3.5 |  4.5 |  10 |  21 |  21 |  20 |  10 |  4.5 |  3.5 |  2

Each experiment was conducted multiple times (from 10 to 15 runs), with
different model initialization and random seeds. The results were averaged in
terms of performance. In all cases, the average standard deviation per label,
at the end of training process, was of order of 0.5%.

The datasets, generated for the experiments, featured non-IID label skew that
was introduced to the local datasets by different label occurrence
probabilities. Example of such probability set is given in Table 1. The
probability was used to generate a subset of labels that were presented on the
“local device”. During experiments with MNIST, only 3 labels were presented in
individual local datasets, while for CIFAR-10 experiments the number of unique
labels was 4.

###  4.2 Experimental results and their analysis

In the first MNIST experiment, the label distribution histogram, illustrating
how many clients had the specific label inside their local dataset, the IID
test accuracy reported from the holdout dataset, and per class accuracy for
the holdout IID dataset, were calculated. The results are shown in Figure 3.
The experiment was conducted with several labels being overrepresented (more
than 50% of clients had these labels in local data), these are labels 3, 4,
and 5. Labels 1, 2, 6, and 7 were mildly underrepresented (15 to 30% of client
devices had them), and labels 0, 8, and 9 were highly underrepresented.

Figure 3: Results of the first MNIST label skew experiment.

The total IID accuracy difference between the AdFL and the baseline FL is
around 1.5%. Although a group of labels (0, 1, 4 - 7, and 9) is predicted with
the same accuracy by both algorithms, two labels have notable difference in
performance: the 3rd label predicted more acurately by 8-11%, and the gain in
accuracy for the 8th label is also around 5-10%.

The results of the second MNIST experiment, in which the distribution of
labels was shifted towards labels 0, 1, and 2, are shown in Figure 4.

Figure 4: Results of the second MNIST label skew experiment.

The difference in total IID test accuracy between federated pipeline without
modifications and the AdFL approach is 1.5-2%. As previously, two
underrepresented labels “performed better”: the accuracy improvement for the
3rd label is around 5 to 9%, and the respective improvement for the 5th label
is in range from 10 to 14%.

The final MNIST experiment featured label distribution being skewed towards
labels 7, 8, and 9, with the rest being less populated on local devices
(Figure 5).

Figure 5: Results of the third MNIST label skew experiment.

As previously, total IID accuracy improvement was around 3% over baseline FL,
and the source of this improvement came from the 1st label, with accuracy
increase in range of 9 to 12%, and the 2nd label, with gain 11-13%.

In summary, for MNIST experiments, the AdFL accuracy improvement came from
mildly underrepresented labels (15-30% of local clients have the label inside
the local dataset) with approximately similar performance for other labels, as
compared to the baseline FL method.

A similar behaviour can be spotted for the CIFAR-10 dataset. Figure 6 shows
the results of one of the conducted experiments. Label distribution in the
experiment is not as severe as in the MNIST experiments – underrepresented
labels are presented in less than 30% of local clients.

Figure 6: CIFAR-10 label skew experiment.

The difference in total IID accuracy between the two compared algorithms is
around 1%, with two labels being predicted better by AdFL: the 1st label
improvement is about 2%, and the 3rd label improvement equals 3.5-4%.

##  5 Concluding remarks

The initial experiments with the Adversarial Federated Learning (AdFL)
algorithm, proposed in this paper, show accuracy improvement on two standard
benchmark sets with non-IID label distribution. The accuracy on individual
labels raises from 2 to 14% depending on the specific distribution scheme.
Obviously, more experiments and a closer look at their results are required.

The main areas of future research include (1) verification of the AdFL
efficacy on more complex datasets. As the experiments presented in this work
were performed with a relatively small number of clients (30), (2) more
populated experiments with more nodes and similar non-IID data partitions are
also needed to confirm the initial conclusions. Using (3) new datasets will
help to understand if the achieved results are data-independent, find cases
where adversarial inputs fail to transfer from one client to another and
explore new behavior patterns of the algorithm.

####  5.0.1 Acknowledgements

Research funded in part by the Centre for Priority Research Area Artificial
Intelligence and Robotics of Warsaw University of Technology within the
Excellence Initiative: Research University (IDUB) programme.

## References

  * [1] Goodfellow, I.J., Shlens, J., Szegedy, C.: Explaining and harnessing adversarial examples (2014). https://doi.org/10.48550/ARXIV.1412.6572, https://arxiv.org/abs/1412.6572
  * [2] Huang, Q., Katsman, I., He, H., Gu, Z., Belongie, S., Lim, S.N.: Enhancing adversarial example transferability with an intermediate level attack (2019). https://doi.org/10.48550/ARXIV.1907.10823, https://arxiv.org/abs/1907.10823
  * [3] Kairouz, P., McMahan, H.B., Avent, B., Bellet, A., Bennis, M., Bhagoji, A.N., Bonawitz, K.A., Charles, Z., Cormode, G., Cummings, R., D’Oliveira, R.G.L., Rouayheb, S.E., Evans, D., Gardner, J., Garrett, Z., Gascón, A., Ghazi, B., Gibbons, P.B., Gruteser, M., Harchaoui, Z., He, C., He, L., Huo, Z., Hutchinson, B., Hsu, J., Jaggi, M., Javidi, T., Joshi, G., Khodak, M., Konečný, J., Korolova, A., Koushanfar, F., Koyejo, S., Lepoint, T., Liu, Y., Mittal, P., Mohri, M., Nock, R., Özgür, A., Pagh, R., Raykova, M., Qi, H., Ramage, D., Raskar, R., Song, D., Song, W., Stich, S.U., Sun, Z., Suresh, A.T., Tramèr, F., Vepakomma, P., Wang, J., Xiong, L., Xu, Z., Yang, Q., Yu, F.X., Yu, H., Zhao, S.: Advances and open problems in federated learning. CoRR abs/1912.04977 (2019), http://arxiv.org/abs/1912.04977
  * [4] Krizhevsky, A.: Learning multiple layers of features from tiny images (2009)
  * [5] LeCun, Y., Cortes, C.: MNIST handwritten digit database (2010), http://yann.lecun.com/exdb/mnist/
  * [6] Li, L., Fan, Y., Tse, M., Lin, K.Y.: A review of applications in federated learning. Computers & Industrial Engineering 149, 106854 (2020). https://doi.org/https://doi.org/10.1016/j.cie.2020.106854, https://www.sciencedirect.com/science/article/pii/S0360835220305532
  * [7] Li, T., Sahu, A.K., Zaheer, M., Sanjabi, M., Talwalkar, A., Smith, V.: Federated optimization in heterogeneous networks (2018). https://doi.org/10.48550/ARXIV.1812.06127, https://arxiv.org/abs/1812.06127
  * [8] Li, Y., Cheng, M., Hsieh, C.J., Lee, T.C.M.: A review of adversarial attack and defense for classification methods. The American Statistician pp. 1–17 (jan 2022). https://doi.org/10.1080/00031305.2021.2006781, https://doi.org/10.1080%2F00031305.2021.2006781
  * [9] Madry, A., Makelov, A., Schmidt, L., Tsipras, D., Vladu, A.: Towards deep learning models resistant to adversarial attacks (2017). https://doi.org/10.48550/ARXIV.1706.06083, https://arxiv.org/abs/1706.06083
  * [10] McMahan, H.B., Moore, E., Ramage, D., Hampson, S., y Arcas, B.A.: Communication-efficient learning of deep networks from decentralized data (2017)
  * [11] Nowroozi, E., Mekdad, Y., Berenjestanaki, M.H., Conti, M., Fergougui, A.E.: Demystifying the transferability of adversarial attacks in computer networks (2021). https://doi.org/10.48550/ARXIV.2110.04488, https://arxiv.org/abs/2110.04488
  * [12] Papernot, N., McDaniel, P., Goodfellow, I., Jha, S., Celik, Z.B., Swami, A.: Practical black-box attacks against machine learning. In: Proceedings of the 2017 ACM on Asia Conference on Computer and Communications Security. p. 506–519. ASIA CCS ’17, Association for Computing Machinery, New York, NY, USA (2017). https://doi.org/10.1145/3052973.3053009, https://doi.org/10.1145/3052973.3053009
  * [13] Szegedy, C., Zaremba, W., Sutskever, I., Bruna, J., Erhan, D., Goodfellow, I., Fergus, R.: Intriguing properties of neural networks (2013). https://doi.org/10.48550/ARXIV.1312.6199, https://arxiv.org/abs/1312.6199
  * [14] Tramèr, F., Papernot, N., Goodfellow, I., Boneh, D., McDaniel, P.: The space of transferable adversarial examples (2017). https://doi.org/10.48550/ARXIV.1704.03453, https://arxiv.org/abs/1704.03453
  * [15] Wang, R., Guo, Y., Yang, R., Wang, Y.: Exploring transferable and robust adversarial perturbation generation from the perspective of network hierarchy (2021). https://doi.org/10.48550/ARXIV.2108.07033, https://arxiv.org/abs/2108.07033
  * [16] Wu, L., Zhu, Z., Tai, C., Ee, W.: Understanding and enhancing the transferability of adversarial examples (02 2018)
  * [17] Xiao, C., Wang, S.: An experimental study of class imbalance in federated learning (2021)
  * [18] Zhao, Y., Li, M., Lai, L., Suda, N., Civin, D., Chandra, V.: Federated learning with non-iid data. CoRR abs/1806.00582 (2018), http://arxiv.org/abs/1806.00582
  * [19] Zhu, H., Xu, J., Liu, S., Jin, Y.: Federated learning on non-iid data: A survey. CoRR abs/2106.06843 (2021), https://arxiv.org/abs/2106.06843


††thanks: These authors contributed equally to this work.††thanks: These
authors contributed equally to this work.

# Origin of magic angles in twisted bilayer graphene: The magic ring

Wei-Chen Wang  Institute of Atomic and Molecular Sciences, Academia Sinica,
Taipei 106319, Taiwan  Department of Physics, National Taiwan University,
Taipei 106319, Taiwan  Feng-Wu Chen  Institute of Atomic and Molecular
Sciences, Academia Sinica, Taipei 106319, Taiwan  Kuan-Sen Lin  Institute of
Atomic and Molecular Sciences, Academia Sinica, Taipei 106319, Taiwan
Department of Physics, National Taiwan University, Taipei 106319, Taiwan
Justin T. Hou  Institute of Atomic and Molecular Sciences, Academia Sinica,
Taipei 106319, Taiwan  Department of Physics, National Taiwan University,
Taipei 106319, Taiwan  Ho-Chun Lin  Institute of Atomic and Molecular
Sciences, Academia Sinica, Taipei 106319, Taiwan  Department of Physics,
National Taiwan University, Taipei 106319, Taiwan  Mei-Yin Chou
<EMAIL_ADDRESS>Institute of Atomic and Molecular Sciences,
Academia Sinica, Taipei 106319, Taiwan  Department of Physics, National Taiwan
University, Taipei 106319, Taiwan

###### Abstract

The unexpected discovery of superconductivity and strong electron correlation
in twisted bilayer graphene (TBG), a system containing only sp electrons, is
considered as one of the most intriguing developments in two-dimensional
materials in recent years. The key feature is the emergent flat energy bands
near the Fermi level, a favorable condition for novel many-body phases, at the
so-called “magic angles”. The physical origin of these interesting flat bands
has been elusive to date, hindering the construction of an effective theory
for the unconventional electron correlation. In this work, we have identified
the importance of charge accumulation in the AA region of the moiré supercell
and the most critical role of the Fermi ring in AA-stacked bilayer graphene.
We show that the magic angles can be predicted by the moiré periodicity
determined by the size of this Fermi ring. The resonant criterion in momentum
space makes it possible to coherently combine states on the Fermi ring through
scattering by the moiré potential, leading to flat bands near the Fermi level.
We thus establish the physical origin of the magic angles in TBG and identify
the characteristics of one-particle states associated with the flat bands for
further many-body investigations.

††preprint: APS/123-QED

##  I Introduction

Two graphene sheets with a small rotation angle $\theta$ between them creates
a special system of twisted bilayer graphene (TBG), in which a moiré pattern
emerges with a spatial periodicity inversely proportional to the twist angle
[1, 2]. It has been shown that the variation of the twist angle in moiré two-
dimensional (2D) materials could modify the electronic properties, giving rise
to an interesting research area of twistronics [3]. Previous calculations on
TBG using different theoretical methods including the low-energy continuum
model [4, 2, 5], the tight-binding method [1, 6, 7, 8], and density functional
theory (DFT) [9] all had similar findings: as $\theta$ decreases, the Fermi
velocity of the linear bands also becomes progressively smaller. This results
from an enhanced interlayer hybridization between the two rotated Dirac cones
of different layers as they get closer in momentum space when $\theta$ gets
smaller. These calculations also revealed that as the bandwidth of the low-
lying bands is reduced at small $\theta$ angles, the corresponding electron
charge becomes localized to the AA-stacked region in the moiré supercell [2,
7, 9]. One would anticipate that the Fermi velocity should decrease smoothly
as $\theta$ approaches zero; however, it was found in these calculations that
at certain so-called “magic angles” the velocity deviates away from the
expected monotonic curve and drops to almost zero [5, 2, 8], with the first of
these magic angles being around $1.1^{\circ}$ . The resulting unusually flat
bands thus have a large ratio between the Coulomb and kinetic energies and
favor the formation of possible many-body phases.

Experimentally, with the “tear-and-stack” technique [10, 11, 12, 13, 14],
researchers have been able to control the rotation angle of TBG precisely. It
was first discovered by Y. Cao et al. [15] that around the first magic angle
there exist insulating gaps when the doping level is tuned to integer
fillings. Moreover, superconducting states appear [16] in the intermediate
filling regions. These surprising reports stimulated further transport
measurements for detailed studies on correlated insulating states at several
integer fillings [17], superconductivity domes in a wide, continuous range of
doping [17, 18, 19, 20, 21], ferromagnetism [17, 22], and the integer quantum
anomalous Hall effect [23] at certain odd fillings. The existence of flat
bands in TBG was confirmed by scanning tunneling spectroscopy (STS) studies
[24, 25, 26] that revealed two pronounced van Hove singularities very close to
each other in the density of states (DOS) near the first magic angle . Using
angle-resolved photoemission spectroscopy with nanoscale resolution (nano-
ARPES), energy bands with little dispersion were directly observed in momentum
space at $\theta\approx 0.96^{\circ}$ [27] and $1.34^{\circ}$ [28]. The charge
localization in the AA region at small $\theta$ was also confirmed using
scanning tunneling microscope (STM) [12, 24, 25, 26]. Furthermore, both
transport and STS measurements in the presence of magnetic field have recently
identified correlated gaps that account for both integer [29, 30, 31, 32, 33,
34, 35] and fractional [31, 36] Chern insulating states, as a result of strong
correlation in combination with the flat band topology.

It is unprecedented that a carbon material with only $sp$ electrons can
exhibit such rich phenomena of electron correlation driven by interaction.
Thus, magic-angle twisted bilayer graphene (MATBG) and related systems have
presented special challenges in our understanding of two-dimensional physics.
In order to explain these exotic properties observed, various theoretical
investigations have been reported. Most of these studies [37, 38, 39, 40]
relied on building a minimum model/basis set consisting of, for example,
Wannier orbitals, which reproduces the single-particle flat band dispersion
while satisfying the original symmetry and topology properties as much as
possible. Many-body interactions were then included in the framework of
(generalized) Hubbard model to shed light on the nature of many-body states
observed experimentally. Some other studies [41, 42, 43] pursued alternative
approaches to avoid possible bias toward the basis set. Nevertheless, no
convincing theory that has been proposed to date can fully explain the
interesting and intriguing experimental findings in MATBG yet.

The existence of multiple magic angles with extremely flat bands was found in
many one-particle calculations. Yet their origin and why a series of flat
bands occur at these specific magic angles are still not understood until now.
Building an effective many-body theory for electron correlation and
superconductivity requires, as the starting point, a good knowledge of the
characteristics of one-particle electronic states at the magic angles.
Therefore, in this theoretical work we aim to uncover the unique features in
the electronic properties of TBG. Since the system is too large for full-scale
DFT calculations, we mostly used the tight-binding method with newly developed
parameters for accurate angle-dependent interlayer interactions [44, 45].
First, we provide an explanation using DFT calculations for the reason why
electrons near the Fermi level become localized in the AA region of TBG when
the twist angle $\theta$ gets small. As a result, the local electronic
structure of AA-stacked bilayer graphene (AABLG) with a characteristic Fermi
ring becomes the key feature to be considered for generating a special effect
at certain twist angles. Second, we demonstrate that when the size of this
Fermi ring matches the reciprocal lattice vectors of the moiré superlattice at
a set of angles, the Dirac points of TBG fall on the Fermi ring. Multiple
states on the ring can then be coupled coherently by the moiré potential,
leading to energy bands with extremely small dispersions. This matching
condition in momentum space generates a series of discrete magic angles in
TBG. Computational results on TBG with different interlayer coupling strengths
when it is under uniaxial pressure are presented to illustrate our theory.
This work provides the physical explanation on the origin of magic angles. The
result also indicates that magic angles are special features in twisted
graphene systems arising from their particular electronic structure.

##  II Results

### Electron accumulation in the AA region

Figure 1: $\mid$ Moiré effect on electronic density of states. a, Moiré
pattern of TBG at $\theta=3.15^{\circ}$ with six representative local atomic
configurations $\mathrm{T1}$ \- $\mathrm{T6}$ in addition to AA and AB(BA).
The black (gray) atoms are in the first (second) layer. The arrow in the
insets shows how the second layer of AA-stacked graphene is shifted to
generate a particular local stacking configuration. The thick hexagon marks
the Wigner-Seitz cell of the moiré supercell. b, Density of states (DOS) of
representative bilayer stacking configurations in a, obtained from density-
functional-theory (DFT) calculations. The uppermost panel covers a larger
energy range, while the lower two panels compare the DOS variations near the
Fermi level along two paths
($\textrm{AA}\rightarrow\mathrm{T1}\rightarrow\mathrm{T2}\rightarrow\mathrm{T3}\rightarrow\textrm{AB}$
and
$\textrm{AA}\rightarrow\mathrm{T4}\rightarrow\mathrm{T5}\rightarrow\mathrm{T6}$).
It is noted that the AA region has significantly more low-energy states around
the Fermi level. c, Dirac cone in momentum space for monolayer graphene. d,
Two Dirac cones shifted in energy for AA-stacked bilayer graphene with the
Fermi ring marked by a dashed circle.

In our tight-binding calculation, we have considered commensurate TBG
configurations with moiré lattice vectors of
$\mathbf{L_{1}}=m\,\mathbf{a_{1}}+(m+1)\,\mathbf{a_{2}}$ and
$\mathbf{L_{2}}=-(m+1)\,\mathbf{a_{1}}+(2m+1)\,\mathbf{a_{2}}$, where
$\mathbf{a_{1}}$ and $\mathbf{a_{2}}$ are lattice vectors for monolayer
graphene, and $m$ is an integer [2]. The rotation starts from AABLG with the
axis passing through a vertical pair of carbon atoms. We neglect the possible
interlayer corrugation and intralayer relaxation in this work. Figure 1a shows
the atomic arrangement of a typical moiré pattern in TBG with a rotation angle
of $\theta=3.15^{\circ}$ ($m=10$). Because of the relative rotation between
two graphene layers, the local atomic arrangements exhibit different stacking
patterns at different locations. Some noticeable ones, such as the AA-stacked
and Bernal-stacked (AB or BA) regions, are labeled in Fig. 1a, as well as six
other representative stacking patterns (T1-T6) in between. Thus, the Bloch
electrons will experience different lattice potentials in different stacking
regions within a single moiré supercell, which is an important feature for
this system. To understand the stacking effect, we have calculated the
electronic DOS for these different lattice potentials with DFT, and the
results are shown in Fig. 1b. It turns out that the AA stacking gives
significantly more low-energy states near the Fermi level, while the AB(BA)
stacking has the fewest states in the same energy range. This can be
understood by examining the band dispersion. Compared with monolayer graphene
with one Dirac cone (Fig. 1c), AABLG has two cones shifted vertically in
energy (Fig. 1d), giving rise to a sizable constant DOS between the two
original band crossings. In contrast, the Bernal stacking gives roughly two
low-energy parabolic bands joined at the Dirac point with a smaller DOS. A
tight-binding model with only the nearest-neighbor interlayer interaction will
give a ratio of 4:1 for the low-energy DOS in these two systems. As the size
of the moiré superlattice increases with decreasing $\theta$, so is the size
of a region with a particular stacking pattern. Therefore, this explains the
observation that, at small $\theta$ values, the electrons near the Fermi level
tend to accumulate around the AA region in the moiré cell, because more states
are available to them there.

Figure 2: $\mid$ Fermi velocity and localized electrons. a-e, Band structure
(upper panel) of twisted bilayer graphene with angles of $9.43^{\circ}$,
$3.89^{\circ}$, $2.13^{\circ}$, $1.12^{\circ}$ and $0.86^{\circ}$,
respectively, calculated by the tight-binding method. The four bands near the
Fermi level (energy zero) are marked in red. Note that the bands become
flatter when the angle decreases, with a minimal Fermi velocity near
$1.12^{\circ}$, but the slope at $0.86^{\circ}$ bounces back after passing
through this magic angle. The lower panel shows the electronic probability
distribution of the zero-energy states at $K$. The dark and light dots mark
the atoms in layer 1 and 2, respectively, and the dashed dark hexagon
represents the Wigner-Seitz cell with the real-space dimension labels in units
of Å. The size of the red dots represents the probability of finding the
electrons at each atom. f, Calculated Fermi velocity (red curve) and electron
localization (blue curve) as a function of the twist angle $\theta$. The Fermi
velocity is normalized with respect to that in monolayer graphene $v_{F}^{0}$.
The electron localization is measured by the total fraction of $K$ electrons
within the the blue dashed circle centered at the AA region that has a radius
of 20% of the moiré lattice constant. The Fermi velocity has a sharp dip at
$1.12^{\circ}$ near the magic angle. In contrast, the electron localization
measure increases smoothly and monotonically as the twist angle decreases.

To illustrate the major change in the electronic structure as $\theta$ varies,
we present in the upper panels of Fig. 2a-e the band structure of TBG by
tight-binding calculations for five different twist angles (9.43∘, 3.89∘,
2.13∘, 1.12∘ and 0.86∘). Note the change of the energy scale in these plots.
The size of the moiré supercell is 1.5 nm, 3.6 nm, 6.6 nm, 13 nm, and 16 nm,
respectively. The four bands in red near the Fermi energy gradually flatten as
the twist angle decreases in Fig. 2a-c, which is expected for an increased
coupling between the Dirac cones of two layers. Then there is a sudden change
near the first magic angle of 1.1∘ where the bands become extremely flat (Fig.
2d). After passing through the first magic angle, the bands regain their
slopes (Fig. 2e). This trend can be seen clearly from Fig. 2f in which the red
curve shows the normalized Fermi velocity of TBG (with respect to the Fermi
velocity $v^{0}_{F}$ of monolayer graphene) as a function of angle $\theta$
and exhibits a clear sharp dip at the first magic angle. The fact that the
extremely flat bands only occur at certain discrete magic angles indicates
that there exists a characteristic condition that breaks the expected
monotonic behavior.

It is interesting to examine the developed electron accumulation in the AA
region. We present the probability of finding electrons at each atomic site in
the lower panels of Fig. 2a-e. The dark and light dots mark the atoms in the
first and second layer, respectively, and the dashed dark hexagon marks the
Wigner-Seitz cell of the moiré superlattice. The electron distribution of the
zero-energy states at ${K}$ is found to be similar to that integrated over the
four red bands in the upper panels. Therefore, we show only the distribution
of the four degenerate states at ${K}$, represented by the size of red dots in
the lower panels of Fig. 2a-e. As the twist angle is reduced and the AA region
clearly develops, the electrons gradually become accumulated in the AA region
at the center of the Wigner-Seitz cell. We measure the degree of localization
by the integrated fraction of electrons within the dashed blue circle of a
radius 0.2$L$, where $L$ is the lattice constant of the moiré supercell at
each angle. This localization measure is found to increase smoothly as the
twist angle decreases, as shown by the blue curve in Fig. 2f. As discussed
above, this behavior can be explained by the fact that the AA region provides
significantly more low-energy states than other stacking patterns.
Surprisingly, this smooth behavior is not affected by the presence of the
magic angle at all, in contrast to the variation of the Fermi velocity in Fig.
2f. It was proposed previously [2, 8] that the series of magic angles are
associated with quantization conditions of the confined states in the AA
region. However, results in Fig. 2f indicate that the flat band dispersion at
the magic angle may not reflect the most localized electronic distribution and
that the emergence of the magic angle requires a mechanism other than electron
confinement alone.

### Connecting AABLG Fermi ring and magic angles in TBG

As the low-energy electrons are accumulated in the AA region at small angles,
apart from the moiré potential they mostly experience the local lattice
potential of AABLG, which has a ring-shaped Fermi surface. This Fermi ring
(Fig. 1d) has a particular radius ($Q_{F}$) in momentum space that is
determined by the interlayer interaction. Next we propose a physical
explanation for the development of flat bands at the magic angles that is
connected to the Fermi ring of AABLG. In particular, we will show that a
special situation exists if the TBG Dirac points (in the extended Brillouin
zones) fall on the Fermi ring of AABLG. We will then demonstrate that the
magic angles occur under these matching conditions and discuss the
characteristics of the states associated with the resulting flat bands.

Table 1: First magic angle $\theta$ and corresponding moiré lattice constant $L$ of TBG obtained from tight-binding calculations for different vertical compression levels . These are determined by a commensurate supercell of $m$ that has the smallest Fermi velocity. compression | 0% | 3% | 5% | 10% | 15% | 20%
---|---|---|---|---|---|---
$m$ | 29 | 26 | 23 | 17 | 13 | 10
$\theta$ (∘) | 1.12 | 1.25 | 1.41 | 1.89 | 2.45 | 3.15
$L$ (nm) | 12.6 | 11.3 | 10.0 | 7.46 | 5.75 | 4.45
Figure 3: $\mid$ Connection between magic-angle moiré periodicity and AABLG
Fermi ring. a-d, Brillouin zone (dashed green hexagon) and Fermi ring (red
circle) of AABLG superimposed on the Brillouin zones (small black hexagons) of
TBG at different angles. For a clear illustration, the Fermi ring for AABLG at
20% compression of the interlayer separation is shown. Both the dashed green
hexagon and the red ring are of the same size in a-d, but the Brillouin zones
of TBG shrink as the moiré size increases with a decreasing twist angle. The
matching condition (see text) is reached at $\theta=3.15^{\circ}$, which is
the magic angle for TBG at 20% compression. e, Correlation between the moiré
size $L$ at the first magic angle ($S=1$) and the inverse of the radius
$Q_{F}$ of the Fermi ring in AABLG with different interlayer compression. The
linear behavior follows the relationship $L=4/\sqrt{3}\times(\pi/Q_{F})$
(dashed line). The horizontal error bars on the data points mark the variation
in the size of the Fermi ring due to trigonal warping, while the vertical
error bars represent the uncertainty in determining the magic angle of TBG
since our calculations used a particular collection of commensurate moiré
supercells. f, Similar linear connection is found for the second
($S=\sqrt{3}$), third ($S=2$), and the fourth ($S=\sqrt{7}$) magic angles for
the TBG system with 10% compression.

The interaction between the layers can be changed by applying uniaxial
pressure, namely, by imposing a compressive vertical stress. With different
reduced interlayer distances, both the magic angle and the size of the Fermi
ring are separately modified. This provides a collection of sample systems
that allow us to examine our proposed mechanism systematically by tight-
binding calculations. The band dispersion of TBG with different compressed
vertical strain values (0%, 3%, 5%, 10%, 15%, and 20%) was evaluated as a
function of the twist angle in order to determine their corresponding magic
angles. The complete results are listed in Table 1. When the interlayer
separation is compressed, the magic angle occurs at a larger angle [45]. For
example, the first magic angle shifts from $\theta$ = 1.12∘ at 0% compression
to 3.15∘ at 20% compression. At the same time, the radius of the Fermi ring in
AABLG also enlarges as the bilayer is compressed. For large compression
values, a trigonal warping may be found, and an approximated Fermi ring is
determined by a circle with the same area. This deviation will be considered
in later discussions.

Since the compressed bilayer graphene has a smaller moiré supercell at the
first magic angle, it is easier to graphically visualize this effect in a
compressed bilayer system. We plot in Fig. 3a-d the 2D Brillouin zone (dashed
green hexagon) and Fermi ring (red circle) of AABLG with 20% vertical
compression. Superimposed in the figures are the Brillouin zones (small black
hexagons) of TBG at different angles, where the $K$ and $K^{\prime}$ Dirac
points are marked by solid and open circles, respectively. Figure 3a shows the
configuration for $\theta$ = 7.34∘, with an enlarged view in Fig. 3b. At this
angle, the Dirac points of TBG do not seem to be influenced by the AABLG Fermi
ring. Note that the Brillouin zone of TBG shrinks as the moiré size increases
with a decreasing twist angle. A special matching condition is found in Fig.
3d, where the Fermi ring of AABLG centered at $K$ reaches the first star of
Dirac $K$ points in the extended Brillouin zones of TBG. We will show below
that this particular matching condition yields the first magic angle in
systems we have investigated.

Figure 4: $\mid$ Higher-order magic angles. a-c, Matching conditions between
the AABLG Fermi ring (red circle) and the first, second, and third stars of
the $K$ points (black filled circles) of TBG in the extended zones (black
hexagons). These give rise to the first, second, and third magic angles,
respectively. The magnitude of $\bm{G}$, the smallest reciprocal lattice
vector of the moiré supercell decreases in a-c. d, Calculated (normalized)
Fermi velocity as a function of the twist angle for TBG with 10% vertical
compression, showing four dips corresponding to the first four magic angles.
e, Enlarged matching condition in a showing the six $K$ points on the Fermi
ring and possible scattering between them via reciprocal lattice vectors of
the moiré supercell.

We plot in Fig. 3e the moiré supercell size $L$ at the first magic angle found
in the calculation (determined by the dip in the calculated Fermi velocity)
versus the inverse of $Q_{F}$ (the radius of the AABLG Fermi ring) for bilayer
systems with various compression levels. These two quantities were
independently obtained for two differernt systems by our tight-binding
calculations. The horizontal error bars on the data points mark the variation
in the size of the Fermi ring due to the trigonal warping: the smallest
(largest) radius in the ${K\mathit{\Gamma}}$ (${KM}$) direction gives the
upper (lower) bound of $\pi/Q_{F}$. The vertical error bars represent the
uncertainty in determining the exact magic angle, because our calculations
used a particular collection of commensurate moiré supercells. If the largest
Fermi velocity dip was found in the calculation for a moiré supercell with a
particular $m$ value, then the moiré size for the exact magic angle is bounded
by those corresponding to $m-1$ and $m+1$ moiré supercells. As one can see
from Fig. 3e, an almost perfect linear behavior is found with a slope of
$4/\sqrt{3}$. In other words, the size of the AABLG Fermi ring determines the
moiré periodicity where the magic angle occurs. And the first magic angle
appears when the radius of the Fermi ring matches $G$, the smallest magnitude
of nonzero reciprocal lattice vectors of the moiré supercell, as shown in Fig.
3d and Fig. 4a. As will be shown below, the critical matching condition is
$Q_{F}=SG$, where $S=1$ corresponds to the first magic angle.

When the radius of the AABLG Fermi ring matches the second and third stars of
the Dirac $K$ points, as shown in Fig. 4b and c, we expect to see the second
($S=\sqrt{3}$) and third ($S=2$) magic angles, respectively. We show in Fig.
4d a more complete plot of the calculated Fermi velocity versus $\theta$ for
TBG with a 10% compression, and data up to the fourth magic angle (0.71∘,
$S=\sqrt{7}$) was obtained. The corresponding moiré size at the first four
magic angles versus $\pi S/Q_{F}$ is plotted in Fig. 3f. Again, a consistent
linear behavior is found that matches the relationship

$L=\frac{4}{\sqrt{3}}\frac{\pi}{Q_{F}}S,$ (1)

where $S=1,\sqrt{3},2,\sqrt{7}$. Given that the minimal length of the
reciprocal lattice vector is $G=4\pi/(\sqrt{3}L)$, $Q_{F}=S\,G$ is then the
criterion in momentum space that determines the occurrence of the series of
magic angles. Possible $S$ values are expected to reflect the radius of the
stars of $K$ points: $S=1,\sqrt{3},2,\sqrt{7},3,2\sqrt{3},...$. Under these
matching conditions, the zero-energy states within AABLG (Fermi ring) and the
zero-energy states within TBG (Dirac points) are aligned.

At small angles, we have $L\approx a/\theta$, where $a$ is the lattice
constant of monolayer graphene. Therefore, the magic angles in TBG can be
expressed as

$\theta_{magic}\approx\frac{Q_{F}}{\sqrt{3}\,\,l_{K}}\frac{1}{S}\,,$ (2)

where $l_{K}=4\pi/(3a)$ is the distance from $\mathit{\Gamma}$ to $K$ for
monolayer graphene, and $S$ is not necessarily an integer. Our DFT value for
$Q_{F}$ is about 0.055 Å-1. Therefore, the predicted magic angles are 1.07∘,
0.62∘, 0.53∘, 0.40∘, $\ldots$ for TBG. If a tight-binding model considers only
the nearest neighbor coupling $t_{\bot}$ for the interlayer interaction, the
AABLG Fermi ring will have a radius of $Q_{F}=t_{\bot}/(\hbar v^{0}_{F})$.

### Wave functions of the Dirac point at magic angles

Away from the magic angle, the energy bands in Fig. 2 can be explained by the
following picture: As the two layers are rotated by an angle in real space,
the Dirac cones of the two layers are also rotated with respect to each other
in momentum space; the interlayer interaction induces coupling between the two
sets of linear bands, giving rise to a reduced slope (Fermi velocity) that is
supposed to vary monotonically with the angle. This simple picture is
understandable, since each moiré supercell has regions of different stacking
patterns (see Fig. 1a) that cannot influence the energy bands in an abrupt
way. On the other hand, when the angle gets small, the relevant low-energy
electrons are accumulated in the AA region, so they should be highly
influenced by the lattice potential of AABLG. When the matching condition
$Q_{F}=S\,G$ is satisfied at the magic angles, it would be reasonable to
construct the wave functions for the Dirac point of TBG from the Bloch states
on the AABLG Fermi ring.

We use the first magic angle ($Q_{F}=G$) as an example. Figure 4e shows the
detailed geometry with six Dirac points on the AABLG Fermi ring; they could be
folded in TBG to the $K$ point at the ring center (marked by a red point). The
Bloch states on the AABLG Fermi ring are doubly degenerate and half occupied
in a neutral system. Using the nearest-neighbor tight-binding model with an
orthogonal basis of $p_{z}$ orbitals on the $A$ and $B$ sublattices in layer 1
and 2: {$A_{1}$, $A_{2}$, $B_{1}$, $B_{2}$}, these two degenerate states can
be chosen as $\left|{k}\right\rangle$ $\propto$ $(1\;0\;0\;{-e^{i\phi}})^{T}$
and $(0\;1\;{-e^{i\phi}}\;0)^{T}$ with $\phi$ =
${\tan^{-1}}\left({{k_{y}}/{k_{x}}}\right)$. These two degenerate states will
not be coupled with each other by a local moiré potential, so in the following
analysis we make this assumption and deal with one set of them at a time. Each
pair of the six $k$ points in Fig. 4e are connected by a reciprocal lattice
vector of the moiré supercell, so they will be coupled with each other via the
scattering by a 2D moiré potential
$U(\bm{r})=\Sigma_{\bm{G^{\prime}}}{U_{\bm{G^{\prime}}}}e^{i{\bm{G^{\prime}}}\cdot{\bm{r}}}$,
where $\bm{G^{\prime}}$ is a reciprocal lattice vector of the moiré
superlattice. With the combination of ${C_{3}}$ (vertical axis) and ${C_{2}}$
(horizontal axes) rotation symmetries and the fact that the moiré potential is
real ($U_{\bm{G^{\prime}}}$ = $U^{*}_{-\bm{G^{\prime}}}$), the number of
independent Fourier coefficients is greatly reduced. We denote that
$U_{\bm{G^{\prime}}=0}$ = $U_{0}$, $U_{\bm{k_{1}}-\bm{k_{2}}}$ =
$U_{\bm{k_{2}}-\bm{k_{3}}}$ = $\ldots$ = $U_{1}$, and
$U_{\bm{k_{1}}-\bm{k_{3}}}$ = $U^{*}_{\bm{k_{2}}-\bm{k_{4}}}$ =
$U_{\bm{k_{3}}-\bm{k_{5}}}$ = $\ldots$ = $U_{2}$. $U_{1}$ is real, and $U_{2}$
can be made real by choosing the 2D origin at an inversion center.

After including the symmetry of the system and assuming that the long-scale
moiré potential varies little over a single graphene unit cell (see Methods),
the perturbation Hamiltonian of the moiré potential for one set of the wave
functions at the six $k$ points ($\bm{k_{1}}$, $\bm{k_{2}}$, $\ldots$,
$\bm{k_{6}}$ in Fig. 4e) has this form:

$H=\left({\begin{array}[]{*{20}{c}}{u_{11}}&{{u_{12}}}&{{u_{13}}}&0&{u_{13}^{*}}&{u_{12}^{*}}\\\
{}\hfil&{u_{11}}&{{u_{12}}}&{{e^{i2\pi/3}}{u^{*}_{13}}}&0&{e^{-i2\pi/3}}{u_{13}}\\\
{}\hfil&{}\hfil&{u_{11}}&{{u_{12}}}&{{u_{13}}}&0\\\
{}\hfil&{}\hfil&{}\hfil&{u_{11}}&{{u_{12}}}&{{e^{i2\pi/3}}{u^{*}_{13}}}\\\
{}\hfil&{\mathrm{c.}\,\mathrm{c.}}&{}\hfil&{}\hfil&{u_{11}}&{{u_{12}}}\\\
{}\hfil&{}\hfil&{}\hfil&{}\hfil&{}\hfil&{u_{11}}\end{array}}\right),$ (3)

where ${u_{ij}}$ = $\left\langle{k_{i}}\right|U\left|{k_{j}}\right\rangle$
with the following values: ${u_{11}}$ = $U_{0}$, ${u_{12}}$ =
$(\sqrt{3}/2)\,{U_{1}}\,{e^{i\pi/6}}$, ${u_{13}}$ =
$(1/2)\,U_{2}\,{e^{i\pi/3}}$, and ${u_{14}}=0$. This perturbation Hamiltonian
turns out to have three doubly degenerate eigenstates. After a constant shift
of $U_{0}-\operatorname{Re}[U_{2}]$ to align the energy of the middle state
with zero, we have the three eigenvalues:

${E_{0}}=0\,;\;\;{E_{\pm}}=\frac{3}{2}\left[\operatorname{Re}[U_{2}]\pm\sqrt{U_{1}^{2}+\left({\operatorname{Im}[U_{2}]}/\sqrt{3}\right)^{2}}\,\right]\,.$
(4)

In the space of the six $k$ states on the Fermi ring of AABLG, the two
degenerate eigenvectors $v_{1}$ and $v_{2}$ for $E_{0}$ are:

${v_{1}}=\left({\begin{array}[]{*{20}{c}}0\\\ 1\\\ 0\\\ {{e^{i2\pi/3}}}\\\
0\\\
{{e^{i4\pi/3}}}\end{array}}\right)\,{\mathrm{and}}\;\;{v_{2}}=\left({\begin{array}[]{*{20}{c}}1\\\
0\\\ {{e^{i2\pi/3}}}\\\ 0\\\ {{e^{i4\pi/3}}}\\\ 0\end{array}}\right)\;,$ (5)

where $v_{1}$ is a linear combination of $k_{1}$, $k_{3}$, and $k_{5}$ states,
while $v_{2}$ is made of $k_{2}$, $k_{4}$, and $k_{6}$ states. Both $v_{1}$
and $v_{2}$ are invariant under the $C_{3}$ rotation and with the same
eigenvalue $\omega$. (The value of $\omega$ could be $e^{-i2\pi/3}$ or 1 with
the rotation axis chosen at a C atom or the center of a six-atom ring.) The
results for the other set of degenerate Bloch states at the six $k$ points are
the same. Therefore, we obtain four zero-energy degenerate states at the $K$
point of TBG that all transform in the same way under the $C_{3}$ rotation.
For higher magic angles, the matrix elements in Eq. (3) will involve higher
Fourier coefficients of the moiré potential. Nevertheless, the features of the
final eigenvectors will be similar to those for $S$ = 1.

Next we show that the bands turn out to be flat going away from the $K$ point
by the $k\cdot p$ perturbation theory. The main $O(k)$ correction in energy
away from the $K$ point comes from the $k\cdot p$ intraband coupling among the
fourfold degenerate states at $E$=0. These intraband couplings, however, turn
out to vanish because all the four eigenstates ($v_{l}$) at $E$=0 have the
same $C_{3}$ rotation symmetry with the same eigenvalue of $\omega$. This can
be seen as follows. First, we rewrite
$k\cdot\hat{p}={{k}_{+}}{{{\hat{p}}}_{-}}+{{k}_{-}}{{{\hat{p}}}_{+}},$ where
${{k}_{\pm}}=({{k}_{x}}\pm i{{k}_{y}})/\sqrt{2}$ and
${{{\hat{p}}}_{\pm}}=({{{\hat{p}}}_{x}}\pm i{{{\hat{p}}}_{y}})/\sqrt{2}$. The
operator ${\hat{p}}_{\pm}$ transforms as
$C_{n}{\hat{p}}_{\pm}C_{n}^{{\dagger}}$ = $e^{\mp i2\pi/n}{\hat{p}}_{\pm}$
under the $C_{n}$ rotation. Then
$\left\langle{{v}_{l}}\right|{{\hat{p}}_{+}}\left|{{v}_{m}}\right\rangle$ =
$\left\langle{{v}_{l}}\right|C_{3}^{{\dagger}}{{C}_{3}}\,{{\hat{p}}_{+}}C_{3}^{{\dagger}}{{C}_{3}}\left|{{v}_{m}}\right\rangle$
=
$\left\langle{{v}_{l}}\right|{\omega^{*}}(e^{-i2\pi/3}{{\hat{p}}_{+}})\omega\left|{{v}_{m}}\right\rangle$
=
$e^{-i2\pi/3}\left\langle{{v}_{l}}\right|{{{\hat{p}}}_{+}}\left|{{v}_{m}}\right\rangle$.
Therefore,
$\left\langle{{v}_{l}}\right|{{\hat{p}}_{+}}\left|{{v}_{m}}\right\rangle$ = 0,
and similarly
$\left\langle{{v}_{l}}\right|{{\hat{p}}_{-}}\left|{{v}_{m}}\right\rangle$ = 0.
As such, the resonant scattering by the moiré potential regulates the $C_{3}$
rotation symmetry of the degenerate states at $E$=0, kills the intraband
$k\cdot p$ coupling, and thus significantly flattens the bands away from $K$.

In addition, the net $O({k^{2}})$ correction from the interband coupling also
tends to be small. If we neglect the long-range interaction and set $U_{2}$=0,
we have $E_{+}$ = ${-E_{-}}$ and $\sigma H\sigma={-H}$ with $\sigma$ =
Diag[1,-1,1,-1,1,-1]. Consequently, contributions from the interband coupling
cancel each other: the ${E_{+}}$ bands tend to bend downward the ${E_{0}}$
band, while the ${E_{-}}$ bands tend to bend it upward, rendering a net zero
correction in energy to the order of $O({k^{2}})$.

##  III Discussion

In this work we have identified the following key factors for the occurrence
of magic angles in TBG, which may suggest the directions for searching for
this interesting phenomenon in other 2D systems.

First, there exists a particular stacking pattern within the the moiré
supercell that attracts low-energy electrons to accumulate. This charge
concentration becomes significant at small twist angles, because the region
associated with the particular stacking pattern becomes substantial over many
bond lengths. This happens in the AA region of TBG that has more low-energy
states available compared with regions with other stacking patterns. If this
condition does not occur and the electrons spread out throughout the moiré
supercell, the averaging effect will only decrease the bandwidth monotonically
as the supercell size increases.

Second, this particular region has its characteristic band structure near the
Fermi level that will respond to the perturbing moiré potential in a
particular way only at certain twist angles, where the resulting electronic
states can have special features. In the case of TBG, the AA region has a
characteristic Fermi ring, and the states on the ring can be coherently
coupled when the size of the ring matches the reciprocal lattice vectors of
the moiré supercell ($Q_{F}=SG$). Therefore, the size of this Fermi ring
determines the moiré periodicity where the series of magic angles occur, and
the corresponding TBG wave functions at the Dirac $K$ point are special linear
combinations of wave functions on the AABLG Fermi ring.

Third, the characteristic wave functions at the magic angles have certain
symmetry properties that are able to induce a flat band dispersion away from
the reference $k$ point. This is possible because with a proper symmetry
transformation property for $\hat{p}_{\pm}$, the first-order correction within
the $k\cdot p$ perturbation theory can vanish. Possible symmetries include
inversion, $C_{n}$ rotations, and reflections. At the magic angles in TBG, the
fourfold degenerate $E$=0 states at $K$ have the same eigenvalue under the
$C_{3}$ rotation. The transformation property
$C_{3}{\hat{p}}_{\pm}C_{3}^{{\dagger}}$ = $e^{\mp i2\pi/3}{\hat{p}}_{\pm}$
sets the intraband $k\cdot p$ coupling to zero, giving rise to flat bands out
of $K$ to the linear order.

The calculations presented in this work did not include the lattice relaxation
effect in TBG. Previous theoretical studies have shown that the relaxation
pattern in the AA region can simply be described as adding an extra small
local rotation around the center of the AA region [46, 47]. For example, for
the case of $\theta=1.05^{\circ}$, the local rotation angle between the two
layers is increased from $1.05^{\circ}$ to $1.63^{\circ}$ upon relaxation
[46], and for very small $\theta$ values, the relaxed local twist angle at the
AA stacking converges to $1.7^{\circ}$ [47]. As a consequence, the AA stacking
remains intact after relaxation, albeit with a reduced extent. Nevertheless,
the remaining size is still of the order of a few nanometers for small twist
angles (see Fig. 5 in [9] for $\theta=0.99^{\circ}$ and Fig. 7 in [46] for
$\theta=1.05^{\circ}$). This is consistent with the STM topographic images at
around $\theta=1.1^{\circ}$ that also showed bright spots of AA stacking
regions with a size of a few nanometers [26, 48]. Since there still exist
sizable AA regions after one takes into account the lattice relaxation, our
theory based on the Fermi ring of the AA stacking should remain valid.

A previous study focusing on an idealized continuum model by switching off the
interlayer coupling parameter $w_{AA}$ for the AA stacking yielded perfectly
flat bands at multiple magic angles [49]. Intriguingly, the wave functions of
these flat bands at the magic angles in this particular model were found to be
reminiscent of quantum Hall wave functions on the torus. Because the condition
$w_{AA}$ = 0 gives rise to perfectly flat bands, it was further speculated
that the fundamental features of TBG are mainly connected with the interlayer
coupling $w_{AB}$ in the AB stacking. Since this idealized model does not
correspond to any realistic graphene systems, it was hoped that one can extend
the results to cases with nonzero $w_{AA}$ by treating $w_{AA}$/$w_{AB}$ as a
perturbation. However, this may be difficult because $w_{AA}$ and $w_{AB}$
should be of the same order of magnitude in realistic TBG systems. In
contrast, the current work focuses on the important effect of the $w_{AA}$
interaction based on the fact that low-energy electrons in TBG are
concentrated in the AA region. Our proposed theory is therefore able to fully
explain the occurrence and positions of magic angles in terms of the local
electronic properties of the AA region.

In conclusion, this work resolves the long-standing question why magic angles
happen in twisted bilayer graphene. The origin of the magic angles turns out
to be the “magic” Fermi ring in AA-stacked bilayer graphene. We have
theoretically derived the series of magic angle values. In addition, we have
identified the composition of the one-particle wave functions associated with
the flat bands for further many-body investigations.

## Methods

### Tight-binding calculations

To calculate the electronic properties of TBG, we used an extended tight-
binding model proposed by S. Fang and E. Kaxiras [44] with parameters
determined from calculations with density functional theory (DFT). This model
includes intralayer couplings up to eighth nearest neighbors, and a strong
angular dependence of interlayer couplings between atoms in two different
layers. The spin-orbit interaction is neglected. The interlayer couplings
under uniaxial pressure were also obtained using similar procedures [45]. The
Hamiltonian was then diagonalized by the FEAST eigenvalue solver [50], and the
eigenvalues and eigenvectors were verified by Mathematica in some cases. The
density of states were calculated by the tetrahedral method [51], and the
slope at the $K$ point (Fermi velocity) was evaluated based on the Hellmann-
Feynman theorem [52] as follows. The velocity operator is given by
$\vec{v}(\vec{k})={\partial_{\vec{k}}}H(\vec{k})$, where $\vec{k}$ is the
crystal momentum, and $H(\vec{k})$ is the tight-binding Hamiltonian of the
moiré structure. Due to the fourfold degeneracy at the Dirac point, we formed
a velocity matrix ${v_{ij}}=\left\langle i\mid{v}(K)\mid j\right\rangle$,
where $i$ and $j$ are labels for the degenerate states at the Dirac point. The
eigenvalues of the velocity matrix give the velocities; the one with the
smallest absolute value was taken as the velocity result. Since the velocity
is almost isotropic, we show the result for the ${K\rightarrow M}$ direction
in most cases. When there was uncertainty, the average between ${K\rightarrow
M}$ and ${K\rightarrow\Gamma}$ was taken.

### Density functional calculations

We also performed DFT calculations for bilayer graphene with different
stacking patterns using the Vienna Ab initio Simulation Package (VASP) [53,
54]. We used the projector augmented wave (PAW) method [55, 56] to treat core
electrons and the optB86b functional to include the van der Waals correction
[57, 58], with a plane-wave cutoff energy of 400 eV. The Monkhorst-Pack
k-point mesh of 21$\times$21$\times$1 was used in the self-consistent
calculations for the 1$\times$1 unit cell. The calculation of the density of
states used a more denser 91$\times$91$\times$1 grid.

### Matrix elements of the moiré potential

The Bloch-state components constructed from the local orbitals at any of the
four sites $\left\\{{{A_{1}},{A_{2}},{B_{1}},{B_{2}}}\right\\}$ in the unit
cell of AABLG have the form:

$\left\langle\bm{r}\right|\left.{{\alpha_{i}}}\right\rangle=\frac{1}{{\sqrt{{N_{\alpha,i}}}}}\sum\limits_{\bm{R_{\alpha,j}}}{e^{i\bm{k}\cdot\bm{R_{\alpha,j}}}}{\zeta_{{p_{z}}}}(\bm{r}-\bm{R_{\alpha,j}}-{\delta_{i,2}}c\hat{z})\;,$

where $\bm{R_{\alpha,j}}$ is the in-plane position vector for sublattice
$\alpha$ in the $j{\rm{th}}$ unit cell; $c$ is the interlayer separation, $i$
refers to the layer index {1,2}, and ${\delta_{i,2}}$ = 1 if $i$=2.
${\zeta_{{p_{z}}}}(\bm{r})$ describes the orthonormal ${p_{z}}$ atomic
orbital, and ${N_{\alpha,i}}$ = $N$ is the total number of atoms pertaining to
sublattice $\alpha$ and layer $i$ in AABLG. For the six $\bm{k}$ vectors on
the Fermi ring in Fig. 4e, there are two degenerate states, and the wave
functions can be chosen as $\left|{k}\right\rangle$ $\propto$
$(1\;0\;0\;{-e^{i\phi}})^{T}$ and $(0\;1\;{-e^{i\phi}}\;0)^{T}$ with $\phi$ =
${\tan^{-1}}\left({{k_{y}}/{k_{x}}}\right)$. These two states will not be
coupled with each other by the local moiré potential, so we only have to deal
with one of them at a time. Both of them are eigenstates of the $C_{3}$
rotation.

The ${C_{3}}$ rotation symmetry exists for both TBG [and thus the moiré
potential $U(\bm{r})$] and AABLG. It connects the matrix elements ${u_{ij}}$ =
$\left\langle{k_{i}}\right|U\left|{k_{j}}\right\rangle$, for example,
$\left\langle{{k_{1}}}\right|U\left|{{k_{2}}}\right\rangle$ =
$\left\langle{{k_{1}}}\right|C_{3}^{{\dagger}}(C_{3}UC_{3}^{{\dagger}}){C_{3}}\left|{{k_{2}}}\right\rangle$
= $\left\langle{{k_{3}}}\right|U\left|{{k_{4}}}\right\rangle$. Therefore, we
have (as shown in Fig. 4e) ${u_{12}}={u_{34}}={u_{56}}$ (navy dashed lines);
${u_{23}}={u_{45}}={u_{61}}$ (blue lines); ${u_{13}}={u_{35}}={u_{51}}$
(orange dashed line); ${u_{24}}={u_{46}}={u_{62}}$; and
${u_{14}}={u_{36}}={u_{52}}$ (brown dashed line). As a consequence, the
$C_{3}$ symmetry reduces the number of variables from
$\left({{6^{2}}-6}\right)/2=15$ to 5, which are
${u_{12}},{u_{23}},{u_{13}},{u_{24}},{\rm{and}}{u_{14}}$.

Since the moiré supercell is much larger than the unit cell of graphene, we
can assume that the effective moiré potential varies little within one unit
cell of AABLG. Therefore, we have the following approximation:

$\displaystyle\int{\mathrm{d^{3}}{r}}\zeta_{{p_{z}}}^{*}(\bm{r}-\bm{R_{\alpha,j}}-{\delta_{i,2}}c\hat{z})\,U(\bm{r})\,$
$\displaystyle{\zeta_{{p_{z}}}}(\bm{r}-\bm{R_{\alpha^{\prime},j^{\prime}}}-{\delta_{i^{\prime},2}}c\hat{z})$
$\displaystyle\approx
U(\bm{R_{\alpha,j}}){\delta_{i,i^{\prime}}}{\delta_{j,j^{\prime}}}{\delta_{\alpha,\alpha^{\prime}}}\,.$

Therefore, the matrix elements of the 2D moiré potential for six
$\left|{k}\right\rangle$ $\propto$ $(1\;0\;0\;{-e^{i\phi}})^{T}$ are:

${}\begin{split}{u_{lm}}&=\left\langle{{k_{l}}}\right|U\left|{{k_{m}}}\right\rangle\\\
&\approx\frac{1}{{2N}}\left[\sum\limits_{{R_{A,j}}}{{e^{i({\bm{k_{m}}-\bm{k_{l}}})\cdot\bm{R_{A,j}}}}}\,U(\bm{R_{A,j}})\right.\\\
&\left.\hskip
39.83368pt+{e^{i({\phi_{m}}-{\phi_{l}})}}\sum\limits_{{R_{B,j}}}{{e^{i(\bm{k_{m}}-\bm{k_{l}})\cdot\bm{R_{B,j}}}}}\,U(\bm{R_{B,j}})\right]\\\
&\approx\frac{1}{{2A}}\left({1+{e^{i({\phi_{m}}-{\phi_{l}})}}}\right)\int{e^{i(\bm{k_{m}}-\bm{k_{l}})\cdot\bm{r}}U(\bm{r})\mathrm{d^{2}}{r}}\\\
&=\frac{1}{2}\left({1+{e^{i({\phi_{m}}-{\phi_{l}})}}}\right){U_{\bm{k_{l}}-\bm{k_{m}}}}\,,\end{split}$

where $A$ denotes the whole area of the TBG. With the combination of ${C_{3}}$
and ${C_{2}}$ rotation symmetries and the fact that the moiré potential is
real ($U_{\bm{G^{\prime}}}$ = $U^{*}_{-\bm{G^{\prime}}}$), we have
$U_{\bm{k_{1}}-\bm{k_{2}}}$ = $U_{\bm{k_{2}}-\bm{k_{3}}}$ = $\ldots$ =
$U_{1}$, and $U_{\bm{k_{1}}-\bm{k_{3}}}$ = $U^{*}_{\bm{k_{2}}-\bm{k_{4}}}$ =
$U_{\bm{k_{3}}-\bm{k_{5}}}$ = $\ldots$ = $U_{2}$. $U_{1}$ is real, and $U_{2}$
is generally complex but can be made real by choosing the 2D origin at an
inversion center. Given that $\phi_{2}$ \- $\phi_{1}$ = $\pi/3$, $\phi_{3}$ \-
$\phi_{1}$ = $2\pi/3$, $\phi_{4}$ \- $\phi_{1}$ = $\pi$, $(1+e^{i\theta})/2$ =
$\cos(\theta/2)e^{i\theta/2}$, and $u_{lm}$ = $u_{ml}^{*}$, we then have
${u_{11}}$ = $U_{G=0}$, ${u_{12}}$ = $(\sqrt{3}/2)\,U_{1}\,{e^{i\pi/6}}$ =
${u_{23}}$, ${u_{13}}$ = $(1/2)\,U_{2}\,{e^{i\pi/3}}$, ${u_{24}}$ =
$e^{i2\pi/3}u_{13}^{*}$ and ${u_{14}}=0$.

###### Acknowledgements.

We gratefully acknowledge Li-Feng Yen and Po-Tung Fang for sharing the tight-
binding codes and for helping with technical matters. Helpful discussions with
Wen-Ying Ruan, Chi-Ruei Pan, Jie-Cheng Chen, Wei-En Tseng, and Martin Callsen
are acknowledged. This work is supported by Academia Sinica, Taiwan.

## References

  * Shallcross  _et al._ [2010] S. Shallcross, S. Sharma, E. Kandelaki, and O. A. Pankratov, Electronic structure of turbostratic graphene, Phys. Rev. B 81, 165105 (2010).
  * Lopes dos Santos  _et al._ [2012] J. M. B. Lopes dos Santos, N. M. R. Peres, and A. H. Castro Neto, Continuum model of the twisted graphene bilayer, Phys. Rev. B 86, 155449 (2012).
  * Carr  _et al._ [2017] S. Carr, D. Massatt, S. Fang, P. Cazeaux, M. Luskin, and E. Kaxiras, Twistronics: Manipulating the electronic properties of two-dimensional layered structures through their twist angle, Phys. Rev. B 95, 075420 (2017).
  * Lopes dos Santos  _et al._ [2007] J. M. B. Lopes dos Santos, N. M. R. Peres, and A. H. Castro Neto, Graphene bilayer with a twist: Electronic structure, Phys. Rev. Lett. 99, 256802 (2007).
  * 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).
  * Suárez Morell  _et al._ [2010] E. Suárez Morell, J. D. Correa, P. Vargas, M. Pacheco, and Z. Barticevic, Flat bands in slightly twisted bilayer graphene: Tight-binding calculations, Phys. Rev. B 82, 121407 (2010).
  * Trambly de Laissardière  _et al._ [2010] G. Trambly de Laissardière, D. Mayou, and L. Magaud, Localization of Dirac Electrons in Rotated Graphene Bilayers, Nano Lett. 10, 804 (2010).
  * Trambly de Laissardière  _et al._ [2012] G. Trambly de Laissardière, D. Mayou, and L. Magaud, Numerical studies of confined states in rotated bilayers of graphene, Phys. Rev. B 86, 125413 (2012).
  * Uchida  _et al._ [2014] K. Uchida, S. Furuya, J.-I. Iwata, and A. Oshiyama, Atomic corrugation and electron localization due to moiré patterns in twisted bilayer graphenes, Phys. Rev. B 90, 155451 (2014).
  * Cao  _et al._ [2016] Y. Cao, J. Y. Luo, V. Fatemi, S. Fang, J. D. Sanchez-Yamagishi, K. Watanabe, T. Taniguchi, E. Kaxiras, and P. Jarillo-Herrero, Superlattice-Induced Insulating States and Valley-Protected Orbits in Twisted Bilayer Graphene, Phys. Rev. Lett. 117, 116804 (2016).
  * Kim  _et al._ [2016] K. Kim, M. Yankowitz, B. Fallahazad, S. Kang, H. C. P. Movva, S. Huang, S. Larentis, C. M. Corbet, T. Taniguchi, K. Watanabe, S. K. Banerjee, B. J. LeRoy, and E. Tutuc, Van der Waals Heterostructures with High Accuracy Rotational Alignment, Nano Lett. 16, 1989 (2016).
  * Kim  _et al._ [2017] K. Kim, A. DaSilva, S. Huang, B. Fallahazad, S. Larentis, T. Taniguchi, K. Watanabe, B. J. LeRoy, A. H. MacDonald, and E. Tutuc, Tunable moiré bands and strong correlations in small-twist-angle bilayer graphene, PNAS 114, 3364 (2017).
  * Frisenda  _et al._ [2018] R. Frisenda, E. Navarro-Moratalla, P. Gant, D. P. D. Lara, P. Jarillo-Herrero, R. V. Gorbachev, and A. Castellanos-Gomez, Recent progress in the assembly of nanodevices and van der Waals heterostructures by deterministic placement of 2D materials, Chem. Soc. Rev. 47, 53 (2018).
  * Ribeiro-Palau  _et al._ [2018] R. Ribeiro-Palau, C. Zhang, K. Watanabe, T. Taniguchi, J. Hone, and C. R. Dean, Twistable electronics with dynamically rotatable heterostructures, Science 361, 690 (2018).
  * 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).
  * 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, Nat. Phys. 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).
  * Cao  _et al._ [2021] Y. Cao, D. Rodan-Legrain, J. M. Park, N. F. Q. Yuan, K. Watanabe, T. Taniguchi, R. M. Fernandes, L. Fu, and P. Jarillo-Herrero, Nematicity and competing orders in superconducting magic-angle graphene, Science 372, 264 (2021).
  * Liu  _et al._ [2021] X. Liu, Z. Wang, K. Watanabe, T. Taniguchi, O. Vafek, and J. I. A. Li, Tuning electron correlation in magic-angle twisted bilayer graphene using Coulomb screening, Science 371, 1261 (2021).
  * Sharpe  _et al._ [2019] A. L. Sharpe, E. J. Fox, A. W. Barnard, J. Finney, K. Watanabe, T. Taniguchi, M. A. Kastner, and D. Goldhaber-Gordon, Emergent ferromagnetism near three-quarters filling in twisted bilayer graphene, Science 365, 605 (2019).
  * Serlin  _et al._ [2020] M. Serlin, C. L. Tschirhart, H. Polshyn, Y. Zhang, J. Zhu, K. Watanabe, T. Taniguchi, L. Balents, and A. F. Young, Intrinsic quantized anomalous Hall effect in a moiré heterostructure, Science 367, 900 (2020).
  * Li  _et al._ [2010] G. Li, A. Luican, J. M. B. Lopes dos Santos, A. H. Castro Neto, A. Reina, J. Kong, and E. Y. Andrei, Observation of van hove singularities in twisted graphene layers, Nature Physics 6, 109 (2010).
  * Brihuega  _et al._ [2012] I. Brihuega, P. Mallet, H. González-Herrero, G. Trambly de Laissardière, M. M. Ugeda, L. Magaud, J. M. Gómez-Rodríguez, F. Ynduráin, and J.-Y. Veuillen, Unraveling the intrinsic and robust nature of van hove singularities in twisted bilayer graphene by scanning tunneling microscopy and theoretical analysis, Phys. Rev. Lett. 109, 196802 (2012).
  * Kerelsky  _et al._ [2019] A. Kerelsky, L. J. McGilly, D. M. Kennes, L. Xian, M. Yankowitz, S. Chen, K. Watanabe, T. Taniguchi, J. Hone, C. Dean, A. Rubio, and A. N. Pasupathy, Maximized electron interactions at the magic angle in twisted bilayer graphene, Nature 572, 95 (2019).
  * 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, Nat. Phys. 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, Nat. Phys. 17, 189 (2021).
  * Nuckolls  _et al._ [2020] K. P. Nuckolls, M. Oh, D. Wong, B. Lian, K. Watanabe, T. Taniguchi, B. A. Bernevig, and A. Yazdani, Strongly correlated Chern insulators in magic-angle twisted bilayer graphene, Nature , 1 (2020).
  * Choi  _et al._ [2021] Y. Choi, H. Kim, Y. Peng, A. Thomson, C. Lewandowski, R. Polski, Y. Zhang, H. S. Arora, K. Watanabe, T. Taniguchi, J. Alicea, and S. Nadj-Perge, Correlation-driven topological phases in magic-angle twisted bilayer graphene, Nature 589, 536 (2021).
  * Wu  _et al._ [2021] S. Wu, Z. Zhang, K. Watanabe, T. Taniguchi, and E. Y. Andrei, Chern insulators, van Hove singularities and topological flat bands in magic-angle twisted bilayer graphene, Nat. Mater. 20, 488 (2021).
  * 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, Nat. Phys. 17, 478 (2021).
  * Das  _et al._ [2021] I. Das, X. Lu, J. Herzog-Arbeitman, Z.-D. Song, K. Watanabe, T. Taniguchi, B. A. Bernevig, and D. K. Efetov, Symmetry-broken Chern insulators and Rashba-like Landau-level crossings in magic-angle bilayer graphene, Nat. Phys. 17, 710 (2021).
  * Park  _et al._ [2021] J. M. Park, Y. Cao, K. Watanabe, T. Taniguchi, and P. Jarillo-Herrero, Flavour Hund’s coupling, Chern gaps and charge diffusivity in moiré graphene, Nature 592, 43 (2021).
  * Pierce  _et al._ [2021] A. T. Pierce, Y. Xie, J. M. Park, E. Khalaf, S. H. Lee, Y. Cao, D. E. Parker, P. R. Forrester, S. Chen, K. Watanabe, T. Taniguchi, A. Vishwanath, P. Jarillo-Herrero, and A. Yacoby, Unconventional sequence of correlated Chern insulators in magic-angle twisted bilayer graphene, Nat. Phys. 17, 1210 (2021).
  * Xie  _et al._ [2021] Y. Xie, A. T. Pierce, J. M. Park, D. E. Parker, E. Khalaf, P. Ledwith, Y. Cao, S. H. Lee, S. Chen, P. R. Forrester, K. Watanabe, T. Taniguchi, A. Vishwanath, P. Jarillo-Herrero, and A. Yacoby, Fractional Chern insulators in magic-angle twisted bilayer graphene, Nature 600, 439 (2021).
  * Yuan and Fu [2018] N. F. Q. Yuan and L. Fu, Model for the metal-insulator transition in graphene superlattices and beyond, Phys. Rev. B 98, 045103 (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).
  * 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).
  * Po  _et al._ [2018] H. C. Po, L. Zou, A. Vishwanath, and T. Senthil, Origin of Mott Insulating Behavior and Superconductivity in Twisted Bilayer Graphene, Phys. Rev. X 8, 031089 (2018).
  * Rademaker and Mellado [2018] L. Rademaker and P. Mellado, Charge-transfer insulation in twisted bilayer graphene, Phys. Rev. B 98, 235158 (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).
  * 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).
  * Fang and Kaxiras [2016] S. Fang and E. Kaxiras, Electronic structure theory of weakly interacting bilayers, Phys. Rev. B 93, 235153 (2016).
  * Carr  _et al._ [2018a] S. Carr, S. Fang, P. Jarillo-Herrero, and E. Kaxiras, Pressure dependence of the magic twist angle in graphene superlattices, Phys. Rev. B 98, 085144 (2018a).
  * Nam and Koshino [2017] N. N. T. Nam and M. Koshino, Lattice relaxation and energy band modulation in twisted bilayer graphene, Phys. Rev. B 96, 075311 (2017).
  * Carr  _et al._ [2018b] S. Carr, D. Massatt, S. B. Torrisi, P. Cazeaux, M. Luskin, and E. Kaxiras, Relaxation and domain formation in incommensurate two-dimensional heterostructures, Phys. Rev. B 98, 224102 (2018b).
  * Tilak  _et al._ [2021] N. Tilak, X. Lai, S. Wu, Z. Zhang, M. Xu, R. d. A. Ribeiro, P. C. Canfield, and E. Y. Andrei, Flat band carrier confinement in magic-angle twisted bilayer graphene, Nature Communications 12, 4180 (2021).
  * 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).
  * Polizzi [2009] E. Polizzi, Density-matrix-based algorithm for solving eigenvalue problems, Phys. Rev. B 79, 115112 (2009).
  * Lehmann and Taut [1972] G. Lehmann and M. Taut, On the numerical calculation of the density of states and related properties, physica status solidi (b) 54, 469 (1972).
  * Alon and Cederbaum [2003] O. E. Alon and L. S. Cederbaum, Hellmann-Feynman theorem at degeneracies, Phys. Rev. B 68, 033105 (2003).
  * Kresse and Furthmüller [1996a] G. Kresse and J. Furthmüller, Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set, Phys. Rev. B 54, 11169 (1996a).
  * Kresse and Furthmüller [1996b] G. Kresse and J. Furthmüller, Efficiency of ab-initio total energy calculations for metals and semiconductors using a plane-wave basis set, Computational Materials Science 6, 15 (1996b).
  * Blöchl [1994] P. E. Blöchl, Projector augmented-wave method, Phys. Rev. B 50, 17953 (1994).
  * Kresse and Joubert [1999] G. Kresse and D. Joubert, From ultrasoft pseudopotentials to the projector augmented-wave method, Phys. Rev. B 59, 1758 (1999).
  * Klimeš  _et al._ [2011] J. Klimeš, D. R. Bowler, and A. Michaelides, Van der Waals density functionals applied to solids, Phys. Rev. B 83, 195131 (2011).
  * Román-Pérez and Soler [2009] G. Román-Pérez and J. M. Soler, Efficient Implementation of a van der Waals Density Functional: Application to Double-Wall Carbon Nanotubes, Phys. Rev. Lett. 103, 096102 (2009).