# Theory for constructing effective models for electrons in generic bilayer
graphene

H. Minh Lam  V. Nam Do<EMAIL_ADDRESS>Department of Basic
Science, Phenikaa Institute for Advanced Study (PIAS), A1 Building, Phenikaa
University, Hanoi 10000, Vietnam

###### Abstract

We present and discuss in detail practical techniques in formulating effective
models to describe the dynamics of low-energy electrons in generic bilayer
graphene. Starting from a tight-binding model using the $p_{z}$ orbital of
carbon atoms as a representation basis set, we reformulate it into the problem
of coupling between Bloch states defined in each graphene layer. This approach
allows transferring the original problem into the determination of Bloch
states in two independent material layers and coupling rules of such states.
We show two schemes to parameterize coupled Bloch state vectors. For the
bilayer graphene configurations of small twist angle in which the long
wavelength approximation is applicable, we show that an effective Hamiltonian
can be written in the canonical form of a kinetic term defined by the momentum
operator and a potential term defined by the position operator. The validity
of effective models of different sophistication levels and their potential
application in treating various physical aspects are numerically discussed.

††preprint: Preprint not for distribution CONFIDENTIAL, Version of

##  I Introduction

In recent years, research on engineering two-dimensional (2D) materials has
been developed intensively and extensively.Sierra  _et al._ (2021); Liaoa  _et
al._ (2019); Sangwan and Hersam (2018); Bao and Hoh (2019) The ultimate aim is
to tailor properties of material platforms to match technical requirements.
The 2D materials have been now artificially created by stacking a number of
atom layers together.Rozhkov  _et al._ (2016) The bilayer graphene system is
such a typical 2D material in which the cohesion between two graphene layers
is governed by the van der Waals (vdW) forces. Working with these 2D vdW
materials, the electronic properties can be tuned by a number of new fashions,
for example, by sliding and/or twisting the two graphene layers. By these
ways, it opens a possibility of engineering the electronic properties of
materials called twistronics.Carr  _et al._ (2017) Theoretically, conventional
approaches to study the electronic properties of materials are based on the
calculation using the density functional theory (DFT), empirical tight-binding
(TB) methods, and also effective models. The first two are the atomic-scale
approaches. They are practically built on the notion of translation invariant
unit cell, which allows exploiting the Bloch theorem to parameterize all
electron eigen-states by a wave vector defined in a Brillouin zone. However,
for a given twist angle it is not always to find a translation invariant unit
cell, except for the case of commensurate alignment between the two graphene
lattices. Therefore, it limits the use of the first two approaches. Besides,
even in the presence of a translation invariant unit cell, the number of atoms
contained in a unit cell can be very large, particularly in the case of small
twist angles. Thus, the computational power required for evaluating
fundamental properties of TBGs is really demanding.Uchida  _et al._ (2014);
Lucignano  _et al._ (2019); de Laissardiere  _et al._ (2012); Morell  _et al._
(2010)

Working with multiple layer vdW materials one question that naturally arises
is: whether can we understand the properties of the complex system if already
knowing basic properties of each individual material layer? This is actually
the idea for the approach of effective models to the twisted bilayer graphene
(TBG) system since 2007.dos Santos  _et al._ (2007); Shallcross  _et al._
(2008); Bistritzer and MacDonald (2011); Shallcross  _et al._ (2013); dos
Santos  _et al._ (2012); Zhang (2019) Though the variation in the scale of the
graphene lattice constant is ignored, effective models provide surprising
results in describing the electronic structure of the TBG configurations of
tiny twist angles in low energy range. Effective models of different
sophistication levels were derived.Rozhkov  _et al._ (2016) For example,
models similar to the Dirac Hamiltonian were established on the basis of
symmetry analysis of the TBG lattice.Mele (2010, 2011); de Gail  _et al._
(2011) Techniques for establishing larger dimension Hamiltonian matrix were
also proposed to obtain better quantitative results.Weckbecker  _et al._
(2016); Tritsaris  _et al._ (2020) Recently, it was found efforts to include
effects of lattice relaxation and lattice vibrations.Koshino and Nam (2020)
Constructing effective models is a rather technical work. However, available
instructions are, in general, presented succinctly and not in detail. It thus
may lead to the difficulty in understanding and/or the lack of background for
the late use of such models to investigate various aspects of the electronic
properties of the TBG systems, and especially, to develop for other vdW
materials. Our aim in this paper is to develop a general formulation for the
construction of effective models for electrons in generic bilayer systems.
Methodologically, we show that effective models can be constructed by
collecting an appropriate sub-set of single layer electronic states to
represent the states of electrons in the complex systems. We will especially
emphasize on the detail of practical techniques as deriving appropriate thumb
rules. Using single-layer Bloch states as a general representation basis set
we show two schemes to parameterize such coupled Bloch states that help to
isolate from the total Hamiltonian matrix a block as an effective Hamiltonian.
We also show that in the limit of long wavelength the procedure of
constructing such an effective Hamiltonian block can be significantly
simplified. It results in a useful continuum effective model that consists of
a kinetic and a potential term canonically defined by the momentum and
position operators, respectively. The solutions to derived effective models
are presented in detail. We also discuss the validity of effective models at
several approximation levels to suggest their application in studying
different aspects of the bilayer systems.

The contents of the paper are arranged as follows. In Sec. II we present the
description of the atomic lattice of the bilayer graphene. In Sec. III we
start by presenting in sub-Sec. III.1 a general formulation for the
Hamiltonian based on the tight-binding description for the $p_{z}$ atomic
orbitals as the vector basis set. The translation symmetry of the graphene
lattice is exploited to transform the atomic orbital basis set to the new one
of Bloch states. A selection rule for the coupling between such Bloch states
in the same material layer is stated in sub-Sec. III.2. Another selection for
the coupling between Bloch states in different layers is presented in sub-Sec.
III.3. We apply these selection rules to build a procedure for collecting a
subset of coupled Bloch vectors that allows to establish an effective
Hamiltonian model. Practical techniques are presented and discussed in detail
in sub-Sec. III.4. In sub-Sec. III.5 we show how the developed procedure is
employed to derive a continuum effective model in the limit of long
wavelength. The case of commensurate TBGs is detailed in sub-Sec. III.6. To
demonstrate the correctness of the calculation procedure we present in Sec. IV
a solution for Bistritzer-MacDonald model. Finally, we present the conclusion
of the paper in Sec. V.

##  II Lattice geometry

Let us consider a system of two flat graphene layers stacked together with the
interlayer distance $d_{GG}=3.35$ Å and relatively twisted and displaced to
each other. The twist is measured by an angle $\theta$ and the displacement is
measured by a vector $\boldsymbol{\tau}$. We assume that the twisted bilayer
graphene system (TBG) is made via three steps: firstly, two graphene sheets
are coincidentally stacked together; secondly, the two graphene layers are
relatively rotated around an axis that goes perpendicularly through the
graphene sheets at the center point of the atomic hexagonal ring; and thirdly,
the two graphene layers are relatively displaced via the vector
$\boldsymbol{\tau}$. Accordingly, when $\theta=0$ and $\boldsymbol{\tau}=0$
the TBG system corresponds to the AA-stacked configuration. By this
construction, the atomic lattice of TBG has the symmetry of the $D_{3}$ point
group if $\boldsymbol{\tau}=0$. However, one should notice the effects of the
displacement and the twist that while the former destroys all point
symmetries, the latter generally breaks the translation symmetry of the
complex lattices. To determine the position of each carbon atom in the TBG
lattice we firstly define two basis vectors of the AA-stacked configuration.
For instance, we can choose:

$\mathbf{a}_{1}=a\left(\frac{\sqrt{3}}{2}\hat{\mathbf{x}}+\frac{1}{2}\hat{\mathbf{y}}\right);\mathbf{a}_{2}=a\left(\frac{\sqrt{3}}{2}\hat{\mathbf{x}}-\frac{1}{2}\hat{\mathbf{y}}\right)$
(1)

where $a=\sqrt{3}a_{CC}$ is the lattice constant and $a_{CC}=1.45$ Å is the
distance between two nearest lattice sites; $\hat{\mathbf{x}}$ and
$\hat{\mathbf{y}}$ denote the unit vectors of the Cartesian coordinate frame.
These two vectors define a unit cell with the area of
$S_{c}=\|\mathbf{a}_{1}\times\mathbf{a}_{2}\|=a^{2}\sin(\pi/3)$. The unit cell
of the AA-stacked configuration contains 4 Carbon atoms whose positions are
given by:

$\displaystyle\mathbf{d}_{A_{1}}$
$\displaystyle=\mathbf{OA}_{1}=-\frac{1}{2}\mathbf{d}_{1},$ (2a)
$\displaystyle\mathbf{d}_{B_{1}}$
$\displaystyle=\mathbf{OB}_{1}=+\frac{1}{2}\mathbf{d}_{1},$ (2b)
$\displaystyle\mathbf{d}_{A_{2}}$
$\displaystyle=\mathbf{OA}_{2}=-\frac{1}{2}\mathbf{d}_{1}+d_{GG}\hat{\mathbf{z}},$
(2c) $\displaystyle\mathbf{d}_{B_{2}}$
$\displaystyle=\mathbf{OB}_{2}=+\frac{1}{2}\mathbf{d}_{1}+d_{GG}\hat{\mathbf{z}},$
(2d)

where $\mathbf{d}_{1}=(\mathbf{a}_{1}+\mathbf{a}_{2})/3$.

We also define two reciprocal lattice vectors
$\mathbf{a}_{1}^{\star},\mathbf{a}_{2}^{\star}$ associating to the vectors
$\mathbf{a}_{1}$ and $\mathbf{a}_{2}$ according to the conditions
$\mathbf{a}_{i}^{\star}\mathbf{a}_{j}=2\pi\delta_{ij}$. For $\mathbf{a}_{1}$
and $\mathbf{a}_{2}$ given by Eq. (1) we obtain:

$\mathbf{a}_{1}^{\star}=\frac{4\pi}{\sqrt{3}a}\left(\frac{1}{2}\hat{\mathbf{x}}+\frac{\sqrt{3}}{2}\hat{\mathbf{y}}\right);\mathbf{a}_{2}^{\star}=\frac{4\pi}{\sqrt{3}a}\left(\frac{1}{2}\hat{\mathbf{x}}-\frac{\sqrt{3}}{2}\hat{\mathbf{y}}\right).$
(3)

These vectors define a Brillouin zone taking the shape of a hexagon with six
corner points, called the $K$ points, determined by:

$\displaystyle\mathbf{K}_{1}$
$\displaystyle=\frac{2}{3}\mathbf{a}_{1}^{\star}+\frac{1}{3}\mathbf{a}_{2}^{\star}=-\mathbf{K}_{4},$
(4a) $\displaystyle\mathbf{K}_{2}$
$\displaystyle=\frac{1}{3}\mathbf{a}_{1}^{\star}-\frac{1}{3}\mathbf{a}_{2}^{\star}=-\mathbf{K}_{5},$
(4b) $\displaystyle\mathbf{K}_{6}$
$\displaystyle=\frac{1}{3}\mathbf{a}_{1}^{\star}+\frac{2}{3}\mathbf{a}_{2}^{\star}=-\mathbf{K}_{3}.$
(4c)

These six points are classified into two equivalent classes
$[\mathbf{K}_{\xi}=\mathbf{K}_{1}]=\\{\mathbf{K}_{1},\mathbf{K}_{3},\mathbf{K}_{5}\\}$
with $\xi=-1$, and
$[\mathbf{K}_{\xi}=\mathbf{K}_{4}]=\\{\mathbf{K}_{2},\mathbf{K}_{4},\mathbf{K}_{6}\\}$
with $\xi=+1$. The Brillouin zone has the area of
$S_{\text{BZ}}=\|\mathbf{a}_{1}^{\star}\times\mathbf{a}_{2}^{\star}\|=(4\pi/\sqrt{3}a)^{2}\sin(2\pi/3)$.

The position of the carbon atoms on the $\ell$ graphene layer ($\ell=1,2$) are
determined by the vectors:

$\displaystyle\mathbf{x}_{\alpha}^{(\ell)}$
$\displaystyle=\mathbf{R}^{(\ell)}+\mathbf{d}_{\alpha}^{(\ell)}+(\ell-1)\boldsymbol{\tau},$
(5)

where $\alpha=A,B$ and

$\mathbf{R}^{(\ell)}=m\mathbf{a}_{1}^{(\ell)}+n\mathbf{a}_{2}^{(\ell)};m,n\in\mathbb{Z}$
(6)

is the vector defining the Bravais lattice $\Gamma^{(\ell)}$ of the $\ell$
graphene layer, and
$\boldsymbol{\tau}=\eta\mathbf{a}_{1}^{(2)}+\zeta\mathbf{a}_{2}^{(2)}$ with
$\eta,\zeta\in[0,1)$ is the displacement vector. The two basis vectors
$\mathbf{a}_{1,2}^{(\ell)}$ of $\Gamma^{(\ell)}$ are obtained by rotating
appropriately the vectors $\mathbf{a}_{1,2}$:

$\mathbf{a}_{1,2}^{(\ell)}=R_{z}^{(\ell)}\cdot\mathbf{a}_{1,2},$ (7)

where $R_{z}^{(\ell)}$ is a matrix representing the rotation operation of an
angle $(-1)^{\ell}\theta/2$ around the twist $Oz$ axis. Denote
$\Gamma^{\star(\ell)}$ the reciprocal lattice associated with the Bravais
lattice $\Gamma^{(\ell)}$. It is a set of vectors given by:

$\mathbf{G}^{(\ell)}=m\mathbf{a}_{1}^{\star(\ell)}+n\mathbf{a}_{2}^{\star(\ell)};m,n\in\mathbb{Z},$
(8)

where $\mathbf{a}_{1,2}^{\star(\ell)}$ are the two basis vectors defined by:

$\mathbf{a}_{1,2}^{\star(\ell)}=R_{z}^{(\ell)}\cdot\mathbf{a}_{1,2}^{\star}.$
(9)

Actually, under the twist, all vectors associated with the graphene layer
$\ell$ are transformed via Eqs.(5) and (9). For instance, vector
$\mathbf{d}_{\alpha}^{(\ell)}$ in Eq. (5) is given by
$\mathbf{d}_{\alpha}^{(\ell)}=R_{z}^{(\ell)}\cdot\mathbf{d}_{\alpha_{\ell}}$.

There is a special case of which the bilayer graphene systems is in the
commensurate stacking between two graphene layers. It means that we can find
two pairs of integer numbers $(m,n)$ and $(p,q)$ such that:

