# Chiral correlation of drag currents inducing optical activity of twisted
bilayer graphene

S. Ta Ho  National University of Civil Engineering (NUCE), 55 Giai Phong road,
Hanoi 10000, Vietnam  V. Nam Do<EMAIL_ADDRESS>Department of
Basic Science, Phenikaa Institute for Advanced Study (PIAS), A1 Building,
Phenikaa University, Hanoi 10000, Vietnam

###### Abstract

The mechanisms of optical activity and quantum transport of twisted bilayer
graphene are studied. The formation of unique electron states in the bilayer
systems is studied using an effective continuum model. Such states are shown
to support the correlation of transverse motions of electrons in two graphene
layers. Because of the chiral structure of the atomic lattices, the
contribution of such drag correlations is incompletely cancelled, thus
resulting in a drag term of the optical conductivity tensor. We show that the
drag term of the conductivity is the manifestation of the spatial dispersion.
We show how to analyze and calculate the components of the conductivity
tensors that governs the optical activity of the systems. The DC conductivity
of the twisted bilayer graphene system is also calculated. It shows the
existence of a quantum conductivity value $\propto e^{2}/h$ at the intrinsic
Fermi energy.

††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. Achievements shows great
potentials of using multiple layer materials to applications in electronics,
optoelectronics, spintronics and valleytronics.Sierra  _et al._ (2021); Liaoa
_et al._ (2019); Sangwan and Hersam (2018); Bao and Hoh (2019) Using multiple
layer materials with the van der Waals (vdW) interlayer coupling brings about
many solutions to tune the electronic structure, e.g., changing the interlayer
distance, changing the lattice alignment by twisting and/or sliding
constitution layers and applying external fields. Physically, such solutions
of material engineering typically cause the breaking of spatial and temporal
symmetries that govern the dynamics of electrons confined inside material
layers. For example, a magnetic field breaks the time reversal symmetry. It
results in many intriguing phenomena such as the quantum Hall and magneto-
optical effects, including the Kerr and Faraday rotation of the polarization
plane of linearly polarized light as well as the circular dichroism.Tse and
MacDonald (2011); Crassee  _et al._ (2011); Nandkishore and Levitov (2011);
Shimano  _et al._ (2013) Such striking optical phenomena have been exploited
to develop devices in many areas of application.Izake (2007); Liu  _et al._
(2005); Solomon  _et al._ (2019) However, using external fields does not allow
designing small-scale devices. Accordingly, searching materials with intrinsic
properties and solutions to break spatial symmetries of an electronic system
to match technical requirements is paramount.

Twisted-bilayer graphene (TBG) is a typical 2D vdW material that has been
attracting an intensive study from research community.Andrei and MacDonald
(2020); Nimbalkar and Kim (2020); Hill (2020) The TBG system is created by
stacking two graphene layers on each other such that the alignment between two
hexagonal lattices is characterized by a rotation, or twist, angle.
Experimental researches showed the ability to tune the twist angle
mechanically.Ribeiro-Palau  _et al._ (2018); Cai and Yu (2021) For the
electron system inside the TBG lattices the time reversal symmetry is not
broken. The TBG lattices only process rotational symmetries of the group
$D_{3}$ or $D_{6}$.Zou  _et al._ (2019); Do  _et al._ (2020) As a consequence,
the response of the TBG systems to an electromagnetic field is characterized
by a diagonal optical conductivity tensor whose elements are
$\sigma_{xx}(\omega)=\sigma_{yy}(\omega)\neq 0$ and $\sigma_{xy}(\omega)=0$.
However, the TBG systems were pointed out to support the circular dichroism
phenomenon.Kim  _et al._ (2016) It was generally supposed to lie in the fact
that the TBG atomic lattices is chiral, i.e., the left-handed and right-handed
twist configurations are not the mirror images of each other. The optical
activity of crystalline solids can be generally explained as a magneto-
electric effect.Condon (1937); Kim  _et al._ (2016); Nandkishore and Levitov
(2011); Stauber  _et al._ (2018) However, microscopic mechanisms for the
optical activity of diverse materials are not unique. It is thus desirable to
clarify a particular mechanism for novel multiple layer vdW materials in
general and for twisted bilayer graphene in particular.

Microscopically, Morell employed an effective continuum model wherein the
dynamics of electrons is described as the coupling of Dirac fermions in each
graphene layer.Morell  _et al._ (2017) The Dirac fermions are characterized by
a helicity. Twisting two graphene layer makes the non-colinear of the two
helicities in two layers. Combined with a symmetry analysis, they discussed
that the origin of the observed CD spectrum is due to the misalignment of the
helicities. Independently, Do et. al. performed a numerical study on the basis
of a tight-binding model to simulate the evolution of electron states in the
whole TBG lattices.Do  _et al._ (2020) Initial states are created as
excitations on different lattice nodes of graphene layers. They showed a
spreading of the electron wave functions not only in the atomic lattice, but
also the interchange between the two layers. It was shown that special excited
states can similarly occur to the picture of neutrino oscillations.Xian  _et
al._ (2013); Do  _et al._ (2019) Especially, due to the misalignment of two
hexagonal lattices, the wave spreading in one graphene layer results in the
spreading in the other layer with a twist of the wavefront shape. Effectively,
a motion of electron in one layer induces a transverse motion in the other
layer, forming a picture similar to that of the classical Hall effect. The
simulation results of the evolution of Gaussian wave packets show the
trajectory of the wave centroid in curly curves, deflecting from the initial
direction of the wave packets. Numerical calculation shows a non-zero Hall
conductivity, demonstrating the existence of correlation between the
longitudinal and transverse components of the current densities in two
graphene layers. This result seems to violate the symmetry properties of
thermodynamics like the Onsager reciprocal relations.de Groot and Mazur (1984)
However, it is totally not because it is the result only from a particular
state excited from one layer. An opposite behavior will occur if one excites
an initial state on the other graphene layer. Eventually, taking average over
an ensemble of possible initial states will exactly result in a total zero
transverse conductivity $\sigma_{xy}(\omega)$. That simulation is consistent
with optical experimental setups wherein incident photons know which graphene
layer they hit first and which one they hit later. The non-zero Hall
conductivity supported by individual electron states in the TBG lattices
therefore should be understood as the drag Hall conductivity. Accordingly, in
order to correctly study the optical response of multiple layer vdW systems
the spatial dispersion along the direction perpendicular to material layers
must be appropriately taken into account.Agranovich and Ginzburg (2013);
Poshakinskiy  _et al._ (2018)

In this paper we present a theory for analyzing the response of multiple layer
vdW systems to external electromagnetic fields by explicitly taking into the
spatial dispersion. This theory naturally leads to the decomposition of the
conductivity tensor into two parts. The first part, called the local one,
defines the local relation of the current density and the electric field in
each layer. The second part, however, defines the linear relation of the
current density in one layer to the value of the electric field in other
layers. The second part therefore encodes the non-local effects and is called
the drag one. We will show how the elements of these two conductivity parts
are calculated using the Kubo formula. We then apply this theory to the TBG
systems and show how the local and drag conductivities play their role in
governing the chiral response of the material systems. For this, we have to
solve the problem of propagation of an electromagnetic wave through a TBG
sheet without ignoring its thickness though small. It means that we explicitly
consider each graphene layer as an interface of a material slap separating two
dielectric mediums. These interfaces, however, support two current densities
which are mutually influenced to each other, different from ordinary
dielectric slaps. The mutual relation is governed by the interlayer coupling.
The hybridization of electronic states in two graphene layers results in
unique states in layer-coupled systems. Interestingly, we will show that these
hybridized states also support a quantum value of the DC conductivity at the
intrinsic Fermi energy.

This article has three main sections. Aside from the introduction, in the next
section, Sec. II, we present the theory for decomposing the total conductivity
tensor into two parts and show how they are calculated via the Kubo formula.
Subsequently, we present a solution to an effective model for electrons in the
TBG lattices. In Sec. III, we present numerical results for the electronic
band structure of several TBG configurations, results for the analysis of the
optical conductivity tensor elements, and the spectra of circular dichroism
with respect to photon frequency. Finally, the conclusion is given in section
IV.

##  II Theory and calculation methods

###  II.1 AC and DC conductivities

The optical properties of bilayer graphene systems are noticeable.Moon and
Koshino (2013); Stauber  _et al._ (2018); Le and Do (2018); Zhang (2019) In
particular, when the layer resolution should be taken into account, effects
involved in the spatial dispersion would emerge. To describe such effects, we
decompose the current responding to the action of an external electric field
of a monochromatic light beam incident perpendicular to the material plane 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).$
(1)

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}).$
(2)

Eq. (1) 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).$
(3)

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),$
(4)

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. (4) 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},$ (11)

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 (11) is identical to results reported in Refs.
Stauber  _et al._ , 2018; 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}).$ (12)

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{i\hbar
e^{2}}{S}\sum_{m,n}\frac{f_{n}-f_{m}}{E_{m}-E_{n}}\frac{O^{(c)mn}_{\alpha\beta}}{\hbar(\omega+i\eta)-(E_{m}-E_{n})},$
(13)

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,$ (14a) $\displaystyle
O_{\alpha\beta}^{(drag)mn}$ $\displaystyle=\langle
m|\hat{v}_{\alpha}^{(1)}|n\rangle\langle n|\hat{v}_{\beta}^{(2)}|m\rangle.$
(14b)

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),$ (15)

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}).$ (16)

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),$
(17)

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

###  II.2 Effective continuum model

Developing effective continuum models for low-energy states of electrons in
TBG systems has been proceeded since 2007 by Lopes et al.dos Santos  _et al._
(2007) However, the model driven by Bistritzer and MacDonald in 2011 was well
known and commonly used. Here, we present our solution to the Bistritzer-
MacDonald model.Bistritzer and MacDonald (2011) We employ the implementation
for the 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.\boldsymbol{e}_{x}$ and
$\boldsymbol{a}_{2}=a\cos\pi/3.\boldsymbol{e}_{x}+a\sin\pi/3.\boldsymbol{e}_{y}$,
where $a$ is the lattice constant of graphene, and
$\boldsymbol{e}_{x},\boldsymbol{e}_{y}$ the unit vectors defining the $x,y$
directions, respectively. The two associated vectors defining the reciprocal
lattice are thus
$\boldsymbol{a}_{1}^{\star}=(4\pi/a)(\cos\pi/6.\boldsymbol{e}_{x}+\sin\pi/6.\boldsymbol{e}_{y})$
and $\boldsymbol{a}_{2}^{\star}=(4\pi/a)\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,\xi=\pm
1$.

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

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

When $|m-n|=1$, the moire reciprocal lattice is defined by two unit vectors
$\boldsymbol{A}_{i}^{\star}=\boldsymbol{a}_{i}^{\star(1)}-\boldsymbol{a}_{i}^{\star(2)}$.
The mini-Brillouin zone is defined by six $\boldsymbol{K}^{M}_{i}$ points
given by:
$\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$.
The other points are found by shifting these points via the reciprocal unit
vectors.

Given a Hamiltonian $H$ for electrons, we first use the layer-resolution
vector basis set $\\{|\ell\rangle\,|\,\ell=1,2\\}$ to represent it. Using the
identity $1=\sum_{\ell=1}^{2}|\ell\rangle\langle\ell|$ we have:

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

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

