# Optical activity and transport in twisted bilayer graphene: the essence of
spatial dispersion effects

S. Ta Ho  Hanoi University of Civil Engineering (HUCE), 55 Giai Phong road,
Hanoi 10000, Vietnam  V. Nam Do<EMAIL_ADDRESS>Department of
Basic Science, Phenikaa Institute for Advanced Studies (PIAS), 25th Floor, A9
Building, Phenikaa University, Hanoi 10000, Vietnam

###### Abstract

This study investigates optical activity and quantum transport in twisted
bilayer graphene (TBG) systems, demonstrating that the former results from
spatial dispersion effects. The transfer matrix method is used to solve the
propagation of electromagnetic waves through two graphene layers that act as
the coupling surfaces of a dielectric slab. The resulting optical conductivity
tensor is decomposed into a local and a drag part, with the drag transverse
conductivity $\sigma_{xy}^{(drag)}$ governing the TBG system’s optical
activity. An effective continuum model is employed to analyze electron state
formation and calculate relevant parts of the optical conductivity tensor.
Correlation of electron motions leads to incomplete cancellation and a finite
$\sigma_{xy}^{(drag)}$ in the chiral TBG lattice. The study also calculates DC
conductivity, showing TBG supports quantum conductivity proportional to
$e^{2}/h$ at the intrinsic Fermi energy.

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

##  I Introduction

The study of two-dimensional (2D) materials has received significant attention
in recent years due to its potential applications in electronics,
optoelectronics, spintronics, and valleytronics.Sierra  _et al._ (2021); Liaoa
_et al._ (2019); Sangwan and Hersam (2018); Bao and Hoh (2019) The use of
multiple-layer materials with van der Waals (vdW) interlayer coupling offers
various means to tune the electronic structure, such as altering the
interlayer distance, changing the lattice alignment through twisting or
sliding the layers, and applying external fields. Physically, these material
engineering solutions often result in breaking the spatial and temporal
symmetries that govern the dynamics of electrons within the material layers.
For instance, a magnetic field breaks the time-reversal symmetry, leading to
exciting phenomena such as the quantum Hall effect, magneto-optical effects,
and the Kerr and Faraday rotations of the polarization plane of linearly
polarized light and circular dichroism.Tse and MacDonald (2011); Crassee  _et
al._ (2011); Nandkishore and Levitov (2011); Shimano  _et al._ (2013) These
intriguing optical phenomena have been utilized in the development of various
devices for various applications.Izake (2007); Liu  _et al._ (2005); Solomon
_et al._ (2019) However, the use of external fields has limitations for small-
scale device design. As a result, it is crucial to search for materials with
intrinsic properties and solutions that break the spatial symmetries of an
electronic system to meet technological demands.

Twisted-bilayer graphene (TBG) is a 2D van der Waals (2D vdW) material that
has garnered significant attention from the research community due to its
interesting electronic properties.Andrei and MacDonald (2020); Nimbalkar and
Kim (2020); Hill (2020) The TBG system is formed by stacking two graphene
layers on top of each other with a rotation, or twist, angle between the two
hexagonal lattices. The TBG system has been shown to be optically active,
meaning that it has the ability to rotate the polarization plane of linearly
polarized light by reflection and transmission and to absorb left- and right-
handed light differently.Kim  _et al._ (2016) This, combined with the ability
to mechanically tune the twist angle,Ribeiro-Palau  _et al._ (2018); Cai and
Yu (2021) makes TBG a potential candidate for advanced optical applications.
The optical activity of a dielectric medium is typically attributed to the
anisotropy of its atomic lattice,Landau and Lifshitz (1984) or to the magneto-
electric effect.Condon (1937) External electric, magnetic, and mechanical
forces can also induce anisotropy in a medium.Kim  _et al._ (2016);
Nandkishore and Levitov (2011); Stauber  _et al._ (2018a) It has been shown
that strained monolayer graphene is optically active due to the presence of
the off-diagonal element $\sigma_{xy}(\omega)$ of the optical conductivity
tensor.Oliva-Leyva and Naumis (2016) The AA- and AB-stacked configurations of
the TBG system are isotropic and characterized by the point groups $D_{6h}$
and $D_{3d}$, respectively.Zou  _et al._ (2019); Do  _et al._ (2020) The
breaking of the lattice symmetry from the $D_{6h}$ and $D_{3d}$ groups to the
$D_{6}$ and $D_{3}$ groups, respectively, has been suggested as the mechanism
causing the optical activity in TBG.Morell  _et al._ (2017) However, even with
this symmetry, the TBG system is characterized by a total diagonal optical
conductivity tensor, i.e., with zero off-diagonal components
$\sigma_{xy}(\omega)=0$. So, the mechanism supporting the optical activity of
TBG remains to be further explored.

In this paper, we analyze the optical activity of twisted bilayer graphene by
treating simultaneously two essential aspects of the field-matter interaction
problem, namely the propagation of the electromagnetic waves through a
material layer and the response of a material layer to the action of an
external field. We show that the optical activity of the TBG system is a
result of spatial dispersion effects. Our theory is based on the following
analysis. Working with 2D van der Waals systems, the physical quantities that
describe them, such as the dielectric tensor and the conductivity tensor, are
often considered as if they were ideal 2D systems, i.e., with zero thickness,
even though the interlayer spacing is taken into account in the description of
interlayer coupling.Moon and Koshino (2013); Tabert and Nicol (2013); Stauber
_et al._ (2013) However, in reality, even for graphene with only one layer of
carbon atoms, an effective thickness can be defined based on the width of the
quantum well potential caused by the carbon atoms.Do  _et al._ (2017) As a
result, 2D vdW systems cannot be considered ideal 2D systems but rather
quasi-2D (Q2D) systems. This means that when describing the interaction
between an electromagnetic field and a Q2D vdW material system at a
macroscopic level, special considerations must be made. Typically, field
equations are determined by averaging a set of microscopic field equations
over infinitesimal volume elements.Landau and Lifshitz (1984) However, this
averaging procedure is only valid when the wavelength of the fields is much
smaller than the dimensionalities of the physical system. For Q2D materials,
averaging cannot be performed along the $Oz$ direction perpendicular to the
$Oxy$ plane of the atomic lattice. As a result, we need to use a hybrid basis
of reciprocal-real spaces to represent physical quantities and field
equations. For example, the equation

$J_{\alpha}(z,\mathbf{q}_{\parallel},\omega)=\int
dz^{\prime}\sigma_{\alpha\beta}(z,z^{\prime},\mathbf{q}_{\parallel},\omega)E_{\beta}(z^{\prime},\mathbf{q}_{\parallel},\omega)$
(1)

describes the relationship between the component
$J_{\alpha}(z,\mathbf{q}_{\parallel},\omega)$ of an electrical current density
induced in the material by an electric field
$E_{\beta}(z,\mathbf{q}_{\parallel},\omega)$, where $\alpha,\beta=x,y$ are the
indices of the components of vector fields, $\mathbf{q}_{\parallel}$ is the
wave vector in the material plane, $\omega$ is the frequency, and
$\sigma_{\alpha\beta}(z,z^{\prime},\mathbf{q}_{\parallel},\omega)$ is a
component of the electrical conductivity tensor of the material layer. In the
optical limit, where $\|\mathbf{q}_{\parallel}\|a\ll 1$ and $a$ is the lattice
constant, we can approximate all quantities at $\mathbf{q}_{\parallel}=0$,
i.e., effectively ignoring the spatial effects in the material plane but not
along the thickness.Nazarov (2015); Agranovich and Ginzburg (2013);
Poshakinskiy  _et al._ (2018) This description highlights the importance of
taking into account the spatial dispersion along the $Oz$ direction, which is
typically the transmission direction of electromagnetic fields used to study
the material’s properties.

Theoretical studies of optical activity in TBG have been discussed in several
references, including Kim  _et al._ , 2016; Morell  _et al._ , 2017; Stauber
_et al._ , 2018a. These studies have revealed the role of the spacing $a$
between two graphene layers and explained optical activity as a magneto-
electric effect. In this approach, a magnetization $\mathbf{m}_{\parallel}$
and a polarization $\mathbf{p}_{\parallel}$ emerge in the atomic plane due to
the difference in current densities in the two graphene layers. To establish a
relationship between these two quantities with electric and magnetic fields, a
symmetry analysis and some assumptions were made,Stauber  _et al._ (2018a)
such as the linear change of the fields in the space between the two graphene
layers and the zero-frequency limit,Stauber  _et al._ (2018b) or the dephasing
of the current operators in the two graphene layers.Morell  _et al._ (2017)
The magnetization and polarization affect the direction and magnitude of the
electric and magnetic field vectors within the material layer. However, the
problem is that the transverse component $\sigma_{xy}^{(drag)}$ of the
conductivity tensor, which plays a role in creating the effect of one layer on
the other when an electromagnetic field is transmitted through the system, is
treated as the total Hall conductivity of the system in this wave propagation
problem. This may be misleading as the system does not support a finite Hall
conductivity due to time-reversal symmetry. Therefore, it is essential to
consider in detail the electromagnetic wave propagation through the system,
where it is understood that the electromagnetic wave must pass through one
layer before reaching the second layer, and the coupling between the two
layers must play a crucial role. In the next section, we will use Eq. (1)
without considering symmetry analysis and the mentioned assumptions to obtain
the most important theoretical results in Ref. (Stauber  _et al._ , 2018a),
see Eq. (2) therein, and clarify the role of $\sigma_{xy}^{(drag)}$ in the
electromagnetic wave propagation through the TBG system.

Our study encompasses both macroscopic and microscopic perspectives. At the
macroscopic level, we use the transfer matrix method to analyze the light
propagation through a 2D van der Waals system composed of multiple material
layers, each treated as a boundary between two vacuum layers of finite
width.Zhan  _et al._ (2013); Széchenyi  _et al._ (2016) The material system’s
conductivity tensor is decomposed into local and drag terms, with the local
term defining the local relationship between current density and electric
field in each layer, and the drag term defining the non-local relationship
between current density in one layer and the electric field in other layers.
Technically, the finite thickness of the material system makes the assumption
of the spatial homogeneity along the thickness invalid. It therefore requires
evaluating expressions that relate physical quantities through integrals over
$z$ and $z^{\prime}$. This results in increased computational complexity and
the need for a microscopic description of the material surfaces, such as the
potential caused by atoms. To simplify the analysis, we assume that each layer
of atoms has zero thickness. It allow us to express relevant quantities as
sums of delta functions:

$J_{\alpha}(z,\omega)=\sum_{\ell}J_{\alpha}^{(\ell)}\delta(z-z_{\ell}),$ (2)

where $z_{\ell}$ is the position of the $\ell$th material layer, and
$J_{\alpha}^{(\ell)}$ is the current density on that layer.

At the microscopic level, we use the Kubo formula to calculate each term of
the conductivity tensor required for the macroscopic level of description. We
apply an effective model to determine the low-energy states of electrons in
the TBG material lattice, and show that the electron eigen-states can be
expressed as a hybridization of single-layer states in the two graphene
layers. If the atomic lattice has mirror symmetry (the AA-stacked
configuration) or glide symmetry (the AB-stacked configuration) in the lattice
plane, the drag transverse conductivity will be vanished due to the cancelling
contribution of the electron states. However, generic twisted bilayer graphene
lattices have a chiral structure and do not possess these symmetries, leading
to a non-zero drag transverse conductivity that determines the optical
activity of the system. In addition to determining the components of the
optical conductivity tensor, the reflection, transmission, absorption
coefficients, and the circular dichroism spectrum, we also calculate the DC
conductivity of the TBG system. We find that the hybridized states support a
quantum value of $2\sigma_{0}^{DC}$, where $\sigma_{0}^{DC}=\frac{4e^{2}}{\pi
h}$, at the intrinsic Fermi energy. This holds true even for the TBG
configuration with a magic twist angle, which has a high density of states
there due to the presence of a flat band.

The structure of this paper is divided into three sections. After the
“Introduction” section, the next section (Sec. II) provides a solution to the
problem of wave transmission through the twisted bilayer graphene (TBG) system
in Sec. II.1. The solution is based on the transfer matrix method, which
illustrates the contribution of the conductivity tensor components to the
transfer matrix elements such as the reflection and transmission coefficients
and the absorption coefficient. In Sec. II.2, the local and drag components of
the optical conductivity tensor are calculated through the use of the Kubo
formula, while the DC conductivity of the TBG system is calculated using the
Kubo-Greenwood formula. In Sec. II.3, an effective model for the electrons in
the TBG lattices is presented. In Sec. III, numerical results for the
electronic band structure of various TBG configurations, the analysis of the
optical conductivity tensor elements, and the spectra of circular dichroism as
a function of photon frequency are presented. The conclusion is presented in
the final section, Sec. IV.

Figure 1:  The energy band structures of electrons in three TBG configurations
with twist angles $\theta=3.890^{\circ}$ (a, d), $\theta=1.890^{\circ}$ (b,
e), and $\theta=1.050^{\circ}$ (c, f) are analyzed. The blue curves in (a),
(b), and (c) represent the dispersion and density of states of graphene
calculated without considering the interlayer coupling. Panels (d), (e), and
(f) show the dispersion curves for the two valleys $\xi=+1$ and $\xi=-1$ in
blue and red, respectively.

##  II Theory and calculation methods

###  II.1 Optical chiral responses

In this section, we investigate the optical response of TBG systems by
formulating a theory for the propagation of electromagnetic waves through a
generic TBG sheet. Unlike typical dielectric slabs with two separate surfaces,
the electric current densities in the two graphene layers of the TBG sheet are
interdependent. The current density in one graphene layer is affected by both
the electric field within that layer and the electric field in the other
layer. This relationship is described by Eq. (1), which we use to solve the
problem at hand.

We set up the system, following Ref. Széchenyi  _et al._ , 2016, with two
graphene layers separating the entire space into three regions, each
characterized by parameters $(\epsilon_{i},\mu_{i})$, where $i=1,2,3$ labels
the mediums. Our goal is to determine the transfer matrix of the system and to
develop expressions for the transmission and reflection matrices of the
bilayer graphene system. As the inhomogeneity of the space only occurs in the
$z$ direction perpendicular to the TBG surface, the Maxwell equations support
a plane wave solution of the form:

$\boldsymbol{F}_{i}(x,y,z)=\left[\boldsymbol{F}_{i}^{+}e^{ik_{iz}z}+\boldsymbol{F}_{i}^{-}e^{-ik_{iz}z}\right]e^{ik_{ix}x}e^{ik_{iy}y},$
(3)

where $\boldsymbol{F}$ represents $\boldsymbol{E}$ (the electric field),
$\boldsymbol{D}$ (the electric induction), $\boldsymbol{H}$ (the magnetic
field), and $\boldsymbol{B}$ (the magnetic induction). The vector coefficients
$\boldsymbol{F}_{i}^{\pm}$ depend on the wave vector $\boldsymbol{k}_{i}$ in
general. The conditions for the continuity of the vector fields at the
interface between two mediums are given by:

$\displaystyle\boldsymbol{n}\times[\boldsymbol{E}_{2}(1)-\boldsymbol{E}_{1}(1)]=0,$
(4a)
$\displaystyle\boldsymbol{n}\times[\boldsymbol{H}_{2}(1)-\boldsymbol{H}_{1}(1)]=\boldsymbol{J}(1),$
(4b)

for the first graphene layer surface, and

$\displaystyle\boldsymbol{n}\times[\boldsymbol{E}_{3}(2)-\boldsymbol{E}_{2}(2)]=0,$
(5a)
$\displaystyle\boldsymbol{n}\times[\boldsymbol{H}_{3}(2)-\boldsymbol{H}_{2}(2)]=\boldsymbol{J}(2).$
(5b)

for the second graphene layer surface. In Eqs. (4b) and (5b), the distribution
of current densities on each graphene layer is linearly related to the
electric field by the conductivity tensor given via Eq. (27). To the
completeness of the setup, we rewrite this equation in the explicit form:

$\displaystyle\boldsymbol{J}(1)$
$\displaystyle=\boldsymbol{\sigma}^{(1)}\boldsymbol{E}_{1}(1)+\boldsymbol{\sigma}^{(drag)}\boldsymbol{E}_{2}(2),$
(6a) $\displaystyle\boldsymbol{J}(2)$
$\displaystyle=\boldsymbol{\sigma}^{(drag)\dagger}\boldsymbol{E}_{1}(1)+\boldsymbol{\sigma}^{(2)}\boldsymbol{E}_{2}(2)$
$\displaystyle=\boldsymbol{\sigma}^{(drag)\dagger}\boldsymbol{E}_{2}(1)+\boldsymbol{\sigma}^{(2)}\boldsymbol{E}_{3}(2).$
(6b)

The continuity conditions of the fields at interfaces allow us to determine
the relationship between the vector fields
$(\boldsymbol{F}_{i}^{+},\boldsymbol{F}_{i}^{-})$ in two different media at
their interface. In particular, the magnetic field $\boldsymbol{H}_{i}$ is
related to the electric field $\boldsymbol{E}_{i}$ through the Maxwell
equation $\text{curl}(\boldsymbol{E})=-\partial_{t}\boldsymbol{B}$, where
$\boldsymbol{B}=\mu\boldsymbol{H}$. Utilizing the plane wave solution (3), we
obtain the expression for $\boldsymbol{H}_{i}^{\pm}$:

$\boldsymbol{H}_{i}^{\pm}=\pm\sqrt{\frac{\epsilon_{i}}{\mu_{i}}}\frac{\boldsymbol{k}_{i}\times\boldsymbol{E}_{i}^{\pm}}{|\boldsymbol{k}_{i}|},$
(7)

where $\epsilon_{i}$ and $\mu_{i}$ represent the absolute permittivity and
absolute permeability of medium $i$, respectively. This equation makes use of
the dispersion relation $\omega=c|\boldsymbol{k}_{i}|/n_{i}$, where $c$ is the
speed of light in vacuum and $n_{i}$ is the refractive index of medium $i$. It
is worth noting that $\sqrt{\epsilon_{i}/\mu_{i}}$ can be expressed through
$n_{i}$ as $\sqrt{\epsilon_{i}/\mu_{i}}\approx n_{i}\sigma_{0}/2\alpha$, where
$\sigma_{0}=e^{2}/h$ is the unit of quantum conductivity and
$\alpha=e^{2}/(4\pi\epsilon_{0}\hbar c)$ is the fine-structure constant. In
this study, the problem is solved for light transmission along the normal
direction to the TBG layer, so $\boldsymbol{k}_{i}=k_{i}\boldsymbol{n}$. The
transfer matrix $\mathbf{M}{31}$ (of size $4\times 4$) is then derived,
relating the components of the electric fields in the first and third mediums
as follows:

$\displaystyle\left(\begin{array}[]{c}\boldsymbol{E}_{3}^{+}(2)\\\
\boldsymbol{E}3^{-}(2)\end{array}\right)$
$\displaystyle=\mathbf{M}{31}\left(\begin{array}[]{c}\boldsymbol{E}_{1}^{+}(1)\\\
\boldsymbol{E}_{1}^{-}(1)\end{array}\right).$ (12)

Here, $\mathbf{M}_{31}=\mathbf{M}_{32}\mathbf{M}^{f}_{2}\mathbf{M}_{21}$,
where:

$\displaystyle\mathbf{M}_{32}$
$\displaystyle=\left(\begin{array}[]{cc}\boldsymbol{1}&\boldsymbol{1}\\\
-\sqrt{\frac{\epsilon_{3}}{\mu_{3}}}\boldsymbol{1}-\boldsymbol{\sigma}^{(2)}&\sqrt{\frac{\epsilon_{3}}{\mu_{3}}}\boldsymbol{1}-\boldsymbol{\sigma}^{(2)}\end{array}\right)^{-1}\left(\begin{array}[]{cc}\boldsymbol{1}&\boldsymbol{1}\\\
-\sqrt{\frac{\epsilon_{2}}{\mu_{2}}}\boldsymbol{1}+\boldsymbol{\sigma}^{(drag)\dagger}e^{-ik_{2z}d}&\sqrt{\frac{\epsilon_{2}}{\mu_{2}}}\boldsymbol{1}+\boldsymbol{\sigma}^{(drag)\dagger}e^{ik_{2z}d}\end{array}\right),$
(13e) $\displaystyle\mathbf{M}_{21}$
$\displaystyle=\left(\begin{array}[]{cc}\boldsymbol{1}&\boldsymbol{1}\\\
-\sqrt{\frac{\epsilon_{2}}{\mu_{2}}}\boldsymbol{1}-\boldsymbol{\sigma}^{(drag)}e^{ik_{2z}d}&\sqrt{\frac{\epsilon_{2}}{\mu_{2}}}\boldsymbol{1}-\boldsymbol{\sigma}^{(drag)}e^{-ik_{2z}d}\end{array}\right)^{-1}\left(\begin{array}[]{cc}\boldsymbol{1}&\boldsymbol{1}\\\
-\sqrt{\frac{\epsilon_{1}}{\mu_{1}}}\boldsymbol{1}+\boldsymbol{\sigma}^{(1)}&\sqrt{\frac{\epsilon_{1}}{\mu_{1}}}\boldsymbol{1}+\boldsymbol{\sigma}^{(1)}\end{array}\right),$
(13j) $\displaystyle\mathbf{M}_{2}^{f}$
$\displaystyle=\left(\begin{array}[]{cc}e^{ik_{2z}d}\boldsymbol{1}&0\\\
0&e^{-ik_{2z}d}\boldsymbol{1}\end{array}\right),$ (13m)

with $\boldsymbol{1}$ the $2\times 2$ identity matrix.

Let $\mathbf{r}$ and $\mathbf{t}$ be the $2\times 2$ matrices that relate the
components of the electric field vector of the incident light to those of the
reflected and transmitted light, i.e.,
$\boldsymbol{E}_{1}^{-}(1)=\mathbf{r}\cdot\boldsymbol{E}_{1}^{+}$ and
$\boldsymbol{E}_{3}^{+}=\mathbf{t}\cdot\boldsymbol{E}_{1}^{+}$. These matrices
are referred to as the reflection and transmission matrices, respectively.
From Eq. (12), we can deduce the following:

$\displaystyle\mathbf{r}$
$\displaystyle=-[\mathbf{M}_{31}^{22}]^{-1}\mathbf{M}_{31}^{21},$ (14a)
$\displaystyle\mathbf{t}$
$\displaystyle=\mathbf{M}_{31}^{11}-\mathbf{M}_{31}^{12}[\mathbf{M}_{31}^{22}]^{-1}\mathbf{M}_{31}^{21}=[\mathbf{M}_{31}^{-1}]^{11}.$
(14b)

From these matrices, the reflectance ($R$) and transmittance ($T$) are
determined as the ratio of energy fluxes, specifically:

$\displaystyle R$
$\displaystyle=\left|\frac{(\mathbf{r}\cdot\boldsymbol{E}_{1}^{+})^{\dagger}\cdot(\mathbf{r}\cdot\boldsymbol{E}_{1}^{+})}{\boldsymbol{E}_{1}^{+\dagger}\cdot\boldsymbol{E}_{1}^{+}}\right|,$
(15a) $\displaystyle T$
$\displaystyle=\left|\frac{n_{3}}{n_{1}}\frac{(\mathbf{t}\cdot\boldsymbol{E}_{1}^{+})^{\dagger}\cdot(\mathbf{t}\cdot\boldsymbol{E}_{1}^{+})}{\boldsymbol{E}_{1}^{+\dagger}\cdot\boldsymbol{E}_{1}^{+}}\right|.$
(15b)

On the base of the conservation of energy flux, the absorptance is then
determined by $A=1-(R+T)$. The values of $R$, $T$, and $A$ are dependent not
only on the energy, but also on the polarization states of light. It is
important to note that the left- and right-handed polarization states of light
are characterized by electric field vectors,
$\boldsymbol{E}_{L,R}^{+}\propto(1,\pm i)^{T}/\sqrt{2}$ (The supersript $T$
here denotes the transpose operation). To measure the dependence of light
absorption on these polarization states, we calculate the spectrum of circular
dichroism (CD), defined as:

$\mathrm{CD}=\frac{A_{\mathrm{L}}-A_{\mathrm{R}}}{A_{\mathrm{L}}+A_{\mathrm{R}}}.$
(16)

In the following section, we present the method for calculating the
conductivity tensors $\boldsymbol{\sigma}^{(\ell)}$ and
$\boldsymbol{\sigma}^{(drag)}$, and we will use the results obtained to
calculate the transmittance, reflectance, absorbance and circular dichroism
via Eqs. (15a), (15b) and (16).

Figure 2:  The longitudinal DC conductivities $\sigma_{xx}^{DC}$ (represented
by solid curves) and $\sigma_{yy}^{DC}$ (represented by dashed curves) for
three TBG configurations with twist angles of $\theta=3.890^{\circ}$
(represented by blue curves), $\theta=1.890^{\circ}$ (represented by green
curves), and $\theta=1.050^{\circ}$ (represented by red curves).

###  II.2 AC and DC conductivities

The optical properties of bilayer graphene systems are notable, as evidenced
by previous research.Moon and Koshino (2013); Stauber  _et al._ (2018a); Le
and Do (2018); Zhang (2019) When taking into account the layer thickness, the
effects of spatial dispersion become evident. To describe these effects, we
decompose the current response to an external electric field caused by an
incident monochromatic light beam perpendicular to the material plane. To
start, we rewrite Equation (1) in the real space representation as follows:

$\displaystyle J_{\alpha}(\boldsymbol{x},z,\omega)$ $\displaystyle=\int
dz^{\prime}\sigma_{\alpha\beta}(\boldsymbol{x},z;\boldsymbol{x},z^{\prime},\omega)E_{\beta}(\boldsymbol{x},z^{\prime},\omega),$
(17)

where $\mathbf{x}=(x,y)$ represents the position vector in the material plane
with $x$ and $y$ as its coordinates. On the assumption of spatial homogeneity
in the material plane, we can write
$\sigma_{\alpha\beta}(\boldsymbol{x},z;\boldsymbol{x},z^{\prime},\omega)=\sigma_{\alpha\beta}(z,z^{\prime},\omega)$,
and:

$\displaystyle\sigma_{\alpha\beta}(z,z^{\prime},\omega)=$
$\displaystyle\sum_{\ell=1}^{2}\sigma_{\alpha\beta}(z,z^{\prime},\omega)\delta(z^{\prime}-z_{\ell}).$
(18)

Eq. (17) therefore becomes:

$\displaystyle
J_{\alpha}(\boldsymbol{x},z,\omega)=\sum_{\ell=1}^{2}\sigma_{\alpha\beta}(z,z_{\ell},\omega)E_{\beta}(\boldsymbol{x},z_{\ell},\omega).$
(19)

More specifically, we resolve the current density induced in each graphene
layer as follows:

$J_{\alpha}^{(\ell)}(\boldsymbol{x},\omega)=\sigma_{\alpha\beta}^{(\ell)}(\omega)E_{\beta}^{(\ell)}(\boldsymbol{x},\omega)+\sigma_{\alpha\beta}^{(drag)}(\omega)E_{\beta}^{(\ell^{\prime})}(\boldsymbol{x},\omega),$
(20)

where the involved quantities are defined by:
$J_{\alpha}^{(\ell)}(\boldsymbol{x},\omega)=J_{\alpha}(\boldsymbol{x},z_{\ell},\omega);E_{\beta}^{(\ell)}(\boldsymbol{x},\omega)=E_{\beta}(\boldsymbol{x},z_{\ell},\omega)$
and
$\sigma_{\alpha\beta}^{(\ell)}(\omega)=\sigma_{\alpha\beta}(z_{\ell},z_{\ell},\omega)$,
and
$\sigma_{\alpha\beta}^{(drag)}(\omega)=\sigma_{\alpha\beta}(z_{\ell},z_{\ell^{\prime}},\omega)$
where $\ell\neq\ell^{\prime}$. Eqs. (20) can be cast in the matrix form:

$\displaystyle\left(\begin{array}[]{c}\boldsymbol{J}^{(1)}\\\
\boldsymbol{J}^{(2)}\end{array}\right)_{\omega}=\left(\begin{array}[]{cc}\boldsymbol{\sigma}^{(1)}&\boldsymbol{\sigma}^{(drag)}\\\
\boldsymbol{\sigma}^{(drag)\dagger}&\boldsymbol{\sigma}^{(2)}\end{array}\right)_{\omega}\left(\begin{array}[]{c}\boldsymbol{E}^{(1)}\\\
\boldsymbol{E}^{(2)}\end{array}\right)_{\omega},$ (27)

where
$\boldsymbol{\sigma}^{(\ell)}(\omega)=\sigma^{(\ell)}(\omega)\boldsymbol{\tau}_{0}$
is called the local part, and
$\boldsymbol{\sigma}^{(drag)}(\omega)=\sigma^{(drag)}(\omega)\boldsymbol{\tau}_{0}-i\sigma_{xy}^{(drag)}(\omega)\boldsymbol{\tau}_{y}=\sigma^{drag}(\omega)-\sigma_{xy}^{drag}\hat{\boldsymbol{z}}\times$
is called the drag part. Here $\boldsymbol{\tau}_{0}$ is the $2\times 2$
identity matrix and $\boldsymbol{\tau}_{2}$ is the conventional second Pauli
matrix. The compact form (27) is identical to results reported in Refs.
Stauber  _et al._ , 2018a; Ochoa and Asenjo-Garcia, 2020. Importantly, it
guarantees the time-reversal symmetry, the rotation symmetry and the layer
interchange symmetry.

The conductivities $\sigma^{(\ell)}(\omega),\sigma^{(drag)}(\omega)$ and
$\sigma^{(drag)}_{xy}(\omega)$ can be determined from the Kubo formula. There
are a number of versions of the Kubo formula for the electrical conductivity
suitable for implementing it in different situations. However, the most
important ingredient we must specify is the velocity operators
$\hat{v}_{\alpha}$. In general, these operators are determined from the
position operator $\hat{x}_{\alpha}$ and the Hamiltonian $\hat{H}$ via the
Heisenberg equation:

$\hat{v}_{\alpha}=\frac{1}{i\hbar}[\hat{x}_{\alpha},\hat{H}]\rightarrow
v_{\alpha}(\boldsymbol{k})=\frac{1}{\hbar}\frac{\partial}{\partial
k_{\alpha}}H(\boldsymbol{k}).$ (28)

According to the linear response theory, we see that the conductivity has the
paramagnetic part only. Using the eigen-vectors of single-particle Hamiltonian
(without the presence of external vector potential
$\boldsymbol{A}(\boldsymbol{x},t)$) as the representation basis,
$\hat{H}|n\rangle=E_{n}|n\rangle$, the elements of the optical conductivity
tensor are given by this formula: Steven G. Louie (2006)

$\displaystyle\sigma_{\alpha\beta}^{(c)}(\omega)$
$\displaystyle=\frac{ie^{2}}{S}\frac{1}{\omega+i\eta}\sum_{m,n}(f_{m}-f_{n})\frac{O^{(c)mn}_{\alpha\beta}}{E_{m}-E_{n}+\hbar(\omega+i\eta)},$
(29)

where $S$ is the area of material sample; $c=\\{\ell,drag\\}$;
$f_{n}=f(E-\mu,k_{B}T)$ is the occupation weight of the energy level $E_{n}$
determined by the Fermi-Dirac function $f$ with $\mu,k_{B}T$ the chemical
potential and thermal energy; $\eta$ is a positive infinitesimal number, and
$O_{\alpha\beta;mn}^{(c)}$ denotes the product of velocity matrix elements:

$\displaystyle O_{\alpha\beta}^{(\ell)mn}$ $\displaystyle=\langle
m|\hat{v}_{\alpha}^{(\ell)}|n\rangle\langle
n|\hat{v}_{\beta}^{(\ell)}|m\rangle,$ (30a) $\displaystyle
O_{\alpha\beta}^{(drag)mn}$ $\displaystyle=\langle
m|\hat{v}_{\alpha}^{(1)}|n\rangle\langle n|\hat{v}_{\beta}^{(2)}|m\rangle.$
(30b)

To the DC conductivity, it might be optimistic to think that the DC
conductivity can be simply obtained from the expression of the optical
conductivity in the limit of zero frequency. However, it is not, at least in
the numerical calculation. This is because for some finite frequency $\omega$,
there is always a finite length scale governing the behavior of electron. This
length scale is $L_{\omega}=2\pi v_{F}/\omega$. Meanwhile, there is no such
length scale for the DC transport. Furthermore, the static transport has the
diffusion nature due to the similarity of the Schrödinger equation and the
diffusion equation. The DC conductivity can be obtained from the linear
response theory. The vector potential is chosen in the form
$\boldsymbol{A}(t)=(-Et,0,0)$, where $E$ is the intensity of a static electric
field. The DC conductivity is calculated using the Kubo-Greenwood formula:Kubo
(1957); Greenwood (1958)

$\sigma_{\alpha\alpha}(\mu,kT)=-\int_{-\infty}^{+\infty}dE\frac{\partial
f(E-\mu,kT)}{\partial E}\sigma_{\alpha\alpha}(E),$ (31)

where

$\sigma_{\alpha\alpha}(E)=\frac{2\pi e^{2}\hbar}{S}\sum_{m,n}|\langle
m|\hat{v}_{\alpha}|n\rangle|^{2}\delta(E-E_{n})\delta(E-E_{m}).$ (32)

This quantity is seen as the DC conductivity at zero temperature. When
performing the above formula, we use the Gaussian function to approximate the
$\delta$-Dirac function:

$\delta(E-E_{n})\approx\frac{1}{\eta\sqrt{\pi}}\exp\left(-\frac{(E-E_{n})^{2}}{\eta^{2}}\right),$
(33)

where $\eta>0$ is a tiny number that is appropriately chosen to smear the
energy levels for the numerical calculation.

###  II.3 Effective continuum model

Effective continuum models for the low-energy states of electrons in TBG
systems have been developed since 2007 by Lopes et al. dos Santos  _et al._
(2007) However, the model proposed by Bistritzer and MacDonald in 2011 is
better known and widely used. In this paper, we present our solution to the
Bistritzer-MacDonald model. Bistritzer and MacDonald (2011) We implement this
model for TBG lattices in which the first layer is clockwise rotated by a
half-twist angle, $\theta_{1}=-\theta/2$, and the second layer is
counterclockwise rotated by $\theta_{2}=+\theta/2$. The unrotated layers are
defined by unit vectors $\boldsymbol{a}_{1}=a\cdot\boldsymbol{e}_{x}$ and
$\boldsymbol{a}_{2}=a\cos\pi/3\cdot\boldsymbol{e}_{x}+a\sin\pi/3\cdot\boldsymbol{e}_{y}$,
where $a$ is the lattice constant of graphene, and $\boldsymbol{e}_{x}$ and
$\boldsymbol{e}_{y}$ are the unit vectors defining the $x$ and $y$ directions,
respectively. The two associated vectors defining the reciprocal lattice are
thus
$\boldsymbol{a}_{1}^{\star}=(4\pi/a)(\cos\pi/6\cdot\boldsymbol{e}_{x}+\sin\pi/6\cdot\boldsymbol{e}_{y})$
and $\boldsymbol{a}_{2}^{\star}=(4\pi/a)\cdot\boldsymbol{e}_{y}$. Under the
twisting, these vectors defining the atomic lattice of the two graphene layers
as well as their reciprocal lattice are determined by applying the rotation
matrix $R_{z}(\theta_{\ell})$, specifically
$\boldsymbol{a}_{i}^{(\ell)}=R_{z}(\theta_{\ell})\cdot\boldsymbol{a}_{i}$, and
$\boldsymbol{a}_{i}^{\star(\ell)}=R_{z}(\theta_{\ell})\cdot\boldsymbol{a}_{i}^{\star}$.
The first Brillouin zone $\text{BZ}^{(\ell)}$ of each graphene layer is
defined by a set of six corner points $\boldsymbol{K}_{i}^{(\ell)}$. For
instance, the point $\boldsymbol{K}_{1}^{(\ell)}$ is determined by the point
$\boldsymbol{K}_{-}^{(\ell)}$ with
$\boldsymbol{K}_{\xi}^{(\ell)}=-\xi(2\boldsymbol{a}_{1}^{\star(\ell)}+\boldsymbol{a}_{2}^{\star(\ell)})/3$,
where $\xi=\pm 1$.

For commensurate twisting, the twist angle $\theta$ is determined by two
integer numbers $m$ and $n$ via the formula:

$\theta=\arctan\left(\frac{|n^{2}-m^{2}|\sin(\pi/3)}{(n^{2}+m^{2})\cos(\pi/3)+2mn}\right).$
(34)

When the absolute value of the difference between $m$ and $n$ is equal to 1,
two unit vectors
$\boldsymbol{A}_{i}^{\star}=\boldsymbol{a}_{i}^{\star(1)}-\boldsymbol{a}_{i}^{\star(2)}$
define the moiré reciprocal lattice. The mini-Brillouin zone comprises six
$\boldsymbol{K}^{M}_{i}$ points.
$\boldsymbol{K}_{1}^{M}=(-\boldsymbol{A}_{1}^{\star}+\boldsymbol{A}_{2}^{\star})/3$,
and
$\boldsymbol{K}^{M}_{6}=(\boldsymbol{A}_{1}^{\star}+2\boldsymbol{A}_{2}^{\star})/3$
define two of these points. The remaining points can be obtained by shifting
these points using the reciprocal unit vectors.

We begin by representing the Hamiltonian $H$ for electrons using the layer-
resolution vector basis set $\\{|\ell\rangle,|,\ell=1,2\\}$. Using the
identity $1=\sum_{\ell=1}^{2}|\ell\rangle\langle\ell|$, we can express $H$ as
follows:

$H=\sum_{\ell,\ell^{\prime}=1}^{2}|\ell\rangle
H_{\ell,\ell^{\prime}}\langle\ell^{\prime}|,$ (35)

where $H_{\ell,\ell^{\prime}}=\langle\ell|H|\ell^{\prime}\rangle$. Next, we
can further specify $H_{\ell,\ell^{\prime}}$ by using the lattice-resolution
vector basis set
${|\boldsymbol{k},m\rangle=|\boldsymbol{k}+\boldsymbol{G}_{m}\rangle,|,\boldsymbol{k}\in
BZ,m\in\mathbb{Z}}$. We can then use the identity $1=\sum{\boldsymbol{k}\in
BZ}\sum_{m}|\boldsymbol{k},m\rangle\langle\boldsymbol{k},m|$ to obtain:

$\displaystyle H$ $\displaystyle=\sum_{\boldsymbol{k}\in
BZ}\sum_{\ell,m}\sum_{\ell^{\prime},n}|\ell,\boldsymbol{k},m\rangle
H_{\ell,m;\ell^{\prime},n}(\boldsymbol{k})\langle\ell^{\prime},\boldsymbol{k},n|,$
(36)

where $|\ell,\boldsymbol{k},m\rangle=|\ell\rangle|\boldsymbol{k},m\rangle$ and
we use the relation
$\langle\ell,\boldsymbol{k},m|H|\ell^{\prime},\boldsymbol{k}^{\prime},n\rangle=H_{\ell,m;\ell^{\prime},n}(\boldsymbol{k})\delta_{\boldsymbol{k},\boldsymbol{k}^{\prime}}$.
In the case of twisted bilayer graphene (TBG) configurations with tiny twist
angles, the low energy states of electrons are distinguished by a quantum
index $\xi=\pm 1$, which corresponds to the two nonequivalent Dirac valleys
$\boldsymbol{K}_{\xi}^{(1,2)}$ of the graphene monolayers. As the twist angle
becomes smaller, the positions of the two points $\boldsymbol{K}_{\xi}^{(1)}$
and $\boldsymbol{K}_{\xi}^{(2)}$ become closer to each other. The Bistritzer-
MacDonald model is represented by a Hamiltonian in real-space representation
as follows:Bistritzer and MacDonald (2011)