$m\mathbf{a}_{1}^{(1)}+n\mathbf{a}_{2}^{(1)}=p\mathbf{a}_{1}^{(2)}+q\mathbf{a}_{2}^{(2)}.$
(10)

It can be shown that this requirement is satisfied not for any twist angle
$\theta$, but only for the ones satisfying the conditions $p=n,q=m$ and:

$\tan\theta=\frac{(n^{2}-m^{2})\sin(\pi/3)}{(n^{2}+m^{2})\cos(\pi/3)+2mn}.$
(11)

Treating Eq. (10) rigorously may lead to the so-called Diophantine problem as
presented in Ref. [Shallcross  _et al._ , 2008]. The commensurate stacking,
characterized by Eq. (10) or Eq. (11), implies that the TBG lattice has a
translation symmetry, which is defined by two basis vectors:

$\displaystyle\mathbf{A}_{1}$
$\displaystyle=m\mathbf{a}_{1}^{(1)}+n\mathbf{a}_{2}^{(1)}=n\mathbf{a}_{1}^{(2)}+m\mathbf{a}_{2}^{(2)},$
(12a) $\displaystyle\mathbf{A}_{2}$
$\displaystyle=n\mathbf{a}_{1}^{(1)}-(m+n)\mathbf{a}_{2}^{(1)}=m\mathbf{a}_{1}^{(2)}-(m+n)\mathbf{a}_{2}^{(2)}.$
(12b)

A super-cell defined by these two vectors has the area of
$S_{sc}=\|\mathbf{A}_{1}\times\mathbf{A}_{2}\|=L_{A}^{2}\sin(2\pi/3)$, where
$L_{A}=\|\mathbf{A}_{1}\|=\|\mathbf{A}_{2}\|=a\sqrt{m^{2}+mn+n^{2}}$. We
denote $\Gamma_{M}$ the Bravais lattice of the commensurate TBG lattice.
$\Gamma_{M}$ therefore is a set of vectors given by
$\mathbf{R}=m\mathbf{A}_{1}+n\mathbf{A}_{2}$ with $m,n\in\mathbb{Z}$. The
associated reciprocal lattice $\Gamma_{M}^{\star}$ is defined by two basis
vectors:

$\displaystyle\mathbf{A}_{1}^{\star}$
$\displaystyle=\frac{a^{2}}{L_{A}^{2}}\left[(m+n)\mathbf{a}_{1}^{\star(1)}+n\mathbf{a}_{2}^{\star(1)}\right]$
(13a)
$\displaystyle=\frac{a^{2}}{L_{A}^{2}}\left[(m+n)\mathbf{a}_{1}^{\star(2)}+m\mathbf{a}_{2}^{\star(2)}\right],$
(13b) $\displaystyle\mathbf{A}_{2}^{\star}$
$\displaystyle=\frac{a^{2}}{L_{A}^{2}}\left[n\mathbf{a}_{1}^{\star(1)}-m\mathbf{a}_{2}^{\star(1)}\right]$
(13c)
$\displaystyle=\frac{a^{2}}{L_{A}^{2}}\left[m\mathbf{a}_{1}^{\star(2)}-n\mathbf{a}_{2}^{\star(2)}\right]$
(13d)

The two vectors $\mathbf{A}_{1,2}^{\star}$ define a mini-Brillouin zone (MBZ)
with the area of
$S_{\text{MBZ}}=\|\mathbf{A}_{1}^{\star}\times\mathbf{A}_{2}^{\star}\|=(4\pi/\sqrt{3}a)^{2}\sin(\pi/3)/(m^{2}+mn+n^{2})=S_{\text{BZ}}/(m^{2}+mn+n^{2})$.
This result shows that the area of the Brillouin zone of monolayer graphene is
an integer multiple of that of the twisted bilayer lattice,
$S_{\text{BZ}}=(m^{2}+mn+n^{2})S_{\text{MBZ}}$. Thus, it suggests that a BZ
can be covered by $(m^{2}+mn+n^{2})$ pieces of MBZ. The six $K$ points of the
MBZ are given by:

$\displaystyle\mathbf{K}_{1}^{M}$
$\displaystyle=\frac{1}{3}(\mathbf{A}_{1}^{\star}+\mathbf{A}_{2}^{\star})=-\mathbf{K}_{4}^{M},$
(14a) $\displaystyle\mathbf{K}_{2}^{M}$
$\displaystyle=\frac{1}{3}(-\mathbf{A}_{1}^{\star}+2\mathbf{A}_{2}^{\star})=-\mathbf{K}_{5}^{M},$
(14b) $\displaystyle\mathbf{K}_{3}^{M}$
$\displaystyle=\frac{1}{3}(-2\mathbf{A}_{1}^{\star}+\mathbf{A}_{2}^{\star})=-\mathbf{K}_{6}^{M}.$
(14c)

For convenience, Eqs. (13a–13d) are inversely rewritten as follows:

$\displaystyle\mathbf{a}_{1}^{\star(1)}$
$\displaystyle=m\mathbf{A}_{1}^{\star}+n\mathbf{A}_{2}^{\star},$ (15a)
$\displaystyle\mathbf{a}_{2}^{\star(1)}$
$\displaystyle=n\mathbf{A}_{1}^{\star}-(m+n)\mathbf{A}_{2}^{\star},$ (15b)
$\displaystyle\mathbf{a}_{1}^{\star(2)}$
$\displaystyle=n\mathbf{A}_{1}^{\star}+m\mathbf{A}_{2}^{\star},$ (15c)
$\displaystyle\mathbf{a}_{2}^{\star(2)}$
$\displaystyle=m\mathbf{A}_{1}^{\star}-(m+n)\mathbf{A}_{2}^{\star}.$ (15d)

These equations remarkably take the form similar to those in Eqs. (12a, 12b).
Using these expressions we can verify this special relationship:

$\Delta\mathbf{K}=\mathbf{K}_{2}^{(1)}-\mathbf{K}_{2}^{(2)}=(m-n)(\mathbf{K}_{1}^{M}-\mathbf{K}_{2}^{M}).$
(16)

It means that the $\Delta\mathbf{K}$ is always a multiple of the MBZ edge.

Figure 1:  The Fourier image $T(\mathbf{q},d_{GG})$ of the model for the
electron hopping between two graphene layers:
$T(\mathbf{r}+z\mathbf{e}_{z})=T(\mathbf{R})=V_{pp\sigma}(R)(\mathbf{R}\cdot\mathbf{e}_{z}/R)^{2}+V_{pp\pi}(R)[1-(\mathbf{R}\cdot\mathbf{e}_{z}/R)^{2}]$
where $R=\|\mathbf{R}\|$ and
$V_{pp\sigma}(R)=V_{pp\sigma}^{0}\exp[-(R-d_{GG})/R_{0}]$,
$V_{pp\pi}(R)=V^{0}_{pp\pi}\exp[-(R-a_{CC})/R_{0}]$. The model parameters are
$V^{0}_{pp\pi}=-2.70$ eV, $V^{0}_{pp\sigma}=0.48$ eV, $a_{CC}=0.145$ nm,
$d_{GG}=0.335$ nm, and $R_{0}=0.148a$. Here $a=\sqrt{3}a_{CC}$ is the lattice
constant of the graphene lattice.

##  III Construction of effective models

###  III.1 Tight-binding description

In the tight-binding formalism the dynamics of an electron in an atomic
lattice is described by a Hamiltonian $H$ which acts on a Hilbert space
$\mathscr{H}$ spanned by a set of state vectors whose corresponding wave
functions are strongly localized in the lattice. For graphene, we particularly
consider a set of $p_{z}$-symmetry atomic orbitals
$\phi_{p_{z},\alpha}(\mathbf{x}-\mathbf{R}^{(\ell)}-\mathbf{d}_{\alpha}^{(\ell)})$
localized around the position $\mathbf{d}^{(\ell)}_{\alpha}$ in the graphene
unit cell $\mathbf{R}^{(\ell)}$ of the material layer $\ell$. For short we
denote $|\ell,\alpha,\mathbf{x}^{(\ell)}_{\alpha}\rangle$, where
$\mathbf{x}^{(\ell)}_{\alpha}=\mathbf{R}^{(\ell)}+\mathbf{d}^{(\ell)}_{\alpha}$,
the state vectors corresponding to
$\phi_{p_{z},\alpha}(\mathbf{x}-\mathbf{x}^{(\ell)}_{\alpha})$. As the
starting point we use the set of these vectors
$\\{|\ell,\alpha,\mathbf{R}^{(\ell)}+\mathbf{d}^{(\ell)}_{\alpha}\rangle\,|\,\ell=1,2;\alpha=A,B;\mathbf{R}^{(\ell)}\in\Gamma^{(\ell)}\\}$
as a representation basis, where $\Gamma^{(\ell)}$ denotes the Bravais lattice
of graphene layer $\ell$. Since each graphene layer possesses the translation
symmetry, we can express $|\ell,\alpha,\mathbf{x}^{(\ell)}_{\alpha}\rangle$ in
terms of the Fourier transform of the Bloch vectors
$|\ell,\alpha,\mathbf{k}^{(\ell)}\rangle$ as follows:

$|\ell,\alpha,\mathbf{x}^{(\ell)}_{\alpha}\rangle=\frac{1}{\sqrt{N}}\sum_{\mathbf{k}^{(\ell)}}^{\text{BZ}_{\ell}}e^{i\mathbf{k}^{(\ell)}\cdot\mathbf{x}^{(\ell)}_{\alpha}}|\ell,\alpha,\mathbf{k}^{(\ell)}\rangle.$
(17)

Inversely, we have the equation:

$|\ell,\alpha,\mathbf{k}^{(\ell)}\rangle=\frac{1}{\sqrt{N}}\sum_{\mathbf{R}^{(\ell)}}^{\Gamma_{\ell}}e^{-i\mathbf{k}^{(\ell)}\cdot\mathbf{x}^{(\ell)}_{\alpha}}|\ell,\alpha,\mathbf{x}^{(\ell)}_{\alpha}\rangle.$
(18)

The single-electron Hamiltonian in the atomic lattice is, in general,
determined as the summation of energies associated with all possible hopping
processes of electron from a state
$|\ell,\alpha,\mathbf{x}^{(\ell)}_{\alpha}\rangle$ to another
$|\ell^{\prime},\beta,\mathbf{x}^{(\ell^{\prime})}_{\beta}\rangle$.
Accordingly, we can write down the general expression for $H$ as follows:

$\displaystyle H$ $\displaystyle=$
$\displaystyle\sum_{\ell,\ell^{\prime}}\sum_{\alpha,\beta}\sum_{\mathbf{R}^{(\ell)},\mathbf{R}^{(\ell^{\prime})}}|\ell,\alpha,\mathbf{x}^{(\ell)}_{\alpha}\rangle
H(\mathbf{x}^{(\ell)}_{\alpha},\mathbf{x}_{\beta}^{(\ell^{\prime})})\langle\ell^{\prime},\beta,\mathbf{x}_{\beta}^{(\ell^{\prime})}|.$

Here we denote
$H(\mathbf{x}^{(\ell)}_{\alpha},\mathbf{x}^{(\ell^{\prime})}_{\beta})=\langle\ell,\alpha,\mathbf{x}^{(\ell)}_{\alpha}|H|\ell^{\prime},\beta,\mathbf{x}^{(\ell^{\prime})}_{\beta}\rangle$
the matrix elements of the Hamiltonian in the chosen basis. These matrix
elements are interpreted as the hopping integrals. All terms in Eq. (III.1)
can be arranged into three parts: $H_{1},H_{2}$ include the terms defining the
coupling between Bloch states in the same graphene layer, and $H_{12}$
includes the terms of the coupling between Bloch states in two different
graphene layers. In the tight-binding description we approximate these
quantities as the two-center integrals, i.e.,
$H(\mathbf{x}^{(\ell)}_{\alpha},\mathbf{x}^{(\ell^{\prime})}_{\beta})=T(\mathbf{x}^{(\ell)}_{\alpha}-\mathbf{x}^{(\ell^{\prime})}_{\beta})\equiv
T(d_{\alpha\beta}^{(\ell\ell^{\prime})})$ where
$d_{\alpha\beta}^{(\ell\ell^{\prime})}=\|\mathbf{d}_{\alpha\beta}^{(\ell\ell^{\prime})}\|$
and
$\mathbf{d}_{\alpha\beta}^{(\ell\ell^{\prime})}=\mathbf{x}^{(\ell)}_{\alpha}-\mathbf{x}_{\beta}^{(\ell^{\prime})}$.
For numerical calculations we use the following Slater-Koster model for the
hopping integrals:

$T(\mathbf{d})=V_{pp\pi}(d)\sin^{2}\theta^{z}+V_{pp\sigma}(d)\cos^{2}\theta^{z},$
(20)

where $\cos\theta^{z}=(\mathbf{d}\cdot\mathbf{e}_{z})/d$ and

$\displaystyle V_{pp\pi}(d)$ $\displaystyle=$ $\displaystyle
V^{0}_{pp\pi}\exp\left(-\frac{d-a_{CC}}{r_{0}}\right),$ (21a) $\displaystyle
V_{pp\sigma}(d)$ $\displaystyle=$ $\displaystyle
V^{0}_{pp\sigma}\exp\left(-\frac{d-d_{GG}}{r_{0}}\right).$ (21b)

In this model, the parameters are commonly set to $V^{0}_{pp\pi}=-2.7$ eV,
$V_{pp\sigma}^{0}=0.48$ eV, and $r_{0}=0.148a$.Moon and Koshino (2013);
Koshino (2015); Le and Do (2018); Le  _et al._ (2019); Do  _et al._ (2019)

Using the expression given by Eq. (17), the Hamiltonian is rewritten in the
form:

$\displaystyle
H=\sum_{\ell,\ell^{\prime}}\sum_{\alpha,\beta}|\ell,\alpha,\mathbf{k}^{(\ell)}\rangle\sum_{\mathbf{k}^{(\ell)},\mathbf{k}^{(\ell^{\prime})}}^{\text{BZ}_{\ell\ell^{\prime}}}H_{\alpha\beta}(\mathbf{k}^{(\ell)},\mathbf{k}^{(\ell^{\prime})})\langle\ell^{\prime},\beta,\mathbf{k}^{(\ell^{\prime})}|,$
(22)