$\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|,$
(20)

where we denote
$|\ell,\boldsymbol{k},m\rangle=|\ell\rangle|\boldsymbol{k},m\rangle$ and 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}}$.
Accordingly, the low energy states of electrons in the TBG configurations of
tiny twist angles are distinguished by a quantum index $\xi=\pm 1$ that
corresponds to the two nonequivalent Dirac valleys
$\boldsymbol{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 $\boldsymbol{K}_{\xi}^{(1)}$ and $\boldsymbol{K}_{\xi}^{(2)}$ are close
to each other. The Bistritzer-MacDonald model is given by a Hamiltonian in the
real-space representation as follows:Bistritzer and MacDonald (2011)

$\displaystyle\hat{H}^{\xi}=\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).$
(23)

Here $\hat{\boldsymbol{p}}$ and $\hat{\boldsymbol{r}}$ are the momentum and
position operators, respectively. The Hamiltonians
$H_{1}^{\xi}(\hat{\boldsymbol{p}})$ and $H_{2}^{\xi}(\hat{\boldsymbol{p}})$ of
the two uncoupled graphene layers are given 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],$
(24)

where $\theta_{\ell}=(-1)^{\ell}\theta/2$ is the rotation angle of layer
$\ell$ with respect to the fixed Cartesian coordinate axes $Oxyz$;
$R^{z}\left(-\theta_{\ell}\right)$ is a matrix to rotate back the relevant
vectors around the $Oz$ axis to keep the canonical form of the Dirac
Hamiltonian; $\boldsymbol{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{\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}}}$
(29) $\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}}},$
(32)

where $\omega=\exp(i2\pi/3)$ and
$\delta\boldsymbol{k}_{2}=\boldsymbol{A}_{1}^{\star}$ and
$\delta\boldsymbol{k}_{3}=\boldsymbol{A}_{1}^{\star}+\boldsymbol{A}_{2}^{\star}$.
We choose the value of two parameters $u=0.0797$ eV and $u^{\prime}=0.0975$ eV
to take into account of effects of lattice reconstruction.Koshino  _et al._
(2018)

Associated with the Hamiltonian (24), the velocity operator is determined by:

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

where $R^{z}_{\alpha}(-\theta_{\ell})$ denotes the first/second column of the
rotation matrix $R^{z}(-\theta_{\ell})$ if $\alpha=x/y$. It should be noticed
that because of the relativistic form of
$H^{\xi}_{\ell}(\hat{\boldsymbol{p}})$ given by Eq. (19), the current operator
$\hat{j}_{\alpha}$ has only one contribution part of drift current,
$\hat{j}_{\alpha}=e\hat{v}_{\alpha}$, i.e., no diffusion and diamagnetic
parts.

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

$\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.$
(36)

Due to the “approximate” periodicity of the moire lattice of TBG lattices,
electron states are determined as the Bloch state vectors
$|\psi_{\xi,\boldsymbol{k}}\rangle$ that are expanded in terms of plane-wave
vectors as follows:

$|\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,$
(37)

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

$\displaystyle\sum_{m}\langle\boldsymbol{k}+\boldsymbol{G}_{n}|\hat{H}^{\xi}|\boldsymbol{k}+\boldsymbol{G}_{m}\rangle\left(\begin{array}[]{c}C_{1,\xi,\boldsymbol{k}}(\boldsymbol{G}_{m})\\\
C_{2,\xi,\boldsymbol{k}}(\boldsymbol{G}_{m})\end{array}\right)$ (40)
$\displaystyle\hskip
91.04872pt=E\left(\begin{array}[]{c}C_{1,\xi,\boldsymbol{k}}(\boldsymbol{G}_{n})\\\
C_{2,\xi,\boldsymbol{k}}(\boldsymbol{G}_{n}).\end{array}\right)$ (43)

With the notice that

$\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}},$
(44a)
$\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}},$
(44b)
$\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}},$
(44c)

we specify Eq. (40) 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.$
(47) $\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).$ (56)

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

$\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),$
(57c)
$\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},$
(57f)
$\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}.$
(57i)

Diagonalize the matrix $H_{\boldsymbol{k}}^{\xi}$, we obtain all possible
eigen-values of $E_{n}^{\xi}(\boldsymbol{k})$. From these data, we can
represent the electronic energy band structure. The issue here is the value of
$N_{\boldsymbol{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}$, then determine $N_{\boldsymbol{G}}$ the number
of $\boldsymbol{G}_{n}$ vectors such that $\|\boldsymbol{G}_{n}\|\leq
E_{c}/\hbar v_{F}$.

Figure 1:  These present the energy band structures of electron in three TBG
configurations with $\theta=3.890^{\circ}$ (a,d), $1.890^{\circ}$ (b,e), and
$1.050^{\circ}$ (c,f). The blue curves in (a,b,c) are for the dispersion and
density of states of graphene, which are obtained by turning off the
interlayer coupling. The blue and red curves in (d,e,f) are the dispersion
curves obtained for the valley $\xi=+1$ and $\xi=-1$.

##  III Results and discussions

###  III.1 Electronic band structure of TBGs

We employed the numerical method to solve Eq. (47). Specifically, for each
value of $\xi=\pm 1$, a Hermitian matrix $H_{\boldsymbol{k}}^{\xi}$ of size
$4N_{G}\times 4N_{G}$ was constructed and then diagonalized for each value of
$\boldsymbol{k}$ given in the mini-Brillouin zone of TBG systems. The sets of
eigen-values $\\{E_{m}(\boldsymbol{k})\\}$ and of eigen-vectors
$\\{\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_{G}}),C^{m}_{2,\xi,\boldsymbol{k}}(\boldsymbol{G}_{N_{G}})]^{T}\\}$
were determined. Obtained results for three TBG configurations with the twist
angles of $\theta=1.05^{\circ},1.89^{\circ}$ and $3.89^{\circ}$ were presented
in Fig. 1. The blue curves are the results obtained as we artificially turn
off the electronic interlayer coupling. We plot both data on the same figure
for the aim of comparison, i.e., in order to highlight the effects of the
interlayer coupling. The appearance of many dispersion curves in the figure is
the result of the folding of the energy band structure of two graphene layers
due to the enlargement of the unit cell of the TBG lattices, hence the shrink
of the Brillouin zone into the mini-Brillouin zone. However, the electronic
band structure of the TBG systems is not simply formed like that. At the
points of crossing bands of monolayer graphene, the interlayer coupling leads
to the hybridization of Bloch states; the energy degeneracy is thus lifted to
form the new bands. Our obtained results are in agreement 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
work, we would like to emphasize the formation of new electronic states which
are unique in the TBG systems because they govern the electronic, optical and
transport properties. Fig. 1 can be seen as the evolution picture of the band
structure with respect to the twist angle. We realize that in the low energy
range around the intrinsic Fermi energy ($E_{F}=0$), the dispersion curves of
the configuration with $\theta=3.89^{\circ}$ qualitatively have the similar
form of the blue curves although they are not identical. This implies that the
created states in the bilayer system much share the essence of states in
monolayer graphene. For smaller twist angles, in particular the one with
$\theta=1.05^{\circ}$, a band around the Fermi energy with a very narrow
bandwidth is formed and separated from the lower and upper bands by a narrow
gap. These dispersion curves are totally different from any blue curves, see
Fig. 1(e). We therefore conclude that the electronic interlayer coupling
creates unique states in the TBG lattices of tiny twist angles. Although the
TBG systems are two-dimensional, geometrical features of the energy surfaces
in the wave vector space are difficult to be directly analyzed. However,
general features should exhibit in the picture of density of states (DOS) as
the van Hove singularity behaviors shown in the right panels of Figs.
1(a,b,c). In these figures, the blue curves are numerically obtained without
taking into account of the electronic interlayer coupling, and thus the DOS of
two independent graphene layers. The red curves show significant peaks due to
the van Hove singularities, reflecting the existence of extremal points and
saddle points of the energy surfaces. Another important point in the TBG
electronic structure is that the energy bands with the index $\xi+1$ (resulted
from the hybridization of graphene Bloch states with the wave vector in the
$K_{1}$ valley of the first layer and those in the $K_{2}$ valley of the
second layer) is different from those with the index $\xi=-1$ (resulted from
hybridization of Bloch states in the $K_{1}^{\prime}$ valley and in the
$K_{2}^{\prime}$). The difference clearly manifests along the ${\it\Gamma}M$
direction as shown in Figs. 1(d,e,f). However, it does not manifest in the
picture of the density of states. We computed the DOS for $\xi=+1$ and
$\xi=-1$ separately and obtained the identical results. This point implies
that the energy surfaces with $\xi=+1$ are simply the rotation of the
corresponding ones with $\xi=-1$ by an angle of $\pi/3$.

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

###  III.2 DC conductivity

The calculation of the electronic structures of TBG configurations with not
too small twist angles ($>2^{\circ}$) shows the appearance of two energy bands
around the Fermi energy $E_{F}=0$ with the dispersion curves very similar to
those of monolayer graphene, except the value of the Fermi velocity
characterizing the linear dispersion law $v_{F}^{TBG}<v_{F}^{MLG}$. Meanwhile,
if continuing to decrease the twist angle, these two bands shall become closer
to each other and form an isolated energy band with a very narrow bandwidth.
Such a “flat-band” at the Fermi energy was supposed not to contribute to the
transport properties of TBG sheets, because of very small value of the Fermi
velocity of carriers. In other words, electrons occupying states corresponding
to these bands should be spatially confined in the atomic lattice. Actually,
these bands are not totally dispersionless. Hence, it is really intriguing to
investigate the transport properties of the TBG systems. In the framework of
non-interacting electrons approximation, we calculated the DC conductivity
using the Kubo-Greenwood formula given by Eqs. (15) and (16). The results we
have obtained for the longitudinal components $\sigma_{xx}^{DC}$ and
$\sigma_{yy}^{DC}$ of three TBG configurations under study in the regime of
low temperature are presented in Fig. 2.