$\displaystyle\hat{H}^{\xi}=\begin{pmatrix}H_{1}^{\xi}(\hat{\boldsymbol{p}})&T^{\xi}(\hat{\boldsymbol{r}})\\\
T^{\xi\dagger}(\hat{\boldsymbol{r}})&H_{2}^{\xi}(\hat{\boldsymbol{p}})\end{pmatrix}.$
(37)

Here, $H_{1}^{\xi}(\hat{\boldsymbol{p}})$ and
$H_{2}^{\xi}(\hat{\boldsymbol{p}})$ are the blocks for the Hamiltonians of the
two monolayers, and $T^{\xi}(\hat{\boldsymbol{r}})$ is the interlayer coupling
block, which are determined via the momentum and position operators,
$\hat{\boldsymbol{p}}$ and $\hat{\boldsymbol{r}}$, respectively. For the range
of energy around the intrinsic Fermi level the two Hamiltonian blocks
$H_{1}^{\xi}(\hat{\boldsymbol{p}})$ and $H_{2}^{\xi}(\hat{\boldsymbol{p}})$
are well approximated by the 2D Dirac Hamiltonian:

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

where $\sigma_{x},\sigma_{y}$ are two conventional Pauli matrices,
$R^{z}\left(-\theta_{\ell}\right)$ is a matrix that rotates the relevant
vectors around the $Oz$ axis to maintain the canonical form of the Dirac
Hamiltonian, $\boldsymbol{K}_{\xi}^{\ell}$ is the corner point of type
(valley) $\xi$ in the first Brillouin zone of layer $\ell$, and $v_{F}$ is the
Fermi velocity. The minus sign in the equation is due to the negative value of
the hopping parameter $V_{pp\pi}=-2.7$ eV. The block term for interlayer
coupling is expressed as follows:Koshino  _et al._ (2018)

$\displaystyle T^{\xi}(\hat{\boldsymbol{r}})=$
$\displaystyle\left(\begin{array}[]{cc}u&u^{\prime}\\\
u^{\prime}&u\end{array}\right)+\left(\begin{array}[]{cc}u&u^{\prime}\omega^{-\xi}\\\
u^{\prime}\omega^{\xi}&u\end{array}\right)e^{i\xi\delta\boldsymbol{k}_{2}\cdot\hat{\boldsymbol{r}}}$
(43) $\displaystyle+\left(\begin{array}[]{cc}u&u^{\prime}\omega^{\xi}\\\
u^{\prime}\omega^{-\xi}&u\end{array}\right)e^{i\xi\delta\boldsymbol{k}_{3}\cdot\hat{\boldsymbol{r}}},$
(46)

The values of $\delta\boldsymbol{k}_{2}$ and $\delta\boldsymbol{k}_{3}$ are
defined as $\boldsymbol{A}_{1}^{\star}$ and
$\boldsymbol{A}_{1}^{\star}+\boldsymbol{A}_{2}^{\star}$, respectively, where
$\boldsymbol{A}_{1}^{\star}$ and $\boldsymbol{A}_{2}^{\star}$ are reciprocal
lattice vectors and $\omega=\exp(i2\pi/3)$. To account for the effects of
lattice reconstruction, the parameters $u$ and $u^{\prime}$, which are set to
be $0.0797$ eV and $0.0975$ eV, respectively, are chosen based on the
description in reference Koshino  _et al._ , 2018.

The velocity operator corresponding to the Hamiltonian in Eq. (38) can be
expressed as follows:

$\hat{v}_{\alpha}=-v_{F}(\xi\sigma_{x},\sigma_{y})\cdot
R^{z}_{\alpha}(-\theta_{\ell}),$ (47)

Here, $R^{z}_{\alpha}(-\theta_{\ell})$ denotes the first or second column of
the rotation matrix $R^{z}(-\theta_{\ell})$ depending on whether $\alpha=x$ or
$\alpha=y$. It is noteworthy that, because of the relativistic nature of
$H^{\xi}_{\ell}(\hat{\boldsymbol{p}})$ in Eq. (35), the current operator
$\hat{j}_{\alpha}$ only comprises a drift current component,
$\hat{j}_{\alpha}=e\hat{v}_{\alpha}$, with no diffusion or diamagnetic terms.

The spectrum of the Hamiltonian $\hat{H}$ is obtained by solving the secular
equation given by:

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

where the electron states in the twisted bilayer graphene (TBG) lattice are
modeled as Bloch state vectors $|\psi_{\xi,\boldsymbol{k}}\rangle$. Because of
the approximate periodicity of the moire lattice in TBG the state vectors
$|\psi_{\xi,\mathbf{k}}\rangle$ can be expanded in terms of the plane-wave
vectors $|\boldsymbol{k}+\boldsymbol{G}_{m}\rangle$ which are the eigenvectors
of the momentum operator $\hat{\boldsymbol{p}}$, i.e.,
$\hat{\boldsymbol{p}}|\boldsymbol{k}+\boldsymbol{G}_{m}\rangle=(\boldsymbol{k}+\boldsymbol{G}_{m})|\boldsymbol{k}+\boldsymbol{G}_{m}\rangle$
and
$\langle\boldsymbol{r}|\boldsymbol{k}+\boldsymbol{G}_{m}\rangle=e^{i(\boldsymbol{k}+\boldsymbol{G}_{m})\cdot\boldsymbol{r}}$.
Specifically, we have:

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

where $C_{1,\xi,\boldsymbol{k}}(\boldsymbol{G}_{m})$ and
$C_{2,\xi,\boldsymbol{k}}(\boldsymbol{G}_{m})$ are 2D vectors of combination
coefficients that need to be found.
$\\{\boldsymbol{G}_{n}\,|\,n=1,2,3,\ldots,N_{\boldsymbol{G}}\\}$ is a set of
$N_{\boldsymbol{G}}$ vectors of the moire reciprocal lattice. Substituting the
expression (51) into Eq. (50) then left-multiplying both sides with
$\langle\boldsymbol{k}+\boldsymbol{G}_{n}|$ we obtain a set of linear
equations in the form:

$\displaystyle\sum_{m}\left[\left(\begin{array}[]{cc}H_{1}^{\xi}(\boldsymbol{k}+\boldsymbol{G}_{n})&T_{1}^{\xi}\\\
T_{1}^{\xi\dagger}&H_{2}^{\xi}(\boldsymbol{k}+\boldsymbol{G}_{n})\end{array}\right)\delta_{\boldsymbol{G}_{n},\boldsymbol{G}_{m}}+\right.$
(54) $\displaystyle\hskip
56.9055pt+\left.\sum_{j=2}^{3}\left(\begin{array}[]{cc}0&T_{j}^{\xi}\\\
0&0\end{array}\right)\delta_{\boldsymbol{G}_{n},\boldsymbol{G}_{m}+\delta\boldsymbol{k}_{j}}+\sum_{j=2}^{3}\left(\begin{array}[]{cc}0&0\\\
T_{j}^{\xi\dagger}&0\end{array}\right)\delta_{\boldsymbol{G}_{n},\boldsymbol{G}_{m}-\delta\boldsymbol{k}_{j}}\right]\left(\begin{array}[]{c}C_{1,\xi,\boldsymbol{k}}(\boldsymbol{G}_{m})\\\
C_{2,\xi,\boldsymbol{k}}(\boldsymbol{G}_{m})\end{array}\right)=E\left(\begin{array}[]{c}C_{1,\xi,\boldsymbol{k}}(\boldsymbol{G}_{n})\\\
C_{2,\xi,\boldsymbol{k}}(\boldsymbol{G}_{n})\end{array}\right),$ (63)

wherein we denote:

$\displaystyle\langle\boldsymbol{k}+\boldsymbol{G}_{n}|H_{\ell}^{\xi}|\boldsymbol{k}+\boldsymbol{G}_{m}\rangle$
$\displaystyle=H_{\ell}^{\xi}(\boldsymbol{k}+\boldsymbol{G}_{n})\delta_{\boldsymbol{G}_{n},\boldsymbol{G}_{m}},$
(64a)
$\displaystyle\langle\boldsymbol{k}+\boldsymbol{G}_{n}|T^{\xi}(\hat{\boldsymbol{r}})|\boldsymbol{k}+\boldsymbol{G}_{m}\rangle$
$\displaystyle=\sum_{j=1}^{3}T_{j}^{\xi}\delta_{\boldsymbol{G}_{n},\boldsymbol{G}_{m}+\delta\boldsymbol{k}_{j}},$
(64b)
$\displaystyle\langle\boldsymbol{k}+\boldsymbol{G}_{n}|T^{\xi\dagger}(\hat{\boldsymbol{r}})|\boldsymbol{k}+\boldsymbol{G}_{m}\rangle$
$\displaystyle=\sum_{j=1}^{3}T_{j}^{\xi\dagger}\delta_{\boldsymbol{G}_{n},\boldsymbol{G}_{m}-\delta\boldsymbol{k}_{j}},$
(64c)

To obtain a numerical solution, we define a $4N_{\boldsymbol{G}}\times
4N_{\boldsymbol{G}}$ Hermitian matrix $H_{\boldsymbol{k}}^{\xi}$ in block form
for each value of $\xi$ and $\boldsymbol{k}\in\text{MBZ}$:

$\displaystyle[H_{\boldsymbol{k}}^{\xi}]_{n,n}=\left(\begin{array}[]{cc}H_{1}^{\xi}(\boldsymbol{k}+\boldsymbol{G}_{n})&T_{1}^{\xi}\\\
T_{1}^{\xi\dagger}&H_{2}^{\xi}(\boldsymbol{k}+\boldsymbol{G}_{n})\end{array}\right),$
(65c)
$\displaystyle[H_{\boldsymbol{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\boldsymbol{G}_{m_{j}}=\boldsymbol{G}_{n}-\delta\boldsymbol{k}_{j},$
(65f)
$\displaystyle[H_{\boldsymbol{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\boldsymbol{G}_{m_{j}^{\prime}}=\boldsymbol{G}_{n}+\delta\boldsymbol{k}_{j}.$
(65i)

The eigenvalues $E_{n}^{\xi}(\boldsymbol{k})$ and corresponding eigenvectors,
obtained by diagonalizing the matrix $H_{\boldsymbol{k}}^{\xi}$, approximate
the states of low-energy electrons in the TBG lattice. The accuracy of the
effective description is numerically dependent on the number of reciprocal
lattice vectors
$\\{\boldsymbol{G}_{n}\,|\,n=1,2,3,\ldots,N_{\boldsymbol{G}}\\}$ considered in
the calculation. To ensure validity in the low-energy range near the intrinsic
Fermi level, where the energy surfaces of the monolayer graphene take the form
of cones, we control the value of $N_{\boldsymbol{G}}$ by using a cutoff
energy $E_{c}$. Specifically, we count only the $\boldsymbol{G}_{n}$ vectors
such that $|\boldsymbol{G}_{n}|\leq E_{c}/\hbar v_{F}$.

Figure 3:  The real and imaginary components of the overall longitudinal
conductivity for TBG configurations with $\theta=3.890^{\circ}$ (a),
$1.890^{\circ}$ (b), and $1.050^{\circ}$ (c). There are no total transverse
conductivities due to the presence of time-reversal symmetry.

Figure 4:  The real and imaginary parts of the diagonal element
$\sigma_{xx}^{(\ell)}$ of the tensor $\boldsymbol{\sigma}^{(\ell)}(\omega)$
($\ell=1,2$) for three TBG configurations with
$\theta=3.890^{\circ},1.890^{\circ}$ and $1.050^{\circ}$.

##  III Results and discussions

###  III.1 Electronic band structure of TBGs

We used numerical methods to solve Eq. (54). Specifically, for each value of
$\xi=\pm 1$, we constructed a Hermitian matrix $H_{\boldsymbol{k}}^{\xi}$ of
size $4N_{\boldsymbol{G}}\times 4N_{\boldsymbol{G}}$ and diagonalized it for
each value of $\boldsymbol{k}$ in the mini-Brillouin zone of the TBG system.
The sets of eigenvalues ${E_{m}(\boldsymbol{k})}$ and eigenvectors
$\\{\mathbf{C}^{m}_{\xi}(\boldsymbol{k})=[C^{m}_{1,\xi,\boldsymbol{k}}(\boldsymbol{G}_{1}),C^{m}_{2,\xi,\boldsymbol{k}}(\boldsymbol{G}_{1}),\ldots,C^{m}_{1,\xi,\boldsymbol{k}}(\boldsymbol{G}_{N_{\boldsymbol{G}}}),C^{m}_{2,\xi,\boldsymbol{k}}(\boldsymbol{G}_{N_{\boldsymbol{G}}})]^{T}\\}$
were calculated. The results for three TBG configurations with twist angles of
$\theta=1.05^{\circ},1.89^{\circ}$ and $3.89^{\circ}$ are shown in Fig. 1. The
blue curves represent the results obtained when the electronic interlayer
coupling is artificially turned off. Both sets of data are plotted on the same
figure for comparison, to highlight the impact of the interlayer coupling. The
presence of multiple dispersion curves in the figure results from the folding
of the energy band structure of two graphene layers due to the enlargement of
the unit cell of the TBG lattice and the corresponding shrinkage of the
Brillouin zone into the mini-Brillouin zone. However, the electronic band
structure of the TBG system is not simply formed in this manner. At the points
where the bands of monolayer graphene cross, the interlayer coupling leads to
the hybridization of Bloch states and the lifting of energy degeneracy to form
new bands. Our results agree with other available data in the literature
Morell  _et al._ (2010); de Laissardiere  _et al._ (2012); Moon and Koshino
(2013); Koshino  _et al._ (2018). However, in this study, we wish to emphasize
the formation of new electronic states, which are special to TBG systems and
play a crucial role in determining its electronic, optical and transport
properties. Fig. 1 can be considered a visual representation of the evolution
of the band structure with respect to the twist angle. We observe that for the
configuration with $\theta=3.89^{\circ}$, the dispersion curves in the low
energy range around the intrinsic Fermi energy ($E_{F}=0$) have a similar form
to the blue curves, although they are not identical. This implies that the
states created in the bilayer system share similarities with those in
monolayer graphene. For smaller twist angles, particularly
$\theta=1.05^{\circ}$, a band with a very narrow bandwidth forms around the
Fermi energy, separated from the lower and upper bands by a narrow gap. These
dispersion curves are completely different from the blue curves, as seen in
Fig. 1(e). We conclude that the electronic interlayer coupling creates special
states in TBG lattices with small twist angles. Although the TBG systems are
two-dimensional, it can be challenging to analyze their geometrical features
of energy surfaces in wave-vector space. However, the density of states (DOS)
generally exhibits key features, such as van Hove singularity behavior, as
shown in the right panels of Figs. 1 (a, b, c). The blue curves in these
figures were obtained numerically without considering the effects of
electronic interlayer coupling, which is equivalent to the DOS of two
independent graphene layers. The red curves show significant peaks due to van
Hove singularities, indicating the presence of extremal and saddle points of
the energy surfaces.

It’s worth noting that the energy bands with index $\xi=+1$ (resulting from
the hybridization of graphene Bloch states in the $K_{1}$ valley of the first
layer and $K_{2}$ valley of the second layer) differ from those with index
$\xi=-1$ (resulting from the hybridization of Bloch states in the
$K_{1}^{\prime}$ valley and $K_{2}^{\prime}$). This difference is evident
along the ${\it\Gamma}M$ direction as shown in Figs. 1 (d, e, f). However, it
does not appear in the DOS picture. We computed the DOS for both $\xi=+1$ and
$\xi=-1$ separately and obtained the same results. This implies that the
energy surfaces with $\xi=+1$ are simply rotations of the corresponding ones
with $\xi=-1$ by an angle of $\pi/3$.

Figure 5:  The real and imaginary parts of the drag component of the
conductivity tensor, $\boldsymbol{\sigma}^{(drag)}(\omega)$, for three TBG
configurations with $\theta=3.890^{\circ},1.890^{\circ}$, and $1.050^{\circ}$.
(a,c,e) represent the longitudinal (diagonal) component
$\sigma_{xx}^{(drag)}(\omega)$, and (b,d,f) represent the transverse (off-
diagonal) component $\sigma_{xy}^{(drag)}(\omega)$.

###  III.2 DC conductivity

The calculation of the electronic structures of TBG configurations with twist
angles greater than 2 degrees shows the presence of two energy bands near the
Fermi energy ($E_{F}=0$), which have dispersion curves similar to those of
monolayer graphene. However, the Fermi velocity of the linear dispersion law
in TBG ($v_{F}^{TBG}$) is smaller than that of monolayer graphene
($v_{F}^{MLG}$). As the twist angle decreases, these two bands become closer
and eventually form an isolated band with a narrow bandwidth. This “flat band”
near the Fermi energy is expected to have limited impact on the transport
properties of TBG sheets because of its small Fermi velocity. In other words,
electrons occupying these bands are expected to be spatially confined in the
atomic lattice. Although these bands are not completely dispersionless, the
study of their transport properties is still of great interest. To this end,
within the framework of the non-interacting electron approximation, the DC
conductivity was calculated using the Kubo-Greenwood formula as given by Eqs.
(31) and (32). The results of the calculation for the longitudinal components
($\sigma_{xx}^{DC}$ and $\sigma_{yy}^{DC}$) of three TBG configurations at low
temperatures are presented in Fig. 2.

The data indicate that the electrical conductivity of TBG systems is
anisotropic, which is reflected in the two different longitudinal elements of
the conductivity tensor, $\sigma_{xx}^{DC}$ (represented by solid curves) and
$\sigma_{yy}^{DC}$ (represented by dashed curves). The conductivity should
reflect the electronic structure with respect to the Fermi energy, and our
calculations for the two cases of $\xi=\pm 1$ yielded the same result. This is
similar to the picture of the density of states, but the difference between
the two longitudinal conductivities highlights the anisotropic nature of the
energy surfaces. For the TBG configuration with a twist angle of
$\theta=3.890^{\circ}$, the electronic structure is similar to that of
monolayer graphene (except for the value of the Fermi velocity), and the
conductivity curves have a characteristic V-shape (or an M-shape when
considered over a wider energy range, as shown by the blue curves). For TBG
configurations with smaller twist angles, the electronic structure and thus
the conductivity picture becomes much more complex. However, our calculations
show that the conductivity at the intrinsic Fermi level $E_{F}=0$ is always
finite and non-zero. Interestingly, the calculated DC conductivity for TBG
configurations with twist angles that are not magic values is about two times
the quantum value of the minimal conductivity of monolayer graphene, which is
$\sigma^{DC}_{0}=4e^{2}/\pi h$. In the case of the TBG configuration with the
magic twist angle $\theta=1.050^{\circ}$, the conductivity curve shows a small
hump at the Fermi energy $E_{F}=0$ with
$\sigma_{\alpha\alpha}^{DC}(0)\approx\sigma^{DC}_{0}$. This quantum value of
the conductivity of monolayer graphene has been intensively discussed in the
past and is rooted in the fact that the valence and conduction bands touch
each other at the $K$ points, resulting in always-available free carriers due
to fluctuations in the electrostatic profile around the Fermi energy
level.Novoselov  _et al._ (2005); Katsnelson (2006); Tworzydlo  _et al._
(2006); Miao  _et al._ (2007); Blake  _et al._ ; Do and Dollfus (2010) In our
calculation, we approximated the delta-Dirac functions in the Kubo-Greenwood
formula with a Gaussian function with a parameter $\eta$ that defines the
finite width of the peak. The value of $\eta$ was taken to be $\eta<5$ meV to
account for both numerical calculation needs and the broadening of energy
levels and finite lifetime of the electric charge-carrying quasi-particles.
Although the observed DC conductivity value in TBG systems is supported by the
characteristic hybridized states in the TBG lattice, we suppose that it has
the same underlying physical explanation as in monolayer graphene.

###  III.3 AC conductivity and optical Hall drag conductivity

In Fig. 3, we present the results of our calculations for the total optical
conductivity of three different twisted bilayer graphene (TBG) configurations.
The calculation was performed using the Kubo formula for the conductivity
tensor, where the velocity operator $\hat{v}_{\alpha}$ was determined from the
Hamiltonian (37) using Eq. (28). The time reversal symmetry and the Onsager
reciprocal relations eliminate the off-diagonal elements $\sigma_{xy}$ and
$\sigma_{yx}$ of the conductivity tensor, which we numerically confirmed
behave as noise with an amplitude of $10^{-3}$. We also verified all of the
symmetry properties of the conductivity tensor according to Eq. (27),
demonstrating the accuracy of our numerical data.

Figure 3 displays the longitudinal conductivity curves $\sigma_{xx}(\omega)$
with significant peaks, which are attributed to dominant interband transition
processes. This correlation is supported by its relationship to the electronic
band structures of the TBG configurations, shown in Fig. 1. For instance, in
Fig. 1(b), at an energy of 0.5681 eV, there is an optical absorption peak that
can be assigned to a transition from the highest valence band to the lowest
conduction band around the $M$ points in the mini-Brillouin zone. As the twist
angle decreases, the electronic structure becomes more complex and the optical
conductivity exhibits more optical absorption peaks.

For TBG configurations with twist angles greater than 2 degrees, there is a
distinctive feature in the low-energy range of the conductivity curve: the
real part of the conductivity is independent of the photon energy
$\hbar\omega$ and takes a value of $2\sigma_{0}$, where $\sigma_{0}=\pi
e^{2}/2h$. This behavior is similar to that of monolayer graphene and results
from the linear dispersion law of electron states. In the low-energy range,
the energy band structure of TBG systems is qualitatively similar to that of
monolayer graphene, but with a smaller Fermi velocity, as previously analyzed.
The value of $\sigma_{xx}(\omega)$ is equal to $2\sigma_{0}$, as found in
previous studies.Le and Do (2018); Do  _et al._ (2020); Moon and Koshino
(2013)

For the case of the “magic angle” $\theta=1.050^{\circ}$, our calculation
shows a clear dominant contribution from optical transitions from the quasi-
flat bands around the Fermi level to the lowest conduction band around the
$\Gamma$ point, as shown in Fig. 3(c) and Figs. 1(c) and 1(f).

Thus far, we would like to emphasize that the use of the total DC conductivity
tensor is not relevant when studying the wave transmission problem in TBG. As
we have presented in Secs. II.1 and II.2, the required quantities to be
utilized are $\boldsymbol{\sigma}^{(\ell)}(\omega)$ and
$\boldsymbol{\sigma}^{(drag)}(\omega)$. In Figs. 4 and 5, we present the
calculation results for the elements of the conductivity tensor blocks
$\boldsymbol{\sigma}^{(\ell)}(\omega)$ and
$\boldsymbol{\sigma}^{(drag)}(\omega)$ as analyzed in Sec. II.2. Our numerical
calculation reveals that
$\boldsymbol{\sigma}^{(1)}(\omega)=\boldsymbol{\sigma}^{(2)}(\omega)=\sigma_{xx}^{(\ell)}(\omega)\boldsymbol{\tau}_{0}$,
verifying the equivalent role of the two graphene layers in the TBG lattice.
Figs. 4(a), (b), and (c) display the real (red) and imaginary (blue) parts of
$\sigma_{xx}^{(\ell)}(\omega)$ for $\ell=1,2$. It is worth noting that
decreasing the twist angle increases the magnitude of $\sigma_{xx}^{(\ell)}$,
particularly when the imaginary part becomes dominant and varies in accordance
with $1/\hbar\omega$, as seen in Fig. 4(c). Fig. 5 shows the real and
imaginary parts of the off-diagonal element
$\boldsymbol{\sigma}^{(drag)}(\omega)$ of the drag part of the conductivity
tensor. For large twist angles in the TBG configuration, there is a notable
difference between $\sigma_{xx}^{(\ell)}(\omega)$ and
$\sigma_{xx}^{(drag)}(\omega)$, see Figs. 4(a) and 5(a). However, for small
twist angles ($<2^{\circ}$), these curves are quite similar but have opposite
signs. As a result, although there is significant variation in the magnitude
of $\sigma_{xx}^{(\ell)}(\omega)$ and $\sigma_{xx}^{(drag)}(\omega)$, they
compensate each other when summed to give the total value
$\sigma_{xx}(\omega)=2[\sigma_{xx}^{(\ell)}(\omega)+\sigma_{xx}^{(drag)}(\omega)]$,
as seen in Fig. 2.

Figures 5(b), (d), and (f) present the calculated results for the transverse
conductivity component $\sigma_{xy}^{(drag)}(\omega)$. The behavior of
$\sigma_{xy}^{(drag)}(\omega)$ is primarily characterized by close-to-zero
values across the energy range, with some sharp peaks and dips appearing at
energy positions with high optical absorption. This finite value of
$\sigma_{xy}^{(drag)}(\omega)$ can be attributed to both the correlation
between the current densities in the two graphene layers and the chiral
structure of the atomic lattices. The former is considered a necessary
condition, representing the result of the coupling between the two graphene
layers, which leads to the formation of composite states like the
superposition of Bloch states between the two layers. This means that when
electrons move in one graphene layer along the $Ox$ direction, they induce a
motion in the second layer along the transverse $Oy$ direction due to
interlayer coupling. The sufficient condition, however, is the contribution of
the electronic states in the system to the resulting conductivity. If the
system had mirror symmetry, these correlations would cancel each other out
completely. However, the TBG lattices lack mirror symmetry, unlike the
hexagonal lattice of graphene, which results in a finite, non-zero
conductivity of $\sigma_{xy}^{(drag)}(\omega)\neq 0$.

To emphasize this point, the hexagonal mini-Brillouin zone of the TBG lattice
was divided into three equal parts, as depicted in the inset of Fig. 6. This
division was based on the assumption of mirror symmetry through the plane
$M_{yz}$, with domain (1) containing the symmetry, and domains (2) and (3)
being its mirror images. Using the Kubo formula, the value of
$\sigma_{xy}^{(drag)}(\omega)$ was calculated by summing over all
$\boldsymbol{k}$ points in each of the three domains. The results, shown in
Fig. 6, indicate that the results in domains (2) and (3) are similar, opposite
in sign, and much larger than the result in domain (1). Upon summing these
results, a strong but incomplete cancellation takes place, leading to the
final value of $\sigma_{xy}^{(drag)}(\omega)$ displayed in Figs. 5(b), (d),
and (f). For the AA- and AB-stacked configurations with mirror symmetries, the
cancellation was found to be complete, resulting in
$\sigma_{xy}^{(drag)}(\omega)=0$. The symmetry analyses in Refs. Morell  _et
al._ , 2010 and Do  _et al._ , 2020 support this theory.

Figure 6:  Contribution of Bloch States to $\sigma_{xy}^{(drag)}(\omega)$: The
red and blue curves represent the contribution of Bloch states with
$\boldsymbol{k}$ in the first third of the Brillouin zone, while the solid
purple and dashed green curves represent the contribution of Bloch states with
$\boldsymbol{k}$ in the second and third thirds of the Brillouin zone,
respectively. Figure 7:  The transmission spectra for three TBG configurations
are shown in (a), (b), and (c) respectively, with twist angles of
$\theta=3.890^{\circ}$, $\theta=1.890^{\circ}$, and $\theta=1.050^{\circ}$.
Figure 8:  The circular dichroism (CD) spectra of three TBG configurations:
with twist angles of $\theta=3.890^{\circ}$ (a), $\theta=1.890^{\circ}$ (b),
and $\theta=1.050^{\circ}$ (c) are displayed.

###  III.4 Optical characteristics

On the basis of the theory established in Sec. II and the data of the
conductivity tensor already computed, we performed the calculation for the
transfer matrix, then the reflection, transmission and absorption spectra of
several TBG configurations. The obtained results for transmittance and CD are
presented in Figs. 7 and 8, respectively. Since the reflectance is negligible,
it is not shown here. It should be noted that while the optical properties of
TBGs have been extensively studied, they are usually analyzed through the
behavior of the real and imaginary parts of the diagonal components of the
total optical conductivity tensor,Nicol and Carbotte (2008); Moon and Koshino
(2013); Zhang (2019) rather than through the transmittance and reflectance,
which are quantities that can be directly measured in experiments. From the
transmission spectra as shown, we see that the transmittance of all TBG
configurations can get the value of 98% and the absorptance of 2% in average
in a large energy range. The spectral curves clearly exhibit the peaks and
dips consistent with the picture of the longitudinal conductivities. The
transmission and absorption spectra for the left- and right-handed lights are
different to be distinguished. However, the CD spectra is clearly manifested.
From Fig. 8, we see that for the TBG configuration with
$\theta=3.890^{\circ}$, the CD curve exhibits essential features of
experimental data reported by Kim et. al.,Kim  _et al._ (2016) for instance,
with the peaks and dips of large width, except several sharp peaks at the
energy position of the absorption peaks. In order to verify the role of the
drag transverse conductivity $\sigma_{xy}^{(drag)}$, we did not include it
into the calculation. As expected, it results in the zero CD for the total
energy range, but mostly does not change the transmission spectrum. We
therefore conclude that this drag transverse conductivity plays the decisive
role in governing the chiral optical behaviors of the TBG systems. Last but
not least, it is important to note that the CD spectrum shown in Fig. 8 is
extracted from the theory presented in Sec. II.1, which may differ from the
formula used in Refs. Morell  _et al._ , 2017; Stauber  _et al._ , 2018b,
where CD is simply proportional to the real part of $\sigma_{xy}$.

##  IV Conclusion

The field of engineering two-dimensional quantum materials for electronic and
optoelectronic applications has been the subject of extensive research. In
this work, we demonstrate that the optical activity observed in the TBG system
results from the spatial dispersion effects of a typical vdW heterostructure
of two graphene layers lacking mirror and glide symmetries. Our analysis
addresses two fundamental aspects of the field-matter interaction problem: the
propagation of electromagnetic waves through a material layer and the material
layer’s response to an external field, both treated simultaneously using the
macroscopic and microscopic descriptions, respectively.

Technically, we utilize an effective continuum model to demonstrate how
electron states in each graphene layer can hybridize, forming electron states
in the bilayer system that allow for the correlation of transverse motions
between the two layers. In terms of electromagnetic wave propagation, we
clarify how the coupling between the two layers can transform the magnitude
and direction of the fields within the material layer. We provide a detailed
solution to the electromagnetic wave propagation problem through a TBG sheet
by considering each graphene layer as a conducting interface between two
dielectric media. The two interfaces are not independent and conduct mutually
influenced current densities. Our solution to the wave propagation problem
clarifies how the roles of the conductivity components characterized by local
and nonlocal spatial effects manifest themselves in the optical behavior of
the material system.

On the basis of explicitly considering the spatial dispersion effects, our
analysis naturally composes the conductivity tensor into the local part and
the drag part without relying on symmetry analysis or assumptions about the
phase shift between current operators in the two layers. We calculate all
elements of these parts using the Kubo formula, and formulate and calculate
the transfer matrix numerically. We demonstrate that the transverse drag
transverse conductivity $\sigma_{xy}^{(drag)}$ plays a critical role in
defining the chiral optical response of the TBG system. This conductivity
arises from the chiral structure of the atomic lattice, which leads to
hybridized states that lack mirror symmetry and support the correlation of
transverse motions between the two graphene layers. We use the results of our
calculations of the components of the optical conductivity tensor to calculate
the transmission and circular dichroism spectra. We show that the overall
transparency of twisted bilayer graphene (TBG) samples, on average, is
generally 98%, but it has the potential to absorb light up to 2%.
Significantly, we establish that the drag transverse conductivity is a
decisive factor in determining the circular dichroism, thereby affirming that
the optical activity of the TBG system is a manifestation of spatial effects.

Furthermore, we found that the DC conductivity of the TBG system exhibits a
quantum conductivity value of $\propto e^{2}/h$ at the intrinsic Fermi energy,
which is determined by the hybridized states in the bilayer system. This
quantum value is obtained from a single-particle approximation. However, a
rigorous discussion of the intrinsic transport properties of TBG
configurations with magic twist angles should be based on a strong correlation
picture of low-energy excited states as Dirac fermions. In summary, our study
emphasizes the importance of considering the finite thickness of twisted
bilayer graphene when analyzing its optical response and transport properties,
owing to the significance of spatial dispersion effects. This theoretical
approach can also be extended to other van der Waals material systems with
multiple layers.

## Acknowledgements

One of the authors, S.T.H., acknowledges the financial support of Hanoi
University of Civil Engineering (HUCE) for his work under grant number
28-2021/KHXD-TD.

## 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).
  * Tse and MacDonald (2011) W.-K. Tse and A. H. MacDonald, Phys. Rev. B 84, 205327 (2011).
  * Crassee  _et al._ (2011) I. Crassee, J. Levallois, A. L. Walter, M. Ostler, A. Bostwick, E. Rotenberg, T. Seyller, D. van der Marel, and A. B. Kuzmenko, Nat. Phys. 7, 48 (2011).
  * Nandkishore and Levitov (2011) R. Nandkishore and L. Levitov, Phys. Rev.Lett. 107, 097402 (2011).
  * Shimano  _et al._ (2013) R. Shimano, G. Yumoto, J. Y. Yoo, R. Matsunaga, S. Tanable, H. Hibino, T. Morimoto, and H. Aoki, Nat. Commun. 4, 1841 (2013).
  * Izake (2007) E. L. Izake, Journal of Pharmaceutical Sciences 96, 1659 (2007).
  * Liu  _et al._ (2005) W. Liu, J. Gan, D. Schlenk, and W. A. Jury, Proceedings of the National Academy of Sciences 102, 701 (2005).
  * Solomon  _et al._ (2019) M. L. Solomon, J. Hu, M. Lawrence, A. García-Etxarri, and J. A. Dionne, ACS Photonics 6, 43 (2019).
  * Andrei and MacDonald (2020) E. Y. Andrei and A. H. MacDonald, Nature Materials 19, 1265 (2020).
  * Nimbalkar and Kim (2020) A. Nimbalkar and H. Kim, Nano-Micro Letters 12, 8712 (2020).
  * Hill (2020) H. M. Hill, Physics Today 73, 18 (2020).
  * Kim  _et al._ (2016) C.-J. Kim, A. Sáchez-Castillo, Z. Ziegler, Y. Ogawa, C. Noguez, and J. Park, Nat. Nanotechnol. 11, 520 (2016).
  * Ribeiro-Palau  _et al._ (2018) R. Ribeiro-Palau, C. Zhang, K. Watanabe, T. Taniguchi, J. Hone, and C. Dean, Science 361, 690 (2018).
  * Cai and Yu (2021) L. Cai and G. Yu, Advanced Materials 33, 2004974 (2021).
  * Landau and Lifshitz (1984) L. D. Landau and E. M. Lifshitz,  _Electrodynamics of continuous media_ , Volume 8 of course of Theoretical Physics (Permagon Press, 1984).
  * Condon (1937) E. U. Condon, Rev. Mod. Phys. 9, 432 (1937).
  * Stauber  _et al._ (2018a) T. Stauber, T. Low, and G. Gómez-Santos, Phys. Rev. Lett. 120, 046801 (2018a).
  * Oliva-Leyva and Naumis (2016) M. Oliva-Leyva and G. G. Naumis, Phys. Rev. B 93, 035439 (2016).
  * Zou  _et al._ (2019) J. Zou, Z. He, and G. Xu, npj Comput. Mater. 5, 96 (2019).
  * Do  _et al._ (2020) V. N. Do, H. A. Le, V. D. Nguyen, and D. Bercioux, Phys. Rev. Research 2, 043281 (2020).
  * Morell  _et al._ (2017) E. S. Morell, L. Chico, and L. Brey, 2D Mater. 4, 035015 (2017).
  * Moon and Koshino (2013) P. Moon and M. Koshino, Phys. Rev. B 87, 205404 (2013).
  * Tabert and Nicol (2013) C. J. Tabert and E. J. Nicol, Phys. Rev. B 87, 121402 (2013).
  * Stauber  _et al._ (2013) T. Stauber, P. San-Jose, and L. Brey, New J. Phys. 15, 113050 (2013).
  * Do  _et al._ (2017) V. N. Do, H. A. Le, and V. T. Vu, Phys. Rev. B 95, 165130 (2017).
  * Nazarov (2015) V. U. Nazarov, New J. Phys. 17, 073018 (2015).
  * Agranovich and Ginzburg (2013) V. M. Agranovich and V. Ginzburg,  _Crystal optics with spatial dispersion, and excitons_, Springer Series in Solid-State Sciences (Springer, 2013).
  * Poshakinskiy  _et al._ (2018) A. V. Poshakinskiy, D. R. Kazanov, T. V. Shubina, and S. A. Tarasenko, Nanophotonics 7, 753 (2018).
  * Stauber  _et al._ (2018b) T. Stauber, T. Low, and G. Gómez-Santos, Phys. Rev. B 98, 195414 (2018b).
  * Zhan  _et al._ (2013) T. Zhan, X. Shi, Y. Dai, X. Liu, and J. Zi, J. Phys.: Condens. Matter 25, 215301 (2013).
  * Széchenyi  _et al._ (2016) G. Széchenyi, M. Vigh, A. Kormányos, and J. Cserti, J. Phys.: Condens. Matter 28, 375802 (2016).
  * Le and Do (2018) H. A. Le and V. N. Do, Phys. Rev. B 97, 125136 (2018).
  * Zhang (2019) M. Zhang, ed.,  _Handbook of graphene_, Volume 3: Graphene-like 2D materials (John Wiley & Sons, 2019) p. 177.
  * Ochoa and Asenjo-Garcia (2020) H. Ochoa and A. Asenjo-Garcia, Phys. Rev. Lett. 125, 037402 (2020).
  * Steven G. Louie (2006) M. L. C. Steven G. Louie, ed.,  _Conceptual Foundations of Materials_, A Standard Model for Ground- and Excited-State Properties, Vol. 2 (2006) Chap. 6.
  * Kubo (1957) R. Kubo, J. Phys. Soc. Japan 12, 570 (1957).
  * Greenwood (1958) D. A. Greenwood, Proc. Phys. Soc. 71, 585 (1958).
  * 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).
  * Bistritzer and MacDonald (2011) R. Bistritzer and A. H. MacDonald, PNAS 108, 12233 (2011).
  * 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).
  * Morell  _et al._ (2010) E. S. Morell, J. D. Correa, P. Vargas, M. Pacheco, and Z. Barticevic, Phys. Rev. B 82, 121407 (2010).
  * de Laissardiere  _et al._ (2012) G. T. de Laissardiere, D. Mayou, and L. Magaud, Phys. Rev. B 86, 125413 (2012).
  * Novoselov  _et al._ (2005) K. S. Novoselov, A. K. Geim, S. V. Morozov, D. J. M. I. Katsnelson, I. V. Grigorieva, S. V. Dubonos, and A. A. Firsov, Nature 438, 197 (2005).
  * Katsnelson (2006) M. I. Katsnelson, Eur. Phys. J. B 51, 51157 (2006).
  * Tworzydlo  _et al._ (2006) J. Tworzydlo, B. Trauzettel, M. Titov, A. Rycerz, and C. W. J. Beenakker, Phys. Rev. Lett. 96, 246802 (2006).
  * Miao  _et al._ (2007) F. Miao, S. Wijeratne, Y. Zhang, U. C. Coskun, W. Bao, and C. N. Lau, Science 317, 1530 (2007).
  * (50) P. Blake, R. Yang, S. V. Morozov, F. Schedin, L. A. Ponomarenko, A. A. Zhukov, R. R. Nair, I. V. Grigorieva, K. S. Novoselov, and A. K. Geim, Solid State Communications 149, 1068.
  * Do and Dollfus (2010) V. N. Do and P. Dollfus, J. Phys.: Condens. Matter 22, 425301 (2010).
  * Nicol and Carbotte (2008) E. J. Nicol and J. P. Carbotte, Phys. Rev. B 77, 155409 (2008).