where we denote
$H_{\alpha\beta}(\mathbf{k}^{(\ell)},\mathbf{k}^{(\ell^{\prime})})=\langle\ell,\alpha,\mathbf{k}^{(\ell)}|H|\ell^{\prime},\beta,\mathbf{k}^{(\ell^{\prime})}\rangle$
the matrix elements of the Hamiltonian in the basis set of Bloch vectors. The
expression of
$H_{\alpha\beta}(\mathbf{k}^{(\ell)},\mathbf{k}^{(\ell^{\prime})})$ reads:

$\displaystyle
H_{\alpha\beta}(\mathbf{k}^{(\ell)},\mathbf{k}^{(\ell^{\prime})})=$
$\displaystyle\frac{1}{N}\sum_{\mathbf{R}^{(\ell)},\mathbf{R}^{(\ell^{\prime})}}T(\mathbf{x}_{\alpha}^{(\ell)}-\mathbf{x}_{\beta}^{(\ell^{\prime})})$
$\displaystyle\hskip 28.45274pt\times
e^{-i\mathbf{k}^{(\ell)}\cdot\mathbf{x}_{\alpha}^{(\ell)}}e^{i\mathbf{k}^{(\ell^{\prime})}\cdot\mathbf{x}_{\beta}^{(\ell^{\prime})}}.$
(23)

In the following we will derive rules for the coupling between Bloch states.
To do so, in the expression of $H$, we distinguish the so-called intra-layer
coupling terms, $H_{\ell}$, including two terms when $\ell=\ell^{\prime}=1$
and 2, and the so-called inter-layer coupling term, denoted by $H_{12}$, when
$\ell=1,\ell^{\prime}=2$ and $\ell=2,\ell^{\prime}=1$.

###  III.2 Rule for the intra-layer coupling

We now try to simplify the expression of
$H_{\alpha\beta}(\mathbf{k}^{(\ell)},\mathbf{k}^{(\ell^{\prime})})$ defined by
Eq. (23). For the $H_{\ell}$ terms with $\ell^{\prime}=\ell$, we replace
$\mathbf{k}^{(\ell^{\prime})}$ by $\mathbf{k}^{(\ell)\prime}$. The summation
over $\mathbf{x}_{\alpha}^{(\ell)}$ and $\mathbf{x}_{\beta}^{(\ell)}$ can be
cast into the product of two independent sums:

$\displaystyle H_{\alpha\beta}(\mathbf{k}^{(\ell)},\mathbf{k}^{(\ell)\prime})$
$\displaystyle=\left(\frac{1}{N}\sum_{\mathbf{R}^{(\ell)}}e^{-i(\mathbf{k}^{(\ell)}-\mathbf{k}^{(\ell)\prime})\cdot\mathbf{x}_{\alpha}^{(\ell)}}\right)$
$\displaystyle\times\left(\sum_{\mathbf{R}^{(\ell)}}T(\mathbf{x}_{\alpha}^{(\ell)}-\mathbf{x}_{\beta}^{(\ell)})e^{-i\mathbf{k}^{(\ell)\prime}\cdot(\mathbf{x}_{\alpha}^{(\ell)}-\mathbf{x}_{\beta}^{(\ell)})}\right).$
(24)

For the first sum, by noticing Eq. (5) we generally have:

$\frac{1}{N}\sum_{\mathbf{R}^{(\ell)}}e^{-i(\mathbf{k}^{(\ell)}-\mathbf{k}^{(\ell)\prime})\cdot\mathbf{x}_{\alpha}^{(\ell)}}=\sum_{\mathbf{G}^{(\ell)}}\delta_{\mathbf{k}^{(\ell)}-\mathbf{k}^{(\ell)\prime},\mathbf{G}^{(\ell)}}e^{i\mathbf{G}^{(\ell)}\cdot\mathbf{d}_{\alpha}^{(\ell)}}.$
(25)

However, since both $\mathbf{k}^{(\ell)}$ and $\mathbf{k}^{(\ell)\prime}$ are
restricted in the first Brillouin zone BZℓ, only the term with
$\mathbf{G}^{(\ell)}=0$ contributes to the sum. It means that:

$\frac{1}{N}\sum_{\mathbf{R}^{(\ell)}}e^{-i(\mathbf{k}^{(\ell)}-\mathbf{k}^{(\ell)\prime})\cdot\mathbf{x}_{\alpha}^{(\ell)}}=\delta_{\mathbf{k}^{(\ell)},\mathbf{k}^{(\ell)\prime}}.$
(26)

The result implies that the dynamics of electron in each individual graphene
layer is governed by only the coupling of Bloch states of the same wave vector
$\mathbf{k}^{(\ell)}$. In other words, given a Bloch state with the wave
vector $\mathbf{k}^{(\ell)}$, it only couples with other Bloch states of the
same wave vector to form the intra-coupling Hamiltonian terms.

For the second sum, by setting
$\mathbf{d}_{\alpha\beta}^{(\ell)}=\mathbf{x}_{\alpha}^{(\ell)}-\mathbf{x}_{\beta}^{(\ell)}$,
it is sufficiently interpreted as the summation over some neighborhoods of the
state located at $\mathbf{x}_{\alpha}^{(\ell)}$. Thanks to the “coupling rule”
just deduced, we directly calculate the second summation by counting $N_{n}$
neighbors at $\\{\mathbf{d}_{\alpha\beta}^{(\ell)j}\\}$ of each
$\mathbf{x}_{\alpha}^{(\ell)}$. For simplicity we set
$H^{(\ell)}_{\alpha\beta}(\mathbf{k})=H_{\alpha\beta}(\mathbf{k}^{(\ell)},\mathbf{k}^{(\ell)})$
where

$\displaystyle H^{(\ell)}_{\alpha\beta}(\mathbf{k}^{(\ell)})$ $\displaystyle=$
$\displaystyle\sum_{j=1}^{N_{n}}T(\mathbf{d}^{(\ell)j}_{\alpha\beta})e^{-i\mathbf{k}^{(\ell)}\cdot\mathbf{d}_{\alpha\beta}^{(\ell)j}}.$
(27)

Technically, the quantities $H_{\alpha\beta}^{(\ell)}(\mathbf{k}^{(\ell)})$
defined by this equation are the Fourier transforms of the hopping integrals
$T(\mathbf{x})$. They comprise the elements of the Hamiltonian matrix in the
representation of Bloch state vectors. Treating the coupling matrix as above
and now substituting expressions (26) and (27) into Eq. (22) we obtain the
expression for the intra-layer coupling Hamiltonian:

$H_{\ell}=\sum_{\mathbf{k}^{(\ell)}}^{\text{BZ}_{\ell}}\sum_{\alpha,\beta}|\ell,\alpha,\mathbf{k}^{(\ell)}\rangle
H^{(\ell)}_{\alpha\beta}(\mathbf{k}^{(\ell)})\langle\ell,\beta,\mathbf{k}^{(\ell)}|.$
(28)

This Hamiltonian allows to state the following rule for the intra-coupling
($\bigstar$): Two Bloch states $|\ell,\alpha,\mathbf{k}^{(\ell)}\rangle$ and
$|\ell,\beta,\mathbf{k}^{(\ell)}\rangle$ always couple to each other with the
strength $H_{\alpha\beta}^{(\ell)}(\mathbf{k}^{(\ell)})$.

###  III.3 Rule for the inter-layer coupling

For the inter-layer coupling terms of the Hamiltonian, the expression of
$H_{\alpha\beta}(\mathbf{k}^{(\ell)},\mathbf{k}^{(\ell^{\prime})})$ given in
Eq. (23), in principle, can be treated in the form of Eq. (III.2). However,
since $\mathbf{k}^{(\ell)}$ and $\mathbf{k}^{(\ell^{\prime})}$ belong to
different Brillouin zones, the first sum should be given by Eq. (25) rather
than by Eq. (26). The expression is thus not further simplified. Importantly,
in order to symmetrize the expression of
$H_{\alpha\beta}(\mathbf{k}^{(\ell)},\mathbf{k}^{(\ell^{\prime})})$ we should
note the localization of the hopping integrals. It allows us to express
$T(\mathbf{x}_{\alpha}^{(\ell)}-\mathbf{x}_{\beta}^{(\ell^{\prime})})$ in the
form of a Fourier series:

$T(\mathbf{x}_{\alpha}^{(\ell)}-\mathbf{x}_{\beta}^{(\ell^{\prime})})=\frac{1}{N}\sum_{\mathbf{q}\in\mathbb{R}^{2}}T(\mathbf{q},d_{GG})e^{i\mathbf{q}\cdot(\mathbf{x}_{\alpha}^{(\ell)}-\mathbf{x}_{\beta}^{(\ell^{\prime})})},$
(29)

where $\mathbf{q}$ is the two-dimensional vector, $d_{GG}=z_{2}-z_{1}$, and
$T(\mathbf{q},d_{GG})$ is inversely determined by

$T(\textbf{q})=\frac{1}{S_{c}}\int_{S}d^{2}\mathbf{x}T(\mathbf{x},z)e^{-i\mathbf{q}\cdot\mathbf{x}},$
(30)

where $S_{c}=S/N$ is the area of the unit cell of the $\ell$ graphene layer.
Using Eq. (29) the expression of
$H_{\alpha\beta}(\mathbf{k}^{(\ell)},\mathbf{k}^{(\ell^{\prime})})$ is now
manipulated as follows:

$\displaystyle
H_{\alpha\beta}(\mathbf{k}^{(\ell)},\mathbf{k}^{(\ell^{\prime})})$
$\displaystyle=\sum_{\mathbf{q}\in\mathbb{R}^{2}}T(\mathbf{q},d_{GG})\frac{1}{N^{2}}\sum_{\mathbf{R}^{(\ell)},\mathbf{R}^{(\ell^{\prime})}}e^{-i\mathbf{k}^{(\ell)}\cdot\mathbf{x}_{\alpha}^{(\ell)}}e^{i\mathbf{q}\cdot(\mathbf{x}_{\alpha}^{(\ell)}-\mathbf{x}_{\beta}^{(\ell^{\prime})})}e^{i\mathbf{k}^{(\ell^{\prime})}\cdot\mathbf{x}_{\beta}^{(\ell^{\prime})}}$
$\displaystyle=\sum_{\mathbf{q}\in\mathbb{R}^{2}}T(\mathbf{q},d_{GG})\left(\frac{1}{N}\sum_{\mathbf{R}^{(\ell)}}e^{i(\mathbf{q}-\mathbf{k}^{(\ell)})\cdot\mathbf{x}_{\alpha}^{(\ell)}}\right)\left(\frac{1}{N}\sum_{\mathbf{R}^{(\ell^{\prime})}}e^{-i(\mathbf{q}-\mathbf{k}^{(\ell^{\prime})})\cdot\mathbf{x}_{\beta}^{(\ell^{\prime})}}\right)$
$\displaystyle=\sum_{\mathbf{q}\in\mathbb{R}^{2}}\sum_{\mathbf{G}^{(\ell)},\mathbf{G}^{(\ell^{\prime})}}T(\mathbf{q},d_{GG})e^{i\mathbf{G}^{(\ell)}\cdot\mathbf{d}_{\alpha}^{(\ell)}}e^{-i\mathbf{G}^{(\ell^{\prime})}\cdot\mathbf{d}_{\beta}^{(\ell^{\prime})}}\delta_{\mathbf{q}-\mathbf{k}^{(\ell)},\mathbf{G}^{(\ell)}}\delta_{\mathbf{q}-\mathbf{k}^{(\ell^{\prime})},\mathbf{G}^{(\ell^{\prime})}}$
(31a)
$\displaystyle=\sum_{\mathbf{G}^{(\ell)},\mathbf{G}^{(\ell^{\prime})}}T(\mathbf{k}^{(\ell)}+\mathbf{G}^{(\ell)},d_{GG})e^{i\mathbf{G}^{(\ell)}\cdot\mathbf{d}_{\alpha}^{(\ell)}}e^{-i\mathbf{G}^{(\ell^{\prime})}\cdot\mathbf{d}_{\beta}^{(\ell^{\prime})}}\delta_{\mathbf{k}^{(\ell)}+\mathbf{G}^{(\ell)},\mathbf{k}^{(\ell^{\prime})}+\mathbf{G}^{(\ell^{\prime})}}$
(31b)

The appearance of the Krönecker symbol in Eqs. (31a) and (31b) is remarkable.
It implies that not any two states $|\ell,\alpha,\mathbf{k}^{(\ell)}\rangle$
and $|\ell^{\prime},\beta,\mathbf{k}^{(\ell^{\prime})}\rangle$ can couple to
each other, but only the ones defined by the vectors $\mathbf{k}^{(\ell)}$ and
$\mathbf{k}^{(\ell^{\prime})}$ related by the condition
$\mathbf{k}^{(\ell)}+\mathbf{G}^{(\ell)}=\mathbf{k}^{(\ell^{\prime})}+\mathbf{G}^{(\ell^{\prime})}$,
where $\mathbf{G}^{(\ell)}\in\Gamma_{\ell}^{\star}$ and
$\mathbf{G}^{(\ell^{\prime})}\in\Gamma_{\ell^{\prime}}^{\star}$ are some
vectors of the reciprocal lattices $\Gamma_{\ell}^{\star}$ and
$\Gamma_{\ell^{\prime}}^{\star}$, respectively. This is the rule for the
coupling between Bloch states in different graphene layers.

For short we denote
$M_{\alpha\beta}(\mathbf{G}^{(\ell)},\mathbf{G}^{(\ell^{\prime})})=e^{i\mathbf{G}^{(\ell)}\cdot\mathbf{d}_{\alpha}^{(\ell)}}e^{-i\mathbf{G}^{(\ell^{\prime})}\cdot\mathbf{d}_{\beta}^{(\ell^{\prime})}}$,
then substitute Eqs. (31a) and (31b) into Eq. (22). We obtain two equivalent
representations for the Hamiltonian term $H_{12}$:

$\displaystyle H_{12}=$
$\displaystyle\sum_{\mathbf{q}\in\mathbb{R}^{2}}\sum_{\alpha,\beta}|1,\alpha,\underbrace{\mathbf{q}-\mathbf{G}^{(1)}}_{\mathbf{k}^{(1)}\in\text{BZ}_{1}}\rangle
T(\mathbf{q},d_{GG})M_{\alpha\beta}({\mathbf{G}^{(1)},\mathbf{G}^{(2)}\in\text{BZ}_{2}})\langle
2,\beta,\underbrace{\mathbf{q}-\mathbf{G}^{(2)}}_{\mathbf{k}^{(2)}\in\text{BZ}_{2}}|+$
$\displaystyle\sum_{\mathbf{q}\in\mathbb{R}^{2}}\sum_{\alpha,\beta}|2,\beta,\underbrace{\mathbf{q}-\mathbf{G}^{(2)}}_{\mathbf{k}^{(2)}\in\text{BZ}_{2}}\rangle
T^{\star}(\mathbf{q},d_{GG})M^{\star}_{\alpha\beta}({\mathbf{G}^{(1)},\mathbf{G}^{(2)}})\langle
1,\alpha,\underbrace{\mathbf{q}-\mathbf{G}^{(1)}}_{\mathbf{k}^{(1)}\in\text{BZ}_{1}}|,$
(32) $\displaystyle H_{12}=$
$\displaystyle\sum_{\mathbf{\mathbf{k}^{(1)}}}^{\text{BZ}_{1}}\sum_{\alpha,\beta}|1,\alpha,\mathbf{k}^{(1)}\rangle\left[\sum_{\mathbf{G}^{(1)}}T(\mathbf{k}^{(1)}+\mathbf{G}^{(1)})M_{\alpha\beta}(\mathbf{G}^{(1)},\mathbf{G}^{(2)})\langle
2,\beta,\underbrace{\mathbf{k}^{(1)}+\Delta\mathbf{G}^{(12)}}_{\mathbf{k}^{(2)}\in\text{BZ}_{2}}|\right]+$
$\displaystyle\sum_{\mathbf{\mathbf{k}^{(2)}}}^{\text{BZ}_{2}}\sum_{\alpha,\beta}|2,\beta,\mathbf{k}^{(2)}\rangle\left[\sum_{\mathbf{G}^{(2)}}T(\mathbf{k}^{(2)}+\mathbf{G}^{(2)})M_{\beta\alpha}(\mathbf{G}^{(2)},\mathbf{G}^{(1)})\langle
1,\alpha,\underbrace{\mathbf{k}^{(2)}+\Delta\mathbf{G}^{(21)}}_{\mathbf{k}^{\prime(1)}\in\text{BZ}_{1}}|\right],$
(33)

wherein $\Delta\mathbf{G}^{(12)}=\mathbf{G}^{(1)}-\mathbf{G}^{(2)}$. Eq. (32)
is obtained by summing up over the vectors $\mathbf{k}^{(1)}$ and
$\mathbf{k}^{(2)}$. The sum over $\mathbf{G}^{(1)}$ and $\mathbf{G}^{(2)}$ can
be neglected because of the constrain that $\mathbf{q}-\mathbf{G}^{(\ell)}$
must be a point in the Brillouin zone $\text{BZ}_{\ell}$. In other words,
$\mathbf{G}^{(1)}$ and $\mathbf{G}^{(2)}$ are two specific vectors
accompanying with each value of $\mathbf{q}$ such that this condition
$\mathbf{q}-\mathbf{G}^{(1,2)}\in\text{BZ}_{1,2}$ is satisfied. This equation
suggests that we always find two Bloch states coupling together by starting
from a point $\mathbf{q}$ in the plane $\mathbb{R}^{2}$ and then projecting it
into the first Brillouin zone of each graphene layer to determine the
corresponding vectors $\mathbf{k}^{(1)}$ and $\mathbf{k}^{(2)}$. In the same
sense, the vector $\mathbf{G}^{(2)}$ in Eq. (33) is a specific vector
accompanying with each value of $\mathbf{k}^{(1)}$ and $\mathbf{G}^{(1)}$ to
map the vector $\mathbf{k}^{(1)}+\mathbf{G}^{(1)}$ back into the first
Brillouin zone $\text{BZ}_{2}$ of the second graphene layer, see Fig. 2(a).
Eq. (33) though does not have the symmetric form, it clearly shows a practical
rule for the coupling between Bloch states in two graphene layers: for each
Bloch state $|\ell,\alpha,\mathbf{k}^{(\ell)}\in\Gamma_{\ell}^{\star}\rangle$
in layer $\ell$ we always find a set of Bloch states
$\\{|\ell^{\prime},\beta,\mathbf{k}^{(\ell)}+\Delta\mathbf{G}^{(\ell\ell^{\prime})}\rangle\,|\,\forall\,\mathbf{G}^{(\ell)}\in\Gamma_{\ell}^{\star};\mathbf{G}^{(\ell^{\prime})}\in\Gamma_{\ell^{\prime}}^{\star};\mathbf{k}^{(\ell)}+\Delta\mathbf{G}^{(\ell\ell^{\prime})}\in\text{BZ}_{\ell^{\prime}}\\}$
in the layer $\ell^{\prime}\neq\ell$ that couple to the given state
($\bigstar\bigstar$).

Figure 2:  (a) Illustration for the coupling of one $\mathbf{k}^{(1)}$ point
(the red dot) in the first Brillouin zone BZ1 (the red hexagon) to three
points $\mathbf{k}^{(2)}$ (the blue dots) in the Brillouin zone BZ2 (the blue
hexagon). The circles imply the domains $\mathscr{O}(\mathbf{K}^{(1)}_{\xi})$
and $\mathscr{O}(\mathbf{K}^{(2)}_{\xi})$. (b) Lattice representation for the
coupling of single-layer Bloch states, see Eq. (III.5). The red/blue
square/round nodes represent the Bloch states of the first/second graphene
layer, respectively. (c) Schema for the map of the $\mathbf{q}$ points in the
first Brillouin zone of the TBG lattice (the small green hexagon in the
center) to the points in the $\mathscr{O}(\mathbf{K}^{(\ell)})$ domains. The
map is defined by
$\mathbf{q}\in\text{MBZ}\mapsto\mathbf{q}+\mathbf{G}^{M}_{m}$. The
parameterization procedure using $\mathbf{q}\in\text{MBZ}$ is appropriate to
the commensurate TBG configurations.

###  III.4 Building effective models

In this sub-section we present technical points in building an effective
Hamiltonian. Equation (32) or (33), together with Eq. (28), in principle,
allow constructing a Hamiltonian matrix to describe the dynamics of electrons
in its whole spectrum. Practically, we are usually interested in the behavior
of electrons in a certain energy range, for example, the energy range around
the Fermi energy level. It is thus expected to extract from the total
Hamiltonian matrix an appropriate block whose spectrum is approximate to the
exact one in the interested energy range. Technically, we need to isolate a
sub-set of vectors from the basis set that dominantly couple together. The
selection rule ($\bigstar\bigstar$) obviously helps to do so. Eq. (32) shows
how to recognize Bloch state vectors that couple together. They are the
vectors $|1,\alpha,\mathbf{q}-\mathbf{G}^{(1)}\rangle$ and
$|2,\beta,\mathbf{q}-\mathbf{G}^{(2)}\rangle$ that are parameterized by a
common vector $\mathbf{q}\in\mathbb{R}^{2}$. However, the matter lies in the
fact that many different points $\mathbf{q}$ in the plane $\mathbb{R}^{2}$ can
be mapped into the same point in the Brillouin zone BZ1,2. Furthermore, though
$\mathbf{q}$ varies smoothly, $\mathbf{q}-\mathbf{G}^{(1,2)}$ does not due to
the discreteness of $\mathbf{G}^{(1,2)}$. Therefore, it can lead to a
situation that many terms involved in Bloch states of really different
energies can be arranged into the same block of the Hamiltonian matrix. It
hence does not help in constructing effective models to describe the dynamics
of electrons in certain narrow ranges of energy.

The second form of $H_{12}$ given by Eq. (33) shows that for a Bloch state of
one graphene layer there is a list of Bloch states of other graphene layer
coupling to the given state. It is useful to see Eq. (33) as the tight-binding
form of the Hamiltonian written for Bloch states centered at discrete
$\mathbf{k}$-points in the reciprocal space. It is thus impractical to
consider the coupling of the state defined at the point $\mathbf{k}$ with the
ones defined at the points $\mathbf{k}+\Delta\mathbf{G}^{(12)}$ far from the
former. So, let us first discuss a technique to determine a set of finite
number of vectors $\mathbf{G}^{(1)}$, and the corresponding
$\mathbf{G}^{(2)}$, that allows approximating well the sum in Eq. (33). In
Fig. 1 we show the variation of the coupling function $T(\mathbf{q},d_{GG})$
given by Eqs. (20,21a,21b) with respect to the length of the vector
$\mathbf{q}$. Accordingly, $T(\mathbf{q},d_{GG})$ decays rapidly. So, it is
possible to define a cutoff value $G_{cutoff}$, and $\mathbf{G}^{(1)}$ are
found as the points $\mathbf{k}^{(1)}+\mathbf{G}^{(1)}$ restricted inside the
circle centered at $\mathbf{k}^{(1)}$ with the radius $G_{cutoff}$. Denote $N$
the number of such $\mathbf{G}^{(1)}$ vectors. Similarly, $N$ vectors of
$\mathbf{G}^{(2)}$ are also found such that
$\mathbf{k}^{(1)}+\Delta\mathbf{G}^{(12)}\in$ BZ2. For convenience, these
$\mathbf{G}^{(\ell)}$ vectors are labelled by the subscript $j$:
$\\{\mathbf{G}^{(\ell)}_{j}\,|\,j=1,2,\ldots,N,\|\mathbf{G}^{(\ell)}_{j}\|\leq
G_{max}\\}$ with $\ell=1,2$. Since the value of
$T(\mathbf{k}^{(1)}+\mathbf{G}^{(1)}_{j})$ is assumed to depend only on the
length of $\mathbf{k}^{(1)}+\mathbf{G}^{(1)}_{j}$, $N$ vectors in the set
$\\{\mathbf{G}^{(\ell)}_{j}\\}$ can be further classified in classes
$\mathscr{G}_{1}^{(1)},\mathscr{G}_{2}^{(1)},\mathscr{G}_{3}^{(1)},\ldots$
whose vectors in each class make the vectors
$\mathbf{k}^{(1)}+\mathbf{G}_{j}^{(1)}$ having the same length. For example,
if $\mathbf{k}^{(1)}$ points to a corner point $\mathbf{K}^{(1)}_{\xi}$ of the
Brillouin zone BZ1, then there are three vectors
$\mathbf{G}_{1}^{(1)}=0,\mathbf{G}_{2}^{(1)}=\xi\mathbf{a}_{1}^{\star(1)}$ and
$\mathbf{G}_{3}^{(1)}=\xi(\mathbf{a}_{1}^{\star(1)}+\mathbf{a}_{2}^{\star(1)})$
that can be classified into the class $\mathscr{G}^{(1)}$ since the three
vectors $\mathbf{K}^{(1)}_{\xi}+\mathbf{G}_{j}^{(1)}$ take the same length of
$\|\mathbf{K}_{\xi}^{(1)}\|=4\pi/3a$. Following are three classes
$\mathscr{G}_{j}^{(1)}$ for the point $\mathbf{K}^{(1)}_{\xi}$. The vectors
$\mathbf{G}_{j}^{(1)}$ are given in the format of a pair of integer numbers
$(p,q)$ standing for
$\mathbf{G}_{j}^{(1)}=p\mathbf{a}_{1}^{\star(1)}+q\mathbf{a}_{2}^{\star(1)}$:

$\displaystyle\mathscr{G}_{1\xi}^{(1)}=$
$\displaystyle(0,0),\xi(1,0),\xi(1,1)$ (34a)
$\displaystyle\mathscr{G}_{2\xi}^{(1)}=$
$\displaystyle\xi(0,1),\xi(0,-1),\xi(2,1)$ (34b)
$\displaystyle\mathscr{G}_{3\xi}^{(1)}=$
$\displaystyle\xi(-1,0),\xi(-1,-1),\xi(1,-1),$
$\displaystyle\xi(2,0),\xi(2,2),\xi(1,2)$ (34c)

Similarly, the set of $N$ vectors $\mathbf{G}_{j}^{(2)}$ is also arranged into
the classes $\mathscr{G}_{j}^{(2)}$ in the correspondence to the classes
$\mathscr{G}_{j}^{(1)}$.

We now discuss in detail techniques in the construction of an effective model
that can work for a certain narrow energy range $(E_{min},E_{max})$ of the
whole spectrum. Denote $\mathscr{O}(\mathbf{K}^{(\ell)})$ a domain of
$\mathbf{k}^{(\ell)}$ centered at some $\mathbf{K}^{(\ell)}$ point such that
there exist single-layer Bloch states of energies $E_{n}(\mathbf{k}^{(\ell)})$
in the range $(E_{min},E_{max})$. We expect that the combination of the state
vectors
$\\{|1,\alpha,\mathbf{k}^{(1)}\rangle\,|\,\alpha=A,B,\mathbf{k}^{(1)}\in\mathscr{O}(\mathbf{K}^{(1)})\\}$
with the ones
$\\{|2,\beta,\mathbf{k}^{(2)}\rangle\,|\,\beta=A,B,\mathbf{k}^{(2)}\in\mathscr{O}(\mathbf{K}^{(2)})\\}$
will form the states of the complex system in the considered energy range. It
is suggested from Eq. (33) that the coupled Bloch state vectors, which can be
parameterized by a continuous vector, are collected by a procedure as follows:
(1) Starting from a state vector
$|1,\alpha,\mathbf{k}^{(1)}=\mathbf{k}\rangle$ with
$\mathbf{k}\in\mathscr{O}(\mathbf{K}^{(1)})$ we can find $N$ states
$\\{|2,\beta,\mathbf{k}^{(2)}=\mathbf{k}+\Delta\mathbf{G}^{(12)}_{j}\rangle,j=1,2,\ldots,N\\}$
that couple to the given one. Next, (2) for each state
$|2,\beta,\mathbf{k}+\Delta\mathbf{G}^{(12)}_{j}\in\mathscr{O}(\mathbf{K}^{(2)})\rangle$
we will find $N$ corresponding coupled states
$\\{|1,\alpha,\mathbf{k}^{(1)}=\mathbf{k}+\Delta\mathbf{G}^{(12)}_{j}-\Delta\mathbf{G}^{(12)}_{i}\rangle,\,i=1,2,\ldots,N\\}$.
Repeat this two-step procedure we will collect a desired subset of vectors
from the total basis set. Since $\mathbf{G}^{(\ell)}_{j}$ can be expanded in
terms of the basis vectors of the reciprocal lattice,
$\\{\mathbf{a}_{1}^{\star(\ell)},\mathbf{a}_{2}^{\star(\ell)}\\}$ we have
$\Delta\mathbf{G}^{(12)}_{j}=m_{1}\Delta\mathbf{a}_{1}^{\star(12)}+m_{2}\Delta\mathbf{a}_{2}^{\star(12)}$where
$m_{1,2}\in\mathbb{Z}$. We therefore realize that the coupled Bloch state
vectors related to Eq. (33) can be written in the much more instructive form
as follows:

$\displaystyle|\mathbf{K}^{(1)},\alpha,\mathbf{q},\mathbf{q}^{(1)}+\Delta\mathbf{G}^{(12)}_{i}\rangle,$
(35a)
$\displaystyle|\mathbf{K}^{(2)},\beta,\mathbf{q},\mathbf{q}^{(2)}+\Delta\mathbf{G}^{(12)}_{j}\rangle.$
(35b)

where
$\mathbf{q}^{(1)}=0;\mathbf{q}^{(2)}=\Delta\mathbf{K}^{(12)}=\mathbf{K}^{(1)}-\mathbf{K}^{(2)}$
and $\mathbf{q}$ is a continuous vector as a small parameter measured from the
$\mathbf{K}^{(\ell)}$ points. This denotation is another writing of
$|\ell,\alpha,\mathbf{k}^{(\ell)}\rangle$, where
$\mathbf{k}^{(\ell)}=\mathbf{K}^{(\ell)}+\mathbf{q}+\mathbf{q}^{(\ell)}+\Delta\mathbf{G}^{(12)}_{i}$,
but it resembles that of localized orbitals in the real-space representation.
The coupled state vectors in the collected subset, which are defined by the
same parameter vector $\mathbf{q}$, can be sorted in the order:

$\displaystyle\\{$ $\displaystyle\ldots,$
$\displaystyle|\mathbf{K}^{(1)},\alpha,\mathbf{q},\mathbf{q}^{(1)}+\Delta\mathbf{G}^{(12)}_{i}\rangle,$
$\displaystyle\ldots$
$\displaystyle|\mathbf{K}^{(1)},\alpha,\mathbf{q},\mathbf{q}^{(1)}+\Delta\mathbf{G}^{(12)}_{2}\rangle,$
$\displaystyle|\mathbf{K}^{(1)},\alpha,\mathbf{q},\mathbf{q}^{(1)}+\Delta\mathbf{G}^{(12)}_{1}\rangle,$
$\displaystyle|\mathbf{K}^{(2)},\beta,\mathbf{q},\mathbf{q}^{(2)}+\Delta\mathbf{G}^{(12)}_{1}\rangle,$
$\displaystyle|\mathbf{K}^{(2)},\beta,\mathbf{q},\mathbf{q}^{(2)}+\Delta\mathbf{G}^{(12)}_{2}\rangle,$
$\displaystyle\ldots$
$\displaystyle|\mathbf{K}^{(2)},\beta,\mathbf{q},\mathbf{q}^{(2)}+\Delta\mathbf{G}^{(12)}_{j}\rangle,$
$\displaystyle\ldots\\}$ (36)

By this sorting the subset of Bloch vectors can be expanded to both sides. The
expansion should be terminated when the vectors
$\mathbf{q}+\mathbf{q}^{(1)}+\Delta\mathbf{G}^{(12)}_{i}$ and
$\mathbf{q}+\mathbf{q}^{(2)}+\Delta\mathbf{G}^{(12)}_{j}$ do no longer lie in
the domains $\mathscr{O}(\mathbf{K}^{(1)})$ and
$\mathscr{O}(\mathbf{K}^{(2)})$, respectively. A subset of the basis vectors
determined by Eq. (III.4) therefore allows us to build a size-small matrix,
denoted by $H[\mathbf{K}^{(1)},\mathbf{K}^{(2)}]$, which is actually a block
of the total Hamiltonian matrix. The effective Hamiltonian, taking only the
coupling of the Bloch states defined in the vicinity of the $\mathbf{K}^{(1)}$
and $\mathbf{K}^{(2)}$ points, is formally written as the summation of three
terms
$H_{eff}[\mathbf{K}^{(1)},\mathbf{K}^{(2)}]=H[\mathbf{K}^{(1)}]+H[\mathbf{K}^{(2)}]+H[\mathbf{K}^{(1)},\mathbf{K}^{(2)}]$,
in which:

$\displaystyle H[\mathbf{K}^{(\ell)}]$
$\displaystyle=\sum_{\mathbf{q}\in\mathscr{O}(0)}\sum_{i}\sum_{\alpha,\beta}|\mathbf{K}^{(\ell)},\alpha,\mathbf{q},\mathbf{q}^{(\ell)}+\Delta\mathbf{G}^{(12)}_{i}\rangle
H_{\alpha\beta}^{(\ell)}(\mathbf{K}^{(\ell)}+\mathbf{q}+\mathbf{q}^{(\ell)}+\Delta\mathbf{G}^{(12)}_{i})\langle\mathbf{K}^{(\ell)},\beta,\mathbf{q},\mathbf{q}^{(\ell)}+\Delta\mathbf{G}^{(12)}_{i}|,$
(37a) $\displaystyle H[\mathbf{K}^{(1)},\mathbf{K}^{(2)}]$
$\displaystyle=\sum_{\mathbf{q}\in\mathscr{O}(0)}\sum_{i,j}\sum_{\alpha,\beta}|\mathbf{K}^{(1)},\alpha,\mathbf{q},\mathbf{q}^{(1)}+\Delta\mathbf{G}^{(12)}_{i}\rangle\left[\sum_{s}T{(\mathbf{K}^{(1)}+\mathbf{G}^{(1)}_{s})}M_{\alpha\beta}(\mathbf{G}^{(1)}_{s},\mathbf{G}^{(2)}_{s})\right.$
$\displaystyle\hskip
170.71652pt\left.\times\delta_{\Delta\mathbf{G}^{(12)}_{j},\Delta\mathbf{G}^{(12)}_{i}+\Delta\mathbf{G}^{(12)}_{s}}\right]\langle\mathbf{K}^{(2)},\beta,\mathbf{q},\mathbf{q}^{(2)}+\Delta\mathbf{G}^{(12)}_{j}|+$
$\displaystyle+\sum_{\mathbf{q}\in\mathscr{O}(0)}\sum_{i,j}\sum_{\alpha,\beta}|\mathbf{K}^{(2)},\beta,\mathbf{q},\mathbf{q}^{(2)}+\Delta\mathbf{G}^{(12)}_{j}\rangle\left[\sum_{s}T{(\mathbf{K}^{(2)}+\mathbf{G}^{(2)}_{s})}M_{\beta\alpha}(\mathbf{G}^{(2)}_{s},\mathbf{G}^{(1)}_{s})\right.$
$\displaystyle\hskip
170.71652pt\left.\times\delta_{\Delta\mathbf{G}^{(12)}_{j},\Delta\mathbf{G}^{(12)}_{i}+\Delta\mathbf{G}^{(12)}_{s}}\right]\langle\mathbf{K}^{(1)},\alpha,\mathbf{q},\mathbf{q}^{(1)}+\Delta\mathbf{G}^{(12)}_{i}|.$
(37b)

This effective Hamiltonian $H_{eff}[\mathbf{K}^{(1)},\mathbf{K}^{(2)}]$
defines a Schrodinger equation whose solution should be found in terms of the
coupled basis vectors given in Eqs. (35a,35b) as follows:

$|\psi\rangle=\sum_{\mathbf{q}\in\mathscr{O}(0)}\sum_{\ell,\alpha,j}C_{\ell\alpha
j}(\mathbf{q})|\mathbf{K}^{(\ell)},\alpha,\mathbf{q},\mathbf{q}^{(\ell)}+\Delta\mathbf{G}^{(12)}_{j}\rangle,$
(38)

where $C_{\ell\alpha j}(\mathbf{q})$ are the linear combination coefficients
needed to be determined. It is expected that the found solutions could
describe correctly the electron states in an interested energy range in the
whole spectrum of the bilayer systems.

The effective model constructed via the procedure above is general for the
arbitrary alignment between two graphene layers. In the case that only the
lattice sliding/displacement is considered, the translation symmetry of the
single graphene layer is preserved in the bilayer lattice, the Brillouin zone
BZ1 is identical to BZ2 and defines the Brillouin zone of the bilayer lattice.
In this case, the vector $\mathbf{G}^{(1)}_{j}\equiv\mathbf{G}^{(2)}_{j}$ and
therefore $\Delta\mathbf{G}^{(12)}_{j}=0$. From the selection rules we see the
state $|1,\alpha,\mathbf{k}^{(1)}\rangle$ and the one
$|2,\beta,\mathbf{k}^{(2)}\rangle$ always couple together when
$\mathbf{k}^{(1)}=\mathbf{k}^{(2)}=\mathbf{k}\in$ BZ. The set of basis state
vectors collected by the procedure above allows to full-fill the whole basis
set of the Bloch state vectors of two graphene layers. The summation over
$\mathbf{G}^{(1)}$ and $\mathbf{G}^{(2)}$ in Eq. (33) therefore becomes:

$H_{\alpha\beta}(\mathbf{k}^{(1)},\mathbf{k}^{(2)})=\delta_{\mathbf{k}^{(1)},\mathbf{k}^{(2)}}\sum_{\mathbf{G}}T(\mathbf{k}+\mathbf{G})e^{i\mathbf{G}\cdot(\mathbf{d}_{\alpha}^{(1)}-\mathbf{d}_{\beta}^{(2)})}.$
(39)

With the notice to Eq. (30) the summation over $\mathbf{G}$ in the above
equation leads to the one identical to Eq. (28). Accordingly, in the case of
sliding alignment the procedure presented here results in exactly the
expression of the Hamiltonian of two graphene layers.

###  III.5 Bistritzer-MacDonald model

In this subsection we consider the bilayer graphene characterized by a tiny
twist angle and a small displacement vector. We are going to derive a
Hamiltonian to describe effectively the behavior of electrons in an energy
range around the Fermi energy level. It is well known that the low energy
eigen-states of electrons in a single monolayer $\ell$ are defined by the
vector $\mathbf{k}^{(\ell)}$ in the vicinity of the two independent corner
points $\mathbf{K}_{\xi}^{(\ell)}$ ($\xi=\pm 1$) of the hexagonal Brillouin
zone BZℓ. Applying the theory presented in the previous sub-section we expect
to derive a desired model as the construction of the linear combination of
single-layer Bloch states defined in the vicinity of the
$\mathbf{K}^{(\ell)}_{\xi}$ points. In the limit of tiny twist angles, the two
points $\mathbf{K}_{\xi}^{(1)}$ and $\mathbf{K}_{\xi}^{(2)}$ are close to each
other. The two domains $\mathscr{O}(\mathbf{K}^{(1)}_{\xi})$ and
$\mathscr{O}(\mathbf{K}^{(2)}_{\xi})$ therefore can be seen as identical. We
consider only the coupling of the Bloch states defined in these
$\mathbf{k}$-domains. Accordingly, we have to choose $\mathbf{G}^{(2)}_{j}$
close to $\mathbf{G}^{(1)}_{j}$ so that $\Delta\mathbf{G}^{(12)}_{j}$ are
small vectors. Since
$\mathbf{K}_{\xi}^{(2)}=R_{z}(\theta)\cdot\mathbf{K}_{\xi}^{(1)}$,
$\mathbf{G}^{(2)}_{j}$ should be chosen as
$\mathbf{G}^{(2)}_{j}=R_{z}(\theta)\cdot\mathbf{G}^{(1)}_{j}$. For simplicity,
we consider only the vectors $\mathbf{G}^{(1)}$ and $\mathbf{G}^{(2)}$ in the
first class $\mathscr{G}_{1}^{(1)}$ and $\mathscr{G}_{1}^{(2)}$ to construct
Eq. (31b). These vectors determine three vectors $\Delta\mathbf{G}^{(12)}_{j}$

$\displaystyle\Delta\mathbf{G}_{1}^{(12)}$ $\displaystyle=0,$ (40a)
$\displaystyle\Delta\mathbf{G}_{2}^{(12)}$
$\displaystyle=\xi\Delta\mathbf{a}_{1}^{\star(12)},$ (40b)
$\displaystyle\Delta\mathbf{G}_{3}^{(12)}$
$\displaystyle=\xi(\Delta\mathbf{a}_{1}^{\star(12)}+\Delta\mathbf{a}_{2}^{\star(12)}),$
(40c)

that define the coupling of a Bloch state in one graphene layer to three Bloch
states in the other graphene layer. The coupling strength is evaluated by
$T(\mathbf{K}_{\xi}^{(1)}+\mathbf{G}_{j}^{(1)},d_{GG})=T(4\pi/3a,d_{GG})$. We
determine the factors $M_{\alpha\beta}^{j}$ and arrange them into the matrix
form as follows:

$M^{1}=\left(\begin{array}[]{cc}1&1\\\
1&1\end{array}\right),M^{2}=\left(\begin{array}[]{cc}1&\omega\\\
\omega^{*}&1\end{array}\right),M^{3}=\left(\begin{array}[]{cc}1&\omega^{*}\\\
\omega&1\end{array}\right),$ (41)

where $\omega=e^{i2\pi/3}$. Notice that the expression of these matrices
depends on the choice of the coordinate frame. By defining three vectors

$\mathbf{q}_{i}=\Delta\mathbf{K}^{(12)}+\Delta\mathbf{G}^{(12)}_{j},$ (42)

the subset of basis Bloch vectors can be obtained from the list as follows:

$\displaystyle\\{$
$\displaystyle|\mathbf{K}^{(1)}_{\xi},\alpha,\mathbf{q},\mathbf{Q}_{1}=0\rangle,$
$\displaystyle|\mathbf{K}^{(2)}_{\xi},\beta,\mathbf{q},\mathbf{Q}_{2}=\mathbf{Q}_{1}+\mathbf{q}_{1}\rangle,$
$\displaystyle|\mathbf{K}^{(2)}_{\xi},\beta,\mathbf{q},\mathbf{Q}_{3}=\mathbf{Q}_{1}+\mathbf{q}_{2}\rangle,$
$\displaystyle|\mathbf{K}^{(2)}_{\xi},\beta,\mathbf{q},\mathbf{Q}_{4}=\mathbf{Q}_{1}+\mathbf{q}_{3}\rangle,$
$\displaystyle|\mathbf{K}^{(1)}_{\xi},\alpha,\mathbf{q},\mathbf{Q}_{5}=\mathbf{Q}_{2}-\mathbf{q}_{2}\rangle,$
$\displaystyle|\mathbf{K}^{(1)}_{\xi},\alpha,\mathbf{q},\mathbf{Q}_{6}=\mathbf{Q}_{2}-\mathbf{q}_{3}\rangle,$
$\displaystyle|\mathbf{K}^{(1)}_{\xi},\alpha,\mathbf{q},\mathbf{Q}_{7}=\mathbf{Q}_{3}-\mathbf{q}_{1}\rangle,$
$\displaystyle|\mathbf{K}^{(1)}_{\xi},\alpha,\mathbf{q},\mathbf{Q}_{8}=\mathbf{Q}_{3}-\mathbf{q}_{3}\rangle,$
$\displaystyle|\mathbf{K}^{(1)}_{\xi},\alpha,\mathbf{q},\mathbf{Q}_{9}=\mathbf{Q}_{4}-\mathbf{q}_{1}\rangle,$
$\displaystyle|\mathbf{K}^{(1)}_{\xi},\alpha,\mathbf{q},\mathbf{Q}_{10}=\mathbf{Q}_{4}-\mathbf{q}_{2}\rangle,$
$\displaystyle|\mathbf{K}^{(2)}_{\xi},\beta,\mathbf{q},\mathbf{Q}_{11}=\mathbf{Q}_{5}+\mathbf{q}_{1}\rangle,$
$\displaystyle|\mathbf{K}^{(2)}_{\xi},\beta,\mathbf{q},\mathbf{Q}_{12}=\mathbf{Q}_{5}+\mathbf{q}_{3}\rangle=|\mathbf{K}^{(2)}_{\xi},\beta,\mathbf{q},\mathbf{Q}_{10}+\mathbf{q}_{1}\rangle,$
$\displaystyle|\mathbf{K}^{(2)}_{\xi},\beta,\mathbf{q},\mathbf{Q}_{13}=\mathbf{Q}_{6}+\mathbf{q}_{1}\rangle,$
$\displaystyle|\mathbf{K}^{(2)}_{\xi},\beta,\mathbf{q},\mathbf{Q}_{14}=\mathbf{Q}_{6}+\mathbf{q}_{2}\rangle=|\mathbf{K}^{(2)}_{\xi},\beta,\mathbf{q},\mathbf{Q}_{8}+\mathbf{q}_{1}\rangle,$
$\displaystyle|\mathbf{K}^{(2)}_{\xi},\beta,\mathbf{q},\mathbf{Q}_{15}=\mathbf{Q}_{8}+\mathbf{q}_{2}\rangle,$
$\displaystyle|\mathbf{K}^{(2)}_{\xi},\beta,\mathbf{q},\mathbf{Q}_{16}=\mathbf{Q}_{7}+\mathbf{q}_{2}\rangle,$
$\displaystyle|\mathbf{K}^{(2)}_{\xi},\beta,\mathbf{q},\mathbf{Q}_{17}=\mathbf{Q}_{7}+\mathbf{q}_{3}\rangle=|\mathbf{K}^{(2)}_{\xi},\beta,\mathbf{q},\mathbf{Q}_{9}+\mathbf{q}_{2}\rangle,$
$\displaystyle|\mathbf{K}^{(2)}_{\xi},\beta,\mathbf{q},\mathbf{Q}_{18}=\mathbf{Q}_{9}+\mathbf{q}_{3}\rangle,$
$\displaystyle|\mathbf{K}^{(2)}_{\xi},\beta,\mathbf{q},\mathbf{Q}_{19}=\mathbf{Q}_{10}+\mathbf{q}_{3}\rangle$
$\displaystyle\ldots\\}$ (43)

It is remarkable to notice that three vectors $\mathbf{q}_{j}$ point to three
$\mathbf{k}$ points that make the corners of an equivalent triangle. The
vectors $\mathbf{Q}_{j}$ appearing in the above list therefore can be arranged
into an hexagonal lattice with the vectors $\Delta\mathbf{G}^{(12)}_{2}$ and
$\Delta\mathbf{G}^{(12)}_{3}$ as the two lattice vectors, see Fig. 2(b). This
observation was first realized by Bistritzer and MacDonald in Ref. [Bistritzer
and MacDonald, 2011]. This lattice representation is useful to determine all
possible coupled Bloch state vectors. The subset of four first state vectors
was used to show the influence of one graphene layer on the electronic
structure of the other graphene layer. The dependence of the Fermi velocity on
the twist angle was demonstrated analytically.Bistritzer and MacDonald (2011)
The subset of ten first state vectors was also used to show the energy band
structure of TBGs with relatively large twist angles.Zhang (2019)

Collecting a list of coupled Bloch state vectors as given in Eq. (III.5) to
construct a finite size matrix of effective Hamiltonian is rather technical.
There is another elegantly physical approach that is based on the long
wavelength approximation. Accordingly, we approximate the Bloch wave functions
by the plane wave functions, i.e.,
$\langle\mathbf{x}|\ell,\alpha,\mathbf{k}^{(\ell)}\rangle\rightarrow
e^{i\mathbf{k}^{(\ell)}\cdot\mathbf{x}}/\sqrt{S}$. In the state vector
denotation, this approximation is written as the decomposition
$|\ell,\alpha,\mathbf{k}^{(\ell)}\rangle\rightarrow|\ell,\alpha\rangle\otimes|\mathbf{k}^{(\ell)}\rangle$.
By this approximation we can define an operator
$H_{12}^{\xi}(\hat{\mathbf{x}})$ such that the elements of the inter-layer
coupling Hamiltonian
$H_{\alpha\beta}^{\xi}(\mathbf{k}^{(1)},\mathbf{k}^{(2)})=\langle
1,\alpha,\mathbf{k}^{(1)}|H|2,\beta,\mathbf{k}^{(2)}\rangle$ are determined by
$\langle\mathbf{k}^{(1)}|\left[\langle
1,\alpha|H_{12}^{\xi}(\hat{\mathbf{x}})|2,\beta\rangle\right]|\mathbf{k}^{(2)}\rangle$.
Here $\hat{\mathbf{x}}$ is the position operator
$\hat{\mathbf{x}}|\mathbf{x}\rangle=\mathbf{x}|\mathbf{x}\rangle$, and
$\langle\mathbf{x}|\mathbf{k}^{(\ell)}\rangle=e^{i\mathbf{k}^{(\ell)}\cdot\mathbf{x}}/\sqrt{S}$.
Indeed, we find that the matrix elements of this operator

$H_{12}^{\xi}(\hat{\mathbf{x}})=T(\mathbf{K}_{\xi}^{(1)},d_{GG})\sum_{j=1}^{3}M^{j}e^{-i\Delta\mathbf{G}_{j}^{(12)}\cdot\hat{\mathbf{x}}}$
(44)

in the basis of vectors
$\\{|\ell,\alpha\rangle\otimes|\mathbf{k}^{(\ell)}\rangle\\}$ are identical to
the expression

$H^{\xi}_{\alpha\beta}(\mathbf{k}^{(1)},\mathbf{k}^{(2)})=T(\mathbf{K}_{\xi}^{(1)},d_{GG})\sum_{j=1}^{3}M^{j}_{\alpha\beta}\delta_{\mathbf{k}^{(1)},\mathbf{k}^{(2)}-\Delta\mathbf{G}^{(12)}_{j}},$
(45)

This observation for the matrix elements
$H_{\alpha\beta}^{\xi}(\mathbf{k}^{(1)},\mathbf{k}^{(2)})$ is important, but
it is not yet the whole story. The rest half must lie in a similar analysis
for the intra-layer coupling matrix elements
$H_{\alpha\beta}^{(\ell)}(\mathbf{k}^{(\ell)})$. By assuming that
$\mathbf{k}^{(\ell)}$ is close to $\mathbf{K}^{(\ell)}_{\xi}$, we can define a
small $\mathbf{q}^{(\ell)}$ such that
$\mathbf{k}^{(\ell)}=\mathbf{K}^{(\ell)}_{\xi}+\mathbf{q}^{(\ell)}$. Now
linearizing the Hamiltonian matrix elements
$H^{(\ell)}_{\alpha\beta}(\mathbf{K}_{\xi}^{(\ell)}+\mathbf{q}^{(\ell)})$ with
respect to $\mathbf{q}$ we have
$H^{(\ell)}_{\alpha\beta}(\mathbf{K}_{\xi}^{(\ell)}+\mathbf{q}^{(\ell)})\approx
H^{(\ell)}_{\alpha\beta}(\mathbf{K}_{\xi}^{(\ell)})+\nabla
H^{(\ell)}_{\alpha\beta}(\mathbf{K}_{\xi})\cdot\mathbf{q}^{(\ell)}$. This
result can be also written as

$\displaystyle
H^{(\ell)}_{\alpha\beta}(\mathbf{k}^{(\ell)})=H^{(\ell)}_{\alpha\beta}(\mathbf{K}_{\xi}^{(\ell)})+$
$\displaystyle\hskip
28.45274pt\langle\mathbf{k}^{(\ell)}|\frac{1}{\hbar}\nabla
H^{(\ell)}_{\alpha\beta}(\mathbf{K}_{\xi}^{(\ell)})\cdot(\hat{\mathbf{p}}-\hbar\mathbf{K}^{(\ell)}_{\xi})|\mathbf{k}^{(\ell)}\rangle,$
(46)

where $\hat{\mathbf{p}}=-i\hbar\nabla$ is nothing rather than the momentum
operator as convention. It means that, we can also define the operators:

$\displaystyle
H_{\ell}^{\xi}(\hat{\mathbf{p}})=H^{(\ell)}(\mathbf{K}^{(\ell)}_{\xi})+\frac{1}{\hbar}\nabla
H^{(\ell)}(\mathbf{K}_{\xi}^{(\ell)})\cdot(\hat{\mathbf{p}}-\hbar\mathbf{K}_{\xi}^{(\ell)})$
(47)

that work in the Hilbert space spanned by plane-wave vectors. Combining Eqs.
(45) and (47) it allows us to write down an effective continuum Hamiltonian
operator defined through the canonical operators $\hat{\mathbf{p}}$ and
$\hat{\mathbf{x}}$ that enter separately into the kinetic and potential-like
terms, i.e.,
$H(\hat{\mathbf{p}},\hat{\mathbf{x}})=K(\hat{\mathbf{p}})+U(\mathbf{\hat{\mathbf{x}}})$,
wherein

$\displaystyle K(\hat{\mathbf{p}})$
$\displaystyle=\left(\begin{array}[]{cc}H_{1}^{\xi}(\hat{\mathbf{p}})&0\\\
0&H_{2}^{\xi}(\hat{\mathbf{p}})\end{array}\right),$ (50) $\displaystyle
U(\hat{\mathbf{x}})$
$\displaystyle=\left(\begin{array}[]{cc}0&H_{12}^{\xi}(\hat{\mathbf{x}})\\\
H_{12}^{\xi\dagger}(\hat{\mathbf{x}})&0\end{array}\right),$ (53)

where the detailed expression of $H^{\xi}_{1,2}(\hat{\mathbf{p}})$ is given by
Eq. (59). This continuum effective Hamiltonian allows to describe the low
energy states of electrons in the bilayer graphene system as the combination
of single-layer Bloch states defined in the valleys
$\mathbf{K}^{(1,2)}_{\xi}$. This obtained Hamiltonian
$H(\hat{\mathbf{p}},\hat{\mathbf{x}})$ is identical to the one driven by
Bistritzer and MacDonald in Ref. [Bistritzer and MacDonald, 2011]. It is worth
noticing that the derivation of the Bistritzer-MacDonald here does not rely
rigidly on the twist only, but on the genetic alignment of two graphene
lattices. The layer displacement effects are, in fact, easily included in the
model by modifying the elements of the matrix $M^{j}$ with the displacement of
the atom position
$\mathbf{d}_{\beta}^{(\ell)}\rightarrow\mathbf{d}_{\beta}^{(\ell)}+(\ell-1)\boldsymbol{\tau}$.

###  III.6 Continuum model for commensurate TBGs

For commensurate TBG configurations, there is another way to parameterize the
coupled single-layer Bloch state vectors using a vector $\mathbf{q}$ defined
in the vicinity of the $\Gamma$ point rather than the $\mathbf{K}^{(\ell)}$
points. Indeed, because of the translation symmetry, the Bloch states of
electron in the complex atomic lattice should be characterized by a vector
$\mathbf{q}$ restricted in a zone, named mini Brillouin zone (MBZ), which is a
mini-zone locating at the centre of the first Brillouin zone of the two
graphene layers. The low energy states of electrons in the TBG lattice are
found as the linear combination of single-layer Bloch states with the wave
vector $\mathbf{k}^{(1,2)}$ defined around the point
$\mathbf{K}_{\xi}^{(1,2)}$. We show in Fig. 2 that the MBZ of the commensurate
TBG lattice can be expanded to cover the BZ of two graphene layers.
Especially, the corner $\mathbf{K}^{M}_{\xi}$ points of the MBZ can be mapped
exactly to the $\mathbf{K}^{(\ell)}_{\xi}$ points of the BZℓ. So, in order to
construct an effective Hamiltonian to determine the states of electrons that
are defined by the parameter vector $\mathbf{q}\in$ MBZ we need to map each
point $\mathbf{q}$ inside the MBZ to as many as possible points in the
vicinity of the $\mathbf{K}^{(\ell)}_{\xi}$ points. This map is realized
thanks to a finite set of the TBG reciprocal lattice vectors
$\\{\mathbf{G}^{M}_{m}\in\Gamma_{M}^{\star}\\}$. By this way, from Eq. (31b)
by replacing $\mathbf{k}^{(\ell)}$ by $\mathbf{q}+\mathbf{G}^{M}_{m}$ and
$\mathbf{k}^{(\ell^{\prime})}$ by $\mathbf{q}^{\prime}+\mathbf{G}^{M}_{n}$ and
the sum over $\mathbf{k}$ by the sums over $\mathbf{q}\in\text{MBZ}$ and over
a set of $\\{\mathbf{G}^{M}_{m}\in\Gamma^{\star}_{M},m=1,2,\ldots\\}$ we
obtain the matrix elements of the Hamiltonian:

$\displaystyle
H^{\ell\ell^{\prime}}_{\alpha\beta;mn}(\mathbf{q})=\langle\ell,\alpha,\mathbf{q}+\mathbf{G}^{M}_{m}|H|\ell^{\prime},\beta,\mathbf{q}+\mathbf{G}^{M}_{n}\rangle$
$\displaystyle\hskip
5.69046pt=\sum_{s}T(\mathbf{K}^{(\ell)}_{\xi}+\mathbf{G}^{(\ell)}_{s})M_{\alpha\beta}(\mathbf{G}^{(\ell)}_{s},\mathbf{G}^{(\ell^{\prime})}_{s})\delta_{\mathbf{G}^{M}_{n},\mathbf{G}^{M}_{m}-\Delta\mathbf{G}^{(\ell\ell^{\prime})}_{s}}.$
(54)