The data show that the ability of electric conducting of the TBG systems is
anisotropic. This is quantified by a conductivity tensor with two different
longitudinal elements $\sigma_{xx}^{DC}$ (the solid curves) and
$\sigma_{yy}^{DC}$ (the dashed curves). In principle, the electronic structure
must exhibit in the picture of the conductivity with respect to the Fermi
energy. We calculated these two quantities for the two cases of $\xi=\pm 1$
and obtained the same value. This point is similar to the DOS picture.
However, the difference of two longitudinal conductivities is the reflection
of the anisotropic nature of the energy surfaces. For the TBG configuration
with $\theta=3.890^{\circ}$, the electronic structure is qualitatively similar
to that of monolayer graphene (except the value of the Fermi velocity),
therefore the conductivity curves have the typical V-shape (or the M-shape if
considered in a larger energy range, see the blue curves). For the
configurations with tiny twist angles, the picture of the electronic
structure, and thus the conductivity picture, is much more complicated.
However, our calculation shows that the conductivity at the intrinsic Fermi
level $E_{F}=0$ always takes a finite non-zero value. Interestingly, this
value, for the TBG configurations whose twist angles are not the magic ones,
is about two times of $\sigma^{DC}_{0}=4e^{2}/\pi h$, a quantum value of the
minimal conductivity of monolayer graphene. For the TBG configuration of the
magic twist angle $\theta=1.050^{\circ}$, the conductivity curve exhibits a
small hump at $E_{F}=0$ with
$\sigma_{\alpha\alpha}^{DC}(0)\approx\sigma^{DC}_{0}$. The quantum value of
graphene conductivity was a topic intensively discussed around 2006.Novoselov
_et al._ (2005); Katsnelson (2006); Tworzydlo  _et al._ (2006); Miao  _et al._
(2007); Blake  _et al._ ; Do and Dollfus (2010) Its origin lies in the fact
that the valence and conduction bands are not separated, but elegantly touch
each other at the $K$ points, so there are always available free carriers in
graphene due to the inevitable fluctuations of the electrostatic profile
around the Fermi energy level. Numerically, in our calculation, we
approximated the delta-Dirac functions in the Kubo-Greenwood formula by a
Gaussian function with the parameter $\eta$ that defines the finite width of
the peak. The value of this parameter was taken to be $\eta<5$ meV. The use of
$\eta$ is due to not only the reason for numerical calculation, but also the
fact of the broadening of energy levels, and thus the finite lifetime of
quasi-particles carrying the electric charge. The quantum value of the DC
conductivity observed here for the TBG systems should have the same reason.
However, it should be emphasized that it is supported by typical hybridized
states in the TBG lattices.

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

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

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

In Fig. 3 we present the calculated results for the total optical conductivity
of three TBG configurations. The calculation was performed on the use of the
Kubo formula for the conductivity tensor wherein the velocity operator
$\hat{v}_{\alpha}$ is determined from the Hamiltonian (23) using Eq. (12). In
principle, the time reversal symmetry as well as the Onsager reciprocal
relations rules out the off-diagonal elements of the conductivity tensor. We
numerically checked this point and observed that the data of
$\sigma_{xy}(\omega)$ behaves as a noise picture with the amplitude of
$10^{-3}$, correctly verified the general principle. Additionally, all the
symmetry properties of the conductivity tensor according to Eq. (11) are
verified. It thus emphasizes the correctness of the numerically obtained data.
Fig. 3 shows the longitudinal conductivity curves $\sigma_{xx}(\omega)$ with
significant peaks, which are assigned to the dominant contribution of
interband transition processes. This analysis is supported by associating it
with the electronic band structures of the TBG configurations shown in Fig. 1.
For example, in Fig. 1(b), at the energy position $E=0.5681$ eV there appears
an optical absorption peak that can be assigned to the transition from the
highest valence band to the lowest conduction band at the value of
$\boldsymbol{k}$ around the $M$ points in the mini-Brillouin zone. When
decreasing the twist angle value, the picture of the electronic structure
becomes more complicated. As a consequence, the optical conductivity shows
more optical absorption peaks. For the TBG configurations with not too small
twist angles ($>2^{\circ}$), there is a noticeable feature of the conductivity
curve in the low energy range: the real part of the conductivity behaves
independently of the photon energy $\hbar\omega$ and takes the value of
$2\sigma_{0}$ where $\sigma_{0}=\pi e^{2}/2h$. This behaviour is similar to
that of monolayer graphene, which is pointed out to be a consequence of the
linear dispersion law of electron states. In the low energy range, the energy
band structure of the TBG systems qualitatively resembles that of monolayer
graphene,Nicol and Carbotte (2008) but with a smaller Fermi velocity as
analyzed previously. The value of $\sigma_{xx}(\omega)$ is, however, equal to
$2\sigma_{0}$.Le and Do (2018); Do  _et al._ (2020); Moon and Koshino (2013)
For the case of magic angle $\theta=1.050^{\circ}$, our calculation shows
clearly the dominant contribution of the optical transitions from the quasi-
flat bands around the Fermi level to the lowest conduction band around the
$\Gamma$ point, see Fig. 3(c) together with Figs. 1(c,f).

In Fig. 4 and Fig. 5, we present the calculation results for the elements of
the blocks $\boldsymbol{\sigma}^{(\ell)}(\omega)$ and
$\boldsymbol{\sigma}^{(drag)}(\omega)$ of the conductivity tensor according to
the analysis given in Sec. II.A. Our numerical calculation obviously shows
$\boldsymbol{\sigma}^{(1)}(\omega)=\boldsymbol{\sigma}^{(2)}(\omega)=\sigma_{xx}^{(\ell)}(\omega)\boldsymbol{\tau}_{0}$.
This result verifies the equivalent role of two graphene layers in the TBG
lattice. Figs. 4(a,b,c) shows the real part (red) and the imaginary part
(blue) of the element $\sigma_{xx}^{(\ell)}(\omega)$ where $\ell=1,2$. One
noticeable point is that when decreasing the value of the twist angle, the
magnitude of $\sigma_{xx}^{(\ell)}$ increases; especially, the imaginary part
becomes dominant and varying via the law of $1/\hbar\omega$, see Fig. 4(c).
Next, 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 the TBG configurations of large twist angle, there is a
significant difference between the curves $\sigma_{xx}^{(\ell)}(\omega)$ and
$\sigma_{xx}^{(drag)}(\omega)$, see Fig. 4(a) and Fig. 5(a). However, for the
configurations with tiny twist angles ($<2^{\circ}$), these curves are mostly
the same but they are opposite in sign. Accordingly, though there is a large
variation of the magnitude of $\sigma_{xx}^{(\ell)}(\omega)$ and
$\sigma_{xx}^{(drag)}(\omega)$, they compensate each other in the summation to
result in the total value
$\sigma_{xx}(\omega)=2[\sigma_{xx}^{(\ell)}(\omega)+\sigma_{xx}^{(drag)}(\omega)]$
as shown in Fig. 2. For the behaviour of $\sigma_{xy}^{(drag)}(\omega)$, it
mainly shows the zero-close value in almost energy range except for some sharp
peaks and dips at the energy positions of strong optical absorption.

The appearance of the transverse conductivity $\sigma_{xy}^{(drag)}(\omega)$
is strikingly important. This is because it governs typical optical properties
of the TBG systems. This conductivity does not result from only the
correlation of current densities distributed on two graphene layers, but also
the chirality of the atomic lattices. The consequence of the latter is that
the TBG atomic lattices have no mirror symmetries, totally different from the
hexagonal lattice of graphene. The studies made by Morell and by Do show that
the hybridized states formed in the bilayer systems always support a finite
value of the correlation of two transverse components of electron motion. In
other words, a motion of electron in one graphene layer along the $Ox$
direction always induces another motion in the second graphene layer along the
transverse direction $Oy$ due to the interlayer coupling. However, the
contribution of such correlations to the total value of the conductivity must
cancel each other completely if the system has a mirror symmetry. Since the
TBG lattices are chiral, such cancellation should be incomplete and thus
results in the finite non-zero conductivity $\sigma_{xy}^{(drag)}(\omega)\neq
0$. In order to emphasize this point, we divided the hexagonal mini-Brillouin
zone of the TBG lattice into three equal parts as shown in the inset of Fig.
6. The division is based on the assumption that there is a mirror symmetry via
the plane $M_{yz}$. Geometrically, the domain (1) would contain the symmetry
$M_{yz}$ meanwhile the domains (2) and (3) would be the mirror images of each
other. Using the Kubo formula we calculated the value of
$\sigma_{xy}^{(drag)}$ by taking the summation over all the $\boldsymbol{k}$
points in these three domains separately. Obtained results are shown in Fig.
6. We clearly see the results obtained in the domains (2) and (3) are very
similar to each other, but opposite in sign, and much larger than that
obtained in the domain (1). Obviously, when taking into account the summation
of the three results, a strong, but not complete, cancellation of the data
occurs and results in the final value of $\sigma_{xy}^{(drag)}(\omega)$ as
shown in Figs. 5(b,d,f). For the two special AA- and AB-stacked
configurations, which processes the mirror symmetries, we verified that the
cancellation takes place completely, and thus
$\sigma_{xy}^{(drag)}(\omega)=0$. We call this result the chiral correlation
of drag current densities distributed in the two graphene layers. The
essential role of the chirality of the TBG atomic lattices is therefore
demonstrated. However, we should emphasize that its essence lies in governing
symmetrical properties of electronic states. This can be seen as the
sufficient condition for $\sigma_{xy}^{(drag)}(\omega)\neq 0$. Meanwhile, the
necessary condition must be the electronic interlayer coupling, which results
in the formation of electronic states that always support the drag
correlation. The Kubo formula allows to calculate each elements of the
conductivity tensor independently. Technically, this can be seen as the result
of decomposing the total velocity operator into appropriate terms, similar to
Morell’s strategy. However, it should be kept in mind that the calculation of
the velocity-velocity correlation is actually realized on the basis of the
electronic states of the TBG systems, not of Bloch states of individual
graphene.

Figure 5:  These present the real and imaginary parts of the drag term
$\boldsymbol{\sigma}^{(drag)}(\omega)$ of the conductivity tensor of three TBG
configurations with $\theta=3.890^{\circ},1.890^{\circ}$ and $1.050^{\circ}$:
(a,c,e) for the longitudinal (diagonal) element
$\sigma_{xx}^{(drag)}(\omega)$, and (b,d,f) for the transverse (off-diagonal)
element $\sigma_{xy}^{(drag)}(\omega)$. Figure 6:  This presents the
contribution of Bloch states to $\sigma_{xy}^{(drag)}(\omega)$: The red and
blue curves are from the Bloch states with $\boldsymbol{k}$ in the first one-
third Brillouin zone; the solid/dashed purpe and green curves are from the
Bloch states with $\boldsymbol{k}$ in the second/third one-third Brillouin
zone. Figure 7:  This presents the transmission spectra of three TBG
configurations: $\theta=3.890^{\circ}$ (a), $\theta=1.890^{\circ}$ (b), and
$\theta=1.050^{\circ}$ (c). Figure 8:  This presents the circular dichroism
(CD) spectra of three TBG configurations: $\theta=3.890^{\circ}$ (a),
$\theta=1.890^{\circ}$ (b), and $\theta=1.050^{\circ}$ (c).

###  III.4 Optical chiral responses

In order to investigate the optical response of the TBG systems, in this
section, we first establish a theory for the problem of electromagnetic wave
propagation through a generic TBG sheet whose thickness $d_{GG}$ is not
ignored. Different from an ordinary dielectric slap with two independent
surfaces, the electric current densities in two graphene layers are mutually
related in the meaning that: a current density in one graphene layer is
induced by the value of electric field not only in that layer, but also in the
other layer. This relation is given by Eq. (4). We used this equation to solve
the problem defined here.

Specifically, we set up the system like the one in Ref. Széchenyi  _et al._ ,
2016 in which two graphene layers separates the whole space into three parts,
which are characterized by the parameters $(\epsilon_{i},\mu_{i})$, where
$i=1,2,3$ is the index labeling the mediums. Our aim is to establish the
transfer matrix of the whose system in order to build the expression for the
transmission and reflection matrices for the system of bilayer graphene. Since
the whole space is only inhomogeneous in the direction perpendicular to the
TBG surface, called the $z$ direction, the Maxwell equations support a plane
wave solution that is given by the expression:

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

