]www.ndsu.edu/pubweb/$\sim$carswagn

# Lattice Gas with Molecular Dynamics Collision Operator

M. Reza Parsa  Alexander J. Wagner<EMAIL_ADDRESS>[ Department of
Physics, North Dakota State University, Fargo, North Dakota 58108, USA
Program in Materials and Nanotechnology, North Dakota State University, Fargo,
North Dakota 58108, USA

###### Abstract

We introduce a lattice gas implementation that is based on coarse-graining a
Molecular Dynamics (MD) simulation. Such a lattice gas is similar to standard
lattice gases, but its collision operator is informed by an underlying MD
simulation. This can be considered an optimal lattice gas implementation
because it allows for the representation of any system that can be simulated
with MD. We show here that equilibrium behavior of the popular lattice
Boltzmann algorithm is consistent with this optimal lattice gas. This
comparison allows us to make a more accurate identification of the expressions
for temperature and pressure in lattice Boltzmann simulations which turn out
to be related not only to the physical temperature and pressure but also to
the lattice discretization. We show that for any spatial discretization we
need to choose a particular temporal discretization to recover the lattice
Boltzmann equilibrium.

Lattice Boltzmann, Molecular Dynamics, kinetic theory

##  I Introduction

Lattice Boltzmann methods are an important computational tool that is most
commonly employed to simulate hydrodynamic systemsSucci (2001), and it has
been adapted to address many complex phenomena from turbulence Teixeira
(1998); Yu  _et al._ (2005) over multi-phase and multi-componet flowSwift  _et
al._ (1995); Shan and Doolen (1995); He and Luo (1997); Briant  _et al._
(2004); Yuan and Schaefer (2006) to pore-scale simulations of porous
mediaTölke (2002); Kang  _et al._ (2006) and simulations of immersed
boundaries Dupin  _et al._ (2003); Jansen and Harting (2011). It derives its
power from an underlying mesoscopic description that ensures exact mass and
momentum conservation. The exact physical meaning of the lattice Boltzmann
densities, however, remains poorly understood.

The lattice Boltzmann method was derived as a theoretical tool for the
analysis of lattice gas methods McNamara and Zanetti (1988). Lattice gas
methods consist of particles moving on a lattice with velocities that connect
neighboring sites. After the particles have moved a stochastic collision step
rearranges the particles. If these collisions conserve both the number of
particles and the total momentum there will be a hydrodynamics limit for mass
and momentum conservation equation. The introduction of a hexagonal instead of
a square lattice by Frisch, Hasslacher and Pomeau Frisch  _et al._ (1986)
recovered the necessary isotropy to allow the momentum equation to be related
to the Navier-Stokes equation.

These lattice gas models had some deficiencies, one unfavorable feature was a
large and essentially uncontrolled amount of noise that required a significant
amount of averaging. To derive the Navier-Stokes equations from the lattice
gas dynamics a theoretical ensemble average was performed, leading to a
lattice Boltzmann representation. Higuera then proposed to simulate the
ensemble averaged lattice Boltzmann evolution equation directly, and thereby
avoid the need to average results of the lattice gas equation Higuera and
Jimenez (1989); Higuera  _et al._ (1989). The collision operation of this
first lattice Boltzmann method could be mapped one-to-one to the lattice gas
and shared some of the positive features of the lattice gas, like the
existence of an H-theorem with unconditional stability, and also some of its
deficiencies like velocity dependent viscosities.

It was then realized that there existed much more freedom in the choice of the
collision operator, and in particular the relaxation towards a local
equilibrium function, often called the Bhatnagar-Gross-Krook (BGK) approach,
allowed the full recovery of the Navier-Stokes equation to second order Qian
_et al._ (1992).

At this time a second approach to derive the lattice Boltzmann equation
directly from the continuous Boltzmann equation with a BGK collision operator
gained popularity He and Luo (1997). Over the years several different local
equilibrium distributions have been proposed, and currently the most popular
method is a standard form of a second order expansion in terms of velocities.

Typically these lattice Boltzmann methods are validated by their ability to
recover the Navier-Stokes equation. Here, however, we want to establish a
relation to an underlying Molecular Dynamics simulation. For any Molecular
Dynamics simulation, we can bin the particles into lattice cells corresponding
to the lattice Boltzmann lattice. We can then observe where the particles in
cell $x$ migrate to after a time $\Delta t$, and associate these particles
with a lattice velocity vector $v_{i}=x(t)-x(t+\Delta t)$. These particles
will collect at their new lattice cells. After another timestep $\Delta t$
these particles are re-distributed to new lattice sites, and can be associated
with new lattice velocities. We call this representation of the MD simulation
Molecular-Dynamics-lattice-gas (MDLG). This redistribution can be understood
to be an effective MDLG-collision operator. In some very fundamental sense
this is the collision operator that the lattice Boltzmann approach is trying
to mimic. The purpose of this paper is to understand the physical meanting of
the lattice Boltzmann densities in terms of this fundamental MDLG
representation.

The paper is organized as follows: we first introduce a general idea of a
lattice gas and then derive a new lattice gas which consists of a coarse-
graining of an underlying MD simulation. We then apply this general idea to a
specific MD simulation of a Lennard Jones gas in two dimensions. We analyse
the equilibrium properties of the associated MDLG method and show that we are
able to predict its mathematical form analytically. We then introduce the
lattice Boltzmann method and compare the equilibrium properties of the MDLG
method to the lattice Boltzmann equilibrium. We show that there are particular
choices for the coarsgraining time and space discretization that lead to
equilibria that are compatible with the lattice Boltzmann results.

##  II Lattice Gas

At its very basis a lattice gas consists of particles, all located on lattice
points, that move with a lattice velocities $v_{i}$. What we mean by lattice
velocity is that if $x$ is a lattice point, so it $x+v_{i}$. There are
$n_{i}(x,t)$ particles at time $t$ at position $x$ moving with velocity
$v_{i}$. The evolution consists of two steps. A collision step that
redistributes the particles at the lattice point $x$ to different velocities:

$n_{i}^{*}(x,t)=n_{i}(x,t)+\Xi_{i}(\\{n_{j}(x,t)\\})$ (1)

where the collision operator $\Xi_{i}$ is a function of all the particles and
their velocities that are located at lattice point $x$ at time $t$. This
collision operator will ensure that none of the locally conserved quantities
are changed in the collision process. These locally conserved quantities will
vary, depending on the desired physical system that one wants to model. In the
majority of cases one will ensure mass and momentum conservation. Early
lattice gases restricted the number of particles to at most one per velocity
$v_{i}$ at a lattice site, and the velocity vectors all had the same length,
ensuring that mass and energy conservation were synonymous. Each conserved
quantity will lead to a corresponding hydrodynamic equation. Most applications
focused on the fluid flow, and the key hydrodynamic equations to recover were
the continuity and Navier-Stokes equations. Energy conservation is often
abandoned in favor of an isothermal condition for many practical applications.
Local mass and momentum densities are defined as

$\displaystyle\rho$ $\displaystyle=\sum_{i}n_{i}$ (2) $\displaystyle\rho
u_{\alpha}$ $\displaystyle=\sum_{i}v_{i\alpha}n_{i}$ (3)

and the conservation of these quantities then implies

$\displaystyle\sum_{i}\Xi_{i}$ $\displaystyle=0,$ (4)
$\displaystyle\sum_{i}v_{i\alpha}\Xi_{i}$ $\displaystyle=0.$ (5)

For the new kind of lattice gas collision operator proposed in section III we
will see that mass conservation of Eq. (4) is indeed fulfilled, but the
momentum conservation of Eq. (5) is not strictly obeyed. Because the new
algorithm is based on MD, however, the algorithm conserves momentum
rigorously. Its representation of momentum through Eq. (3), however, is
inexact.

This collision is then followed by a streaming step

$n_{i}(x+v_{i},t+1)=n_{i}^{*}(x,t)$ (6)

where the particles move to the lattice site indicated by the velocity index
$i$, i.e. they move from $x$ to $x+v_{i}$. The full evolution equation for
these densities can then be written as

$n_{i}(x+v_{i},t+1)=n_{i}(x,t)+\Xi_{i}.$ (7)

Of course to make this description complete we need to define the collision
operator. Originally lattice gases were defined such that there could be at
most one particle for each $n_{i}$ Frisch  _et al._ (1986). For the purpose of
this paper, however, we will make no such restriction. Instead we investigate
a collision operator that is defined by an underlying molecular dynamic
simulation.

##  III Lattice Gas with Molecular Dynamics collision operator

In principle most systems of interest for a lattice gas (LG) simulation could
be simulated using a Molecular Dynamics (MD) approach as well. MD is a
standard tool which follows classical particle trajectories for particles
interacting with a pair-potential by numerically integrating Newton’s equation
of motion.

Figure 1: Sketch of the MDLG algorithm: a lattice (blue line) is superimposed
on the domain of the MD simulation. Particles in the reciprocal lattice cells
(indicated by the red boundaries), are associated with the corresponding
lattice point. Particles then get associated with the $n_{i}$ for the $v_{i}$
which corresponds to their lattice displacement in the time-interval $\Delta
t$.

To construct a lattice gas method from a molecular dynamics simulation we
overlay a lattice onto our MD simulation. The number of particles in each
reciprocal lattice cell around the lattice position $x$ then corresponds to
the lattice gas density $\rho(x)$, as shown in Fig. 1. If we then choose a
time-step $\Delta t$ we can observe where particles ending up in cell $x$ came
from. The number of particles moving from cell $x-v_{i}$ to cell $x$ then
corresponds to the lattice gas occupation number $n_{i}(x,t)$. What is
important to note is that this resulting lattice gas model is fundamentally
correct in the sense that will obey the continuity and Navier-Stokes equations
simply because the molecular dynamics simulation does so.

The initial condition for a set of $N$ particles in a finite container with
periodic boundary conditions is given by their initial positions $x_{i}(0)$
and velocities $v_{i}(0)$. These particles then interact through an
interaction pair-potential that we take here to only depend on the distance
between the two particles: $V_{ij}=V(|x_{i}(t)-x_{j}(t)|)$. The MD simulation
then provides (too good accuracy) the trajectories $x_{i}(t)$ which solve
Newton’s second law:

$\displaystyle\frac{dx_{i}(t)}{dt}$ $\displaystyle=v_{i}(t)$ (8)
$\displaystyle\frac{dv_{i}(t)}{dt}$ $\displaystyle=-\frac{\partial}{\partial
x_{i}}\left(\frac{1}{2}\sum_{j\neq i}V_{ij}\right)$ (9)

We now superimpose a lattice onto the computational domain of the MD
simulation. For simplicity we can imagine a cubic lattice of lattice spacing
$\Delta x$, although any lattice will do here. Let us define a function that
determines if a particle resides in a specific cell of the reciprocal lattice
associated with a lattice point $x$:

$\Delta_{x}(x^{\prime})=\left\\{\begin{array}[]{lcl}1&&x_{\alpha}<x^{\prime}_{i\alpha}(t)\leq
x_{\alpha}+\Delta x\;\forall\alpha\in\\{x,y,z\\}\\\
0&&\mbox{otherwise}\end{array}\right.$ (10)

Next we pick a time step $\Delta t$. We can now determine the lattice
displacement each particle originally residing in a lattice point $x$
experiences. The set of all such displacements makes up the minimal set of
lattice velocities $v_{i}$ for our lattice gas method, and the number of
particles associated with this displacement makes up the lattice gas densities
$n_{i}(x,t)$. We define

$n_{i}(x,t)=\sum_{j}\Delta_{x}(x_{j}(t))\Delta_{x-v_{i}}(x_{j}(t-\Delta t)).$
(11)

This definition ensures that the particle numbers $n_{i}(x,t)$ will undergo a
streaming step given by Eq. (6). For any given MD simulation we then know all
$n_{i}(x,t)$. From Eq. (7) we see that the collision operator is then given by

$\Xi_{i}=n_{i}(x+v_{i},t+1)-n_{i}(x,t).$ (12)

This fully defines the MDLG algorithm, a lattice gas with a collision operator
that is defined through an underlying MD simulation. In some sense this is an
ideal lattice gas model that can handle even the most complex situations, i.e.
anything that can be addressed by MD, correctly. The key question is whether
this collsion operator can be reduced to some stochastic collision operator
that only depends on the local $n_{i}(x,t)$. Clearly this will only be the
case for very simple systems since the MDLB collision operator contains
information about temporal and spatial correllations of the underlying MD
algorithm and can, in principle, deal with many complex phenomena like liquid-
gas-solid coexistence, large varieties of transport parameters, including
phenomena at high Knudsen, high Mach, and/or high Reynolds numbers, which we
don’t expect to be accessible to a simple lattice gas algorithm of Eq. (7)
with a local collision operator. Such extensions will be subject of future
research, but are outside the scope of the current paper.

The local number of particles in lattice cell $x$ at time $t$ is give by

$N(x,t)=\sum_{j}\Delta_{x}(x_{j}(t)).$ (13)

This is consistent with the lattice gas definition of the local density
because

$\displaystyle N(x,t)$ $\displaystyle=\sum_{i}n_{i}(x,t)$ (14)
$\displaystyle=\sum_{i}\sum_{j}\Delta_{x}(x_{j}(t))\Delta_{x-v_{i}}(x_{j}(t-\Delta
t))$ (15) $\displaystyle=\sum_{j}\Delta_{x}(x_{j}(t)).$ (16)

The last equality follows because

$\sum_{i}\Delta_{x-v_{i}}(x_{j}(t-\Delta t))=1,$ (17)

i.e. every particle will be found somewhere on the lattice. Note that we have
not yet restricted the velocity set. We will use as many velocities as needed.
Mass conservation of Eq. (4) is clearly fulfilled since

$\displaystyle\sum_{i}\Xi_{i}$ (18) $\displaystyle=$
$\displaystyle\sum_{i}[n_{i}(x+v_{i},t+1)-n_{i}(x,t)]$ (19) $\displaystyle=$
$\displaystyle\sum_{i}\left\\{\sum_{j}\Delta_{x+v_{i}}[x_{j}(t+\Delta
t)]\Delta_{x}(x_{j}(t))\right.$
$\displaystyle\left.-\sum_{k}\Delta_{x}(x_{k}(t))\Delta_{x-v_{i}}(x_{k}(t-\Delta
t))]\right\\}$ (20) $\displaystyle=$ $\displaystyle\rho(x,t)-\rho(x,t)$ (21)
$\displaystyle=$ $\displaystyle 0$ (22)

| x=1 | x=2 | x=3
---|---|---|---
t=1 | $\mathbf{f_{0}=1}$ | $f_{0}=0$ | $f_{0}=0$
| $f_{1}=0$ | $f_{1}=0$ | $f_{1}=0$
t=2 | $\mathbf{f_{0}=1}$ | $f_{0}=0$ | $f_{0}=0$
| $f_{1}=0$ | $f_{1}=0$ | $f_{1}=0$
t=3 | $f_{0}=0$ | $f_{0}=0$ | $f_{0}=0$
| $f_{1}=0$ | $\mathbf{f_{1}=1}$ | $f_{1}=0$
t=4 | $f_{0}=0$ | $\mathbf{f_{0}=1}$ | $f_{0}=0$
| $f_{1}=0$ | $f_{1}=0$ | $f_{1}=0$
t=5 | $f_{0}=0$ | $\mathbf{f_{0}=1}$ | $f_{0}=0$
| $f_{1}=0$ | $f_{1}=0$ | $f_{1}=0$
t=6 | $f_{0}=0$ | $f_{0}=0$ | $f_{0}=0$
| $f_{1}=0$ | $f_{1}=0$ | $\mathbf{f_{1}=1}$

Figure 2: Simple thought experiment for lattice gas representation of a
particle moving with a constant velocity. Clearly the lattice gas definition
of momentum will fluctuate as a function of time even though the underlying MD
momentum is conserved. Averaging over all lattice placements, however, will
recover the correct momentum.

The definition of momentum in the lattice gas sense is typically defined as:

$N(x,t)u(x,t)=\sum_{i}n_{i}(x,t)v_{i}$ (23)

However, relating this to the underlying momentum of the MD simulation is not
exact, as can be seen in the example of a single MD particle moving with a
velocity less that is not a lattice velocity shown in Fig. 2. The
correspondence could be made exact if we were to introduced an average over
all possible placements of the lattice. Such an average would make no
difference to the global equilibrium distribution, which is the main focus of
the remaining paper. We therefore avoid this additional complication for the
current paper.

Similarly, momentum conservation of Eq. (5) is only exact if we introduce an
average over lattice placements:

$\displaystyle\sum_{i}v_{i\alpha}\Xi_{i}\neq 0\textbf{ in general}$ (24)

Of course this does not mean that there is a problem with momentum
conservation. Instead the problem arises due to the definition of momentum
through measured mass transfer between sites for a fixed lattice.

Despite the apparent lack of momentum conservation the MDLG collision rules
are still correct, even without the averaging, since the underlying MD
simulation respects momentum conservation. As such the apparent violation of
momentum conservation of the MDLG model is benign. We reserve a closer
examination of this averaged lattice gas implementation for a followup paper.

The key question is then whether the collision operator (Eq. (12)) can take
the form of Eq. (1), i.e. a stochastic collision operator that only depends on
the current local occupation numbers $n_{j}(x,t)$. Since the is a whole
ensemble of MD simulations that is consistent with a set of $n_{j}(x,t)$, and
these different MD simulations will lead to different collision terms, it is
clear that there can be no exact mapping. However, it is reasonable to hope
that we will be able to construct a stochastic lattice gas collision operator
that is statistically equivalent to the collision operators for the ensemble
of corresponding molecular dynamics simulations. Establishing this is not a
trivial task, and we will focus on the easier problem of showing that these
collision operators are consistent with the equilibrium behavior of the
lattice gas. In the next section we will present the lattice Boltzmann method
which conceptually represents the ensemble average of a lattice gas method.
Given the complexity of the task we focus in this paper on examining for
which, if any, discretizations the MDLG and the standard lattice Boltzmann
method give an equivalent equilibrium behavior.

##  IV MDLG for an two-dimensional Lennard-Jones gas

As a test case we use for our underlying MD simulations particles interacting
with the standard Lennard-Jones interaction potential, which is given by

$V(x)=4\epsilon\left[\left(\frac{\sigma}{x}\right)^{12}-\left(\frac{\sigma}{x}\right)^{6}\right].$
(25)

The interaction strength is controlled by $\epsilon$ and the spatial scaling
by $\sigma$. If $m$ is the mass of a particle we can construct the time-scale

$\tau=\sqrt{\frac{m\sigma^{2}}{\epsilon}}.$ (26)

We performed a molecular dynamics simulation using LAMMPS for $100\,000$
particles in a two-dimensional box with the length of $1000\sigma$,
corresponding to a nominal volume fraction of 0.0785 if we approximate the
particles as circles with diameter $\sigma$. The system was initialized with a
homogeneous distribution of particles with a kinetic energy corresponding to a
temperature of 50 in the LJ units defined above. This corresponds to a gas at
a high temperature, dense enough so that there are a significant number of
collisions. The temperature is well above the critical temperature for a
liquid-gas coexistence of $T_{c}=1.3120(7)$Potoff and Panagiotopoulos (1998).
We ran an simulation of an equilibrium system with a time-step of $0.0001\tau$
for $3\,000\,000$ time-steps, i.e. up to $\tau=300$. Early time data of
$1\,000\,000$ timesteps was discarded, to ensure that we were only probing the
equilibrium dynamics.

Figure 3: The numbering convention for the velocities $v_{i}$ in two
dimensions. The central point is $0$ and corresponds to velocity
$v_{0}=(0,0)$, and the other velocities are given by the connecting vector
between the central point and the lattice point in question.

We then analyse the resulting MD trajectories to obtain the resulting averaged
MDLG occupation numbers $n_{i}(x,t)$ from Eq. (11). We should note here that
results for different mean velocities $U$ can be obtained simply be using a
lattice displaced by $-U\;\Delta t$ for the MDLG analysis. It is therefore not
necessary to re-run the MD simulations to examine different mean velocities.

The first information to be gleaned from this is the resulting velocity set
for the $v_{i}$. For small times $\Delta t$, only the nearest neighbors have
non-negligible contributions, but as $\Delta t$ is increased more densities
get populated. We identify the velocities $v_{i}$ using the numbering scheme
shown in Fig. 3. So for $\Delta t\rightarrow 0$ only $v_{0}-v_{8}$ will have
contributions. These velocities form a complete shell around the central point
$v_{0}$. Most standard lattice Boltzmann methods work hard to make due with
this minimal velocity set. This comes at some cost, the most important one is
that only one temperature, $\theta=1/3$ in lattice unites, is allowable in
lattice units to recover the correct viscous stress tensor (see Eq. (57)). For
larger $\Delta t$ particle will travel further and a larger velocity set is
required. The average occupation numbers are given by the global equilibrium
distribution. The next subsection discusses how these equilibrium distribution
can be obtained analytically.

###  IV.1 Global Equilibrium Distribution

For a system in thermal equilibrium, sufficient averaging will give an
equilibrium distribution

$\displaystyle f_{i}^{eq}=$ $\displaystyle\langle n_{i}\rangle$ (27)
$\displaystyle=$
$\displaystyle\sum_{j}\langle\Delta_{x}(x_{j}(t))\Delta_{x-v_{i}}(x_{j}(t-\Delta
t))\rangle.$ (28)

We can numerically approximate this equilibrium density by averaging the
values of the $n_{i}$ from Eq. (11) over the whole lattice and for the
duration of the simulation. For a given MD simulation the results will depend
both on the lattice spacing $\Delta x$ and on the time step $\Delta t$. As
mentioned above the MD simulations considered in this paper deal with fairly
hot gases, that should be reasibably well approximated by an ideal gas.

To theoretically calculate this expectation value we assume that the particles
are uniformly distributed, so only the displacement during the time-interval
$\Delta t$ will enter the averaging:

$\delta x=x(t)-x(t-\Delta t).$ (29)

The key for our averaging will then be the probability of finding such a
displacement $P(\delta x)$, and this allows us to write the average as

$\displaystyle f_{i}^{eq}=$ $\displaystyle\frac{\rho^{eq}}{(\Delta x)^{d}}\int
dx\int d(\delta x)\Delta_{x}(x)\Delta_{x-v_{i}}(x-\delta x)P(\delta x)$
$\displaystyle=$ $\displaystyle\frac{\rho^{eq}}{(\Delta x)^{d}}\int dx\int
d(\delta x)\Delta_{x}(x)\Delta_{x}(x+v_{i}-\delta x)P(\delta x)$
$\displaystyle=$ $\displaystyle\rho^{eq}\int d(\delta x)P(\delta
x)W(v_{i}-\delta x)$ $\displaystyle=$ $\displaystyle\rho^{eq}\int d(\delta
x)P(\delta x+v_{i})W(\delta x)$ (30)

where $W(x)$ is the d-dimensional wedge function defined as

$\displaystyle W(x)=$ $\displaystyle\prod_{\alpha=1}^{d}W_{\alpha}(x)$ (31)
$\displaystyle=$
$\displaystyle\prod_{\alpha=1}^{d}\left(1-\frac{|x_{\alpha}|}{\Delta
x}\right)\Theta\left(1-\frac{|x_{\alpha}|}{\Delta x}\right)$ (32)

where $\Theta$ is the Heaviside function and $\alpha$ denotes cartesian
coordinate index.

For very short times $\Delta t$, i.e. times shorter than the mean free time
between two collisions, particles simply undergo ballistic motion. The
velocity distribution of the particles is given by the Maxwell-Boltzmann
distribution

$P(v)=\frac{1}{(2\pi
k_{B}T)^{d/2}}\exp\left(-\frac{(v-u)^{2}}{2k_{B}T}\right).$ (33)

With this, and neglecting any collisions between the particles, we get for the
mean squared displacement in one dimension

$\langle(\delta x_{\alpha})^{2}\rangle^{bal}=k_{B}T(\Delta t)^{2}.$ (34)

The probability distribution for the displacement is then given by

$P^{bal}(\delta x)=\frac{1}{(2\pi k_{B}T)^{d/2}(\Delta
t)^{d}}\exp\left(-\frac{\left(\delta x-u\Delta t\right)^{2}}{2k_{B}T(\Delta
t)^{2}}\right).$ (35)

For times much longer than the mean free time particles undergo multiple
collisions and instead of following a ballistic motion they will diffuse. If
we call the self-diffusion constant $D$ we when have

$\langle(\delta x_{\alpha})^{2}\rangle^{dif}=2D(\Delta t).$ (36)

This implies that the probability of the discplacement is given by

$P^{diff}(\delta x)=\frac{1}{(4\pi(\Delta
t)D)^{d/2}}\exp\left(-\frac{\left(\delta x-u\Delta t\right)^{2}}{4D(\Delta
t)}\right).$ (37)

Now since both limiting displacements are given by Gaussian distributions is
it reasonable to expect that the intermediate probabilities are also well
approximated by a Gaussian and, if we know the mean squared displacement in
one dimension $\langle(\delta x_{\alpha})^{2}\rangle$ (and assume isotropy),
we get for the probability

$P(\delta x)=\frac{1}{(2\pi\langle(\delta
x_{\alpha})^{2}\rangle)^{d/2}}\exp\left(-\frac{\left(\delta x-u\Delta
t\right)^{2}}{2\langle(\delta x_{\alpha})^{2}\rangle}\right).$ (38)

In all of these cases the probability distribution factorizes

$P(\delta x)=\prod_{\alpha=1}^{d}P_{\alpha}(\delta x)$ (39)

where

$P_{\alpha}(\delta x)=\frac{1}{\sqrt{2\pi\langle(\delta
x_{\alpha})^{2}\rangle}}\exp\left(-\frac{(\delta x_{\alpha}-u_{\alpha}\Delta
t)^{2}}{2\langle(\delta x_{\alpha})^{2}\rangle}\right)$ (40)

and we can write Eq. (30) as a product of Gaussian integrals:

$f_{i}^{eq}=\rho^{eq}\prod_{\alpha=1}^{d}\int d(\delta x_{\alpha})P(\delta
x_{\alpha}+v_{i\alpha})W_{\alpha}\left(\delta x\right).$ (41)

The solution is given by

$\frac{f_{i}^{eq}}{\rho^{eq}}=\prod_{\alpha=1}^{d}f_{i,\alpha}^{eq}$ (42)

where

$\displaystyle f_{i,\alpha}^{eq}=$ $\displaystyle
N\left(e^{-\frac{(u_{i,\alpha}-1)^{2}}{2a^{2}}}-2e^{-\frac{u_{i,\alpha}^{2}}{2a^{2}}}+e^{-\frac{(u_{i,\alpha}+1)^{2}}{2a^{2}}}\right)$
$\displaystyle+\frac{u_{i,\alpha}-1}{2}\left[\operatorname{erf}(\frac{u_{i,\alpha}-1}{\sqrt{2}a})-\operatorname{erf}(\frac{u_{i,\alpha}}{\sqrt{2}a})\right]$
$\displaystyle+\frac{u_{i,\alpha}+1}{2}\left[\operatorname{erf}(\frac{u_{i,\alpha}+1}{\sqrt{2}a})-\operatorname{erf}(\frac{u_{i,\alpha}}{\sqrt{2}a})\right]$
(43)

where

$\displaystyle a^{2}$ $\displaystyle=\frac{\langle(\delta
x)^{2}\rangle}{(\Delta x)^{2}}$ (44) $\displaystyle N$
$\displaystyle=\frac{a}{\sqrt{2\pi}}$ (45) $\displaystyle u_{i,\alpha}$
$\displaystyle=v_{i\alpha}-u_{\alpha}$ (46)

This is a lattice equilibrium distribution function derived from first
principles. At first glance it looks different than other lattice equilibrium
distributions, and we will examine its relation to know equilibrium
distribution functions below.

(a)

(b)

Figure 4: (a) Measured velocity correlation function from MD simulation data
compared to the exponential fit. (b) Measured mean square displacement from MD
simulation data compared to the predicted value according to the Eq. (47).

First we need to fully define the lattice equilibrium distribution. To do so
we need to obtain the mean square displacement $\langle(\delta
x_{\alpha})^{2}\rangle$. In general the mean square displacement can be
measured in our MD simulations, but this would require us to consider a whole
function of $\Delta t$ as in input parameter. For our simple ideal gas system
we can obtain a simpler dependence on a single parameter by expressing it in
terms of the velocity corellation function as

$\langle(\delta
x_{\alpha})^{2}\rangle=2\int_{0}^{t}dt^{\prime}\;(t-t^{\prime})\langle
v_{\alpha}(t^{\prime})v_{\alpha}(0)\rangle.$ (47)

For gases this velocity correlation function is typically well approximated by
a simple exponential decay. There is also a long range $1/t$ contribution to
the velocity correlation function for our two dimensional system, but for the
times $\Delta t$ that are of interest here this divergent contribution does
not yet contribute noticeably.

$\langle
v_{\alpha}(t)v_{\alpha}(0)\rangle=k_{B}T\exp\left(-\frac{t}{\tau}\right)$ (48)

For our system we compare this prediction of an exponentially decaying
velocity correlation function to the measured corellation function in Fig.
4(a). We see that for early times we see good agreement with this prediction
for $\tau=0.5(0)$. We also see the long time-tale typcical for a two
dimensional system, which we ignore here. This is justified below Uhlenbeck
and Ornstein (1930); Tsang and Tang (1977); Green (1952, 1954); Kubo (1957);
Weitz  _et al._ (1989).

Then the mean squared displacement can be predicted according to Eq. (47) as

$\displaystyle\langle(\delta
x_{\alpha})^{2}\rangle=2k_{B}T\tau^{2}\bigg{(}e^{-\frac{t}{\tau}}+\frac{t}{\tau}-1\bigg{)}.$
(49)

We show that this prediction recovers the measured mean squared displacement
well in Fig. 4(b). Deviations resulting from the long time tails of the
velocity corellation function only show up for later times and larger
displacements considered in this paper, which justifies our ignoring these
long time tails here.

This fully completes the definition of the MDLG equilibrium function in the
case of gases. To verify our results we compare a numerically measured
equilibrium distribution with the theoretically predicted one for different
discretizations. The results are shown in Fig. 5. The agreement between our
theoretical results and the experimental ones is excellent.

In the next section we introduce the lattice Boltzmann method and then examine
the relation of this MDLG equilibrium function with existing lattice
equilibrium distribution functions derived for lattice Boltzmann methods.

##  V Lattice Boltzmann

Lattice Boltzmann methods were derived as ensemble averages of lattice
Boltzmann methods. The variables in a lattice Boltzmann method are
distribution functions

$f_{i}=\langle n_{i}\rangle_{neq}$ (50)

where the $\langle\cdots\rangle_{neq}$ represents a non-equilibrium ensemble
average over microscopic lattice gas states. Taking the same ensemble average,
the evolution equation for these lattice Boltzmann densities derives from the
underlying lattice gas evolution Eq. (7)

$f_{i}(x+v_{i},t+1)=f_{i}(x,t)+\Omega_{i}$ (51)

where the collision operator $\Omega_{i}=\langle\Xi_{i}\rangle$ is a
deterministic function of all the densities at lattice point $x$. We will
investigate later if this collision operator can be, at least approximately
and for some suitable discretization, be cast in the standard BGK form
typically employed for lattice Boltzmann simulations. This question is a
crucial first step if one wants to relate lattice Boltzmann to an explicit
molecular system which will be represented by our MDLB algorithm.

This standard LB collision operator is a first order BGK approximation and can
be written as

$\Omega_{i}=\sum_{j}\Lambda_{ij}(f_{j}^{0}-f_{j})$ (52)

where $f_{j}^{0}$ is a local equilibrium distribution which depends only on
the locally conserved quantities

$\displaystyle\rho$ $\displaystyle=\sum_{i}f_{i}$ (53) $\displaystyle\rho
u_{\alpha}$ $\displaystyle=\sum_{i}v_{i\alpha}f_{i}.$ (54)

although other collision operators are also being used Ansumali and Karlin
(2002); Ansumali  _et al._ (2003); Ansumali and Karlin (2005); Geier  _et al._
(2015) and it is a longer term goal of the MDLG method to help identify which
of these collision operators are most realistic.

To ensure that the lattice Boltzmann equation reproduces the continuity and
Navier-Stokes equations in the hydrodynamic limit it is necessary that the
equilibrium distribution matches the first four (apart sometimes from a
$u^{3}$ term) velocity moments of the Maxwell-Boltzmann distribution:

$\displaystyle\sum_{i}f_{i}^{0}$ $\displaystyle=\rho$ (55)
$\displaystyle\sum_{i}(v_{i\alpha}-u_{\alpha})f_{i}^{0}$ $\displaystyle=0$
(56)
$\displaystyle\sum_{i}(v_{i\alpha}-u_{\alpha})(v_{i\beta}-u_{\beta})f_{i}^{0}$
$\displaystyle=\rho\theta\delta_{\alpha\beta}$ (57)
$\displaystyle\sum_{i}(v_{i\alpha}-u_{\alpha})(v_{i\beta}-u_{\beta})(v_{i\gamma}-u_{\gamma})f_{i}^{0}$
$\displaystyle=Q_{\alpha\beta\gamma}$ (58)

where $Q_{\alpha\beta\gamma}$ should be zero. For velocity sets including only
one shell we have $v_{i\alpha}\in\\{-1,0,1\\}$. For these velocity sets these
moments overconstrain the equilibrium distributions. In particular we have
$v_{i\alpha}^{3}=v_{i\alpha}$, which couples the first and the third moment.
This is a key source of Galilean invariance violations in lattice
BoltzmannWagner and Li (2006). These moments can only be reconciled for the
special choice of

$\theta=1/3$ (59)

and the third order tensor $Q_{\alpha\beta\gamma}=\rho
u_{\alpha}u_{\beta}u_{\gamma}\ll 1$ which is assumed to be small because
$u<0.1$ in typical situations. The equilibrium distribution is typically given
in terms of an expansion in terms of the local velocity $u$ up to second order
as

$f_{i}^{0}=\rho
w_{i}\left(1+\frac{v_{i\alpha}u_{\alpha}}{\theta}+\frac{1}{2}\frac{v_{i\alpha}u_{\alpha}v_{i\beta}u_{\beta}}{\theta^{2}}-\frac{1}{2}\frac{u_{\alpha}u_{\alpha}}{\theta}\right)$
(60)

where the weights $w_{i}$ depend on the velocity set and summation over
repeated Greek indices is implied. In this article we focus on the question
whether this form of an equilibrium distribution is compatible with a concrete
MDLB implementation.

This collision operator together with the local equilibrium distribution
implies mass and momentum conservation

$\displaystyle\sum_{i}\Omega_{i}$ $\displaystyle=0,$ (61)
$\displaystyle\sum_{i}v_{i\alpha}\Omega_{i}$ $\displaystyle=0.$ (62)

which is consistent with the typical conditions for lattice gases of Eqs. (4)
and (5).

In the following we will examine the MDLG method for the example of a hot,
dilute gas. For this lattice gas we examine the resulting distribution
functions and see under which circumstances this lattice gas can reproduce (to
some approximation) the lattice Boltzmann method equilibrium distribution Eq.
(60).

##  VI Relation of MDLG equilibrium functions to lattice Boltzmann equilibria

We are now in a position to predict for which set of parameters $\Delta
x,\Delta t$, if any, we can recover the traditional form of the lattice
Boltzmann equilibrium from our MDLB algorithm. Most lattice Boltzmann methods
use a limited velocity set that corresponds to a single shell in Fig. 3. For
our two dimensions this corresponds to 9 velocities. The corresponding
equilibrium distribution is typically given as a second order polynomial in
the velocities, as we have presented earlier in Eq. (60). For the two
dimensional D2Q9 lattice Boltzmann method we consider here the weights $w_{i}$
in Eq. (60) are given by

$\displaystyle w_{0}$ $\displaystyle=4/9$ (63) $\displaystyle w_{1-4}$
$\displaystyle=1/9$ (64) $\displaystyle w_{5-8}$ $\displaystyle=1/36$ (65)

where the velocity indices correspond to the numbering of Fig. 3.

Figure 5: Measured equilibrium distributions $f_{i}^{eq}$ as a function of the
mean-squared displacement measure a from Eq. (44). They are compared to the
analytic solution given by Eq. (42). We find excellent agreement between the
predicted and measured equilibrium distributions. The horizontal lines
indicate the value of D2Q9 lattice Boltzmann weights, and the green vertical
line indicates the value of $a^{2}=1/6$ for which these weights agree with the
MDLG results.

In Fig. 5 we show the lattice Boltzmann weights $w_{i}$ as horizontal lines.
To match the MDLG and LB equilibria we require

$f_{i}^{eq}(a^{2})/\rho^{eq}=w_{i}.$ (66)

This is a over-determined system of equations.

Fortuitously the solutions for the three distinct $w_{i}$ for a D2Q9 lattice
Boltzmann give the same value for $a^{2}\approx 1/6$ to very good
approximation. This is shown as the green vertical line in Fig. 5.

Figure 6: Pairs of spatial and temporal discretizations $(\Delta x,\Delta t)$
that lead to equilibrium distributions of MDLG that are consistent with the
lattice Boltzmann equilibrium distribution $f_{i}^{0}(\rho,0)=w_{i}$. Notice
the two scaling regimes of $\Delta t\propto\Delta x$ for the ballistic regime
for small times and the diffusive regime $(\Delta t)^{2}\propto\Delta x$ for
large times.

This suggests that matching lattice Boltzmann and MD simulations would likely
benefit from using the conditions where the $f_{i}^{eq}$ match up, and the
methodology explained above would give guidance on the appropriate time step
$\Delta t$ for a given spatial discretization $\Delta x$. Given a $\Delta x$
can numerically solve Eq. (66) for $\Delta t$ for a system with zero mean
momentum. This is shown in Fig. 6. We find that there is close agreement
between the solutions for different velocities $v_{i}$. Corresponding to the
transition from ballistic to diffusive regime around $\Delta t=\tau$ that we
saw in Fig. 4(b) we also see a transition here from a $\Delta t\propto t$
regime for $\Delta t\ll 1$ to a $\Delta t\propto t^{2}$ regime for $\Delta
t\gg 1$. We expect this relation that gives $\Delta t$ in terms of $\Delta x$
to be valuable when one tries to generate a coarse-graining transition between
an MD and LG region in a multi-scale numerical method.

Figure 7: Dependence of the equilibrium distribution function on an imposed
velocity $U$ for $\Delta x=100\;\sigma$ and $\Delta t=34.16\;\tau$,
corresponding to the green line in Fig. 5b. These results are compared to the
lattice Boltzmann equilibrium distribution for D2Q9. We find good agreement
for velocities below about 0.2.

So far we have only matched the equilibrium distribution at zero velocity. The
theory contains the mean velocity as $u$ in Eq. (43). For the measurements we
could set up different simulations for mean velocities using the algorithm
described above in Eq. (11). It is more practical, however, to move the Grid
instead or equivalently use Galilean transformed particle positions of
$\hat{x}_{i}(t)=x_{i}(t)+ut$ instead. Using this approach to find the
equilibrium distribution for different mean velocities $u$ we show our
comparison between the measured discrete equilibria (Eq. (11)), their
theoretical prediction (Eq. (42)), and the lattice Boltzmann equilibrium
distribution (Eq. (60)) in Fig. 7. For small velocities $|u|<0.1$ we find good
agreement between all three densities. This is the relevant range, as lattice
Boltzmann is only considered reliable for small enough velocities. The
agreement between the measured and predicted MDLG equilibrium distributions
continues to be excellent throughout the whole regime. Note, that the
agreement would continue to be excellent for larger velocities. We would only
need to adapt the velocity set we consider as velocities with magnitude larger
than 0.5 are the same as velocities with magnitude larger than 0.5 plus an
additional integer lattice displacement.

###  VI.1 Moments of the equilibrium distribution

The key property of the equilibrium distribution in kinetic theory are the
velocity moments. For the derivation of the Navier-Stokes equation moments up
to third order are required. It is therefore helpful to examine the moments of
the discrete MDLG equilibrium distribution, compare them to the expected
moment for a lattice Boltzmann equilibrium distribution and examine how these
moments relate to the continuous velocity distribution function for the MD
simulation.

Let us spend a moment considering these different concepts, since they are not
usually clearly separated in a LB derivation. Firstly we have the velocity
distribution of the MD simulation, given by Eq. (33). The moments of this
velocity distribution are usually used as a rational for constructing LB
equilibrium distributions such that the relevant moments of the discrete LB
equilibrium distribution function match those of the continuous Maxwell
Boltzmann distribution. This is the rational behind the moment Eqs. (55–58)
for the lattice Boltzmann equilibrium function.

The moments of the discrete MDLG method are not a priori constrained to obey
such a constraint, and indeed we don’t expect such an agreement for two
reasons. First the underlying displacement probabilities (Eq. (38)) are in
general different from the displacement probabilities (Eq. (35)) directly
related to the Maxwell-Boltzmann distribution (Eq. (33)). Second there is no
reason to believe that the averaging procedure of Eq. (41) will preserve the
moments in general.

Now let us consider the first three velocity moments specifically. The zeroth
moment relates to the total mass. Since the algorithm conserves mass exactly
we expect all three approaches to agree on this moment. Indeed, as we saw in
Eq. (22) mass is clearly conserved, and consequently this moment will agree
for all of the above approach.

The first moment relates to the local momentum. Even if the MDLG approach does
not locally conserve the momentum, the averaged momentum of the equilibrium
distribution remains exact. This simply follows from the fact that this
discrete moment corresponds to the net mass flow through the lattice. Even
though this flow can be inexact at any instance in time (since particles may
not cross a boundary despite the fact that they are moving) on average the
count of particles crossing boundaries has to give the exact mass current.

Figure 8: The normalized second moment of Eq. (67) as a function of the second
moment $a^{2}$ from Eq. (44) of the mean squared displacement probability. The
value of $\Psi$ is shown for different 11 different values of $u_{x}$ between
0 and 0.5. The values for $u_{x}=0$ correspond to the bottom points of the
graph. For a Galilean-Invariant discrete equilibrium distribution moment the
value of $\Psi$ should be independent of $u$. This is compared to a standard
lattice Boltzmann second moment with $\theta=1/3$.

Let us next consider the second moment. This second moment in the lattice
Boltzmann approach (Eq. (57)) is related ideal gas equation of state
$p=\rho\theta$. We calculate the second moment of our discrete MDLG
equilibrium distribution

$\Psi_{\alpha\beta}(a,u)=\frac{\sum_{i}f_{i}^{eq}(v_{i\alpha}-u_{\alpha})(v_{i\beta}-u_{\beta})}{\rho^{eq}a^{2}}.$
(67)

For an equilibrium distribution that obeys the lattice Boltzmann moment (Eq.
(57)) with a temperature $\theta=a^{2}$ this expression would give exactly
one. For the MDLG equilibrium distribution this second moment is shown in Fig.
8. This shows that the discrete second moment does only agree with the MD
temperature for large $a^{2}\gg 1$. As we see in Fig. 5 this corresponds to a
situation where the populated set of velocities encompasses several shells in
Fig. 3. For lower values of $a^{2}$ we find that the second moment diverges.
The reason lies in the way we define the discrete equilibrium distribution.
Even for very small $\langle(\delta x)^{2}\rangle$, corresponding to very
small $\Delta t$, we will identify a fraction of particles that happen to
cross from one lattice point to the next and are therefore assigned a lattice
velocity $v_{i}$ of order one. This appearance of apparent large displacements
causes the divergence of the discrete second moment. This effect is
significantly enhanced by imposed velocites $u$.

For $a^{2}<0.1$ we see that we get significantly diverging values of $\Psi$
for different velocities $u$. This implies that this discrete second moment is
not Galilean invariant. It is important to note that despite such violations
of Galilean invariance of the discrete moments the full MDLG algorithm does
not suffer from a Galilean invariance problem. Instead this is an indication
that the collision operator Eq. (12) must exactly compensate this apparent
Galilean invariance violation. For lattice Boltzmann methods one typically
tries to avoid Galilean invariance violations by ensuring that both the local
equilibrium distributions and the collision operator independently obey
Galilean invariance.

For the lattice Boltzmann method we expect this second moment to be
$\Psi_{\alpha\beta}=\theta/a^{2}\delta_{\alpha\beta}$. As we saw above [see
Eq. (59)]lattice Boltzmann methods which have a velocity set consisting of a
single shell in velocity space require $\theta=1/3$. This value of $\theta$ is
consistent with the moment of the MDLG equilibrium for a value of
$a^{2}\approx 1/6$. We find $\theta\approx 2a^{2}=\langle(\delta
x)^{2}\rangle/(\Delta x)^{2}$. This corresponds to about the lowest value for
$a^{2}$ where $\Psi$ does not strongly depend on $u$ and would therefore
violate Galilean invariance.

Figure 9: Third velocity moment given by Eq. (68) for an unrestricted velocity
set and a velocity set restricted to only a single shell of velocities for 10
equally spaced velocities $u_{x}$ between 0 and 0.001.

The condition that $\theta=1/3$ came out of a consideration for the third
moment for a minimal velocity set with $v_{ix}=v_{ix}^{3}$. We can define a
third moment as

$\Phi_{\alpha\beta\gamma}=\frac{\sum_{i}f_{i}^{eq}(v_{i\alpha}-u_{\alpha})(v_{i\beta}-u_{\beta})(v_{i\gamma}-u_{\gamma})}{\rho^{eq}a^{3}}$
(68)

We examine the behavior of $\Phi_{xxx}$ in Fig. 9. In fact this third moment
should be zero, and for sufficiently large $a$ it converges to zero
exponentially. However, if we artificially restrict our velocity set to a
single shell, neglecting the small densities for discrete velocities outside
the first shell, we find that there is a collapse of $\Phi_{xxx}$ to zero for
the same $a^{2}\approx 1/6$ that we found in Fig. 8. For the full velocity
set, however, nothing special occurs at $a^{2}=1/6$. This initially surprised
us because we see in Fig. 5 that the densities associated with the second
shell are still a factor of 30 smaller than the densities of the first shell
and might therefore be taken to be negligible. The second shell velocities
that are about a factor 2 larger than the first shell velocities, so these
densities get multiplied by a factor of $2^{3}$ which allows these densities
to contribute enough that we now have

$\sum_{i}f_{i}^{eq}v_{ix}\neq\sum_{i}f_{i}^{eq}v_{ix}^{3}$ (69)

and they are different enough that the cancellation of terms that is supposed
to lead to for $\Phi_{xxx}=0$ at $a^{2}=1/6$ no longer exists. Instead
$\Phi_{xxx}$ monotonously approaches zero.

Figure 10: Two manifestations of the temperature $\theta$ from the second and
third moments, given by Eqs. (70) and (71). We show the moments for the full
velocity set and moments for a velocity set restricted to the first shell
(thinner lines). For the restricted velocity set we only see agreement between
the two moments for $\theta=1/3$ corresponding to $a^{2}=1/6$.

The second and third moments of the distribution functions both relate to the
lattice Boltzmann temperature of Eq. (57). We can derive two expressions for
the temperature:

$\displaystyle\theta_{2}$ $\displaystyle=a^{2}\Psi_{xx}$ (70)
$\displaystyle\theta_{3}$
$\displaystyle=\frac{1}{3\rho^{eq}}\partial_{u_{x}}\sum_{i}f_{i}^{eq}(\rho^{eq},u)(v_{ix})^{3}$
(71)

The dependence of these two quantities on the mean-square displacement is
shown in Fig. 10. We see that these two definitions only agree with each other
for large $a^{2}$, where agreement is facilitated by utilizing a large set
velocity for the discretization of $f_{i}^{eq}$. For small $a^{2}$, and
therefore a minimal velocity set, we get $\theta_{3}=1/3$, for the same reason
that was discussed before, i.e. $v_{ix}=v_{ix}^{3}$.

If we artificially restrict the velocity set to a single shell (e.g. 9
velocities in our two-dimensional example) we affect both definitions of the
LB temperature. In this case the agreement for large $a^{2}$ disappears, and
we have exactly one point for $\theta=1/3$ corresponding to $a^{2}=1/6$ for
which both expressions for the temperature agree. This corresponds to the
special point in Fig. 9 where the results for the one-shell velocity set
become independent of $u$ (for small enough $u$).

This suggests that there is a serendipitous agreement between the MDLG
equilibrium distribution and the standard LB equilibrium distribution for one-
shell velocity sets. It allows us to recover velocity moments up to order 3,
which is exactly the order required by kinetic theory to recover the
continuity and Navier-Stokes equations. As seen in Fig. 8 and 9, this value is
just large enough to avoid apparent Galilean-invariance violations in the
moments of the MDLG equilibrium distribution, and just small enough so that
the values of the $f_{i}$ for velocities with $v_{i\alpha}>1$ are small enough
not to contribute considerably to the moments of the equilibrium
distributions. These effects can only be noticed for the third moment, where
they conspire to induce agreement between the measures of temperature implied
by the second and third moments, as seen in Fig. 10.

##  VII Consistent discretizations

Up to now we have seen how the moments of the equilibrium distribution depend
on $a^{2}$ which is a measure of the mean squared displacement. We saw that
the discrete moments of $f_{i}^{eq}$ differ from the continuous moments of the
Maxwell-Boltzmann distribution even in the ballistic regime. The underlying
reason for this disagreements result from two conspiring effects. Firstly we
only know the position of our particles to lie somewhere within their assigned
lattice cell. This uncertainty enters the definition of the $f_{i}^{eq}$ in
Eq. (30). Secondly we use a discrete second moment. Here we show how both of
these give an offset of $1/12$ in $a^{2}$ giving rise to a total shift of
$1/6$ observed in the previous numerical results.

For simplicity let us consider large times away from the ballistic regime.
This occurs without loss of generality if our assumption of a Gaussian
distribution in Eq. (38) remains correct. We can then assume the motion of the
particles is entirely diffusive with some diffusion constant $D$:

$\partial_{t}\rho(x,t)=D\partial_{x}^{2}\rho(x,t)$ (72)

If we had known the position of the particle initially, i.e.
$\rho_{c}(x,0)=\delta(x)$, then the particle probability density would evolve
as

$\rho_{c}(x,t)=\frac{1}{\sqrt{4\pi Dt}}e^{-\frac{x^{2}}{4Dt}}$ (73)

which has a second moment of

$\frac{1}{\sqrt{4\pi
Dt}}\int_{-\infty}^{\infty}x^{2}e^{-\frac{x^{2}}{4Dt}}dx=2Dt$ (74)

If we only know that at time $t=0$ the particle is inside a lattice cell
centered around the origin with width $\Delta x$, then the density $\rho$ as a
function of time is given by

$\rho_{d}(x,t)=\frac{1}{2\Delta x}\left[\mbox{erf}\left(\frac{x+\Delta
x/2}{2\sqrt{Dt}}\right)-\mbox{erf}\left(\frac{x-\Delta
x/2}{2\sqrt{Dt}}\right)\right]$ (75)

which for $t\rightarrow 0$ gives a density that is $1/\Delta x$ inside the
interval $[-\Delta x,\Delta x]$ and zero outside it. This probability
distribution then spreads out and approaches a Gaussian at late times. The
second moment of the position is then given by

$\int_{-\infty}^{\infty}x^{2}\rho_{d}(x,t)dx=2Dt+\frac{(\Delta x)^{2}}{12}$
(76)

i.e. a simple offsett of $(\Delta x)^{2}/12$ that does not depend on time. We
can identify

$a^{2}=\frac{2Dt}{(\Delta x)^{2}}$ (77)

from the definition of Eq. (44). We then we see that the results are expected
to be shifted by $1/12$. However, in Fig. 10 we clearly see that this only
represents half of the observed shift $1/6$.

Figure 11: The difference of the discrete second moment $m_{2}$ of Eq. (78)
minus $a^{2}$ converges quickly to a constant $1/6$.

The second effect relates to discretizing the position into displacement bins.
We now calculate the discrete second moment (normalized by $(\Delta x)^{2}$)as

$m_{2}=\sum_{i=-\infty}^{\infty}i^{2}\int_{i\Delta x-\frac{\Delta
x}{2}}^{i\Delta x+\frac{\Delta x}{2}}\rho_{d}(x,t)\;dx$ (78)

In Fig. 11 we show that this discrete moment quickly converges to $a^{2}+1/6$.
We clearly see that the missing additional offset of $1/12$ that we observed
in Fig. 10 is the result of taking the discrete moment.

This shows that there are two effects of discretization. Firstly the initially
broader distribution of particles confined to a lattice site rather than a
point shifts the second moment by exactly 1/12. The second effect of
discretization is more complicated, particularly for small $a^{2}$. But for
$a^{2}>0.2$ this effect quickly converges to another offset in $a^{2}$ of
$1/12$. Together they make up the offset of $1/6$, seen repeatedly in our
numerical results of Figs. 5, 8, 9, and 10.

##  VIII Outlook

In this paper we have introduced a new tool for comparing the results of
Molecular Dynamics simulations with those of coarse-grained lattice gas or
lattice Boltzmann methods. It consists of re-interpreting the MD results as a
lattice gas which we call the MDLG. The dynamics of this special kind of
lattice gas is entirely given by the MD simulation, and therefore will be able
to give a coarse-grained picture for any results that are obtainable with MD
simulations.

Here we focus entirely on the averaged equilibrium behavior and show that
there exists a close connection between the equilibria of lattice Boltzmann
methods and the equilibrium for the MDLG method, when applied to a hot dilute
gas. We were able to determine this equilibrium distribution analytically and
were able to verify this analytical solution with the results of the MDLG
method. Importantly there is a surprisingly good agreement between our
equilibrium distribution and the standard lattice Boltzmann result for
carefully chosen (and analytically known) pairs of time and space
discretizations $\Delta t$ and $\Delta x$. We were able to understand the
observed offset of $1/6$ in the dimensionless measure of the mean squared
displacement $a^{2}$ in terms of our discretization procedure.

This opens the way for a more careful analysis of the fundamental
underpinnings of lattice gas and lattice Boltzmann methods. We intend to
utilize our MDLG method to investigate the fundamental properties of the
collision operator, including its fluctuating properties. Further down the
line we hope to investigate how the behavior of liquids alters the behavior of
the MDLG method and examine if MDLG can also be matched with lattice Boltzmann
methods. We anticipate that this method will also be instrumental in putting
lattice Boltzmann methods for non-ideal and multi-component systems on a
firmer footing.

## References

  * Succi (2001) S. Succi,  _The lattice Boltzmann equation: for fluid dynamics and beyond_ (Oxford university press, 2001).
  * Teixeira (1998) C. M. Teixeira, International Journal of Modern Physics C 9, 1159 (1998).
  * Yu  _et al._ (2005) H. Yu, S. S. Girimaji, and L.-S. Luo, Journal of Computational Physics 209, 599 (2005).
  * Swift  _et al._ (1995) M. R. Swift, W. Osborn, and J. Yeomans, Physical Review Letters 75, 830 (1995).
  * Shan and Doolen (1995) X. Shan and G. Doolen, Journal of Statistical Physics 81, 379 (1995).
  * He and Luo (1997) X. He and L.-S. Luo, Physical Review E 56, 6811 (1997).
  * Briant  _et al._ (2004) A. Briant, A. Wagner, and J. Yeomans, Physical Review E 69, 031602 (2004).
  * Yuan and Schaefer (2006) P. Yuan and L. Schaefer, Physics of Fluids 18, 042101 (2006).
  * Tölke (2002) J. Tölke, Philosophical Transactions of the Royal Society of London A: Mathematical, Physical and Engineering Sciences 360, 535 (2002).
  * Kang  _et al._ (2006) Q. Kang, P. C. Lichtner, and D. Zhang, Journal of Geophysical Research: Solid Earth 111 (2006).
  * Dupin  _et al._ (2003) M. Dupin, I. Halliday, and C. Care, Journal of Physics A: Mathematical and General 36, 8517 (2003).
  * Jansen and Harting (2011) F. Jansen and J. Harting, Physical Review E 83, 046707 (2011).
  * McNamara and Zanetti (1988) G. R. McNamara and G. Zanetti, Physical review letters 61, 2332 (1988).
  * Frisch  _et al._ (1986) U. Frisch, B. Hasslacher, and Y. Pomeau, Physical review letters 56, 1505 (1986).
  * Higuera and Jimenez (1989) F. Higuera and J. Jimenez, EPL (Europhysics Letters) 9, 663 (1989).
  * Higuera  _et al._ (1989) F. Higuera, S. Succi, and R. Benzi, EPL (Europhysics Letters) 9, 345 (1989).
  * Qian  _et al._ (1992) Y. Qian, D. d’Humières, and P. Lallemand, EPL (Europhysics Letters) 17, 479 (1992).
  * Potoff and Panagiotopoulos (1998) J. J. Potoff and A. Z. Panagiotopoulos, The Journal of chemical physics 109, 10914 (1998).
  * Uhlenbeck and Ornstein (1930) G. E. Uhlenbeck and L. S. Ornstein, Physical review 36, 823 (1930).
  * Tsang and Tang (1977) T. Tsang and H. Tang, Physical Review A 15, 1696 (1977).
  * Green (1952) M. S. Green, The Journal of Chemical Physics 20, 1281 (1952).
  * Green (1954) M. S. Green, The Journal of Chemical Physics 22, 398 (1954).
  * Kubo (1957) R. Kubo, Journal of the Physical Society of Japan 12, 570 (1957).
  * Weitz  _et al._ (1989) D. Weitz, D. Pine, P. Pusey, and R. Tough, Physical review letters 63, 1747 (1989).
  * Ansumali and Karlin (2002) S. Ansumali and I. V. Karlin, Physical Review E 65, 056312 (2002).
  * Ansumali  _et al._ (2003) S. Ansumali, I. V. Karlin, and H. C. Öttinger, EPL (Europhysics Letters) 63, 798 (2003).
  * Ansumali and Karlin (2005) S. Ansumali and I. V. Karlin, Physical review letters 95, 260605 (2005).
  * Geier  _et al._ (2015) M. Geier, M. Schönherr, A. Pasquali, and M. Krafczyk, Computers & Mathematics with Applications 70, 507 (2015).
  * Wagner and Li (2006) A. Wagner and Q. Li, Physica A: Statistical Mechanics and its Applications 362, 105 (2006).


††thanks: These two authors contribute equally.††thanks: These two authors
contribute equally.††thanks: Correspondence to<EMAIL_ADDRESS><EMAIL_ADDRESS>
# Moiré Semiconductors on Twisted Bilayer Dice Lattice

Di Ma1 Yu-Ge Chen2† Yue Yu1 Xi Luo3 1\. Department of Physics, Fudan
University, Shanghai 200433,China
2\. Beijing National Laboratory for Condensed Matter Physics and Institute of
Physics, Chinese Academy of Sciences, Beijing 100190, China
3\. College of Science, University of Shanghai for Science and Technology,
Shanghai 200093, China

###### Abstract

We propose an effective lattice model for the Moiré structure of twisted
bilayer dice lattice. We find that there are flat bands near zero energy level
at any twist angle besides the magic ones. The flat bands contain both bands
with zero Chern number which are originated from the destructive interference
of the dice lattice and the topological non-trivial ones at the magic angle.
The existence of the flat bands can be detected from the peak-splitting fine
structure of the optical conductance at all angles, while the transition peaks
do not split and only occur at magic angles in twisted bilayer graphene.

##  I  Introduction.

The search for flat-band system has become one of the new trends over the last
few decades. Due to the large effective mass of quasi-particles, the density
of states (DOS) is large, and the kinetic energy of the carriers is strongly
quenched. Therefore, the flat-band system is a good candidate for studying
strongly correlated electronic states induced by strong Coulomb interaction,
such as ferromagnetism JOP1991 ; PRA2010 ; PRL1992 , heavy fermion nature2020
; np2019 , fractional Chern insulator IJMP2013 ; PRL2011 , Wigner crystal
PRL2007 ; PRB2008 , and unconventional superconductivity RMP1990 ; PCS2007 ;
IOS2020 .

Traditionally, the nearly flat-band system can be achieved by invoking fine-
tuned nearest neighbor hoppings, long-ranged hoppings, or by breaking time-
reversal symmetry. Several lattice models have been proposed along this line
in kagomé PTP1951 , Lieb PRL1989 , and dice lattices PRB1986 . The existence
of the flat-band is guaranteed by destructive interference of the Wannier
functions of lattice structure, and the flat-band states are identified as
compact localized states PRL2007 ; PRB2008 ; PRB1986 ; PRL2013 ; CPB2014 .
This destructive interference protection can also be generalized to lattices
with mirror symmetry chen2022 . Usually, this flat-band has zero Chern number
in the lattice model with nearest neighbor hopping only. On the dice lattice,
the flat-band can acquire a non-zero Chern number by invoking Rashba spin-
orbital coupling, and exhibits anomalous quantum Hall effect by adding onsite
Hubbard interactions, which could be realized in the transition-metal oxide
SrTiO3/SrIrO3/SrTiO3 trilayer heterostructure by growing along the (111)
direction wang2011 . Materials with flat band structures along this line are
also reported in Cu(111) confined by CO molecules np2017 , optical lattice and
cold atom systems shen2010 ; njp2014 ; prl2015a ; prl2015b ; sa2015 ; ol2016 .
In 3D, a famous example is the Kane semimetal, where the flat band structure
is associated with the triplet degenerate nodes and can be described by a 3D
Lieb lattice model Luo2018 . The low energy quasi-particle, the Kane fermion,
can be viewed as a fermionic photon, i.e., a spin-1 fermion, and the flat band
of Kane fermion corresponds to the ”longitudinal mode” of the photon.
Experimentally, the Kane fermion has been reported in Hg1-x CdxTe hgcdte  and
Cd3As2 cdas . The existence of the flat band is shown in the optical
conductance by the large peaks near zero frequency hgcdte ; Luo2019 .

Recently, a new mechanism for generating flat band is found in the twisted
bilayer graphene (TBG) at a magic twist angle $\theta\sim 1.08^{\degree}$
Bistritzer-MacDonald2011 ; mac2020 . Unlike the destructive interference
induced flat band, the flat band structure of TBG originates form the
extremely large band folding of the Moiré structure, and TBG becomes a
strongly-correlated electron system. Very soon, superconductivity is reported
in TBG cao2018 ; bernevig2018 . And the flat band in TBG has a non-trivial
Chern number bernevig2020 , which can be explained by the zeroth chiral Landau
levels of Dirac/Weyl fermions Liu2019Pseudo . Away from the half-filling, the
fractional Chern insulator phase is also proposed and reported in TBG
Tarnopolsky2019 ; ashvin2020 ; ashvin2021 , and the twisted bilayer MoTe2
xu2023 ; xu2023b ; shan2023 .

In this paper, we consider the combination of Moiré structure and destructive
interference induced flat band, namely, a twisted bilayer dice lattice (TBD).
We construct a lattice model for TBD in the reciprocal space. We find that
there are flat bands near zero energy at any twist angle besides the magic
ones. The flat bands are contributed from the ones with zero Chern number as
in the Dice lattice and the non-trivial ones at the magic angle, therefore,
TBD is a playground for studying the interplay between the zero Chern number
flat-bands and non-trivial ones. We further confirm this scenario by
considering the chiral limit of the pseudo Landau level description
Liu2019Pseudo .

The paper is organized as follows. In Sec. II, we provide the detailed
construction of the lattice model of TBD. In Sec. III, we calculate the Bloch
band structure numerically, and compare it with that of TBG. There are flat
bands in TBD at all angles other than the magic ones. In Sec. IV, we use the
Landau level language to describe the physics of the flat bands in TBD, where
the effective magnetic field is caused by the interlayer hopping. One can
directly find that, besides the topological zeroth Landau level, the higher
Landau levels which are topological trivial also contribute to the flat bands.
In Sec. V, we use the degenerate perturbation method to calculate the optical
conductance of TBD, and we find that the flat bands contribute a peak
splitting fine structure, which is a smoking-gun experimental prediction for
the existence of the flat bands. This phenomenon also exists at all angles
besides the magic ones when comparing with TBG. And the final section is
devoted to conclusions and discussions.

##  II  Effective model of TBD

The dice lattice can be viewed as two honeycomb lattices (A-B and B-C) sharing
a same sublattice site (B), See Fig. 1. For simplicity, we consider a spinless
lattice model with nearest neighbor hopping only. Due to this similarity, the
TBD has the same Moiré structure as TBG. There is a flat band $E=0$ in the
lattice spectrum. The low energy behavior near the Dirac point can be captured
by $H_{eff}={\bf k\cdot S}$, where ${\bf k}=(k_{1},k_{2},0)$ is the lattice
momentum, and ${\bf S}$ is the spin-1 generalization of Pauli matrix that acts
on sublattice space.

$\displaystyle S_{1}=\begin{pmatrix}0&1&0\\\ 1&0&1\\\
0&1&0\end{pmatrix},S_{2}=\begin{pmatrix}0&-i&0\\\ i&0&-i\\\
0&i&0\end{pmatrix},S_{3}=\begin{pmatrix}1&0&0\\\ 0&0&0\\\
0&0&-1\end{pmatrix}.$ (A1)

(a)

(b)

Figure 1: (color online) (a): Diagrams show the structure of dice lattice,
blue dot indicate $A$ atoms, red dot $B$ indicate atoms and $C$ atoms for
black dot. The lattice vectors described as
$\bm{a_{1}}=a(\frac{1}{2},\frac{\sqrt{3}}{2}),\bm{a_{1}}=a(-\frac{1}{2},\frac{\sqrt{3}}{2})$,
and corresponding reciprocal vectors
$\bm{b_{1}}=\frac{4\pi}{\sqrt{3}a}(\frac{\sqrt{3}}{2},\frac{1}{2}),\bm{b_{2}}=\frac{4\pi}{\sqrt{3}a}(-\frac{\sqrt{3}}{2},\frac{1}{2})$,
where $a=\sqrt{3}d$ is lattice constant. (b): Bilayer dice lattice in A-B
stacking (Bernal). The hopping between A and C atoms in different layer are
forbidden.

Since dice lattice can be realized in SrTiO3/SrIrO3/SrTiO3 trilayer
heterostructure by growing along the (111) direction wang2011 , the TBD is
possible to be realized in this material by twisting when growing.

For TBG, Bistritzer and MacDonald have proposed a low energy effective
continuum Dirac model of the Moiré structure for small twist angle
$\theta<10^{\degree}$. The effective model consists two isolated graphene
layers and hopping terms between them. With this model they reveal flat Bloch
bands in the electric structure at magic twist angles which give rise to large
DOS Bistritzer-MacDonald2011 . Similar to the case of TBG, we follow Ref.
Bistritzer-MacDonald2011  to constructed a continuum model for TBD. In TBD, we
keep the top layer (Layer 1) fixed, and rotate the bottom layer (Layer 2) by
$\theta$ with respect to Layer 1. The effective TBD Hamiltonian contains
intralayer and interlayer parts. The low energy intra Hamiltonian reads

$\displaystyle H=v_{f}\sum_{k}\psi_{1,\mathbf{k}}^{\dagger}\mathbf{S\cdot
k}\psi_{1,\mathbf{k}}+v_{f}\sum_{k}\psi_{2,\mathbf{k}}^{\dagger}\mathbf{S_{\theta}\cdot
k}\psi_{2,\mathbf{k}},$ (A2)

where $S_{\theta}=e^{(i\theta/2)S_{3}}(S_{1},S_{2})e^{(-i\theta/2)S_{3}}$.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

Figure 2: (color online) (a), (c), (e) The band structure of TBG lattices with
hopping parameters $W_{AA^{\prime}}=0$ at twist angle
$\theta=0.5^{\degree},1.08^{\degree},5^{\degree}$ respectively. The flat band
only occurs at the magic angle $\theta=1.08^{\degree}$. (b), (d), (f) The band
structure of TBD lattice as we set $W_{AA^{\prime}}=0$, totally flat band with
high level degeneracy appears at zero energy at all angles. (g) and (h) are
the band structures of TBG and TBD at twist angle $\theta=5^{\degree}$,
$W_{AA^{\prime}}/W_{AB^{\prime}}=0.5$ respectively. In TBD, with the change of
$W_{AA^{\prime}}$ from 0 to 1, the flat bands gradually spread. It is
important to notice that we already set $W_{AC^{\prime}}=0$, which will break
particl-hole symmetry at align bilayer dice lattice case. The momentum labels
$1\rightarrow 2\rightarrow 3\rightarrow 4\rightarrow 1$ is defined in Fig.
3(b).

For the interlayer term $H_{\bot}$, we consider nearest neighbor hopping from
layer 1 in $\alpha(\alpha=A,B,C)$ sublattice to closest sublattice
$\beta(\beta=A^{\prime},B^{\prime},C^{\prime})$ in layer 2 (see Fig. 1).

$\displaystyle
H_{\bot}=\sum_{\mathbf{r}_{1},\mathbf{r}_{2},\alpha,\beta}c_{1,\alpha}^{\dagger}(\mathbf{r}_{1})t^{\alpha,\beta}_{1,2}(\mathbf{r}_{1},\mathbf{r}_{2})c_{2,\beta}(\mathbf{r}_{2})+h.c,$
(A3)

where
$t_{1,2}^{\alpha,\beta}=\bra{1,\mathbf{r}_{1},\alpha}H_{\bot}\ket{2,\mathbf{r}_{2},\beta}$.
The band structure of aligned bilayer graphene is changed by the type of
stacking, which is known as $AA$-stacked and $AB$-stacked (Bernal). For TBD,
we set the aligned configuration coordinate as
$\tau_{1,A}=\tau_{2,B^{\prime}}=(0,0),\tau_{1,B}=\tau_{2,C^{\prime}}=(0,d),\tau_{1,C}=\tau_{2,A^{\prime}}=(0,2d)$(see
Fig. 1). To compare the results with TBG, we choose $1.42\text{\AA}$, the
lattice constant of graphene.

The transition formula becomes

$\displaystyle
T^{\alpha,\beta}_{1,2}(\mathbf{k_{1}},\mathbf{k_{2}})=\frac{1}{A_{u.c}}\sum_{\mathbf{G}_{1},\mathbf{G}_{2}}e^{-i\mathbf{G}_{1}\cdot\tau_{1,\alpha}}t(K_{1}+\mathbf{q}_{1}+\mathbf{G}_{1})\cdot$
$\displaystyle
e^{-i\mathbf{G}_{2}\cdot\tau_{1,\beta}}\delta_{K_{1}+\mathbf{q}_{1}+\mathbf{G}_{1},K_{2}+\mathbf{q}_{2}+\mathbf{G}_{2}}$
(A4)

where $A_{u.c}$ is unit cell area, $K_{1,2}$ represent Dirac point for each
layer which satisfy $K_{2}=R_{\theta}\cdot K_{1}$ where $R_{\theta}$ is the
rotation matrix. Here we assume that interlayer hopping obeys
$t^{\alpha,\beta}_{1,2}(\mathbf{r}_{1},\mathbf{r}_{2})=t^{\alpha,\beta}_{1,2}(\mathbf{r}_{1}+\tau_{1,\alpha}-\mathbf{r}_{2}-\tau_{2,\beta})$,
which decays rapidly if $k$ in reciprocal space exceeds the Dirac
pointBistritzer-MacDonald2011 . Considering this property we can only choose
three vectors ${\bf 0}$, and
$\mathbf{G}_{l}=\mathbf{g}_{l,1},\mathbf{g}_{l,2},\mathbf{g}_{l,3}$ with
$\mathbf{g}_{l,1}=0,\mathbf{g}_{l,2}=\mathbf{b}_{l,2},\mathbf{g}_{l,3}=-\mathbf{b}_{l,1}$
being reciprocal lattice vectors and $l=1,2$ being layer index. Then we have
three hopping matrices reads

$\displaystyle T_{1}=W\begin{pmatrix}1&1&1\\\ 1&1&1\\\
1&1&1\end{pmatrix},T_{2}=We^{-i\mathbf{g}_{1,2}\cdot\mathbf{\tau}_{0}}\begin{pmatrix}e^{i\phi}&1&e^{-i\phi}\\\
e^{-i\phi}&e^{i\phi}&1\\\ 1&e^{-i\phi}&e^{i\phi}\end{pmatrix},$ $\displaystyle
T_{3}=We^{-i\mathbf{g}_{1,3}\cdot\mathbf{\tau}_{0}}\begin{pmatrix}e^{-i\phi}&1&e^{i\phi}\\\
e^{i\phi}&e^{-i\phi}&1\\\ 1&e^{i\phi}&e^{-i\phi}\end{pmatrix},$ (A5)

where $W=\frac{t(|K|)}{A_{u.c}}=110meV$ with $|K|=4\pi/(3\sqrt{3}d)$ as in
graphene Bistritzer-MacDonald2011 , and $\phi=2\pi/3$.

(a)

(b)

Figure 3: (color online) (a): The Moiré period structure of dice lattice which
unit cell described by its Moiré lattice vectors reads
$\bm{a_{1}^{m}}=L_{s}(\frac{\sqrt{3}}{2},-\frac{1}{2}),\bm{a_{2}^{m}}=L_{s}(\frac{\sqrt{3}}{2},\frac{1}{2})$,
where $L_{s}=\frac{a}{2\sin(\theta/2)}$ is size of Moiré unit cell. (b):
Reciprocal space of TBD with reciprocal vectors reads
$\bm{b^{m}_{1}}=\frac{4\pi}{\sqrt{3}L_{s}}(\frac{1}{2},-\frac{\sqrt{3}}{2}),\bm{b^{m}_{2}}=\frac{4\pi}{\sqrt{3}L_{s}}(\frac{1}{2},\frac{\sqrt{3}}{2})$,
and red dash line mark the first Brillouin zone with green dash line mark the
reciprocal path of TBD band structure. We use Bravais lattice vectors mapping
the reciprocal space of TBD like
$\bm{b}=l_{1}\bm{b^{m}_{1}}+l_{2}\bm{b^{m}_{2}}$, the first number represent
the layer index which also represent by red point(layer 1) and blue
point(layer2).

##  III  Comparing the band structures of TBG and TBD

After building effective model of TBD, we now can compare the band structures
of TBG and TBD. In our previous calculation, we assume the hopping term is
isotropic, namely, the hopping amplitude of different sublattice equal to each
other, which induces $W_{AA^{\prime}}=W_{AB^{\prime}}$. In TBG, the absolutely
flat band at magic angle $\theta\approx 1.08^{\degree}$ can be obtained in the
chiral limit, namely, by choosing $W_{AA^{\prime}}=0$ Tarnopolsky2019 . In the
chiral limit of TBD, highly degenerate flat bands at zero energy also exist.
In our model, the chiral limit corresponds to setting
$W_{AA^{\prime}}=W_{AC^{\prime}}=0$. $W_{AC^{\prime}}$ is a new ingredient in
TBD when comparing with TBG, and the effects of non-zero $W_{AC^{\prime}}$
will be discussed at the end of this section. We set parameter $W=110mev$, and
truncate Hamiltonian at order of $l_{1}=l_{2}=3$, meanwhile we choose three
moiré angle $\theta=0.5^{\degree},1.08^{\degree},5^{\degree}$, and final
results of the band structure are summarized in Fig.2.

Unlike TBG which has flat bands near magic angles, TBD has flat bands at all
angles, which is the manifestation of flat band of single layer dice lattice
model (see Fig. 2). These TBD flat bands are highly degenerate in the chiral
limit. When $W_{AA^{\prime}}$ and $W_{AC^{\prime}}$ are nonzero, the
degeneracy is lifted and behaves as Landau levels (see Fig. 5b).

When $W_{AA^{\prime}}$ is non-zero, the degeneracy of the central flat bands
is lifted. As $W_{AA^{\prime}}$ increases, we can observe from Fig. 2h that
the range of the central flat band gradually expands, where we have set
$W_{AC^{\prime}}=0$. Interestingly, in the TBD system, even if we set the
parameter $W_{AA^{\prime}}$ = 0, changing the value of $W_{AC^{\prime}}$ can
also lead to the lifting of degeneracy in
$\theta=0.5^{\degree},1.08^{\degree},5^{\degree}$. Therefore, to obtain the
band structure of the central flat band, both parameters, $W_{AA^{\prime}}$
and $W_{AC^{\prime}}$, must be set to zero. Additionally, when the two layers
are aligned, the parameter $W_{AC^{\prime}}$ affects particle-hole symmetry.
Meanwhile, in TBD under the condition both $W_{AA^{\prime}}\neq 0$ and
$W_{AC^{\prime}}\neq 0$, the particle-hole symmetry will also be broken.

After discussing the effects of parameters $W_{AA^{\prime}}$ and
$W_{AC^{\prime}}$ on the flat band structure of TBD system, we now turn to
their effects on the band structure in align bilayer case. In both perfect
aligned stacking and A-B stacking bilayer dice lattice, the parameter
$W_{AC^{\prime}}$ will significantly affect the particle-hole symmetry of the
band. First we have the Hamiltonian of A-B align bilayer dice lattice reads

$\displaystyle H=\begin{pmatrix}0&h(\bm{k})&0&0&t_{\bot}&0\\\
h^{*}(\bm{k})&0&f(\bm{k})&0&0&t_{\bot}\\\ 0&f^{*}(\bm{k})&0&0&0&0\\\
0&0&0&0&h(\bm{k})&0\\\ t_{\bot}&0&0&h^{*}(\bm{k})&0&f(\bm{k})\\\
0&t_{\bot}&0&0&f^{*}(\bm{k})&0\end{pmatrix},$ (A6)

where
$h(\bm{k})=-\tau(e^{i\bm{k}\cdot\tau_{B}}+e^{i\bm{k}\cdot(\tau_{B}-a_{1})}+e^{i\bm{k}\cdot(\tau_{B}-a_{2})})$,
and corresponding
$f(k)=-\tau(e^{i\bm{k}\cdot(\tau_{B}-\tau_{C})}+e^{i\bm{k}\cdot(\tau_{B}-a_{1}-\tau_{C})}+e^{i\bm{k}\cdot(\tau_{B}-a_{2}-\tau_{C})})$
which $\tau=\sqrt{2}v_{f}\hbar/(3d)$. The hopping amplitude of A-B sublattice
same with graphene and the basis of bilayer align dice lattices is
$\Psi^{\dagger}=[c^{\dagger}_{1,A},c^{\dagger}_{1,B},c^{\dagger}_{1,C},c^{\dagger}_{2,A^{\prime}},c^{\dagger}_{2,B^{\prime}},c^{\dagger}_{2,C^{\prime}}]$.
The position of atoms B and C in one unit cell are
$\tau_{B}=d(0,1),\tau_{C}=d(0,2)$, and $d$ is lattice constant. It is worth
noticing that we set hopping parameter $W_{AC^{\prime}}=0$ in order to
maintain particle-hole symmetry of system. The band structure of two different
ways of stacking shown in Fig.4

(a)

(b)

Figure 4: (color online) Here we set reciprocal path same with
graphenecatarina2019twisted , lattices vector and coordinate of atoms in unit
cell same with Fig.1. (a) The band structure of A-A stacking(perfect align),
here we set hopping amplitude
$t_{A,A^{\prime}}=t_{B,B^{\prime}}=t_{C,C^{\prime}}=0.33eV$. (b) The band
structure of A-B stacking(Bernal-stack), in this case hopping between $A$ and
$C^{\prime}$ sublattices is forbidden because of particle-hole symmetry,
$t_{A,B^{\prime}}=t_{B,C^{\prime}}=0.33eV$2016PhR…648….1R .

##  IV  $SU(2)$ symmetry and Landau level structure

We have assumed chiral limit in the previous calculations of the band
structure, while in real TBG system at small twist angle $\theta\leq
10^{\degree}$, the lattice period will be dependent on different atomic
stacking and cause atomic layer deformation PhysRevB.90.155451 , which
eventually drives the system away from the chiral limit. This effect can be
described by introducing interlayer couplings which can be represented by an
$SU(2)$ gauge potential Liu2019Pseudo ; San-Jose2012Non-Abelian ; ren2021 .
Similar to TBG, one can write the low energy Hamiltonian as

$\displaystyle H=-\hbar
v_{f}(\mathbf{k}-\frac{e}{\hbar}\mathbf{A}\tau_{2})\cdot S+3u_{0}\tau_{1},$
(A7)

where $3u_{0}$ represents the hopping parameter $W_{AA^{\prime}}$ which
transforms into a scalar potential Liu2019Pseudo , and $\tau_{1,2}$ are the
Pauli matrices that act on the layer index. The gauge potential reads
$\mathbf{A}=(2\pi u^{\prime}_{0})/(L_{s}ev_{f})(y,-x)$, $u^{\prime}_{0}$
denote hopping parameters $W_{AB^{\prime}},W_{BC^{\prime}}$ in our system and
$u_{0}$ denote $W_{AA^{\prime}}$. The effective magnetic field ${\bf
B}=\nabla\times{\bf A}\approx 120T$ for TBG at $\theta\approx 1.08^{\degree}$
Liu2019Pseudo . The effective magnetic field has opposite directions on each
layer, therefore time-reversal symmetry is preserved. In the following
numerical calculations, we set $u^{\prime}_{0}=0.01eV,u_{0}=0.1eV$. A
comparison between the effective continuous Hamiltonian and the Landau level
description is presented in Appendix A, which also supports the usage of the
Landau level method.

To discuss the Landau level structure, we transform the Hamiltonian (A7) into
basis of $\tau_{2}$

$\displaystyle{\scriptsize
H(\mathbf{k})=\begin{pmatrix}0&\pi_{x}+i\pi_{y}&0&3iu_{0}&0&0\\\
\pi_{x}-i\pi_{y}&0&\pi_{x}+i\pi_{y}&0&3iu_{0}&0\\\
0&\pi_{x}-i\pi_{y}&0&0&0&3iu_{0}\\\
-3iu_{0}&0&0&0&\pi^{\prime}_{x}+i\pi^{\prime}_{y}&0\\\
0&-3iu_{0}&0&\pi^{\prime}_{x}-i\pi^{\prime}_{y}&0&\pi^{\prime}_{x}+i\pi^{\prime}_{y}\\\
0&0&-3iu_{0}&0&\pi^{\prime}_{x}-i\pi^{\prime}_{y}&0\end{pmatrix}},$ (A8)

where $\pi_{x}=\hbar v_{f}k_{x}-ev_{f}A_{x}$, $\pi_{y}=-\hbar
v_{f}k_{y}+ev_{f}A_{y}$, $\pi^{\prime}_{x}=\hbar v_{f}k_{x}+ev_{f}A_{x}$, and
$\pi^{\prime}_{y}=-\hbar v_{f}k_{y}-ev_{f}A_{y}$. Since $u_{0}$ is small, we
can treat it as a perturbation. Then for the unperturbed Hamiltonian $H_{0}$,
we define the Landau level creation operators of the two layers,
$b^{\dagger}=\sqrt{\frac{L_{s}}{8\pi u^{\prime}_{0}\hbar
v_{f}}}(\pi_{x}-i\pi_{y})$, and $a^{\dagger}=\sqrt{\frac{L_{s}}{8\pi
u^{\prime}_{0}\hbar v_{f}}}(\pi^{\prime}_{x}+i\pi^{\prime}_{y})$. And the
nonzero commutators are $[b,b^{\dagger}]=[a,a^{\dagger}]=1$. Denoting the
unperturbed eigenfunction as $\ket{\Psi}$,

$\displaystyle H_{0}(\mathbf{k})\ket{\Psi}=\hbar
w\begin{pmatrix}0&b&0&0&0&0\\\ b^{\dagger}&0&b&0&0&0\\\
0&b^{\dagger}&0&0&0&0\\\ 0&0&0&0&a^{\dagger}&0\\\ 0&0&0&a&0&a^{\dagger}\\\
0&0&0&0&a&0\end{pmatrix}\begin{pmatrix}A^{1}_{n}\ket{n-1}\\\
A^{2}_{n}\ket{n}\\\ A^{3}_{n}\ket{n+1}\\\ B^{1}_{m}\ket{m+1}\\\
B^{2}_{m}\ket{m}\\\ B^{3}_{m}\ket{m-1}\end{pmatrix},$ (A9)

where $\hbar w=\sqrt{\frac{8\pi u^{\prime}_{0}\hbar v_{f}}{L_{s}}}$,
$\ket{n},\ket{m}$ denote Landau levels with
$b\ket{n}=\sqrt{n}\ket{n-1},b^{\dagger}\ket{n}=\sqrt{n+1}\ket{n+1},a^{\dagger}\ket{m}=\sqrt{m+1}\ket{m+1},a\ket{m}=\sqrt{m}\ket{m-1}$.

The spectra of $H_{0}(\mathbf{k})$ For $n,m>0$ reads

$\displaystyle E_{\pm}^{n,m}=\pm\hbar w\sqrt{2(n,m)+1},\quad E_{0}^{n,m}=0.$
(A10)

Due to the quasi spin-1 nature of the sublattice structure of Dice model, For
$n,m=0$,

$\displaystyle E^{0}_{\pm}=\pm\hbar w,$ (A11)

Which correspond to the 1st Landau levels in the upper and lower layers. And
the $E=0$ for the 0th Landau level with $n,m=-1$. As mentioned before, $n(m)$
labels the states in Layer 1(2), therefore we use $1(2)$ to denote the
corresponding wave functions.

The wavefunctions for $m,n=1,2,3...$ are

$\displaystyle\ket{\Psi}^{1}_{\pm}=\begin{pmatrix}\sqrt{\frac{n}{4n+2}}\ket{n-1}\\\
\pm\sqrt{\frac{2n+1}{4n+2}}\ket{n}\\\ \sqrt{\frac{n+1}{4n+2}}\ket{n+1}\\\ 0\\\
0\\\ 0\end{pmatrix},\ket{\Psi}^{2}_{\pm}=\begin{pmatrix}0\\\ 0\\\ 0\\\
\sqrt{\frac{m+1}{4m+2}}\ket{m+1}\\\ \pm\sqrt{\frac{2m+1}{4m+2}}\ket{m}\\\
\sqrt{\frac{m}{4m+2}}\ket{m-1}\\\ \end{pmatrix},$ (A12)
$\displaystyle\ket{\Psi}^{1}_{0}=\begin{pmatrix}-\sqrt{\frac{n+1}{2n+1}}\ket{n-1}\\\
0\ket{n}\\\ \sqrt{\frac{n}{2n+1}}\ket{n+1}\\\ 0\\\ 0\\\
0\end{pmatrix},\ket{\Psi}^{2}_{0}=\begin{pmatrix}0\\\ 0\\\ 0\\\
-\sqrt{\frac{m}{2m+1}}\ket{m+1}\\\ 0\ket{m}\\\
\sqrt{\frac{m+1}{2m+1}}\ket{m-1}\\\ \end{pmatrix}.$ (A13)

For $n,m=0$, we have

$\displaystyle\ket{\Psi}^{1}_{\pm}=\sqrt{\frac{1}{2}}\begin{pmatrix}0\\\
\pm\ket{0}\\\ \ket{1}\\\ 0\\\ 0\\\
0\end{pmatrix},\quad\ket{\Psi}^{2}_{\pm}=\sqrt{\frac{1}{2}}\begin{pmatrix}0\\\
0\\\ 0\\\ \ket{1}\\\ \pm\ket{0}\\\ 0\end{pmatrix},$ (A14)

and for $n,m=-1$

$\displaystyle\ket{\Psi}^{1}_{0}=\begin{pmatrix}0\\\ 0\\\ \ket{0}\\\ 0\\\ 0\\\
0\end{pmatrix},\quad\ket{\Psi}^{2}_{0}=\begin{pmatrix}0\\\ 0\\\ 0\\\
\ket{0}\\\ 0\\\ 0\end{pmatrix}.$ (A15)

The band structure is illustrated in Fig.5, which is similar to the band
structure in Fig. 2. And the lifting of the degeneracy of flat bands is caused
by the perturbation $u_{0}$, the interlayer hopping. After perturbation, a
double degeneracy remains which is related to the $S_{3}$ symmetry (the
pseudo-spin in the $z$ direction) Liu2019Pseudo .

##  V  Optical conductance

In reality, materials that rotates the polarization plane of linearly
polarized light have lots of applications in various devices, which is usually
achieved by the magneto-optical effects, quantum Hall effect, and the Kerr and
Faraday rotations by invoking external magnetic fields oc1 ; oc2 ; oc3 ; oc4 ;
oc5 ; oc6 ; oc7 , which also limits the applications in small-scale devices.
Therefore, searching for materials with intrinsic properties that rotates
light is urgent for recent applications. The TBG is such a candidate for
advanced optical applications due to the tunable twist angles occ1 ; occ2 .
Here we study the optical conductance of TBD, which could be another candidate
for such applications.

We utilize Kubo formula in Landau level basis to calculate the optical
conductance of TBD system PhysRevB.94.125435 ,

$\displaystyle\sigma_{\alpha,\beta}=\frac{ig}{2\pi\hbar
l^{2}_{B}}\sum_{LL_{s}}\frac{f-f^{\prime}}{\varepsilon^{\prime}-\varepsilon}\frac{\braket{\Psi}{j_{\alpha}}{\Psi^{\prime}}\braket{\Psi^{\prime}}{j_{\beta}}{\Psi}}{w-(\varepsilon^{\prime}-\varepsilon)+i\Gamma},$
(A16)

where $l_{B}=\sqrt{(L_{s}hcv_{f})/4\pi u^{\prime}_{0}}$ is magnetic length,
without considering spin degree of freedom we can simply set $g=2$, $f$ is
Fermi distribution meanwhile $w$ is photon energy. Here we shall emphasize
that, although we use the Landau level description, no magnetic field is
applied. The effective magnetic field is induced by the hopping of
$W_{AB^{\prime}}$ and $W_{BC^{\prime}}$, and the effective magnetic field has
opposite directions on top and bottom layers of TBD. Therefore time-reversal
symmetry is preserved.

We take into account the contribution of all Landau levels, however the
remaining double degeneracy of the band structure will cause divergence of
conductance. For doubly-degenerate band, we fix the divergent part of Kubo
formula (A16) into

$\displaystyle\lim_{\epsilon^{\prime}-\epsilon\to
0}\frac{f(\epsilon)-f^{\prime}(\epsilon^{\prime})}{\epsilon^{\prime}-\epsilon}=\frac{e^{(\epsilon_{0}-\mu)/kT}}{(e^{(\epsilon_{0}-\mu)/kT}+1)^{2}}\cdot\frac{1}{kT}.$
(A17)

From Hamiltonian (A8), the current is obtained by
$j_{\alpha,\beta}=\frac{\delta H_{0}(\mathbf{k})}{\delta A_{\alpha,\beta}}$,

$\displaystyle j_{x}=ev_{f}\begin{pmatrix}0&-1&0&0&0&0\\\ -1&0&-1&0&0&0\\\
0&-1&0&0&0&0\\\ 0&0&0&0&1&0\\\ \ 0&0&0&1&0&1\\\ 0&0&0&0&1&0\end{pmatrix},$
(A18) $\displaystyle j_{y}=ev_{f}\begin{pmatrix}0&i&0&0&0&0\\\ -i&0&i&0&0&0\\\
0&-i&0&0&0&0\\\ 0&0&0&0&-i&0\\\ \ 0&0&0&i&0&-i\\\ 0&0&0&0&i&0\end{pmatrix}.$
(A19)

(a)

(b)

(c)

(d)

Figure 5: (color online) (a) The Landau level of TBG lattice in the vicinity
of $K_{1}$ and $K_{2}$ with chemical potential $\mu=0.05eV$ (the red dotted
line). (b) The Landau level of TBD with same chemical potential of TBD. (c)
Optical conductance of TBG with same truncation, $\pm$ represent the two flat
band near zero energy, for the band bigger(smaller) we denote as
$L=0^{\pm},\pm 1,\pm 2,...$. For example, $T_{\pm,0^{+}}$ represents the
contribution from two flat bands to first upper branch Landau level denoted by
$L=0^{+}$. (d) Optical conductance of TBD with truncate $n=m=6$, we denote the
Landau levels(without flat band near zero energy) with integer $L=0^{\pm},\pm
1,\pm 2,...$, $0^{+}$ denoting upper branch and $0^{-}$ lower branch. For flat
band we use their energy value to denote, for example $T_{-0.012,0^{+}}$
respects contribution from flat band with energy $E=-0.012eV$ to first upper
branch Landau level denoted by $L=0^{+}$.

The numerical result for $\mu=0.05eV$ and $T=300K$ of the conductance of TBD
lattice is shown in Fig. 5(c). The double peak structure is caused by the
splitting of double degeneracy, and the width of splitting depends on
$u^{\prime}_{0}$. Due to the choice of chemical potential, Meanwhile, because
of the high DOS in the near zero band of Landau level, all significant
absorption peaks are originated from transitions between the near zero band
and the positive energy band, or between the bottom valence band and the near
zero band. By comparing with the result of TBG in Fig. 5(d), the fine
structure of TBD provides a smoking-gun experiment prediction for showing the
existence of these near zero energy bands.

Although the Landau level description for TBG fails when the twist angle is
away from the magic ones, the Landau level structure remains in TBD (see Figs.
2g and 2h). Therefore in the optical conductance at angles other than the
magic ones, the peak splitting remains in TBD while the peaks raising from the
transitions form the flat bands disappear in TBG, which is also a key
experimental difference between TBG and TBD.

The transitions between these zero energy levels are forbidden in TBD, which
meanings that there are no peaks at low frequency near $\omega\sim 0$, which
is a key difference to the three-dimensional Kane fermion Luo2019 . The reason
of this phenomenon attributed to the structure of current and wavefunction in
2D, that is,
$\braket{\Psi}{j_{\alpha}}{\Psi^{\prime}}\braket{\Psi^{\prime}}{j_{\beta}}{\Psi}=0$
between those near zero band, and their contribution to the optical
conductance is zero, while in three-dimensional, the $k_{3}S_{3}$ part will
have a nontrivial contribution, and the transition between different Landau
levels has a non zero $k_{3}$, which causes the peaks near zero frequency
Luo2019 .

##  VI  Conclusions

In this paper, we constructed a lattice model for TBD. It has flat bands at
all twisted angles besides the magic ones, which could be conformed by the
peak splitting fine structure of the optical conductance near the magic
angles. And away from the magic angles, the peak splitting remains in TBD
while these peaks disappear in TBG due to the nonexistence of the flat bands.
The flat bands in TBD are composed of zero Chern number ones by destructive
interference of the dice lattice as well as the topological nontrivial ones by
the Moiré structure at the magic angles. In this model we have neglected the
spin degrees of freedom of electrons. If the spin orbital coupling interaction
is added, the bands with zero Chern number may become non-trivial wang2011 .
The TBD is possible to be realized in the transition-metal oxide
SrTiO3/SrIrO3/SrTiO3 trilayer heterostructure by growing and twisting along
the (111) direction. As a semiconductor, due to the large DOS of the flat
bands of TBD, it may have potential applications on temperature-sensitive,
photosensitive manipulations. With interactions, TBD may also be a good
candidate for fractional Chern insulator.

###### Acknowledgements.

This work is supported by the National Natural Science Foundation of China
with No. 12174067.

##  Appendix A  Density of states of two different system

In our previous discussions, we follow Ref. Bistritzer-MacDonald2011  to
construct a tight-binding model for TBG and TBD lattice in reciprocal space.
To discuss the influence of parameters $W_{AC^{\prime}}$ and $W_{AA^{\prime}}$
on the spectrum of these models, we also consider the chiral
limitLiu2019Pseudo ; San-Jose2012Non-Abelian  for a comparison. In order to
compute the optical conductance of both models, we use the pseudo-Landau level
description (A7) of the effective theory. In order to verify applicability of
the Landau level scenario, we construct a reciprocal space tight-binding model
in Landau basis, and compare the density of states (DOS) between this model
and the one generated from Ref. Bistritzer-MacDonald2011 . We find that they
have similar behaviors (see Fig. 6), especially, they have a high DOS at the
zero energy near the magic angle of $1.08^{\degree}$, which supports the usage
of the effective description of pseudo-Landau levels.

The low energy Dirac equation reads

$\displaystyle H=-\hbar
v_{f}(\mathbf{k}-\frac{e}{\hbar}\mathbf{A}\tau_{2})\cdot\bm{\sigma}+3u_{0}\tau_{1},$
(A1)

where hopping matrices reads

$\displaystyle
T=\begin{pmatrix}\psi_{1,A}^{\dagger}&\psi_{1,B}^{\dagger}\end{pmatrix}\begin{pmatrix}3iu_{0}&0\\\
0&3iu_{0}\\\ \end{pmatrix}\begin{pmatrix}\psi_{2,A^{\prime}}\\\
\psi_{2,B^{\prime}}\end{pmatrix},$ (A2)

(a)

(b)

Figure 6: (color online) (a) DOS of TBG mode which calculated using Landau
level model with $\theta=1.08^{\degree}$ and $u_{0}/u_{0}^{\prime}=0.2$. (b)
DOS of TBG is calculated by using continuum model with
$W_{AA^{\prime}}=W_{AB^{\prime}}$ and $\theta=1.08^{\degree}$. We calculate
this result in order to compare with Ref. Bistritzer-MacDonald2011 . For more
details see Appendix A.

after transforming into the Landau basis, the Hamiltonian of intralayer reads

$\displaystyle H_{0}(\mathbf{k})=\hbar w\begin{pmatrix}0&\sqrt{n}&0&0\\\
\sqrt{n}&0&0&0\\\ 0&0&0&\sqrt{m}\\\ 0&0&\sqrt{m}&0\end{pmatrix}.$ (A3)

The Hamiltonian of low energy continuum TBG in the first honeycomb(see Fig.
3(b)) reads

$\displaystyle H(q)=\begin{pmatrix}H_{1}^{K_{1}}(q)&T_{1}&T_{2}&T_{3}\\\
T_{1}^{\dagger}&H_{2}^{K_{2}}(q+\mathbf{q_{b}})&0&0\\\
T_{2}^{\dagger}&0&H_{2}^{K_{2}}(q+\mathbf{q_{tr}})&0\\\
T_{3}^{\dagger}&0&0&H_{2}^{K_{2}}(q+\mathbf{q_{tl}})\end{pmatrix},$ (A4)

where $H_{1}^{K_{1}}$ denote Hamiltonian of layer 1 in the vicinity of Dirac
point $K_{1}$, and $1,2$ is layer index. The three hopping direction can be
illustrated as

$\displaystyle\mathbf{q_{b}}=\Delta K,$ (A5)
$\displaystyle\mathbf{q_{tr}}=\Delta K+\mathbf{g}_{1,2}-\mathbf{g}_{2,2},$
(A6) $\displaystyle\mathbf{q_{tl}}=\Delta
K+\mathbf{g}_{1,3}-\mathbf{g}_{2,3}.$ (A7)

Thus, we can use reciprocal vectors $b^{m}_{1}$ and $b^{m}_{2}$ mapping
numerous reciprocal points(see Fig. 3(b)) to construct a continuum model. We
numerically solve the continuum model contains near 250 reciprocal points,
which allow us to calculate DOS of the system in different method. Result
display in Fig. 6, two methods with same rotation angle
$\theta=1.08^{\degree}$ and truncation $l_{1},l_{2}=5$. From Fig. 6 we can see
that the DOS result of two methods show similar structure.

## References

  * (1) A. Mielke, Journal of Physics A: Mathematical and General 24, L73 (1991).
  * (2) H. Tasaki, Physical Review Letters 69, 1608 (1992).
  * (3) S. Zhang, H.-H. Hung, and C. Wu, Physical Review A 82, 053618 (2010).
  * (4) Jun-Won Rhim, Kyoo Kim, and Bohm-Jung Yang, Nature 584, 59-63 (2020).
  * (5) J.-X. Yin, S. S. Zhang, G. Chang, Q. Wang, S. S. Tsirkin, Z. Guguchia, B. Lian, H. Zhou, K. Jiang, I. Belopolski, N. Shumiya, D. Multer, M. Litskevich, T. A. Cochran, H. Lin, Z. Wang, T. Neupert, S. Jia, H. Lei, and M. Z. Hasan, Nature Physics 15, 443-448 (2019).
  * (6) E. Tang, J.-W. Mei, and X.-G. Wen, Physical Review Letters 106, 236802 (2011).
  * (7) E. J. Bergholtz, and Z. Liu, International Journal of Modern Physics B 27, 1330017 (2013).
  * (8) C. Wu, D. Bergman, L. Balents, and S. Das Sarma, Physical Review Letters 99, 070401 (2007).
  * (9) C. Wu and S. Das Sarma, Physical Review B 77, 235107 (2008).
  * (10) R. Micnas, J. Ranninger, and S. Robaszkiewicz, Reviews of Modern Physics 62, 113 (1990).
  * (11) S. Miyahara, S. Kusuta, and N. Furukawa, Physica C: Superconductivity 460, 1145 (2007).
  * (12) H. Aoki, Journal of Superconductivity and Novel Magnetism 33, 2341 (2020).
  * (13) I. Syôzi, Progress of Theoretical Physics 6, 306 (1951).
  * (14) E. H. Lieb, Phys. Rev. Lett. 62, 1201 (1989).
  * (15) B. Sutherland, Phys. Rev. B 34, 5208 (1986).
  * (16) Z. Liu, Z.-F. Wang, J.-W. Mei, Y.-S. Wu, and F. Liu, Phys. Rev. Lett. 110, 106804 (2013).
  * (17) Z. Liu, F. Liu, Y.-S. Wu, Chin. Phys. B 23, 077308, (2014).
  * (18) Y.-G. Chen, J.-T. Huang, K. Jiang, J.-P. Hu, arXiv:2212.13526.
  * (19) F. Wang, and Y. Ran, Phys. Rev. B 84, 241103(R) (2011).
  * (20) M. R. Slot, T. S. Gardenier, P. H. Jacobse, G. C.P. van Miert, S. N. Kempkes, S. J.M. Zevenhuizen, C. M. Smith, D. Vanmaekelbergh, I. Swart, Nature Physics 13, 672-676 (2017).
  * (21) R. Shen, L. B. Shao, B. Wang, and D. Y. Xing. Phys. Rev. B 81, 041410(R) (2010).
  * (22) D. Guzmán-Silva, C. Mejía-Cortés, M. A. Bandres, M. C. Rechtsman, S. Weimann, S. Nolte, M. Segev, A. Szameit, and R. A. Vicencio. New J. Phys. 16, 063061 (2014).
  * (23) S. Mukherjee, , A. Spracklen, D. Choudhury, N. Goldman, P. Ohberg, E. Andersson, and R. R. Thomson. Phys. Rev. Lett. 114, 245504 (2015).
  * (24) R. A. Vicencio, C. Cantillano, L. Morales-Inostroza, B. Real, C. Mejía-Cortés, S. Weimann, A. Szameit, and M. I. Molina. Phys. Rev. Lett. 114, 245503 (2015).
  * (25) S. Taie, H. Ozawa, T. Ichinose, T. Nishio, S. Nakajima, and Y. Takahashi. Sci. Adv. 1, 1500854 (2015).
  * (26) S. Xia, Y. Hu, D. Song, Y. Zong, L. Tang, and Z. Chen. Opt. Lett. 41, 1435 (2016).
  * (27) X. Luo, F.-Y. Li, Y. Yu, New J. Phys. 20, 083036 (2018).
  * (28) M. Orlita, D. M. Basko, M. S. Zholudev, F. Teppe, W. Knap, V. I. Gavrilenko, N. N. Mikhailov, S. A. Dvoretskii, P. Neugebauer, C. Faugeras, A. -L. Barra, G. Martinez, and M. Potemski, Nature Physics 10, 233 (2014).
  * (29) A. Akrap, M. Hakl, S. Tchoumakov, I. Crassee, J. Kuba, M. O. Goerbig, C. C. Homes, O. Caha, J. Novak, F. Teppe, W. Desrat, S. Koohpayeh, L. Wu, N. P. Armitage, A. Nateprov, E. Arushanov, Q. D. Gibson, R. J. Cava, D. van der Marel, B. A. Piot, C. Faugeras, G. Martinez, M. Potemski, and M. Orlita, Phys. Rev. Lett. 117, 136401 (2016).
  * (30) X. Luo, Y.-G. Chen, and Y. Yu, New J. Phys. 21, 083010 (2019).
  * (31) R. Bistritzer and A. H. MacDonald. Proc. Natl. Acad. Sci. U.S.A. 108, 12233–12237 (2011).
  * (32) E. Y. Andrei, and A. H. MacDonald, Nature Materials 19, pages1265-1275 (2020).
  * (33) Y. Cao, V. Fatemi, S. Fang, K. Watanabe, T. Taniguchi, E. Kaxiras, and P. Jarillo-Herrero. Nature 556, 43-50 (2018).
  * (34) B. Lian, Z. Wang, and B. A. Bernevig, Phys. Rev. Lett. 122, 257002 (2019).
  * (35) K. P. Nuckolls, M. Oh, D. Wong, B. Lian, K. Watanabe, T. Taniguchi, B. A. Bernevig, A. Yazdani, Nature 588, 610-615 (2020)
  * (36) J. Liu, J. Liu, and X. Dai. Phys. Rev. B 99, 155415 (2019).
  * (37) G. Tarnopolsky, A.J. Kruchkov, and A. Vishwanath. Phys. Rev. Lett. 122, 106405 (2019).
  * (38) P. J. Ledwith, G. Tarnopolsky, E. Khalaf, A. Vishwanath, Phys. Rev. Research 2, 023237 (2020).
  * (39) Y. Xie, A. T. Pierce, J. Min Park, D. E. Parker, E. Khalaf, P. Ledwith, Y. Cao, S. H. Lee, S. Chen, P. R. Forrester, K. Watanabe, T. Taniguchi, A. Vishwanath, P. Jarillo-Herrero, and A. Yacoby, Nature volume 600, pages439-443 (2021).
  * (40) J. Cai, E. Anderson, C. Wang, X. Zhang, X. Liu, W. Holtzmann, Y. Zhang, F. Fan, T. Taniguchi, K. Watanabe, Y. Ran, T. Cao, L. Fu, D. Xiao, W. Yao, and X. Xu. Nature (2023).
  * (41) H. Park, J. Cai, E. Anderson, Y. Zhang, J. Zhu, X. Liu, C. Wang, W. Holtzmann, C. Hu, Z. Liu, T. Taniguchi, K. Watanabe, J. Chu, T. Cao, L. Fu, W. Yao, C. Chang, D. Cobden, D. Xiao, and X. Xu. Nature (2023).
  * (42) Y. Zeng, Z. Xia, K. Kang, J. Zhu, P. Knüppel, C. Vaswani, K. Watanabe, T. Taniguchi, K. F. Mak, and J. Shan. arXiv preprint arXiv:2305.00973.
  * (43) G. Catarina, B. Amorim, E. V. Castro, J. M. V. P. Lopes, and N. Peres, in Handbook of Graphene Set (John Wiley & Sons, Ltd, 2019), Chap. 6, pp. 177-231.
  * (44) K. Uchida, S. Furuya, J.-I. Iwata, and A. Oshiyama, Phys. Rev. B 90, 155451 (2014).
  * (45) A. V. Rozhkov, A. O. Sboychakov, A. L. Rakhmanov, and Franco Nori, Phys. Rep. 648, 1-104 (2016).
  * (46) P. San-Jose, J. Gonzalez, and F. Guinea, Phys. Rev. Lett. 108, 216802 (2012).
  * (47) Y. Ren, Q. Gao, A.H. MacDonald, and Q. Niu Phys. Rev. Lett. 126, 016404 (2021).
  * (48) W.-K. Tse and A. H. MacDonald, Phys. Rev. B 84, 205327 (2011).
  * (49) I. Crassee, J. Levallois, A. L. Walter, M. Ostler, A. Bost- wick, E. Rotenberg, T. Seyller, D. van der Marel, and A. B. Kuzmenko, Nat. Phys. 7, 48 (2011).
  * (50) R. Nandkishore and L. Levitov, Phys. Rev.Lett. 107, 097402 (2011).
  * (51) R. Shimano, G. Yumoto, J. Y. Yoo, R. Matsunaga, S. Tan- able, H. Hibino, T. Morimoto, and H. Aoki, Nat. Commun. 4, 1841 (2013).
  * (52) E. L. Izake, Journal of Pharmaceutical Sciences 96, 1659 (2007).
  * (53) W. Liu, J. Gan, D. Schlenk, and W. A. Jury, Proceedings of the National Academy of Sciences 102, 701 (2005).
  * (54) M. L. Solomon, J. Hu, M. Lawrence, A. García-Etxarri, and J. A. Dionne, ACS Photonics 6, 43 (2019).
  * (55) C.-J. Kim, A. Sáchez-Castillo, Z. Ziegler, Y. Ogawa, C. Noguez, and J. Park, Nat. Nanotechnol. 11, 520 (2016).
  * (56) S. T. Ho and V. N. Do, arXiv:2303.07330.
  * (57) E. Illes and E. J. Nicol, Phys. Rev. B 94, 125435 (2016).