The eigen-vectors of the obtained effective Hamiltonian matrix are
characterized by the vector $\mathbf{q}\in\text{MBZ}$. They should be found in
the form:

$|\psi(\mathbf{q})\rangle=\sum_{\ell}\sum_{\alpha}\sum_{m}C_{\ell,\alpha,m}(\mathbf{q})|\ell,\alpha,\mathbf{q}+\mathbf{G}^{M}_{m}\rangle.$
(55)

Here we should notice that $\mathbf{q}\in$ MBZ and
$\mathbf{q}+\mathbf{G}^{M}_{m}\in\mathscr{O}(\mathbf{K}_{\xi}^{(\ell)})$. In
the long wavelength approximation, the Bloch state vectors can be
approximately decomposed into
$|\ell,\alpha,\mathbf{q}+\mathbf{G}^{M}_{m}\rangle\rightarrow|\ell,\alpha\rangle\otimes|\mathbf{q}+\mathbf{G}^{M}_{m}\rangle$,
where $|\mathbf{q}+\mathbf{G}^{M}_{m}\rangle$ are the eigen-vectors of the
momentum operator $\hat{\mathbf{p}}$, and hence
$|\mathbf{q}+\mathbf{G}^{M}_{m}\rangle=\exp\left(-i\mathbf{G}^{M}_{m}\cdot\hat{\mathbf{x}}\right)|\mathbf{q}\rangle$.
Eq. (54) is thus converted to the form of Eq. (44). Accordingly, it is worth
understanding that the Bistritzer-MacDonald model is applicable to both
incommensurate and commensurate TBG systems. In the next section we will
present in detail the solution to the Bistritzer-MacDonald model for the
commensurate TBGs. We will show how the plane-wave expansion method can be
employed to solve this model.

Figure 3:  Energy band structures of electrons in three TBG configurations
with $\theta=3.89^{\circ}$ (a), $1.89^{\circ}$ (b), and $1.05^{\circ}$ (c).
The black, green and blue curves are obtained using the effective models based
on the number of coupled Bloch states as the basis set of $N_{\mathbf{Q}}=4$,
10 and 19, respectively. The red curves are obtained from the model using the
second scheme of parameterization with $\mathbf{q}\in$ MBZ. The data are
presented for the variation of the $\mathbf{q}$ vector in the corresponding
$\mathbf{K}^{(1)}_{\xi}$ valley with $\xi=+1$. Figure 4:  (a) Illustration for
anisotropy of the energy band structure of the TBG system. The red-solid and
blue-dashed curves are resulted for $\xi=+1$ and $\xi=-1$, respectively. (b)
The total density of states.

##  IV Electronic band structure

###  IV.1 Solution of the Bistritzer-MacDonald model

Developing effective continuum models for low-energy states of electrons in
TBG systems was proceeded since 2007 by Lopes dos Santos et al. However, the
model driven by Bistritzer and MacDonald in 2011 was well known and commonly
used. In the previous subsections, we present the derivation of this model,
and explain the key points in the two approaches. In this section, we now
present our solution to the Bistritzer-MacDonald model. Accordingly, the low
energy states of electrons in the TBG configurations of tiny twist angles are
distinguished by the quantum indices $\xi=\pm 1$ that correspond to the two
nonequivalent Dirac valleys $\mathbf{K}_{\xi}^{(1,2)}$ of the graphene mono-
layers.Bistritzer and MacDonald (2011) Notice that as the tiny twist angles,
the position of the two points $\mathbf{K}_{\xi}^{(1)}$ and
$\mathbf{K}_{\xi}^{(2)}$ are close to each other. The Bistritzer-MacDonald
model is given by a Hamiltonian in the real-space presentation as follows:

$\displaystyle\hat{H}^{\xi}=\left(\begin{array}[]{cc}H_{1}^{\xi}(\hat{\mathbf{p}})&T^{\xi}(\hat{\mathbf{x}})\\\
T^{\xi\dagger}(\hat{\mathbf{x}})&H_{2}^{\xi}(\hat{\mathbf{p}})\end{array}\right)$
(58)

Here $\hat{\mathbf{p}}$ and $\hat{\mathbf{x}}$ are the momentum operator and
the position operator, respectively. The Hamiltonian
$H_{1}^{\xi}(\hat{\mathbf{p}})$ and $H_{2}^{\xi}(\hat{\mathbf{p}})$ of the two
uncoupled graphene layers are given by the 2D Dirac Hamiltonian:

$H_{\ell}^{\xi}(\hat{\mathbf{p}})=-v_{F}(\xi\sigma_{x},\sigma_{y})\cdot\left[R^{z}\left(-\theta_{\ell}\right)\cdot(\hat{\mathbf{p}}-\hbar\mathbf{K}_{\xi}^{(\ell)})\right],$
(59)

where $\theta_{\ell}=(-1)^{\ell}\theta/2$ is the rotation angle of layer
$\ell$ with respect to the fix coordinate axes $Oxy$;
$R^{z}\left(-\theta_{\ell}\right)$ is a matrix to rotate back the relevant
vectors to the $Oxy$ axes to keep the canonical form of the Dirac Hamiltonian;
$\mathbf{K}_{\xi}^{(\ell)}$ is the corner point of type (valley) $\xi$ of the
first Brillouin zone of the layer $\ell$; $\sigma_{x},\sigma_{y}$ are two
conventional Pauli matrices; and $v_{F}$ is the Fermi velocity (The minus sign
is due to the negative value of the hopping parameter, $V_{pp\pi}=-2.7$ eV).
The interlayer coupling block term is given by:Koshino  _et al._ (2018)

$\displaystyle T^{\xi}(\hat{\mathbf{x}})=$
$\displaystyle\left(\begin{array}[]{cc}u&u^{\prime}\\\
u^{\prime}&u\end{array}\right)e^{-i\Delta\mathbf{G}^{(12)}_{\xi
1}\cdot\hat{\mathbf{x}}}+$ (62)
$\displaystyle\left(\begin{array}[]{cc}u&u^{\prime}\omega^{-\xi}\\\
u^{\prime}\omega^{\xi}&u\end{array}\right)e^{-i\Delta\mathbf{G}^{(12)}_{\xi
2}\cdot\hat{\mathbf{x}}}+$ (65)
$\displaystyle\left(\begin{array}[]{cc}u&u^{\prime}\omega^{\xi}\\\
u^{\prime}\omega^{-\xi}&u\end{array}\right)e^{-i\Delta\mathbf{G}^{(12)}_{\xi
3}\cdot\hat{\mathbf{x}}},$ (68)

where $\omega=\exp(i2\pi/3)$ and $\Delta\mathbf{G}^{(12)}_{\xi
1}=0,\Delta\mathbf{G}^{(12)}_{\xi
2}=\xi(\mathbf{a}_{1}^{\star(1)}-\mathbf{a}_{1}^{\star(2)})$ and
$\Delta\mathbf{G}^{(12)}_{\xi
3}=\xi(\mathbf{a}_{1}^{\star(1)}-\mathbf{a}_{1}^{\star(2)})+\xi(\mathbf{a}_{2}^{\star(1)}-\mathbf{a}_{2}^{\star(2)})$.
For commensurate TBG lattices we determine further that
$\Delta\mathbf{G}^{(12)}_{\xi}=0,\Delta\mathbf{G}^{(12)}_{\xi
2}=\mathbf{A}_{1}^{\star}$ and $\Delta\mathbf{G}^{(12)}_{\xi
3}=\xi(\mathbf{A}_{1}^{\star}+\mathbf{A}_{2}^{\star})$, where
$\mathbf{A}_{1,2}^{\star}$ are the basis vectors of the TBG reciprocal lattice
(see Sec. I). Here we use the values proposed by Koshino for the two
parameters $u=0.0797$ eV and $u^{\prime}=0.0975$ eV. The difference between
these two parameters was discussed to capture effects of the lattice
relaxation.Koshino and Nam (2020)

We find the spectrum of the Hamiltonian $\hat{H}$ by solving this secular
equation:

$\displaystyle\left(\begin{array}[]{cc}H_{1}^{\xi}(\hat{\mathbf{p}})&T^{\xi}(\hat{\mathbf{x}})\\\
T^{\xi\dagger}(\hat{\mathbf{x}})&H_{2}^{\xi}(\hat{\mathbf{p}})\end{array}\right)|\psi_{\mathbf{\xi,k}}\rangle=E|\psi_{\xi,\mathbf{k}}\rangle.$
(71)

Due to the periodicity of the moire lattice of TBG lattices, electron states
are determined as the Bloch state vectors $|\psi_{\xi,\mathbf{k}}\rangle$,
where $\mathbf{k}\in$ MBZ. These state vectors are expanded in terms of plane-
wave vectors as follows:

$|\psi_{\xi,\mathbf{k}}\rangle=\sum_{m}\left(\begin{array}[]{c}C_{1,\xi,\mathbf{k}}(\mathbf{G}_{m}^{M})\\\
C_{2,\xi,\mathbf{k}}(\mathbf{G}_{m}^{M})\end{array}\right)|\mathbf{k}+\mathbf{G}_{m}^{M}\rangle$
(72)

where
$\hat{\mathbf{p}}|\mathbf{k}+\mathbf{G}_{m}^{M}\rangle=(\mathbf{k}+\mathbf{G}_{n}^{M})|\mathbf{k}+\mathbf{G}_{n}^{M}\rangle$,
$\langle\mathbf{x}|\mathbf{k}+\mathbf{G}_{m}^{M}\rangle=e^{i(\mathbf{k}+\mathbf{G}_{m}^{M})\cdot\mathbf{x}}$
and $C_{1,\xi,\mathbf{k}}(\mathbf{G}_{m}^{M})$ and
$C_{2,\xi,\mathbf{k}}(\mathbf{G}_{m}^{M})$ are the 2D vectors of combination
coefficients that are needed to be found. Here
$\\{\mathbf{G}_{n}^{M}\,|\,n=1,\ldots,N_{\mathbf{G}}\\}$ is a set of
$N_{\mathbf{G}}$ vectors of the moire reciprocal lattice. Substituting this
trial expression into Eq. (71) and left-multiplying both sides with
$\langle\mathbf{k}+\mathbf{G}_{n}^{M}|$ we have:

$\displaystyle\sum_{m}\langle\mathbf{k}+\mathbf{G}_{n}^{M}|\hat{H}^{\xi}|\mathbf{k}+\mathbf{G}_{m}^{M}\rangle\left(\begin{array}[]{c}C_{1,\xi,\mathbf{k}}(\mathbf{G}_{m}^{M})\\\
C_{2,\xi,\mathbf{k}}(\mathbf{G}_{m}^{M})\end{array}\right)$ (75)
$\displaystyle\hskip
93.89418pt=E\left(\begin{array}[]{c}C_{1,\xi,\mathbf{k}}(\mathbf{G}_{n}^{M})\\\
C_{2,\xi,\mathbf{k}}(\mathbf{G}_{n}^{M})\end{array}\right)$ (78)

With the notice that

$\langle\mathbf{k}+\mathbf{G}_{n}^{M}|H_{\ell}^{\xi}|\mathbf{k}+\mathbf{G}_{m}^{M}\rangle=H_{\ell}^{\xi}(\mathbf{k}+\mathbf{G}_{n}^{M})\delta_{\mathbf{G}_{n}^{M},\mathbf{G}_{m}^{M}}$
(79)

and

$\displaystyle\langle\mathbf{k}+\mathbf{G}_{n}^{M}|T^{\xi}(\hat{\mathbf{x}})|\mathbf{k}+\mathbf{G}_{m}^{M}\rangle$
$\displaystyle=\sum_{j=1}^{3}T_{j}^{\xi}\delta_{\mathbf{G}_{n}^{M},\mathbf{G}_{m}^{M}-\Delta\mathbf{G}^{(12)}_{\xi
j}},$ (80)
$\displaystyle\langle\mathbf{k}+\mathbf{G}_{n}^{M}|T^{\xi\dagger}(\hat{\mathbf{x}})|\mathbf{k}+\mathbf{G}_{m}^{M}\rangle$
$\displaystyle=\sum_{j=1}^{3}T_{j}^{\xi\dagger}\delta_{\mathbf{G}_{n}^{M},\mathbf{G}_{m}^{M}+\Delta\mathbf{G}^{(12)}_{\xi
j}},$ (81)

we specify Eq. (75) in the following form:

$\displaystyle\sum_{m}\left[\left(\begin{array}[]{cc}H_{1}^{\xi}(\mathbf{k}+\mathbf{G}_{n}^{M})&T_{1}^{\xi}\\\
T_{1}^{\xi\dagger}&H_{2}^{\xi}(\mathbf{k}+\mathbf{G}_{n}^{M})\end{array}\right)\delta_{\mathbf{G}_{m}^{M},\mathbf{G}_{n}^{M}}+\right.$
(84)
$\displaystyle+\left.\sum_{j=2}^{3}\left(\begin{array}[]{cc}0&T_{j}^{\xi}\\\
0&0\end{array}\right)\delta_{\mathbf{G}_{m}^{M},\mathbf{G}_{n}^{M}+\Delta\mathbf{G}^{(12)}_{j}}+\sum_{j=2}^{3}\left(\begin{array}[]{cc}0&0\\\
T_{j}^{\xi\dagger}&0\end{array}\right)\delta_{\mathbf{G}_{m}^{M},\mathbf{G}_{n}^{M}-\Delta\mathbf{G}^{(12)}_{j}}\right]\left(\begin{array}[]{c}C_{1,\xi,\mathbf{k}}(\mathbf{G}_{m}^{M})\\\
C_{2,\xi,\mathbf{k}}(\mathbf{G}_{m}^{M})\end{array}\right)=E\left(\begin{array}[]{c}C_{1,\xi,\mathbf{k}}(\mathbf{G}_{n}^{M})\\\
C_{2,\xi,\mathbf{k}}(\mathbf{G}_{n}^{M})\end{array}\right).$ (93)

With a set of reciprocal lattice vectors
$\\{\mathbf{G}_{n}^{M}\,|\,n=1,\ldots,N_{\mathbf{G}}\\}$, the above equation
is the representative of a set of $N_{G}$ linear equations for the
coefficients
$(C_{1,\xi,\mathbf{k}}(\mathbf{G}_{n}^{M}),C_{2,\xi,\mathbf{k}}(\mathbf{G}_{n}^{M}))^{T}$.
Numerically, for each value of $\xi$ and each value of
$\mathbf{k}\in\text{MBZ}$ we define a $4N_{G}\times 4N_{G}$ Hermitian matrix
$H_{\mathbf{k}}^{\xi}$ that is given in the blocks as follows:

$\displaystyle[H_{\mathbf{k}}^{\xi}]_{n,n}=\left(\begin{array}[]{cc}H_{1}^{\xi}(\mathbf{k}+\mathbf{G}_{n}^{M})&T_{1}^{\xi}\\\
T_{1}^{\xi\dagger}&H_{2}^{\xi}(\mathbf{k}+\mathbf{G}_{n}^{M})\end{array}\right),$
(96)
$\displaystyle[H_{\mathbf{k}}^{\xi}]_{n,m_{j}}=\left(\begin{array}[]{cc}0&T_{j}^{\xi}\\\
0&0\end{array}\right)\hskip 14.22636pt\text{if}\hskip
14.22636pt\mathbf{G}^{M}_{m_{j}}=\mathbf{G}^{M}_{n}+\Delta\mathbf{G}^{(12)}_{\xi
j},$ (99)
$\displaystyle[H_{\mathbf{k}}^{\xi}]_{n,m_{j}^{\prime}}=\left(\begin{array}[]{cc}0&0\\\
T_{j}^{\xi\dagger}&0\end{array}\right)\hskip 11.38092pt\text{if}\hskip
14.22636pt\mathbf{G}^{M}_{m_{j}^{\prime}}=\mathbf{G}^{M}_{n}-\Delta\mathbf{G}^{(12)}_{\xi
j}.$ (102)

Diagonalize the matrix $H_{\mathbf{k}}^{\xi}$ we obtain all possible eigen-
values of $E_{n}^{\xi}(\mathbf{k})$. From these data we can display the
electronic energy band structure. The issue here is the value of
$N_{\mathbf{G}}$. Since the model is valid for low energy range in which the
energy surfaces of the monolayer graphene take the cone geometry. We thus
define a cutoff energy $E_{c}$ and then determine $N_{\mathbf{G}}$ the number
of $\mathbf{G}^{M}_{n}$ vectors such that $\|\mathbf{G}^{M}_{n}\|\leq
E_{c}/\hbar v_{F}$.

###  IV.2 Numerical results

We employed the numerical method to solve Eq. (84) and similar equations.
Concretely, for the parameterization of the coupled Bloch state vectors using
the small vector $\mathbf{q}$ measured from th $\mathbf{K^{(\ell)}_{\xi}}$
points we used the approximations with the number of $\mathbf{Q}$-points
$N_{\mathbf{Q}}=4,10$ and 19 to build the matrix for the effective
Hamiltonian. The purpose of this is to evaluate the ability of describing
electronic states of the TBG lattices from the simple to complicated levels of
effective models. For commensurate TBG configurations we used the second
parameterization scheme using the vector $\mathbf{q}$ defined in the mini
Brillouin zone of the reciprocal lattice. Since the number of basis state
vectors can be easily extended in the latter scheme, it allows to describe
precisely all the energy dispersion curves in the energy range containing the
Fermi energy level, i.e., ranging from -3 eV up to 3 eV. We hence use this
scheme as a benchmark to validate effective models using the first scheme of
parameterization of coupled Bloch states.

When building Hamiltonian matrices we distinguish the contribution of the
Bloch states defined in the nonequivalent $\mathbf{K}$ valleys with $\xi=\pm
1$ as denoted in Fig. 2. Technically, for the first scheme of parameterization
the Hamiltonian matrices have the size of $2N_{\mathbf{Q}}\times
2N_{\mathbf{Q}}$ with $N_{\mathbf{Q}}=4$, 10 and 19, corresponding to three
simple truncations of the list of coupled Bloch vectors. Meanwhile, it is
$4N_{\mathbf{G^{M}}}\times 4N_{\mathbf{G^{M}}}$ for the second scheme of
parameterization, where $N_{\mathbf{G^{M}}}$ is the number of the reciprocal
lattice vectors of the TBG lattice to map the mini Brillouin zone MBZ to the
domains $\mathscr{O}(\mathbf{K}^{(\ell)})$, see Fig. 2c. The Hamiltonian
matrices are Hermitian and depend on the parameter vector $\mathbf{q}$. They
are numerically diagonalized for each value of $\mathbf{q}$ to display the
energy band structure. We investigated the band structure of the TBG
configurations with the twist angle varying in a large range, from
$1.05^{\circ}$ to $30^{\circ}$, but show in Fig. 3 data for three
configurations with $\theta=3.89^{\circ}$ ($m=9$), $1.89^{\circ}$ ($m=18$) and
$1.05^{\circ}$ ($m=32$). Our obtained results are in agreement with other
available data in literature.Morell  _et al._ (2010); de Laissardiere  _et
al._ (2012); Moon and Koshino (2013); Koshino  _et al._ (2018)

We realize that, for the TBG configurations with the twist angles not too
small, $\theta>2.5^{\circ}$, the model with $N_{\mathbf{Q}}=4$ describes well
two energy bands closest to the Fermi level, see the black curves in Fig. 3.
However, it does not determine fully the number of bands in the energy range
of $(-1,1)$ eV because there are only 8 Bloch state vectors used to represent
the electronic states of the complex system. By extending $N_{\mathbf{Q}}$ to
10, the obtained model allows to produce correctly the number of energy bands
in a narrow energy range about the Fermi level, see the green curves in Fig.
3. However, there are green curves that are not identical to the red curves,
which are obtained by using the second parameterization scheme. Extending
$N_{\mathbf{Q}}$ to 19 allows to improve quantitatively the blue energy
dispersion curves to coincide with the red curves. However, deviations are
still observed for the configurations with small twist angles, see Figs. 3(a)
and 3(b), in the higher energy ranges, for instance, $E>0.6$ eV and $E<-0.6$
eV for $\theta=3.89^{\circ}$, and $E>0.3$ eV and $E<-0.3$ eV for
$\theta=1.89^{\circ}$. For the special TBG configuration with
$\theta=1.05^{\circ}$ we clearly realize the failure of the effective model
with $N_{\mathbf{Q}}=4$ in describing the “flat bands” about the Fermi level.
It is surprised that the model with $N_{\mathbf{Q}}=10$ reproduces well these
special bands, but the model with $N_{\mathbf{Q}}=19$ as expected does not.
Notice that in the two models with $N_{\mathbf{Q}}=4$ and 19 there is an
unbalance between the number of Bloch states in two graphene layers
contributing to the electronic states of the complex system, but it is not the
case in the model with $N_{\mathbf{Q}}=10$. It therefore suggests that in
order to describe well the electronic structure of the TBG system, the
effective models should be constructed on the basis set of state vectors with
the balance of the number of states in each graphene layer contributing to the
real system.

We sum up our investigation as follows: The model with $N_{\mathbf{Q}}=4$ is
simple and easy to be constructed, but it allows to describe rather well the
two energy bands closest to the Fermi level. This model was actually used to
show analytically the dependence of the Fermi velocity on the twist angle.dos
Santos  _et al._ (2007); Bistritzer and MacDonald (2011) Because the $8\times
8$ matrix can be arranged into a $4\times 4$ block matrix whose diagonal
blocks can be approximated by the 2D Dirac model, this model can be
transformed into the real-space representation. It is therefore suitable for
the investigation of properties involving spatial effects of the TBG system.
Extending the number of coupled Bloch states to represent the electron states
in the complex systems in general allows to improve quantitatively the
electronic band structure in a large energy range about the Fermi level. The
models with $N_{\mathbf{Q}}=10$ and $19$ involve in small size matrices (of 20
and 38 dimensions, respectively), so they are really efficient in the
numerical calculation viewpoint. These models are useful to describe the
electronic structure and optical properties of incommensurate TBG
configurations. For commensurate TBG configurations, the second scheme of
parameterization is really useful. It leads to a compact effective model that
can be formally separated into a kinetic and a potential part. The so-called
Bistritzer-MacDonald Hamiltonian takes this beauty, validating in the long
wavelength approximation. The model is thus suitable to describe efficiently
not only optical transition processes but also transport properties of
electrons taking place in the bilayer lattices.

##  V Conclusion

We present and discuss in detail practical techniques in building effective
models to describe the dynamics of electrons in certain energy ranges of the
whole electronic spectrum of the generic bilayer graphene lattices. The
electronic states of the bilayer system are determined as the linear
combinations of single-layer Bloch states. We symmetrize the expression of the
interlayer coupling Hamiltonian to figure out a selection rule that allows to
determine coupled Bloch states. Concretely, two Bloch states can couple to
each other if the difference between two wave vectors must be equal to the
difference of two reciprocal lattice vectors of the two layers. On the basis
of this selection rule, we present a procedure to collect a subset of coupled
Bloch states. This subset of states allows to isolate a block from the total
Hamiltonian matrix to describe the dynamics of electrons in a certain narrow
energy range. In the approximation of long wavelength, we show that when the
Bloch functions in the collected subset are replaced by the plane wave
functions, an effective continuum model can be established. Interestingly,
this model is explicitly defined through the momentum operator and the
position operator, which do not mix together, but enter into two independent
terms. That allows to write the effective Hamiltonian as the sum of the
kinetic- and potential-like terms. Applying the established procedure to the
twisted bilayer graphene of tiny twist angles, we recover the model
established by Bistritzer-MacDonald in 2011 by considering only the strongest
coupling of Bloch states defined in the $\mathbf{K}$-valleys of two graphene
layers. We present in detail a plane wave expansion solution to this model and
numerical results for some commensurate TBG configurations. The obtained data
for the electronic structure are in good agreement with those resulted from
the exhausted tight-binding and DFT calculations. The practical rules and
technical discussions presented in this work are expected to provide useful
methodological knowledge as the background to exploit the effective continuum
Bistritzer-MacDonald model for studying various physical aspects of the
bilayer graphene system as well as to develop for other multiple layer van der
Waals material systems.

## References

  * Sierra  _et al._ (2021) J. F. Sierra, J. Fabian, R. K. Kawakami, S. Roche, and S. O. Valenzuela, Nat. Nanotechnol. 16 (2021), doi.org/10.1038/s41565-021-00936-x.
  * Liaoa  _et al._ (2019) W. Liaoa, Y. Huang, H. Wang, and H. Zhang, Applied Materials Today 16, 435 (2019).
  * Sangwan and Hersam (2018) V. K. Sangwan and M. C. Hersam, Annu. Rev. Phys. Chem 69, 12 (2018).
  * Bao and Hoh (2019) Q. Bao and H. Hoh,  _2D Materials for Photonic and Optoelectronic Applications_, Woodhead Publishing Series in Electronic and Optical Materials (Elsevier Science, 2019).
  * Rozhkov  _et al._ (2016) A. V. Rozhkov, A.O.Sboychakov, A.L.Rakhmanov, and F. Nori, Phys. Rep. 648, 1 (2016).
  * Carr  _et al._ (2017) S. Carr, D. Massatt, S. Fang, P. Cazeaux, M. Luskin, , and E. Kaxiras, Phys. Rev. B 95, 075420 (2017).
  * Uchida  _et al._ (2014) K. Uchida, S. Furyya, J. I. Iwata, and A. Oshiyama, Phys. Rev. B 90, 155451 (2014).
  * Lucignano  _et al._ (2019) P. Lucignano, D. Alfe, V. Cataudella, D. Ninno, and G. Cantele, Phys. Rev. B 99, 195419 (2019).
  * de Laissardiere  _et al._ (2012) G. T. de Laissardiere, D. Mayou, and L. Magaud, Phys. Rev. B 86, 125413 (2012).
  * Morell  _et al._ (2010) E. S. Morell, J. D. Correa, P. Vargas, M. Pacheco, and Z. Barticevic, Phys. Rev. B 82, 121407 (2010).
  * dos Santos  _et al._ (2007) J. M. B. L. dos Santos, N. M. R. Peres, and A. H. C. Neto, Phys. Rev. Lett. 99, 256802 (2007).
  * Shallcross  _et al._ (2008) S. Shallcross, S. Sharma, and O. A. Pankratov, Phys. Rev. Lett. 101, 056803 (2008).
  * Bistritzer and MacDonald (2011) R. Bistritzer and A. H. MacDonald, PNAS 108, 12233 (2011).
  * Shallcross  _et al._ (2013) S. Shallcross, S. Sharma, , and O. Pankratov, Phys. Rev. B 87, 245403 (2013).
  * dos Santos  _et al._ (2012) J. M. B. L. dos Santos, N. M. R. Peres, and A. H. C. Neto, Phys. Rev. B 86, 155449 (2012).
  * Zhang (2019) M. Zhang, ed.,  _Handbook of graphene_, Volume 3: Graphene-like 2D materials (John Wiley & Sons, 2019) p. 177.
  * Mele (2010) E. J. Mele, Phys. Rev. B 81, 161405 (2010).
  * Mele (2011) E. J. Mele, Phys. Rev. B 84, 235439 (2011).
  * de Gail  _et al._ (2011) R. de Gail, M. O. Goerbig, F. Guinea, G. Montambaux, and A. H. C. Neto, Phys. Rev. B 84, 045436 (2011).
  * Weckbecker  _et al._ (2016) D. Weckbecker, S. Shallcross, M. Fleischmann, N. Ray, S. Sharma, and O. Pankratov, Phys. Rev. B 93, 035452 (2016).
  * Tritsaris  _et al._ (2020) G. A. Tritsaris, S. Carr, Z. Zhu, Y. Xie, S. B. Torrisi, J. Tang, M. Mattheakis, D. T. Larson, and E. Kaxiras, 2D Materials 7, 035028 (2020).
  * Koshino and Nam (2020) M. Koshino and N. N. T. Nam, Phys. Rev. B 101, 195425 (2020).
  * Moon and Koshino (2013) P. Moon and M. Koshino, Phys. Rev. B 87, 205404 (2013).
  * Koshino (2015) M. Koshino, New. J. Phys. 17, 015014 (2015).
  * Le and Do (2018) H. A. Le and V. N. Do, Physical Review B 97, 125136 (2018).
  * Le  _et al._ (2019) H. A. Le, V. T. Nguyen, V. D. Nguyen, V. N. Do, and S. T. Ho, Commun. Phys. 29, 455 (2019).
  * Do  _et al._ (2019) V. N. Do, H. A. Le, and D. Bercioux, Phys. Rev. B 99, 165127 (2019).
  * Koshino  _et al._ (2018) M. Koshino, N. F. Q. Yuan, T. Koretsune, M. Ochi, K. Kuroki, and L. Fu, Phys. Rev. X 8, 031087 (2018).