where $\boldsymbol{F}$ means $\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. Conditions required 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,$
(59a)
$\displaystyle\boldsymbol{n}\times[\boldsymbol{H}_{2}(1)-\boldsymbol{H}_{1}(1)]=\boldsymbol{J}(1),$
(59b)

for the surface of the first graphene layer, and

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

for the surface of the second graphene layer. In Eqs. (59b) and (60b), the
distribution of current densities on each graphene layer is linearly related
to the electric field by the conductivity tensor given via Eq. (11). 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),$
(61a) $\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).$
(61b)

The equations of continuity conditions allow us to establish the relation
between the vector fields $(\boldsymbol{F}_{i}^{+},\boldsymbol{F}_{i}^{-})$ in
two mediums at their interface. Particularly, we first notice that the
magnetic field $\boldsymbol{H}_{i}$ is related to the electric field
$\boldsymbol{E}_{i}$ via the Maxwell equation
$\text{curl}(\boldsymbol{E})=-\partial_{t}\boldsymbol{B}$ wherein
$\boldsymbol{B}=\mu\boldsymbol{H}$. Using the plane wave solution (58) we
obtain:

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

where $\epsilon_{i}$ and $\mu_{i}$ are the absolute permittivity and the
absolute permeability of medium $i$. In order to obtain the above equation,
the dispersion relation $\omega=c\|\boldsymbol{k}_{i}\|/n_{i}$ is already
used, where $c$ is the light speed in vacuum and $n_{i}$ is the refractive
index of medium $i$. Notice that the value of $\sqrt{\epsilon_{i}/\mu_{i}}$
can be expressed through $n_{i}$, i.e., $\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. We solve the problem for the transmission of light along
the normal direction to the TBG layer, so
$\boldsymbol{k}_{i}=k_{i}\boldsymbol{n}$. It results in the transfer matrix
$\mathbf{M}_{31}$ (of size $4\times 4$) relating the components of the
electric fields in the first and third mediums

$\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).$ (67)

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

$\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),$
(68e) $\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),$
(68j) $\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),$ (68m)

where $\boldsymbol{1}$ is the $2\times 2$ identity matrix. Let $\mathbf{r}$
and $\mathbf{t}$ be the $2\times 2$ matrices relating 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 called the reflection and transmission matrices, respectively. From Eq.
(67) we deduce:

$\displaystyle\mathbf{r}$
$\displaystyle=-[\mathbf{M}_{31}^{22}]^{-1}\mathbf{M}_{31}^{21},$ (69a)
$\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}.$
(69b)

From these matrices, the physical quantities as the reflectance ($R$),
transmittance ($T$) are determined as the ratio of the energy fluxes.
Concretely,

$\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|,$
(70a) $\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|.$
(70b)

From the conservation of the energy flux, the absorptance is determined by
$A=1-(R+T)$. The quantities $R,T$ and $A$ depend on not only the energy but
also the polarization states of light. Remember that the left/right-handed
polarization states of light are characterized by the electric field vectors
$\boldsymbol{E}_{L,R}^{+}\propto(1,\pm i)^{T}/\sqrt{2}$. We calculated the
spectrum of circular dichroism (CD), a quantity is defined to measure the
dependence of the light absorption of the left-handed (L) and right-handed (R)
polarization by the formula:

$\mathrm{CD}=\frac{A_{\mathrm{L}}-A_{\mathrm{R}}}{A_{\mathrm{L}}+A_{\mathrm{R}}},$
(71)

On the basis of the theory established above 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. Obtained results are presented in Fig. 7 for the transmittance
and in Fig. 8 for the CD. The reflectance is very small, so we ignore to show
it here. From the transmission spectra, 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 Hall 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 Hall conductivity plays the decisive role in governing
the chiral optical behaviors of the TBG systems.

##  IV Conclusion

Engineering two-dimensional quantum materials to electronic and opto-
electronic applications is the research direction that has been extensively
developing. In this work, we have presented a study of the mechanism of
optical activity of twisted bilayer graphene, a typical two-dimensional van
der Waals material. The important feature of this material system is that its
atomic lattice processes only rotational symmetries, but not the mirror
symmetry due to the chiral structure. We used an effective continuum model to
show the hybridization of electron states in each graphene layer to form
unique electron states in the bilayer system. Such hybridized states support
the correlation of transverse motions in two graphene layers. We established a
solution to the problem of electromagnetic wave propagation through a TBG
sheet wherein each graphene layer is considered as the conducting interface
separating two dielectric mediums. The two interfaces are not independent, but
they conduct the mutually influenced current densities. By this, we considered
explicitly the spatial dispersion and presented a theory analyzing the
conductivity tensor into two parts, the local and the drag ones. We calculated
all the elements of these parts using the Kubo formula. We then formulated and
numerically calculated the transfer matrix and demonstrated that the
transverse drag Hall conductivity $\sigma_{xy}^{(drag)}$ plays the decisive
role in defining the chiral optical response of the TBG systems. We showed the
existence of $\sigma_{xy}^{(drag)}$ due to the chiral structure of the atomic
lattice. This results in the hybridized states of lacking the mirror symmetry.
Such states not only support the correlation of transverse motions in two
graphene layers, but also leads to the incomplete cancellation of such drag
correlation to form this drag Hall optical conductivity. Another important
result in this work is what we calculated the DC conductivity of the TBG
systems and showed the existence of a quantum conductivity value $\propto
e^{2}/h$ at the intrinsic Fermi energy. This quantum conductivity value is
governed by the hybridized states in the bilayer systems. Our study therefore
emphasises that in order to correctly analyze the optical response and
transport properties, the finite thickness of the twisted bilayer graphene
must be taken into account because of the importance of spatial dispersion
effects. The theoretical approach we developed here is suitable to utilize for
other van der Waals material systems of multiple layers.

## Acknowledgements

The authors are grateful to Shingo Takeuchi for reading carefully the
manuscript before it was submitted. One of the authors, S.T.H., acknowledges
Hanoi University of Civil Engineering (HUCE) for financially funding his work
under the grant numbered 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).
  * 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).
  * 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).
  * 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).
  * Condon (1937) E. U. Condon, Rev. Mod. Phys. 9, 432 (1937).
  * Stauber  _et al._ (2018) T. Stauber, T. Low, and G. Gómez-Santos, Phys. Rev. Lett. 120, 046801 (2018).
  * Morell  _et al._ (2017) E. S. Morell, L. Chico, and L. Brey, 2D Mater. 4, 035015 (2017).
  * Xian  _et al._ (2013) L. Xian, Z. F. Wang, and M. Y. Chou, Nano Lett. 13, 5159 (2013).
  * Do  _et al._ (2019) V. N. Do, H. A. Le, and D. Bercioux, Phys. Rev. B 99, 165127 (2019).
  * de Groot and Mazur (1984) S. R. de Groot and P. Mazur,  _Non-equilibirum thermodynamics_ (Dover Publications, 1984).
  * 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).
  * Moon and Koshino (2013) P. Moon and M. Koshino, Phys. Rev. B 87, 205404 (2013).
  * 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).
  * (44) 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).
  * Széchenyi  _et al._ (2016) G. Széchenyi, M. Vigh, A. Kormányos, and J. Cserti, J. Phys.: Condens. Matter 28, 375802 (2016).


$\displaystyle\mathbb{P}\Big{(}\|Y^{\ell}-\hat{Y}^{k,\ell}\|_{\Lambda_{T}}>(1+T)^{1/2}\Big{(}\rho_{b}(\varepsilon_{m_{k}}^{2})^{1/2}+R_{m_{k}}\frac{\rho_{\sigma}(\varepsilon_{m_{k}}^{2})}{\varepsilon_{m_{k}}}\Big{)}^{1/4}\varepsilon^{1/2}_{m_{k}}\Big{)}$
$\displaystyle\leq\Big{\\{}648+C_{0}e^{C_{0}T}(1+\|y\|^{2}_{\mathcal{H}_{\mu}})\Big{\\}}\Big{(}\rho_{b}(\varepsilon_{m_{k}}^{2})^{1/2}+R_{m_{k}}\frac{\rho_{\sigma}(\varepsilon_{m_{k}}^{2})}{\varepsilon_{m_{k}}}\Big{)}^{1/4}+\mbox{\rm{1}}\mbox{\rm{l}}_{\\{J_{k}<T\\}}.$

Recall the well-known fact that the Lévy–Prokhorov metric between two
probability measures $\nu_{1}$ and $\nu_{2}$ on a Polish space $(E,d)$ is
represented in terms of couplings between $\nu_{1}$ and $\nu_{2}$:

$d_{\mathrm{LP}}(\nu_{1},\nu_{2})=\inf\left\\{\varepsilon>0\mathrel{}\middle|\mathrel{}\inf_{\pi\in\mathscr{C}(\nu_{1},\nu_{2})}\int_{E\times
E}\mbox{\rm{1}}\mbox{\rm{l}}_{\\{d(x_{1},x_{2})>\varepsilon\\}}\,\mathrm{d}\pi(x_{1},x_{2})<\varepsilon\right\\},$

where $\mathscr{C}(\nu_{1},\nu_{2})$ denotes the set of probability measures
on $E\times E$ with $\nu_{1}$ and $\nu_{2}$ as respective first and second
marginals (cf. [12]). By the above observations, we get

$\displaystyle
d_{\mathrm{LP}}\big{(}\mathrm{Law}_{\mathbb{P}}(Y^{\ell}|_{[0,T]}),\mathrm{Law}_{\mathbb{P}}(\hat{Y}^{k,\ell}|_{[0,T]})\big{)}$
$\displaystyle\leq\inf\left\\{\varepsilon>0\mathrel{}\middle|\mathrel{}\mathbb{P}\Big{(}\|Y^{\ell}-\hat{Y}^{k,\ell}\|_{\Lambda_{T}}>\varepsilon\Big{)}<\varepsilon\right\\}$
$\displaystyle\leq\Big{\\{}648+C_{0}e^{C_{0}T}(1+\|y\|^{2}_{\mathcal{H}_{\mu}})\Big{\\}}\Big{(}\rho_{b}(\varepsilon_{m_{k}}^{2})^{1/2}+R_{m_{k}}\frac{\rho_{\sigma}(\varepsilon_{m_{k}}^{2})}{\varepsilon_{m_{k}}}\Big{)}^{1/4}+\mbox{\rm{1}}\mbox{\rm{l}}_{\\{J_{k}<T\\}}.$
(4.19)

Observe that

$\displaystyle
d_{\mathrm{LP}}\big{(}\mathrm{Law}_{\mathbb{P}}(Y^{\ell}|_{[0,T]}),\mathrm{Law}_{\mathbb{P}}(Y^{k}|_{[0,T]})\big{)}$
$\displaystyle\leq
d_{\mathrm{LP}}\big{(}\mathrm{Law}_{\mathbb{P}}(Y^{\ell}|_{[0,T]}),\mathrm{Law}_{\mathbb{P}}(\hat{Y}^{k,\ell}|_{[0,T]})\big{)}+d_{\mathrm{LP}}(\mathrm{Law}_{\mathbb{P}}(\hat{Y}^{k,\ell}|_{[0,T]}),\mathrm{Law}_{\mathbb{P}}(Y^{k}|_{[0,T]})\big{)}$
$\displaystyle\leq
d_{\mathrm{LP}}\big{(}\mathrm{Law}_{\mathbb{P}}(Y^{\ell}|_{[0,T]}),\mathrm{Law}_{\mathbb{P}}(\hat{Y}^{k,\ell}|_{[0,T]})\big{)}+d_{\mathrm{TV}}(\mathrm{Law}_{\mathbb{P}}(\hat{Y}^{k,\ell}|_{[0,T]}),\mathrm{Law}_{\mathbb{P}}(Y^{k}|_{[0,T]})\big{)}$
$\displaystyle\leq
d_{\mathrm{LP}}\big{(}\mathrm{Law}_{\mathbb{P}}(Y^{\ell}|_{[0,T]}),\mathrm{Law}_{\mathbb{P}}(\hat{Y}^{k,\ell}|_{[0,T]})\big{)}+d_{\mathrm{TV}}(\mathrm{Law}_{\mathbb{P}}(\hat{Y}^{k,\ell}),\mathrm{Law}_{\mathbb{P}}(Y^{k})\big{)},$

where the second inequality follows from the well-known inequality
$d_{\mathrm{LP}}(\nu_{1},\nu_{2})\leq d_{\mathrm{TV}}(\nu_{1},\nu_{2})$ for
any probability measures $\nu_{1}$ and $\nu_{2}$ on a Polish space. Therefore,
by (4.19) and (4.17), we obtain

$\displaystyle
d_{\mathrm{LP}}\big{(}\mathrm{Law}_{\mathbb{P}}(Y^{\ell}|_{[0,T]}),\mathrm{Law}_{\mathbb{P}}(Y^{k}|_{[0,T]})\big{)}$
$\displaystyle\leq\Big{\\{}648+C_{0}e^{C_{0}T}(1+\|y\|^{2}_{\mathcal{H}_{\mu}})\Big{\\}}\Big{(}\rho_{b}(\varepsilon_{m_{k}}^{2})^{1/2}+R_{m_{k}}\frac{\rho_{\sigma}(\varepsilon_{m_{k}}^{2})}{\varepsilon_{m_{k}}}\Big{)}^{1/4}+\mbox{\rm{1}}\mbox{\rm{l}}_{\\{J_{k}<T\\}}$
$\displaystyle\hskip
56.9055pt+\frac{1}{2c_{\mathrm{UE}}^{1/2}}\Big{\\{}\varepsilon_{m_{k}}+\Big{(}\rho_{b}(\varepsilon_{m_{k}}^{2})^{1/2}+R_{m_{k}}\frac{\rho_{\sigma}(\varepsilon_{m_{k}}^{2})}{\varepsilon_{m_{k}}}\Big{)}^{1/2}\Big{\\}}.$

Note that the right-hand side tends to zero as $k\to\infty$. Consequently,
$\\{\mathrm{Law}_{\mathbb{P}}(Y^{k}|_{[0,T]})\\}_{k\in\mathbb{N}}$ is a Cauchy
sequence in $(\mathcal{P}(\Lambda_{T}),d_{\mathrm{LP}})$ for any $T>0$. In
particular, $\\{\mathrm{Law}_{\mathbb{P}}(Y^{k})\\}_{k\in\mathbb{N}}$ is
relatively weakly compact in $\mathcal{P}(\Lambda)$. By Prokhorov’s theorem,
we see that $\\{Y^{k}\\}_{k\in\mathbb{N}}$ is tight on the Polish space
$\Lambda$. Then, by 4.4, weak existence holds for the SEE (1.3).

####  4.3.2 Proof of 2: Uniqueness in law, weak continuity and the Markov
property

Next, we prove 2. This can be proved by a similar manner as in the proof of 1.
Let $y\in\mathcal{H}_{\mu}$ be given, and suppose that
$(Y,W,\Omega,\mathcal{F},\mathbb{F},\mathbb{P})$ is a weak solution of the
original SEE (1.3) with the initial condition $y$. For each $k\in\mathbb{N}$,
again we denote by $Y^{k}$ the mild solution of the SEE (4.15) with the
initial condition $y$, which is now defined for the Brownian motion $W$ on the
filtered probability space $(\Omega,\mathcal{F},\mathbb{F},\mathbb{P})$
specified above. Fix $k\in\mathbb{N}$. Here, we apply 4.5 to the given
coefficients $b,\sigma$ and

$\bar{b}=b_{k},\ \bar{\sigma}=\sigma_{k},\ m=m_{k},\ M=\infty,\
\bar{M}=M_{k},\ \lambda=\lambda_{k},\ J=J_{k},$ $\Delta_{0}=\Delta_{0,k},\
\Delta_{1}=\Delta_{1,k},\ \Delta_{2}=\Delta_{2,k},\ \Delta_{3}=\Delta_{3,k}.$

By A.1, there exists a unique mild solution $\hat{Y}^{k}$ of the controlled
SEE

$\begin{dcases}\mathrm{d}\hat{Y}^{k}_{t}(\theta)=-\theta\hat{Y}^{k}_{t}(\theta)\,\mathrm{d}t+b_{k}(\mu_{M_{k}}[\hat{Y}^{k}_{t}])\,\mathrm{d}t+\sigma_{k}(\mu_{M_{k}}[\hat{Y}^{k}_{t}])\,\mathrm{d}W_{t}+\lambda_{k}\mu_{m_{k}}[Y_{t}-\hat{Y}^{k}_{t}]\mbox{\rm{1}}\mbox{\rm{l}}_{[0,\tau_{k}]}(t)\,\mathrm{d}t,\\\
\hskip 199.16928pt\theta\in{\mathrm{supp}\,\mu},\ t>0,\\\
\hat{Y}^{k}_{0}(\theta)=y(\theta),\ \
\theta\in{\mathrm{supp}\,\mu},\end{dcases}$

with the stopping time $\tau_{k}$ defined by

$\tau_{k}:=\inf\left\\{t\geq 0\mathrel{}\middle|\mathrel{}\begin{aligned}
\int^{t}_{0}\big{|}\mu_{m_{k}}[Y_{s}-\hat{Y}^{k}_{s}]\big{|}^{2}\,\mathrm{d}s\geq\Delta_{1,k},\
\int^{t}_{0}\int_{\mathrm{supp}\,\mu}\theta
r_{m_{k}}(\theta)|Y_{s}(\theta)-\hat{Y}^{k}_{s}(\theta)|^{2}\,\mu(\mathrm{d}\theta)\,\mathrm{d}s\geq\Delta_{2,k}\\\
\text{or}\
\|Y_{t}-\hat{Y}^{k}_{t}\|^{2}_{m_{k}}\geq\Delta_{3,k}\end{aligned}\right\\}.$

Define

$\zeta_{k}:=\inf\left\\{t\geq
0\mathrel{}\middle|\mathrel{}\int^{t}_{0}\|Y_{s}\|^{2}_{\mathcal{V}_{\mu}}\,\mathrm{d}s\geq
J_{k}\right\\}\wedge J_{k}$

and

$\hat{\Omega}_{k}=\left\\{\begin{aligned}
&\|Y_{t\wedge\tau_{k}\wedge\zeta_{k}}-\hat{Y}^{k}_{t\wedge\tau_{k}\wedge\zeta_{k}}\|^{2}_{m_{k}}+2\int^{t\wedge\tau_{k}\wedge\zeta_{k}}_{0}\int_{\mathrm{supp}\,\mu}\theta
r_{m_{k}}(\theta)|Y_{s}(\theta)-\hat{Y}^{k}_{s}(\theta)|^{2}\,\mu(\mathrm{d}\theta)\,\mathrm{d}s\\\
&\hskip
142.26378pt+2(\lambda_{k}-1)\int^{t\wedge\tau_{k}\wedge\zeta_{k}}_{0}\big{|}\mu_{m_{k}}[Y_{s}-\hat{Y}^{k}_{s}]\big{|}^{2}\,\mathrm{d}s\\\
&\leq\frac{1}{2}\Big{(}\rho_{b}(\varepsilon_{m_{k}}^{2})^{1/2}+R_{m_{k}}\frac{\rho_{\sigma}(\varepsilon_{m_{k}}^{2})}{\varepsilon_{m_{k}}}\Big{)}^{1/2}\varepsilon_{m_{k}}\hskip
28.45274pt\text{for any $t\geq 0$}\end{aligned}\right\\}.$

Then, by 4.5 and the estimates (4.12), (4.13) and (4.14), we obtain

$\mathbb{P}(\Omega\setminus\hat{\Omega}_{k})\leq
648\Big{(}\rho_{b}(\varepsilon_{m_{k}}^{2})^{1/2}+R_{m_{k}}\frac{\rho_{\sigma}(\varepsilon_{m_{k}}^{2})}{\varepsilon_{m_{k}}}\Big{)}^{1/4}$

and

$d_{\mathrm{TV}}\big{(}\mathrm{Law}_{\mathbb{P}}(\hat{Y}^{k}),\mathrm{Law}_{\mathbb{P}}(Y^{k})\big{)}\leq\frac{1}{2c_{\mathrm{UE}}^{1/2}}\Big{\\{}\varepsilon_{m_{k}}+\Big{(}\rho_{b}(\varepsilon_{m_{k}}^{2})^{1/2}+R_{m_{k}}\frac{\rho_{\sigma}(\varepsilon_{m_{k}}^{2})}{\varepsilon_{m_{k}}}\Big{)}^{1/2}\Big{\\}}.$

Furthermore, by the same arguments as in the proof of 1, we get

$\hat{\Omega}_{k}\subset\\{\zeta_{k}<\tau_{k}\\}.$ (4.20)

Also, by 4.2, we have

$\mathbb{E}\Big{[}\int^{T}_{0}\|Y_{t}\|^{2}_{\mathcal{V}_{\mu}}\,\mathrm{d}t\Big{]}\leq
C_{0}e^{C_{0}T}(1+\|y\|^{2}_{\mathcal{H}_{\mu}})$

for any $T>0$, where
$C_{0}=C_{0}(\mu,|b(0)|,|\sigma(0)|,\rho_{b}(1),\rho_{\sigma}(1))\geq 1$ is a
constant which does not depend on $T$, $y$ or the choice of the weak solution
$(Y,W,\Omega,\mathcal{F},\mathbb{F},\mathbb{P})$. Then, by the same arguments
as in the proof of 1, we can show that

$\displaystyle
d_{\mathrm{LP}}\big{(}\mathrm{Law}_{\mathbb{P}}(Y|_{[0,T]}),\mathrm{Law}_{\mathbb{P}}(Y^{k}|_{[0,T]})\big{)}$
$\displaystyle\leq\Big{\\{}648+C_{0}e^{C_{0}T}(1+\|y\|^{2}_{\mathcal{H}_{\mu}})\Big{\\}}\Big{(}\rho_{b}(\varepsilon_{m_{k}}^{2})^{1/2}+R_{m_{k}}\frac{\rho_{\sigma}(\varepsilon_{m_{k}}^{2})}{\varepsilon_{m_{k}}}\Big{)}^{1/4}+\mbox{\rm{1}}\mbox{\rm{l}}_{\\{J_{k}<T\\}}$
$\displaystyle\hskip
56.9055pt+\frac{1}{2c_{\mathrm{UE}}^{1/2}}\Big{\\{}\varepsilon_{m_{k}}+\Big{(}\rho_{b}(\varepsilon_{m_{k}}^{2})^{1/2}+R_{m_{k}}\frac{\rho_{\sigma}(\varepsilon_{m_{k}}^{2})}{\varepsilon_{m_{k}}}\Big{)}^{1/2}\Big{\\}}$

for any $k\in\mathbb{N}$, $T>0$ and $y\in\mathcal{H}_{\mu}$. Therefore, for
any $T>0$, we have $\lim_{k\to\infty}Y^{k}|_{[0,T]}=Y|_{[0,T]}$ weakly on
$\Lambda_{T}$, and the convergence is uniform with respect to the initial
condition $y$ on each bounded subset of $\mathcal{H}_{\mu}$. This implies that
$\lim_{k\to\infty}Y^{k}=Y$ weakly on $\Lambda$ uniformly in $y$ on each
bounded subset of $\mathcal{H}_{\mu}$. This means that the law of $Y$ is
specified as the weak limit of the sequence $\\{Y^{k}\\}_{k\in\mathbb{N}}$ of
the mild solutions of the approximating SEEs (4.15) which satisfy strong
existence and pathwise uniqueness (and in particular uniqueness in law).
Consequently, uniqueness in law holds for the original SEE (1.3). Furthermore,
by the weak continuity of the map $y\mapsto Y^{k}$, together with the local
uniformity in $y\in\mathcal{H}_{\mu}$ of the weak convergence
$\lim_{k\to\infty}Y^{k}=Y$, we see that $y\mapsto\mathrm{Law}_{\mathbb{P}}(Y)$
is continuous as a map from
$(\mathcal{H}_{\mu},\|\cdot\|_{\mathcal{H}_{\mu}})$ to
$(\mathcal{P}(\Lambda),d_{\mathrm{LP}})$. The Markov property of $Y$ follows
from that of $Y^{k}$ and a standard approximation argument; see the last part
of the proof of 4.3. The Feller property follows from the weak continuity of
$y\mapsto\mathrm{Law}_{\mathbb{P}}(Y)$. This completes the proof of 3.1 in the
“singular case” (i).

###  4.4 Some remarks on the proof of the main theorem

Let us make some remarks on the proof of 3.1. We will give a sketch of the
proof of the “regular case” and discuss about the “balance condition”.

####  4.4.1 The regular case

Now we consider the “regular case” (ii) of 3.1. Namely, assume that the kernel
$K$ is regular (i.e., $\lim_{t\downarrow 0}K(t)<\infty$, or equivalently
$\mu([0,\infty))<\infty$) and that the following holds:

$\liminf_{\delta\downarrow 0}\frac{\rho_{\sigma}(\delta^{2})}{\delta}=0.$

By the regularity of $K$, we can take $r\equiv 1$ in 1. In this case, for any
$m\in[1,\infty)$, we have $r_{m}(\theta)=1$ for any $\theta\in[0,\infty)$, and
thus

$\displaystyle
R_{m}=\int_{\mathrm{supp}\,\mu}r_{m}(\theta)\,\mu(\mathrm{d}\theta)=\mu([0,\infty))<\infty,$
$\displaystyle\varepsilon_{m}=\int_{\mathrm{supp}\,\mu}(1-r_{m}(\theta))^{2}\theta^{-1}r_{m}(\theta)^{-1}\,\mu(\mathrm{d}\theta)=0.$

Take a decreasing sequence $\\{\delta_{k}\\}_{k\in\mathbb{N}}\subset(0,1]$
such that $\lim_{k\to\infty}\delta_{k}=0$ and
$\lim_{k\to\infty}\frac{\rho_{\sigma}(\delta^{2}_{k})}{\delta_{k}}=0$. For
each $k\in\mathbb{N}$, let $m_{k}=M_{k}=1$, and define
$\Delta_{0,k},\Delta_{1,k},\Delta_{2,k},\Delta_{3,k}\in(0,\infty)$,
$\lambda_{k}\in(1,\infty)$ and $J_{k}\in[1,\infty)$ by

$\Delta_{0,k}=\rho_{b}(\delta^{2}_{k})^{1/2}\delta_{k}+\rho_{\sigma}(\delta^{2}_{k}),\
\Delta_{1,k}=\delta^{2}_{k},\ \Delta_{2,k}=1,\
\Delta_{3,k}=\Big{(}\rho_{b}(\delta^{2}_{k})^{1/2}+\frac{\rho_{\sigma}(\delta^{2}_{k})}{\delta_{k}}\Big{)}^{1/2}\delta_{k},$

and

$\lambda_{k}=1+\Big{(}\rho_{b}(\delta^{2}_{k})^{1/2}+\frac{\rho_{\sigma}(\delta^{2}_{k})}{\delta_{k}}\Big{)}^{1/2}\delta^{-1}_{k},\
J_{k}=\Big{(}\rho_{b}(\delta^{2}_{k})^{1/2}+\frac{\rho_{\sigma}(\delta^{2}_{k})}{\delta_{k}}\Big{)}^{-1/4}.$

Construct suitable maps $b_{k}:\mathbb{R}^{n}\to\mathbb{R}^{n}$ and
$\sigma_{k}:\mathbb{R}^{n}\to\mathbb{R}^{n\times d}$ by taking convolutions of
$b$ and $\sigma$ with smooth mollifiers. Then, by the same arguments as in
Section 4.3, we can show that the conclusions of 3.1 hold in the regular case
(ii).

####  4.4.2 The balance condition

Let us remark on the “balance condition”

$\liminf_{m\to\infty}R_{m}\frac{\rho_{\sigma}(\varepsilon^{2}_{m})}{\varepsilon_{m}}=0$
(4.21)

appearing in the singular case (i) of 3.1. Note that, in the singular case,
$\varepsilon_{m}>0$ for any $m\in[1,\infty)$, $\lim_{m\to\infty}R_{m}=\infty$
and $\lim_{m\to\infty}R_{m}\varepsilon_{m}=0$. By means of 4.5, the key step
in Section 4.3 is to construct sequences
$\\{m_{k}\\}_{k\in\mathbb{N}}\subset[1,\infty)$,
$\\{\lambda_{k}\\}_{k\in\mathbb{N}}\subset(1,\infty)$ and
$\\{\Delta_{i,k}\\}_{k\in\mathbb{N}}\subset(0,\infty)$, $i=1,2,3$, such that
$\lim_{k\to\infty}m_{k}=\infty$ and that

$\begin{dcases}\Delta_{3,k}\leq C(\lambda_{k}-1)\Delta_{1,k}\ \text{for any
$k\in\mathbb{N}$ for some constant $C\in(0,\infty)$},\\\ \Delta_{3,k}\leq
C\Delta_{2,k}\ \text{for any $k\in\mathbb{N}$ for some constant
$C\in(0,\infty)$},\\\ \lim_{k\to\infty}\Delta_{3,k}=0,\\\
\lim_{k\to\infty}R_{m_{k}}\Delta^{-1}_{3,k}\rho_{\sigma}(\Delta_{1,k}+\varepsilon_{m_{k}}\Delta_{2,k})=0,\\\
\lim_{k\to\infty}\lambda_{k}\Delta^{1/2}_{1,k}=0.\end{dcases}$ (4.22)

Here, in view of the C-n-R strategy (4.1), the first four lines of (4.22)
correspond to the “control-step” (4.7), and the last line corresponds to the
“reimbursement-step” (4.8). The first two lines are required in order to
ensure that the stopping time $\tau$ is sufficiently large on the event
$\hat{\Omega}$; see (4.18) and (4.20).

We have shown in Section 4.3 that the balance condition (4.21) makes it
possible to construct sequences satisfying all the requirements in (4.22).
Conversely, if the diffusion coefficient $\sigma$ is of “Hölder-type”, i.e.,
$0<\liminf_{t\downarrow
0}\frac{\rho_{\sigma}(t)}{t^{\gamma}}\leq\limsup_{t\downarrow
0}\frac{\rho_{\sigma}(t)}{t^{\gamma}}<\infty$ for some constant
$\gamma\in(0,1]$, then the balance condition (4.21) is necessary for the
existence of such sequences. Indeed, by easy calculations, we can show that
(4.22) implies that

$\lim_{k\to\infty}\Delta_{1,k}=0,\
\lim_{k\to\infty}\varepsilon_{m_{k}}\Delta_{2,k}=0,\
\lim_{k\to\infty}\Delta_{3,k}\Delta^{-1/2}_{1,k}=0,$

and that

$\lim_{k\to\infty}R_{m_{k}}\Delta^{\gamma-1/2}_{1,k}=0,\
\lim_{k\to\infty}R_{m_{k}}\varepsilon^{\gamma}_{m_{k}}\Delta^{\gamma-1}_{2,k}=0,\
\lim_{k\to\infty}R_{m_{k}}\Delta^{\gamma}_{1,k}\Delta^{-1}_{2,k}=0,\
\lim_{k\to\infty}R_{m_{k}}\varepsilon^{\gamma}_{m_{k}}\Delta^{\gamma}_{2,k}\Delta^{-1/2}_{1,k}=0.$
(4.23)

Since $\lim_{k\to\infty}R_{m_{k}}=\infty$ and
$\lim_{k\to\infty}\Delta_{1,k}=0$, the first equality in (4.23) implies that
$\gamma>\frac{1}{2}$. Now assume that the balance condition (4.21) does not
hold; $\liminf_{m\to\infty}R_{m}\varepsilon^{2\gamma-1}_{m}>0$. Then, noting
that $\gamma\in(\frac{1}{2},1)$, the first, second and fourth equalities in
(4.23) imply that

$\lim_{k\to\infty}\frac{\Delta_{1,k}}{\varepsilon^{2}_{m_{k}}}=0,\
\lim_{k\to\infty}\frac{\varepsilon_{m_{k}}}{\Delta_{2,k}}=0\ \text{and}\
\lim_{k\to\infty}\Big{(}\frac{\Delta_{2,k}}{\varepsilon_{m_{k}}}\Big{)}^{\gamma}\Big{(}\frac{\varepsilon^{2}_{m_{k}}}{\Delta_{1,k}}\Big{)}^{1/2}=0,$

respectively, and we get a contradiction.

The above discussions indicate that the balance condition (4.21) is essential
for the C-n-R strategy to prove 3.1 based on the estimates appearing in 4.5.

##  Appendix A Appendix

In this appendix, we prove existence and uniqueness of the mild solutions of
SEEs of the form (4.6) involving stopping times which depend on the solutions.

######  Proposition A.1.

Suppose that 1 holds. Let a filtered probability space
$(\Omega,\mathcal{F},\mathbb{F},\mathbb{P})$ and a $d$-dimensional
$\mathbb{F}$-Brownian motion $W$ be fixed. Let
$b_{1},b_{2}:\Omega\times[0,\infty)\times\mathcal{H}_{\mu}\to\mathbb{R}^{n}$
and
$\sigma_{1},\sigma_{2}:\Omega\times[0,\infty)\times\mathcal{H}_{\mu}\to\mathbb{R}^{n\times
d}$ be measurable maps satisfying the following conditions:

  * •

For $\varphi=b_{1},b_{2},\sigma_{1},\sigma_{2}$, the process
$(\varphi(t,y))_{t\geq 0}$ is predictable for any $y\in\mathcal{H}_{\mu}$.

  * •

There exist a constant $c_{\mathrm{LG}}>0$ and a nonnegative predictable
process $\xi$ satisfying
$\mathbb{E}[\int^{T}_{0}\xi^{2}_{t}\,\mathrm{d}t]<\infty$ for any $T>0$ such
that, for $\varphi=b_{1},b_{2},\sigma_{1},\sigma_{2}$,

$|\varphi(t,y)|\leq c_{\mathrm{LG}}(\xi_{t}+\|y\|_{\mathcal{H}_{\mu}})$

for any $t\geq 0$ and any $y\in\mathcal{H}_{\mu}$ a.s.

  * •

There exist $L_{k}>0$, $k\in\mathbb{N}$, such that, for
$\varphi=b_{1},b_{2},\sigma_{1},\sigma_{2}$,

$|\varphi(t,y)-\varphi(t,\bar{y})|\leq L_{k}\|y-\bar{y}\|_{\mathcal{H}_{\mu}}$

for any $t\geq 0$ and any $y,\bar{y}\in\mathcal{H}_{\mu}$ satisfying
$\|y\|_{\mathcal{H}_{\mu}}\vee\|\bar{y}\|_{\mathcal{H}_{\mu}}\leq k$ a.s., for
any $k\in\mathbb{N}$.

Also, let
$f:\Omega\times[0,\infty)\times\mathcal{V}_{\mu}\to\mathbb{R}^{d_{f}}$ and
$g:\Omega\times[0,\infty)\times\mathcal{H}_{\mu}\to\mathbb{R}^{d_{g}}$ with
$d_{f},d_{g}\in\mathbb{N}$ be measurable maps satisfying the following
conditions:

  * •

The process $(f(t,y))_{t\geq 0}$ is predictable for any
$y\in\mathcal{V}_{\mu}$, and the process $(g(t,y))_{t\geq 0}$ is predictable
for any $y\in\mathcal{H}_{\mu}$.

  * •

There exist a constant $c_{f}>0$ and a nonnegative predictable process $\eta$
satisfying $\int^{T}_{0}\eta_{t}\,\mathrm{d}t<\infty$ a.s. for any $T>0$ such
that

$|f(t,y)|\leq c_{f}(\eta_{t}+\|y\|^{2}_{\mathcal{V}_{\mu}})$

for any $t\geq 0$ and any $y\in\mathcal{V}_{\mu}$ a.s.

  * •

For $\mathbb{P}$-almost every $\omega\in\Omega$, the map
$[0,\infty)\times\mathcal{H}_{\mu}\ni(t,y)\mapsto
g(\omega,t,y)\in\mathbb{R}^{d_{g}}$ is continuous.

Let $O_{f}\subset\mathbb{R}^{d_{f}}$ and $O_{g}\subset\mathbb{R}^{d_{g}}$ be
open sets. Then, for any initial condition $y\in\mathcal{H}_{\mu}$, there
exists a unique mild solution $Y$ to the following SEE:

$\begin{dcases}\mathrm{d}Y_{t}(\theta)=-\theta
Y_{t}(\theta)\,\mathrm{d}t+\big{\\{}b_{1}(t,Y_{t})+b_{2}(t,Y_{t})\mbox{\rm{1}}\mbox{\rm{l}}_{[0,\tau]}(t)\big{\\}}\,\mathrm{d}t+\big{\\{}\sigma_{1}(t,Y_{t})+\sigma_{2}(t,Y_{t})\mbox{\rm{1}}\mbox{\rm{l}}_{[0,\tau]}(t)\big{\\}}\,\mathrm{d}W_{t},\\\
\hskip 199.16928pt\theta\in{\mathrm{supp}\,\mu},\ t>0,\\\
Y_{0}(\theta)=y(\theta),\ \ \theta\in{\mathrm{supp}\,\mu}.\end{dcases}$ (A.1)

with the stopping time $\tau=\tau(Y)$ given by

$\tau=\inf\left\\{t\geq
0\mathrel{}\middle|\mathrel{}\int^{t}_{0}f(s,Y_{s})\,\mathrm{d}s\notin O_{f}\
\text{or}\ g(t,Y_{t})\notin O_{g}\right\\}.$

###### Proof.

By [16, Theorem 2.18], there exists a unique mild solution $\hat{Y}$ to the
SEE

$\begin{dcases}\mathrm{d}\hat{Y}_{t}(\theta)=-\theta\hat{Y}_{t}(\theta)\,\mathrm{d}t+\big{\\{}b_{1}(t,\hat{Y}_{t})+b_{2}(t,\hat{Y}_{t})\big{\\}}\,\mathrm{d}t+\big{\\{}\sigma_{1}(t,\hat{Y}_{t})+\sigma_{2}(t,\hat{Y}_{t})\big{\\}}\,\mathrm{d}W_{t},\
\ \theta\in{\mathrm{supp}\,\mu},\ t>0,\\\ \hat{Y}_{0}(\theta)=y(\theta),\ \
\theta\in{\mathrm{supp}\,\mu}.\end{dcases}$

Define a stopping time $\hat{\tau}=\hat{\tau}(\hat{Y})$ by

$\hat{\tau}:=\inf\left\\{t\geq
0\mathrel{}\middle|\mathrel{}\int^{t}_{0}f(s,\hat{Y}_{s})\,\mathrm{d}s\notin
O_{f}\ \text{or}\ g(t,\hat{Y}_{t})\notin O_{g}\right\\}.$

Again by [16, Theorem 2.18], there exists a unique mild solution $Y$ to the
SEE

$\begin{dcases}\mathrm{d}Y_{t}(\theta)=-\theta
Y_{t}(\theta)\,\mathrm{d}t+\big{\\{}b_{1}(t,Y_{t})+b_{2}(t,Y_{t})\mbox{\rm{1}}\mbox{\rm{l}}_{[0,\hat{\tau}]}(t)\big{\\}}\,\mathrm{d}t+\big{\\{}\sigma_{1}(t,Y_{t})+\sigma_{2}(t,Y_{t})\mbox{\rm{1}}\mbox{\rm{l}}_{[0,\hat{\tau}]}(t)\big{\\}}\,\mathrm{d}W_{t},\\\
\hskip 199.16928pt\theta\in{\mathrm{supp}\,\mu},\ t>0,\\\
Y_{0}(\theta)=y(\theta),\ \ \theta\in{\mathrm{supp}\,\mu}.\end{dcases}$

Also, for any $k\in\mathbb{N}$, there exists a unique mild solution $Y^{k}$ to
the SEE

$\begin{dcases}\mathrm{d}Y^{k}_{t}(\theta)=-\theta
Y^{k}_{t}(\theta)\,\mathrm{d}t+\big{\\{}b^{k}_{1}(t,Y^{k}_{t})+b^{k}_{2}(t,Y^{k}_{t})\big{\\}}\,\mathrm{d}t+\big{\\{}\sigma^{k}_{1}(t,Y^{k}_{t})+\sigma^{k}_{2}(t,Y_{t})\big{\\}}\,\mathrm{d}W_{t},\
\ \theta\in{\mathrm{supp}\,\mu},\ t>0,\\\ Y^{k}_{0}(\theta)=y(\theta),\ \
\theta\in{\mathrm{supp}\,\mu},\end{dcases}$

where
$b^{k}_{1},b^{k}_{2}:\Omega\times[0,\infty)\times\mathcal{H}_{\mu}\to\mathbb{R}^{n}$
and
$\sigma^{k}_{1},\sigma^{k}_{2}:\Omega\times[0,\infty)\times\mathcal{H}_{\mu}\to\mathbb{R}^{n\times
d}$ are defined by

$\varphi^{k}(t,y):=\begin{dcases}\varphi(t,y)\ &\text{if
$\|y\|_{\mathcal{H}_{\mu}}<k$},\\\
\varphi\Big{(}t,\frac{k}{\|y\|_{\mathcal{H}_{\mu}}}y\Big{)}\ &\text{if
$\|y\|_{\mathcal{H}_{\mu}}\geq k$},\end{dcases}$

for $\varphi=b_{1},b_{2},\sigma_{1},\sigma_{2}$. Note that, for
$\varphi=b_{1},b_{2},\sigma_{1},\sigma_{2}$,

$|\varphi^{k}(t,y)-\varphi^{k}(t,\bar{y})|\leq
L_{k}\|y-\bar{y}\|_{\mathcal{H}_{\mu}}$

for any $t\geq 0$ and any $y,\bar{y}\in\mathcal{H}_{\mu}$ a.s. Thus, by [16,
Theorem 2.17], there exists a constant $C_{k}>0$ such that, for any $T>0$,

$\displaystyle\mathbb{E}\Big{[}\sup_{t\in[0,T]}\|Y^{k}_{t\wedge\hat{\tau}\wedge\zeta_{k}}-Y_{t\wedge\hat{\tau}\wedge\zeta_{k}}\|^{2}_{\mathcal{H}_{\mu}}+\int^{T\wedge\hat{\tau}\wedge\zeta_{k}}_{0}\|Y^{k}_{t}-Y_{t}\|^{2}_{\mathcal{V}_{\mu}}\,\mathrm{d}t\Big{]}$
$\displaystyle\leq
C_{k}e^{C_{k}T}\mathbb{E}\Big{[}\int^{T\wedge\hat{\tau}\wedge\zeta_{k}}_{0}\Big{\\{}|b^{k}_{1}(t,Y_{t})+b^{k}_{2}(t,Y_{t})-b_{1}(t,Y_{t})-b_{2}(t,Y_{t})|^{2}$
$\displaystyle\hskip
142.26378pt+|\sigma^{k}_{1}(t,Y_{t})+\sigma^{k}_{2}(t,Y_{t})-\sigma_{1}(t,Y_{t})-\sigma_{2}(t,Y_{t})|^{2}\Big{\\}}\,\mathrm{d}t\Big{]}$

and

$\displaystyle\mathbb{E}\Big{[}\sup_{t\in[0,T]}\|Y^{k}_{t\wedge\hat{\zeta}_{k}}-\hat{Y}_{t\wedge\hat{\zeta}_{k}}\|^{2}_{\mathcal{H}_{\mu}}+\int^{T\wedge\hat{\zeta}_{k}}_{0}\|Y^{k}_{t}-\hat{Y}_{t}\|^{2}_{\mathcal{V}_{\mu}}\,\mathrm{d}t\Big{]}$
$\displaystyle\leq
C_{k}e^{C_{k}T}\mathbb{E}\Big{[}\int^{T\wedge\hat{\zeta}_{k}}_{0}\Big{\\{}|b^{k}_{1}(t,\hat{Y}_{t})+b^{k}_{2}(t,\hat{Y}_{t})-b_{1}(t,\hat{Y}_{t})-b_{2}(t,\hat{Y}_{t})|^{2}$
$\displaystyle\hskip
142.26378pt+|\sigma^{k}_{1}(t,\hat{Y}_{t})+\sigma^{k}_{2}(t,\hat{Y}_{t})-\sigma(t,\hat{Y}_{t})-\sigma_{2}(t,\hat{Y}_{t})|^{2}\Big{\\}}\,\mathrm{d}t\Big{]},$

where $\zeta_{k}$ and $\hat{\zeta}_{k}$ are stopping times defined by

$\zeta_{k}:=\inf\\{t\geq 0\,|\,\|Y_{t}\|_{\mathcal{H}_{\mu}}\geq k\\}\
\text{and}\ \hat{\zeta}_{k}:=\inf\\{t\geq
0\,|\,\|\hat{Y}_{t}\|_{\mathcal{H}_{\mu}}\geq k\\}.$

By the definitions of $b^{k}_{1},b^{k}_{2}$, $\sigma^{k}_{1},\sigma^{k}_{2}$,
$\zeta_{k}$ and $\hat{\zeta}_{k}$, we see that, for any $T>0$,

$Y^{k}_{t}=Y_{t}\ \ \text{for any $t\in[0,T\wedge\hat{\tau}\wedge\zeta_{k}]$
a.s.}$

and

$Y^{k}_{t}=\hat{Y}_{t}\ \ \text{for any $t\in[0,T\wedge\hat{\zeta}_{k}]$
a.s.}$

This implies that

$\mathbb{P}\big{(}Y_{t}=\hat{Y}_{t}\ \text{for any
$t\in[0,T\wedge\hat{\tau}]$}\big{)}\geq\mathbb{P}\big{(}\zeta_{k}>T\
\text{and}\ \hat{\zeta}_{k}>T\big{)}$

for any $T>0$ and any $k\in\mathbb{N}$. By the $\mathcal{H}_{\mu}$-continuity
of the processes $Y$ and $\hat{Y}$, the right-hand side tends to $1$ as
$k\to\infty$. Thus, we see that $Y_{t}=\hat{Y}_{t}$ for any
$t\in[0,T\wedge\hat{\tau}]$ a.s. for any $T>0$, which implies that
$Y_{t}=\hat{Y}_{t}$ for any $t\in[0,\hat{\tau}]$ a.s. Now we define

$\tau:=\inf\left\\{t\geq
0\mathrel{}\middle|\mathrel{}\int^{t}_{0}f(s,Y_{s})\,\mathrm{d}s\notin O_{f}\
\text{or}\ g(t,Y_{t})\notin O_{g}\right\\}.$

Noting that the processes $\int^{\cdot}_{0}f(s,Y_{s})\,\mathrm{d}s$,
$g(\cdot,Y_{\cdot})$, $\int^{\cdot}_{0}f(s,\hat{Y}_{s})\,\mathrm{d}s$ and
$g(\cdot,\hat{Y}_{\cdot})$ are continuous and that
$O_{f}\subset\mathbb{R}^{d_{f}}$ and $O_{g}\subset\mathbb{R}^{d_{g}}$ are open
sets, we have

$\displaystyle\mathbb{P}\big{(}\hat{\tau}<\tau\big{)}$
$\displaystyle=\mathbb{P}\Big{(}\int^{\hat{\tau}}_{0}f(s,Y_{s})\,\mathrm{d}s\in
O_{f},\ g(\hat{\tau},Y_{\hat{\tau}})\in O_{g}\ \text{and}\
\hat{\tau}<\tau\Big{)}$
$\displaystyle\leq\mathbb{P}\Big{(}\int^{\hat{\tau}}_{0}f(s,\hat{Y}_{s})\,\mathrm{d}s\in
O_{f},\ g(\hat{\tau},\hat{Y}_{\hat{\tau}})\in O_{g}\Big{)}=0$

and

$\displaystyle\mathbb{P}\big{(}\tau<\hat{\tau}\big{)}$
$\displaystyle=\mathbb{P}\Big{(}\int^{\tau}_{0}f(s,\hat{Y}_{s})\,\mathrm{d}s\in
O_{f},\ g(\tau,\hat{Y}_{\tau})\in O_{g}\ \text{and}\ \tau<\hat{\tau}\Big{)}$
$\displaystyle\leq\mathbb{P}\Big{(}\int^{\tau}_{0}f(s,Y_{s})\,\mathrm{d}s\in
O_{f},\ g(\tau,Y_{\tau})\in O_{g}\Big{)}=0.$

Thus, we see that $\tau=\hat{\tau}$ a.s. By the construction, the process $Y$
is a mild solution of the SEE (A.1).

Now we prove the uniqueness. Assume that there is another mild solution
$\bar{Y}$ to the SEE (A.1) with the stopping time $\tau$ replaced by

$\bar{\tau}=\inf\left\\{t\geq
0\mathrel{}\middle|\mathrel{}\int^{t}_{0}f(s,\bar{Y}_{s})\,\mathrm{d}s\notin
O_{f}\ \text{or}\ g(t,\bar{Y}_{t})\notin O_{g}\right\\}.$

By the same arguments as above, we see that $Y_{t}=\bar{Y}_{t}$ for any
$t\in[0,\tau\wedge\bar{\tau}]$ a.s. This implies that $\tau=\bar{\tau}$ a.s.
Therefore, $Y$ and $\bar{Y}$ solve the same SEE. By the uniqueness of the mild
solution, we get $Y_{t}=\bar{Y}_{t}$ for any $t\geq 0$ a.s. This completes the
proof. ∎

## References

  * [1] E. Abi Jaber, Weak existence and uniqueness for affine stochastic Volterra equations with $L^{1}$-kernels, Bernoulli, 27(3), 1583–1615, 2021\.
  * [2] E. Abi Jaber, C. Cuchiero, M. Larsson, and S. Pulido, A weak solution theory for stochastic Volterra equations of convolution type, Ann. Appl. Probab., 31(6), 2924–2952, 2021\.
  * [3] E. Abi Jaber and O. El Euch, Multifactor approximation of rough volatility models, SIAM J. Financ. Math., 10(2), 309–349, 2019\.
  * [4] E. Abi Jaber, M. Larsson, S. Pulido, Affine Volterra processes, Ann. Appl. Probab., 29(5), 3155–3200, 2019\.
  * [5] O. Barndorff-Nielsen, F. Benth, and A. Veraart, Ambit processes and stochastic partial differential equations, In Advanced Mathematical Methods for Finance, eds G. Di Nunno and B. Øksendal, 35–74, Springer, Berlin, Heidelberg, 2011\.
  * [6] M.A. Berger and V.J. Mizel, Volterra equations with Itô integrals—I, J. Integral Equations, 2(3), 187–245, 1980\.
  * [7] M.A. Berger and V.J. Mizel, Volterra equations with Itô integrals—II, J. Integral Equations, 2(4), 319–337, 1980\.
  * [8] O. Butkovsky, A. Kulik, and M. Scheutzow, Generalized couplings and ergodic rates for SPDEs and other Markov models, Ann. Appl. Probab., 30, 1–39, 2020\.
  * [9] A.S. Cherny, On the uniqueness in law and the pathwise uniqueness for stochastic differential equations, Theory Probab. Appl., 46(3), 406–419, 2003\.
  * [10] C. Cuchiero and J. Teichmann, Markovian lifts of positive semidefinite affine Volterra-type processes, Decis. Econ. Finance, 42, 407–448, 2019\.
  * [11] C. Cuchiero and J. Teichmann, Generalized Feller processes and Markovian lifts of stochastic Volterra processes: the affine case, J. Evol. Equ., 20, 1301–1348, 2020\.
  * [12] R.M. Dudley, Real Analysis and Probability, Cambridge University Press, 2002\.
  * [13] O. El Euch, M. Fukasawa, and M. Rosenbaum, The microstructural foundations of leverage effect and rough volatility, Finance Stoch., 22, 241–280, 2018\.
  * [14] F. Flandoli, Random Perturbation of PDEs and Fluid Dynamic Models, Springer, Berlin, 2011\.
  * [15] G. Gripenberg, S.O. Londen, and O. Staffans, Volterra Integral and Functional Equations, Cambridge University Press, 1990\.
  * [16] Y. Hamaguchi, Markovian lifting and asymptotic log-Harnack inequality for stochastic Volterra integral equations, preprint, arXiv:2304.06683, 2023\.
  * [17] J. Jacod and A.N. Shiryaev, Limit Theorems for Stochastic Processes, 2nd edition, Springer, Berlin, 2003\.
  * [18] T. Jaisson and M. Rosenbaum, Rough fractional diffusions as scaling limits of nearly unstable heavy tailed Hawkes processes, Ann. Appl. Probab., 26(5), 2860–2882, 2016\.
  * [19] I. Karatzas and S.E. Shreve, Brownian Motion and Stochastic Calculus, 2nd edition, Springer, New York, 1991\.
  * [20] A. Kulik, Ergodic Behavior of Markov Processes: With Applications to Limit Theorems, De Gruyter, 2018\.
  * [21] A. Kulik and M. Scheutzow, Well-posedness, stability and sensitivities for stochastic delay equations: a generalized coupling approach, Ann. Probab., 48(6), 3041–3076, 2020\.
  * [22] T. Kurtz, Weak and strong solutions of general stochastic models, Electron. Commun. Probab., 19(58), 1–16, 2014\.
  * [23] L. Mytnik and T.S. Salisbury, Uniqueness of Volterra-type stochastic integral equations, preprint, arXiv:1502.05513, 2015\.
  * [24] D.J. Prömel and D. Scheffels, On the existence of weak solutions to stochastic Volterra equations, preprint, arXiv:2207.01367, 2022\.
  * [25] D.J. Prömel and D. Scheffels, Pathwise uniqueness for singular stochastic Volterra equations with Hölder coefficients, preprint, arXiv:2212.08029, 2022\.
  * [26] D.J. Prömel and D. Scheffels, Stochastic Volterra equations with Hölder diffusion coefficients, Stochastic Process. Appl., 161, 291–315, 2023\.
  * [27] S.G. Samko, A.A. Kilbas, and O.I. Marichev, Fractional Integrals and Derivatives, Theory and Applications, Gordon and Breach Science Publishers, Yverdon, Switzerland, 1987\.
  * [28] R.L. Schilling, R. Song, and Z. Vondraček, Bernstein Functions: Theory and Applications, Walter de Gruyter, 2009\.
  * [29] D.W. Stroock and S.R.S. Varadhan, Multidimensional Diffusion Processes, Springer, Berlin, 1979\.
  * [30] M. Veraar, The stochastic Fubini theorem revisited, Stochastics, 84(4), 543–551, 2012\.
  * [31] Z. Wang, Existence and uniqueness of solutions to stochastic Volterra equations with singular kernels and non-Lipschitz coefficients, Statist. Probab. Lett., 78, 1062–1071, 2008\.
  * [32] T. Yamada and S. Watanabe, On the uniqueness of solutions of stochastic differential equations, J. Math. Kyoto Univ., 11, 155–167, 1971\.
  * [33] X. Zhang, Stochastic Volterra equations in Banach spaces and stochastic partial differential equation, J. Funct. Anal., 258, 1361–1425, 2010\.