# Analytic modelling of Quantum Capacitance and Carrier Concentration for
$\beta_{12}$ \- Borophene FET based Gas Sensor

Nimisha<EMAIL_ADDRESS>Reeta<EMAIL_ADDRESS>Arindam<EMAIL_ADDRESS>and Saumen
<EMAIL_ADDRESS>Department of Physics, Dibrugarh
University, Dibrugarh 786 004, Assam, India,

###### Abstract

In this work, we investigate the physical and electronic properties of
$\beta_{12}$ \- borophene FET-based gas sensor using a theoretical quantum
capacitance model based on tight - binding approach. We study the impact of
adsorbed NH3, NO, NO2 and CO gas molecule on its density of states, carrier
concentration, quantum capacitance and I-V characteristics. We found a
remarkable variation in the energy band structure and the density of states
(DOS) of the $\beta_{12}$ \- borophene in the presence of the adsorbed gas
molecule. The appearance of non-identical Van-Hove singularities in the DOS in
the presence of adsorbed gas molecules strongly indicates the high sensitivity
of $\beta_{12}$ \- borophene. We found a significant increase in the carrier
concentration for NH3 gas while it decreases for all other gases. Moreover, a
drastic change in quantum capacitance and current-voltage relation is also
observed in adsorbed gases. The different properties of the given gas
molecules are compared with the pristine borophene and found to exhibit
distinct wrinkles in each case, thereby indicating the strong selectivity of
our proposed gas sensor. Though $\beta_{12}$ \- borophene is found to be
highly sensitive for all studied gases, the NO gas is found to be most
sensitive compared to the others.

###### pacs:

72.80.Vp, 84.37.+q, 73.63.-b

##  I Introduction

During the last two decades, two-dimensional (2D) material has garnered
enormous interest due to their unique physical properties, making them a
promising candidate for next-generation electronics and energy conversion
devices novoselove ; lemme ; liu ; xu ; quellmalz ; illarionov ; passi ;
fuechsle ; acharjee . As the first 2D material with significantly different
electronic balandin ; lee ; zhang ; sang , thermal sang ; pop , mechanical
papageorgiou  and optical properties papageorgiou , graphene-based devices
have been widely explored. However, the absence of band gap, low on/off ratio
and extremely high carrier mobility in graphene restrict its use as a
semiconducting system neto ; mccann ; ando . So, there is a need for new
alternative materials with similar but more advanced electronic properties
than graphene. Thus a class of new 2D materials, like phosphorene akhtar ,
silicene vogt , stanene zhu2 , antimonene ji , bismuthene drozdov , germanine
davila , molybdenum disulfide splendiani  etc., have been synthesized as
alternate candidates for graphene in recent times. Depending on their stacking
and composition, these 2D materials can possess distinctive physicochemical
properties, such as electrical conductivity, thermal conductivity, and band
structure around the Fermi level, making them suitable for various
applications akhtar ; vogt ; zhu2 ; ji ; drozdov ; davila ; splendiani .

Alongside the other promising candidates, Borophene, a single element 2D sheet
of boron, has also been synthesized on Ag (111) substrate in recent times
mannix ; zhang2 ; feng . Unlike graphene, borophene has a triangular and
hexagonal array of atoms exhibiting in-plane elasticity and can be more
flexible in some configurations. Apart from that, borophene has high mobility,
electrical and thermal conductivity and also possesses high mechanical
strength, with Young’s modulus much higher than its rivals zhang3 ; mozvashi ;
yuan . Thus borophene has excellent potential in advanced electrical
information, sensing and other optical, mechanical and thermal applications
wang . Moreover, previous works indicate that the borophene-based systems are
active to hazardous gases, thus making them suitable for gas sensing
applications huang ; hou ; shukla . Significant attention also has been given
to the optimization of the borophene/MoS2 heterostructure using density
functional theory (DFT) implemented in the VASP package shen . Although there
exist several boron clusters with fascinating properties, it was found that
the allotrope $\beta_{12}$ \- borophene is thermodynamically most stable as
compared to the other members of the borophene family peng ; feng . It is
observed that $\beta_{12}$ \- borophene can display semiconducting and semi-
metallic properties in the presence of an applied electric field and charged
impurity li2 . Furthermore, this boron cluster can also display anisotropic
Kubo conductivity nobahari  in the presence of an applied electric field,
making them unique from other members. Though several attempts have been made
to understand the adsorption properties in different boron-based clusters both
theoretical as well as experimentally huang ; hou ; shukla ; opoku , the study
of gas sensing in $\beta_{12}$ \- borophene is limited. So, in this work, we
have focused on $\beta_{12}$ phase of borophene to investigate the adsorption
of the gas molecules.

Recently, quantum capacitance has gained much attention in studying 2D
electron systems to reveal interesting many-body effects tsu ; xia ;
hassanzadazar . It also carries information regarding the ground state of the
system, revealing the effect of electron-electron interaction and quantum
correlation. The Quantum capacitance of a system carries essential information
regarding the density of states. Recently quantum capacitance model has been
used to measure the Van-Hove singularities and Luttinger parameters in a one-
dimensional system li4 ; yu3 . Moreover, the quantum capacitance measurements
have also shown the linear density of states of topological insulator Bi2Se3
inhofer  and monolayer graphene cheremisin ; pourasl ; dutta . Thus, it is
necessary to have knowledge about the quantum capacitance for understanding
the system properly. Though quantum capacitance has been widely studied
theoretically and experimentally in graphene, such studies have been limited
to borophene until now. An attempt has been made to explore the quantum
capacitance effect in graphene chemical sensors and $\delta_{6}$ \- borophene
kolavada ; john  to induce changes in the capacitance values upon the
adsorption of molecules on the surface. But quantum capacitance of
$\beta_{12}$ \- borophene has not been explored so far. Therefore, in this
work, we investigated the quantum capacitance of $\beta_{12}$ \- borophene in
the presence of adsorbed gas molecules.

The organisation of the paper is as follows: we present a theoretical
framework based on tight-binding (TB) Hamiltonian and also present a quantum
capacitance model of $\beta_{12}$ \- borophene-based FET based gas sensor in
Section II. In Section III, we study the density of states, carrier
concentration, quantum capacitance and I-V characteristics to understand the
gas sensing property of $\beta_{12}$ \- borophene-based FET. A summary of our
work is presented in Section IV.

##  II Theoretical framework

The schematic illustration of $\beta_{12}$ \- borophene field effect
transistor-based gas sensor is shown in Fig. 1(a). The honeycomb lattice
arrangement and the unit cell of $\beta_{12}$ \- borophene is shown in Fig.
1(b). The unit cell comprises of five boron atoms labelled a, b, c, d and e.
The $\beta_{12}$ \- borophene lattice can be generated using the lattice basis
vectors ($\vec{a},\vec{b}$) $=(\sqrt{3}a_{0}\hat{e}_{x},3a_{0}\hat{e}_{y})$
where $a_{0}=1.74\AA$ is the boron-boron atom distance xia . It is to be noted
that borophene does not allow the formation of in-plane $\sigma$ bonds due to
the lack of one electron in the boron atom compared to the carbon atom. So,
the C-atoms can be treated as a perfect donor in filling the in-plane
hexagonal $\sigma$-bonds. Borophene has the energy levels involving $s$,
$p_{x}$, $p_{y}$ and $p_{z}$ orbitals which are $sp^{2}$ \- hybridized.
However, only $p_{z}$ \- orbitals can contribute to the carrier dynamics. It
is due to the reason that the wave function at the Fermi energy has a
vanishing amplitude at site c resulting in phase cancellation at the six-fold
B-atoms. The band structure of $\beta_{12}$ \- borophene in the presence of
adsorbed gas molecules can be obtained using the tight binding approach
considering only the nearest neighbour approximation. To model the molecular
adsorption effect, we consider a $\beta_{12}$ \- borophene sheet consisting of
N sites. The gas adsorption in the $\beta_{12}$ \- borophene lattice and the
$n^{\text{th}}$ unit cell with its nearest neighbours is shown in Figs. 1(c)
and 1(d).

Figure 1: (a) Schematic illustration of $\beta_{12}$ \- borophene field effect
transistor (FET) based gas sensor. (b) Top view of the geometry structure of
$\beta_{12}$ \- borophene. The unit cell (shaded green region) is one-half of
the honeycomb lattice and consists of five boron atoms. The lattice parameters
of $\beta_{12}$ \- borophene are $a=|\vec{a}|=\sqrt{3}a_{0}$ and
$b=|\vec{b}|=3a_{0}$ where $a_{0}$ is the boron-boron atom distance. (c) Top
view of the adsorption of the gas molecule on $\beta_{12}$ \- borophene
geometry. (d) 3D arrangement of the adsorption gas molecule on $a$ site of
$\beta_{12}$ \- borophene unit cell.

The matrix equation for the $n^{\text{th}}$ unit cell according to tight
binding model can be written as dutta ; pourasl

$\sum_{m}\left[\mathcal{H}_{nm}\right]\\{\Psi_{m}\\}=E\\{\Psi_{n}\\}$ (1)

where $\\{\Psi_{n}\\}$ is a ($b\times 1$) column matrix corresponds to the
wavefunction in the unit cell $n$. It is to be noted that the Hamiltonian
matrix for pristine $\beta_{12}$ \- borophene with five atoms in its unit cell
will be a ($5\times 5$) matrix. To characterize the effect of gas adsorption,
we consider that each gas molecule is adsorbed only in each A site of the unit
cell of $\beta_{12}$ \- borophene as shown in Fig. 1(d). Thus, the Hamiltonian
of the $\beta_{12}$ \- borophene in the presence of a gas molecule will be a
matrix of ($6\times 6$) dimension.

We consider a plane waveform of the wave function, i.e.,
${\Psi_{n}}=\Psi_{0}e^{i\vec{k}.\vec{d}_{n}}$, with $\vec{k}$ and
$\vec{d}_{n}$ being the plane wave vector and distance between a source with
the $n^{\text{th}}$ unit-cell respectively. The band structure of the
$\beta_{12}$ \- borophene can be calculated by solving a matrix eigenvalue
equation of the form dutta ; pourasl

$[h(k)]=\sum_{nm}\left[\mathcal{H}_{nm}\right]e^{i\vec{k}.(\vec{d}_{n}-\vec{d}_{m})}$
(2)

where $\\{\Psi_{n}\\}$ is a ($b\times 1$) column matrix corresponding to the
wavefunction in the unit cell $n$. It is to be noted that the Hamiltonian
matrix for pristine $\beta_{12}$ \- borophene with five atoms in its unit cell
will be a ($5\times 5$) matrix. To characterize the effect of gas adsorption,
we consider that each gas molecule is adsorbed only in each A site of the unit
cell of $\beta_{12}$ \- borophene as shown in Fig. 1(d). Thus, the Hamiltonian
of the $\beta_{12}$ \- borophene in the presence of a gas molecule will be a
matrix of ($6\times 6$) dimension.

$[h(k)]_{nm}=\left(\begin{array}[]{cccccc}0&t_{\text{ab}}\delta_{k}&0&0&0&0\\\
0&0&0&0&0&0\\\ 0&t_{\text{cb}}\delta_{k}&0&t_{\text{cd}}\delta_{k}&0&0\\\
0&0&0&0&0&0\\\ 0&0&0&t_{\text{ed}}\delta_{k}&0&0\\\ 0&0&0&0&0&0\\\
\end{array}\right)+\left(\begin{array}[]{cccccc}0&0&0&0&0&0\\\ 0&0&0&0&0&0\\\
0&0&0&0&0&0\\\ 0&0&0&0&0&0\\\ t_{\text{ea}}&0&0&0&0&0\\\ 0&0&0&0&0&0\\\
\end{array}\right)+\left(\begin{array}[]{cccccc}0&0&0&0&0&0\\\
t_{\text{ba}}\delta_{k}^{\ast}&0&t_{\text{bc}}\delta_{k}^{\ast}&0&0&0\\\
0&0&0&0&0&0\\\
0&0&t_{\text{dc}}\delta_{k}^{\ast}&0&t_{\text{de}}\delta_{k}^{\ast}&0\\\
0&0&0&0&0&0\\\ 0&0&0&0&0&0\\\
\end{array}\right)+\left(\begin{array}[]{cccccc}0&0&0&0&t_{\text{ae}}&0\\\
0&0&0&0&0&0\\\ 0&0&0&0&0&0\\\ 0&0&0&0&0&0\\\ 0&0&0&0&0&0\\\ 0&0&0&0&0&0\\\
\end{array}\right)$ (3)

where, $m=1,2,3,4$ correspond to the four nearest unit cells of the
$n^{\text{th}}$ unit cell. Here, $t_{\text{ij}}$ are the hopping energy
parameter of $i^{\text{th}}$ and $j^{\text{th}}$ boron atoms. The hopping
parameters for the homogeneous model are $t_{\text{ij}}=-2$ eV huang , while
in the inversion non-symmetric (INS) model, the hopping parameters are:
$t_{\text{ab}}=t_{\text{de}}=-2.04$ eV, $t_{\text{ac}}=t_{\text{ce}}=-1.79$
eV, $t_{\text{bc}}=t_{\text{cd}}=-1.84$ eV, $t_{\text{bd}}=-1.91$ eV,
$t_{\text{ad}}=0$ eV and $t_{\text{ae}}=-2.12$ eV huang . Here, we define
$\delta_{k}\equiv\exp\left[-\frac{ia_{0}k_{x}}{2}\right]$. It is to be noted
that one may use the lattice parameters ($\vec{a},\vec{b}$) for the
calculation of distance $\vec{d}_{n}$. However, in this work, the distance
between the lattice points $(\vec{d}_{m}-\vec{d}_{n})$ are obtained by
converting the lattice parameters in terms of bond length $a_{0}$ for
simplicity of the calculations.

In a similar way, the interaction of the $n^{\text{th}}$ unit cell with the
adsorbed gas molecule can be written as

$[h(k)]_{nn}=\left(\begin{array}[]{cccccc}\varepsilon_{\text{a}}&t_{\text{ab}}&t_{\text{ac}}&0&0&t^{\prime}\\\
t_{\text{ba}}&\varepsilon_{\text{b}}&t_{\text{bc}}&t_{\text{bd}}&0&0\\\
t_{\text{ca}}&t_{\text{cb}}&\varepsilon_{\text{c}}&t_{\text{cd}}&t_{\text{ce}}&0\\\
0&t_{\text{db}}&t_{\text{dc}}&\varepsilon_{\text{d}}&t_{\text{de}}&0\\\
0&0&t_{\text{ec}}&t_{\text{ed}}&\varepsilon_{\text{e}}&0\\\
t^{\prime}&0&0&0&0&\varepsilon^{\prime}\\\ \end{array}\right)$ (4)

where, $\varepsilon_{i}$ are the corresponding on-site energy of the
$i^{\text{th}}$ atom. The on-site energies in homogeneous model are
$\varepsilon_{i}=0$ eV while in the INS model on-site energies are:
$\varepsilon_{\text{a}}=\varepsilon_{\text{d}}=0.196$ eV,
$\varepsilon_{\text{b}}=\varepsilon_{\text{e}}=-0.058$ eV and
$\varepsilon_{\text{c}}=-0.845$ eV huang . The parameters $t^{\prime}$ and
$\varepsilon^{\prime}$ characterize the hopping energy between the Boron atom
and the adsorbed molecule and the on-site energies, respectively.

Figure 2: Band structure of $\beta_{12}$ \- borophene using (a) homogenous and
(b) the INS model. The plot (c) represents the band structure of $\beta_{12}$
\- borophene in presence of NH3 using the INS model.

Thus, the total tight binding Hamiltonian of the $\beta_{12}$-borophene system
in presence of adsorbed gas in view of Eq. (3) and Eq. (4) can be expressed as

$\left[h(k)\right]=\left(\begin{array}[]{cccccc}\varepsilon_{\text{a}}&f^{\ast}_{k}t_{\text{ab}}&t_{\text{ac}}&0&t_{\text{ae}}&t^{\prime}\\\
f_{k}t_{\text{ab}}&\varepsilon_{\text{b}}&f_{k}t_{\text{bc}}&t_{\text{bd}}&0&0\\\
t_{\text{ac}}&f^{\ast}_{k}t_{\text{bc}}&\varepsilon_{\text{c}}&f^{\ast}_{k}t_{\text{cd}}&t_{\text{ce}}&0\\\
0&t_{\text{bd}}&f_{k}t_{\text{cd}}&\varepsilon_{\text{d}}&t_{\text{de}}f_{k}&0\\\
t_{\text{ae}}&0&t_{\text{ce}}&f^{\ast}_{k}t_{\text{de}}&\varepsilon_{\text{e}}&0\\\
t^{\prime}&0&0&0&0&\varepsilon^{\prime}\\\ \end{array}\right)$ (5)

where, $f_{k}\equiv 1+\exp\left[\frac{ia_{0}k_{x}}{2}\right]$. The band
structure and electronic dispersion of pristine $\beta_{12}$ \- borophene can
be obtained by diagonalizing the Hamiltonian.

Figure 3: DOS for $\beta_{12}$ \- borophene using INS model in absence and in
presence of (a) NH3, (b) NO (c) NO2 and (d) CO gases.

For non-trivial solutions, we can write,

$det[h(k)-E_{k}\hat{I}]=0,\indent\lambda\in\\{1,....5\\}$ (6)

where, $E_{k}$ are the energy eigen values of the $\left[h(k)\right]$. The
determinant of the matrix can be expressed for pristine $\beta_{12}$ \-
borophene as

$E_{k}^{5}+\mathcal{Q}_{4}E_{k}^{4}+\mathcal{Q}_{3}E_{k}^{3}+\mathcal{Q}_{2}E_{k}^{2}+\mathcal{Q}_{1}E_{k}+\mathcal{R}=0$
(7)

where $\mathcal{Q}_{1}$, $\mathcal{Q}_{2}$, $\mathcal{Q}_{3}$,
$\mathcal{Q}_{4}$ and $\mathcal{R}$ are the $k$ \- dependent variable. An
extensive form of the same is displayed in Appendix A. The energy band and the
brillouin zones of $\beta_{12}$ \- borophene can be obtained by solving the
dispersion relation Eq. (7). Fig. 2, presents the band structure of
$\beta_{12}$ \- borophene in both homogeneous and the INS model. It is evident
from Fig. 2 that $\beta_{12}$ \- borophene has three conduction and two
valence bands. So, $\beta_{12}$ \- borophene has five bands. It is to be noted
that the band edges touch at the high symmetry $K^{\prime}$ and $K$ points
with coordinates $(-2\pi/3a_{0},0)$ and $(2\pi/3a_{0},0)$ indicating the
presence of Dirac fermions at $E_{k}=0$. Also, triplet fermions at X and M
points appear at $E_{k}\neq 0$. Moreover, a direct energy gap exists between
$\sim 0.25$ eV at $K^{\prime}$ and $K$ points for the INS model of
$\beta_{12}$ \- borophene. To understand the molecular gas adsorption effect
on borophene energy, we have plotted the band structure of $\beta_{12}$ \-
borophene in the presence of NH3 gas in Fig. 1 (c). It is evident that the
band gap at $K^{\prime}$, $K$ and $X$ points increases in the presence of
adsorbed NH3 molecule. Moreover, it is to be noted that the widening of the
band gap is dependent on the hopping energy $t^{\prime}$ of the boron atom
with the adsorbed gas molecule.

The hopping energies can be calculated by using the relation pourasl

$\displaystyle t_{xy}=t\left(\frac{a_{0}}{d_{xy}}\right)$ (8)

where, $t_{xy}$ and $d_{xy}$ are the hopping energy and the distance between
$\beta_{12}$ \- borophene surface and the gas molecule respectively.

The density of states (DOS) has an energy dependence which signifies the
number of allowed states per unit area at a given energy range $E$ and $E+dE$
and can be obtained by using the relation dutta .

$D(E)=\frac{\Delta n}{A\Delta E}$ (9)

where, A is the area of the $\beta_{12}$ \- borophene surface and $n$ is the
carrier concentration of the electrons. Although we obtain an analytic
expression for the DOS using Eq. (7) and Eq. (9), in presence of adsorbed gas
molecules it is too large to represent.

The carrier concentration in the $\beta_{12}$ \- borophene based gas sensor
can be obtained by using the standard relation dutta

$n(E)=\int_{0}^{\infty}D(E)f(E)dE$ (10)

where,
$f(E)=\frac{1}{1+\exp\left(\frac{E-E_{\text{F}}}{k_{\text{B}}T}\right)}$ is
the Fermi-Dirac distribution function, $E_{\text{F}}$ is the Fermi energy and
$k_{\text{B}}$ is the Boltzmann constant.

Table 1: Hopping parameter for gas molecules adsorbed on $\beta_{12}$ \- Borophene surface Adsorbed | Adsorbate distance | Hopping energy | On-site Energy
---|---|---|---
Molecule | from Borophene |  $t_{xy}$ (eV) |  $\varepsilon^{\prime}$ (eV)
| surface $d_{xy}$ (Å) |  |
NH3 | 1.63 | -2.113 | 1.11
NO | 1.38 | -2.496 | 0.95
NO2 | 1.57 | -2.194 | 1.75
CO | 1.48 | -2.167 | 1.19

Figure 4: (a) Variation of carrier concentration with $V_{\text{ds}}$ of
$\beta_{12}$ \- borophene in absence and in presence for different adsorbed
gases. (b) Quantum conductance of $\beta_{12}$ \- borophene as a function of
$V_{\text{gs}}$ in absence and in presence of adsorbed gas molecules.

##  III Results and Discussions

###  III.1 Density of states

The density of states (DOS) of $\beta_{12}$ \- borophene in presence of
adsorbed gas molecules is studied using Eq. (10) in Fig. 3. The hopping
parameters for different gases are calculated using Eq. (6) and given in Tab.
1. We observe a significant change in the DOS spectra due to the adsorption of
different gases. It is to be noted that the peak position in DOS spectra
correspond to the flatter band energies in the $E-k$ diagram. Moreover, each
extremum point can also be visualized as a Van-Hove singularity (VHS) in the
DOS. At the same time, the smooth zones in the band structure depict the
presence of localized electrons. As observed earlier from Fig. 2(c) that the
presence of adsorbed NH3 molecule drastically changes the band structure
resulting in a significant change in DOS also. The change in the peak
positions in DOS indicates the high sensitive gas sensing performance of
$\beta_{12}$ \- borophene. A similar characteristic is also observed for the
DOS for pristine $\beta_{12}$ \- borophene in the presence of adsorbed gases
as depicted in Fig. 3. We consider the INS model for all our analysis, and the
effect of the adsorbed gas molecule on DOS can be understood via the hopping
energy parameter $t_{xy}$ from Eq. (8). It is to be noted that the presence of
more number of bands in the $E-k$ diagram indicates the presence of more
degenerate bands in the DOS spectra. A similar characteristic in DOS is also
for $\beta_{12}$ \- borophene in the presence of NO molecule. However, in this
case, the VHS at $E\sim 6.5$ eV shifted towards a higher energy region, as
seen from Fig. 3(b). However, the DOS is found to be quite identical with
$\beta_{12}$ \- borophene for $E<0$ regions. The adsorption of NO2 molecule
significantly changes the band structure, resulting in a significant change in
the DOS. In this case, the VHS at $E\sim 3.5$ eV disappears completely, as
seen in Fig. 3(c). The DOS for $\beta_{12}$ \- borophene in the presence of CO
molecule is similar to that for $\beta_{12}$ \- borophene with adsorbed NO
molecule. However, the peak positions shifted towards higher energy for
adsorption of CO molecule, as seen from Figs. 3(d).

###  III.2 Carrier Concentration and Quantum Capacitance

To obtain an expression for quantum capacitance we consider the device is in
quasi-equilibrium and the carrier distribution shifted by the local
electrostatic potential. Then the charge density $Q$ of the electrons can be
written as dutta

$Q=q\int_{0}^{\infty}D(E)f\left(E+E_{\text{g}}+qV_{\text{a}}\right)dE$ (11)

where, $q$ is the magnitude of the electronic charge and $V_{\text{a}}$ is the
local electrostatic potential and $E_{\text{g}}$ is the band gap. The quantum
capacitance $C_{\text{Q}}$ of the system is defined as hassanzadazar ; yu3 ;
li4

$\displaystyle C_{\text{Q}}$ $\displaystyle=\frac{\partial Q}{\partial
V_{\text{a}}}$
$\displaystyle=\frac{q}{4k_{\text{B}}T}\int_{0}^{\infty}D(E)\text{sech}\left(\frac{E+E_{\text{g}}+qV_{\text{a}}}{2k_{\text{B}}T}\right)dE$
(12)

In Fig. 4(a), we plot the quantum capacitance ($C_{\text{Q}}$) of $\beta_{12}$
\- borophene FET against the gate-source voltage ($V_{\text{gs}}$) for various
adsorbed gas molecules, considering their different hopping energy parameters.
It is observed that at zero gate voltage the $C_{\text{Q}}$ is minimum while
it increases gradually with increase or decrease in $V_{\text{gs}}$. The
charge transfer from the adsorbed molecules results in the modification of the
DOS of $\beta_{12}$ \- borophene resulting in the change in the number of
charge carrier concentration. This change in charge carrier in turn can change
the $C_{\text{Q}}$ between the gate electrode and the $\beta_{12}$ \-
borophene surface. It is seen from the Fig. 4(b) that the $C_{\text{Q}}$
slightly increased when the adsorbed molecule was NH3 while it decreased for
NO, NO2 and CO. This change in $C_{\text{Q}}$ will result in the variations in
I-V characteristics of the $\beta_{12}$ \- borophene.

###  III.3 I-V characteristics

The I-V characteristics of the $\beta_{12}$ \- borophene FET based sensor can
be understood from quantum capacitance using the expression hassanzadazar ;
yu3 ; li4

$I=\mu C_{\text{Q}}V_{\text{ds}}E$ (13)

where, $\mu$ is the mobility of the electrons and $V_{\text{ds}}$ is the
drain-source voltage. The FET based gas sensors respond to gas adsorption by
adjusting the gate bias which controls the carrier concentration in the
channel thus tuning the amount of charge carriers in exchange with the
absorbed gas molecules. In order to confirm the influence of carrier
concentration on the gas response, the response curves of FET-based sensors
with $V_{\text{ds}}$ is shown in Fig 5. It has been clearly seen that the NO
response on the surface, increases with the increase of $V_{\text{ds}}$.

Figure 5: I-V characteristics of the $\beta_{12}$ \- borophene FET in absence
and in presence of different adsorbed gas molecules.

The current – voltage (I-V) characteristics are calculated using our developed
formalism, the results being shown in Fig 5. It reflects the effect of
molecular adsorption of NH3, NO, NO2 and CO on the quantum capacitance of
borophene. The sensitivities of the proposed FET-based sensor for NH3, NO, NO2
and CO increases with the increase in $V_{\text{ds}}$. Among the four gas
molecules, the sensitivity of NH3 is found to be highest for given
$V_{\text{ds}}$, showing better sensing performance. The change in current has
been observed after the adsorption of the gas molecules. The adsorption of the
gas molecules on the $\beta_{12}$ \- borophene surface could also change the
energy bandgap of the $\beta_{12}$ \- borophene which in turn would change the
conductivity of the $\beta_{12}$ \- borophene sensor. On the other hand, the
modulation of the concentration of the charge carriers on the borophene
surface occurs due to the charge transfer between borophene and gas molecules.
The I-V characteristics are calculated for pristine $\beta_{12}$ \- borophene
and in presence of adsorbed gas molecules using Eq. (13) and the results are
shown in Fig. 5. It reflects the effect of molecular adsorption of NH3, NO,
NO2, CO on the quantum capacitance of borophene. The sensitivities of the
proposed FET-based sensor for NH3, NO, NO2, CO increase with the increase in
$V_{\text{ds}}$. Among the four gas molecules, the sensitivity of NO is found
to be highest for given $V_{\text{ds}}$, showing better sensing performance.
The change in current has been observed after the adsorption of the gas
molecules. The adsorption of the gas molecules on the borophene surface could
also change the energy bandgap of the borophene which in turn would change the
conductivity of the borophene sensor. On the other hand, the modulation of the
concentration of the charge carriers on the borophene surface occurs due to
the charge transfer between borophene and gas molecules.

##  IV Conclusions

In this work, we have studied the effects of molecular adsorption on the
physical and electrical properties of $\beta_{12}$ \- borophene. We propose a
theoretical model based on tight binding technique to study the impact of
various adsorbed gases on the number concentration, quantum capacitance and
I-V characteristics of the FET based $\beta_{12}$ \- borophene gas sensor. The
study reveals the variation of borophene energy band gap near Fermi energy
influenced by the adsorption of gas molecules. Furthermore, we study the
density of states in the presence of different gases through the carrier
concentration and developed the quantum capacitance model of FET based
$\beta_{12}$ \- borophene gas sensor. The existence and shifting of the Van-
Hove singularities at different energy sites in presence of gas molecules
indicate the sensing performance of ability of $\beta_{12}$ \- borophene. The
study was mainly focused on the adsorption of NH3, NO2, NO and CO gas
molecules and the effects on the quantum capacitance and the I-V
characteristics of the borophene sensor was investigated. The present study
specifies a significant variation in the band gap and quantum capacitance
after gas adsorption contributing a change in the conductance, carrier
concentration and current in the borophene gas sensor. This work paves the way
for future studies of borophene as gas sensor due to its attractive gas
sensing properties.

##  Appendix A Extended form of $\mathcal{Q}$’s and $\mathcal{R}$

In this section we present an extended form of $\mathcal{Q}$’s and
$\mathcal{R}$ appearing in Eq. (7). For a pristine $\beta_{12}$ \- borophene
i.e., for $t^{\prime}=0$ and $\varepsilon^{\prime}=0$, the $\mathcal{Q}$’s and
$\mathcal{R}$ are defined as

$\mathcal{Q}_{1}=\\{16t_{\text{ab}}^{2}(t_{\text{cd}}^{2}+t_{\text{de}}^{2})+16t_{\text{bc}}^{2}t_{\text{de}}^{2}\\}\cos^{4}\left(\frac{ka}{4}\right)+[8t_{\text{ab}}t_{\text{ac}}\\{t_{\text{bc}}(\varepsilon_{\text{d}}+\varepsilon_{\text{e}})-t_{\text{bd}}t_{\text{cd}}\\}-4\varepsilon_{\text{a}}\varepsilon_{\text{b}}t_{\text{cd}}^{2}-4\varepsilon_{\text{a}}\varepsilon_{\text{b}}t_{\text{de}}^{2}+8\varepsilon_{\text{a}}t_{\text{bc}}t_{\text{bd}}t_{\text{cd}}-4t_{\text{bc}}^{2}\varepsilon_{\text{e}}(\varepsilon_{\text{a}}\\\
+\varepsilon_{\text{d}})-4\varepsilon_{\text{a}}t_{\text{bc}}^{2}\varepsilon_{\text{d}}-4\varepsilon_{\text{a}}\varepsilon_{\text{c}}t_{\text{de}}^{2}+8\varepsilon_{\text{a}}t_{\text{cd}}t_{\text{ce}}t_{\text{de}}-4\varepsilon_{\text{a}}t_{\text{cd}}^{2}\varepsilon_{\text{e}}-8t_{\text{ab}}t_{\text{ae}}(t_{\text{bc}}t_{\text{ce}}+t_{\text{bd}}t_{\text{de}})-4t_{\text{ab}}^{2}\varepsilon_{\text{e}}(\varepsilon_{\text{c}}+\varepsilon_{\text{d}})-4t_{\text{ab}}^{2}\varepsilon_{\text{c}}\varepsilon_{\text{d}}+4t_{\text{ab}}^{2}t_{\text{ce}}^{2}-8t_{\text{ac}}t_{\text{ae}}t_{\text{cd}}t_{\text{de}}+4t_{\text{ac}}^{2}t_{\text{de}}^{2}\\\
+4t_{\text{ae}}^{2}(t_{\text{bc}}^{2}+t_{\text{cd}}^{2})-4\varepsilon_{\text{b}}\varepsilon_{\text{c}}t_{\text{de}}^{2}+8\varepsilon_{\text{b}}t_{\text{cd}}t_{\text{ce}}t_{\text{de}}-4\varepsilon_{\text{b}}t_{\text{cd}}^{2}\varepsilon_{\text{e}}+8t_{\text{bc}}t_{\text{bd}}t_{\text{cd}}\varepsilon_{\text{e}}-8t_{\text{bc}}t_{\text{bd}}t_{\text{ce}}t_{\text{de}}]\cos^{2}\left(\frac{ka}{4}\right)-[\varepsilon_{\text{d}}\varepsilon_{\text{e}}\\{\varepsilon_{\text{c}}(\varepsilon_{\text{a}}+\varepsilon_{\text{b}})+\varepsilon_{\text{a}}\varepsilon_{\text{b}}\\}+\varepsilon_{\text{a}}\varepsilon_{\text{b}}\varepsilon_{\text{c}}\varepsilon_{\text{d}}\\\
+\varepsilon_{\text{a}}\varepsilon_{\text{b}}\varepsilon_{\text{c}}\varepsilon_{\text{e}}+\varepsilon_{\text{a}}\varepsilon_{\text{b}}t_{\text{ce}}^{2}-t_{\text{bd}}^{2}\varepsilon_{\text{e}}(\varepsilon_{\text{a}}+\varepsilon_{\text{c}})-\varepsilon_{\text{a}}t_{\text{bd}}^{2}\varepsilon_{\text{c}}-\varepsilon_{\text{a}}t_{\text{ce}}^{2}\varepsilon_{\text{d}}+2t_{\text{ac}}t_{\text{ae}}t_{\text{ce}}(\varepsilon_{\text{b}}+\varepsilon_{\text{d}})-t_{\text{ac}}^{2}t_{\text{ac}}^{2}t_{\text{bd}}^{2}\varepsilon_{\text{e}}-(\varepsilon_{\text{b}}+\varepsilon_{\text{d}})-t_{\text{ac}}^{2}\varepsilon_{\text{b}}\varepsilon_{\text{d}}-t_{\text{ae}}^{2}\varepsilon_{\text{d}}(\varepsilon_{\text{b}}+\varepsilon_{\text{c}})-t_{\text{ae}}^{2}\varepsilon_{\text{b}}\varepsilon_{\text{c}}\\\
+t_{\text{ae}}^{2}t_{\text{bd}}^{2}-\varepsilon_{\text{b}}t_{\text{ce}}^{2}\varepsilon_{\text{d}}+t_{\text{bd}}^{2}t_{\text{ce}}^{2}]$
(14)
$\mathcal{Q}_{2}=\\{4t_{\text{bc}}^{2}(\varepsilon_{\text{a}}+\varepsilon_{\text{d}}+\varepsilon_{\text{e}})+4\varepsilon_{\text{a}}t_{\text{cd}}^{2}+4\varepsilon_{\text{a}}t_{\text{de}}^{2}-8t_{\text{ab}}t_{\text{ac}}t_{\text{bc}}+4t_{\text{ab}}^{2}(\varepsilon_{\text{c}}+\varepsilon_{\text{d}}+\varepsilon_{\text{e}})+4\varepsilon_{\text{b}}t_{\text{cd}}^{2}+4\varepsilon_{\text{b}}t_{\text{de}}^{2}-8t_{\text{bc}}t_{\text{bd}}t_{\text{cd}}+4\varepsilon_{\text{c}}t_{\text{de}}^{2}-8t_{\text{cd}}t_{\text{ce}}t_{\text{de}}\\\
+4t_{\text{cd}}^{2}\varepsilon_{\text{e}}\\}\cos^{2}\left(\frac{ka}{4}\right)-\varepsilon_{\text{a}}\varepsilon_{\text{e}}(\varepsilon_{\text{b}}+\varepsilon_{\text{c}}+\varepsilon_{\text{d}})-\varepsilon_{\text{a}}\varepsilon_{\text{b}}\varepsilon_{\text{c}}-\varepsilon_{\text{a}}\varepsilon_{\text{b}}\varepsilon_{\text{d}}+t_{\text{bd}}^{2}(\varepsilon_{\text{a}}+\varepsilon_{\text{c}}+\varepsilon_{\text{e}})-\varepsilon_{\text{a}}\varepsilon_{\text{c}}\varepsilon_{\text{d}}+\varepsilon_{\text{a}}t_{\text{ce}}^{2}-2t_{\text{ac}}t_{\text{ae}}t_{\text{ce}}+t_{\text{ac}}^{2}(\varepsilon_{\text{b}}+\varepsilon_{\text{d}}+\varepsilon_{\text{e}})+t_{\text{ae}}^{2}(\varepsilon_{\text{b}}\\\
+\varepsilon_{\text{c}}+\varepsilon_{\text{d}})-\varepsilon_{\text{d}}\varepsilon_{\text{e}}(\varepsilon_{\text{b}}+\varepsilon_{\text{c}})-\varepsilon_{\text{b}}\varepsilon_{\text{c}}\varepsilon_{\text{d}}-\varepsilon_{\text{b}}\varepsilon_{\text{c}}\varepsilon_{\text{e}}+\varepsilon_{\text{b}}t_{\text{ce}}^{2}+t_{\text{ce}}^{2}\varepsilon_{\text{d}}$
(15)
$\mathcal{Q}_{3}=-4\cos^{2}\left(\frac{ka}{4}\right)\left(t_{\text{ab}}^{2}+t_{\text{bc}}^{2}+t_{\text{cd}}^{2}+t_{\text{de}}^{2}\right)-t_{\text{ac}}^{2}-t_{\text{ae}}^{2}-t_{\text{bd}}^{2}-t_{\text{ce}}^{2}+\varepsilon_{\text{a}}\varepsilon_{\text{b}}+\varepsilon_{\text{a}}\varepsilon_{\text{c}}+\varepsilon_{\text{a}}\varepsilon_{\text{d}}+\varepsilon_{\text{a}}\varepsilon_{\text{e}}+\varepsilon_{\text{b}}\varepsilon_{\text{c}}+\varepsilon_{\text{b}}\varepsilon_{\text{d}}+\varepsilon_{\text{b}}\varepsilon_{\text{e}}+\varepsilon_{\text{c}}\varepsilon_{\text{d}}+\varepsilon_{\text{c}}\varepsilon_{\text{e}}\\\
+\varepsilon_{\text{d}}\varepsilon_{\text{e}}$ (16)
$\mathcal{Q}_{4}=-\left(\varepsilon_{\text{a}}+\varepsilon_{\text{b}}+\varepsilon_{\text{c}}+\varepsilon_{\text{d}}+\varepsilon_{\text{e}}\right)$
(17)
$\mathcal{R}=\cos^{2}\left(\frac{ka}{4}\right)\\{4\varepsilon_{\text{a}}\varepsilon_{\text{b}}\varepsilon_{\text{c}}t_{\text{de}}^{2}-8\varepsilon_{\text{a}}\varepsilon_{\text{b}}t_{\text{cd}}t_{\text{ce}}t_{\text{de}}+4\varepsilon_{\text{a}}\varepsilon_{\text{b}}t_{\text{cd}}^{2}\varepsilon_{\text{e}}-8\varepsilon_{\text{a}}t_{\text{bc}}t_{\text{bd}}t_{\text{cd}}\varepsilon_{\text{e}}+8\varepsilon_{\text{a}}t_{\text{bc}}t_{\text{bd}}t_{\text{ce}}t_{\text{de}}+4\varepsilon_{\text{a}}t_{\text{bc}}^{2}\varepsilon_{\text{d}}\varepsilon_{\text{e}}-8t_{\text{ab}}t_{\text{ac}}t_{\text{bc}}\varepsilon_{\text{d}}\varepsilon_{\text{e}}\\\
+8t_{\text{ab}}t_{\text{ac}}t_{\text{bd}}(t_{\text{cd}}\varepsilon_{\text{e}}-t_{\text{ce}}t_{\text{de}})+8t_{\text{ab}}t_{\text{ae}}t_{\text{bc}}t_{\text{ce}}\varepsilon_{\text{d}}+8t_{\text{ab}}t_{\text{ae}}t_{\text{bd}}(\varepsilon_{\text{c}}t_{\text{de}}-t_{\text{cd}}t_{\text{ce}})-4t_{\text{ab}}^{2}\varepsilon_{\text{d}}(t_{\text{ce}}^{2}-\varepsilon_{\text{c}}\varepsilon_{\text{e}})+8t_{\text{ac}}t_{\text{ae}}\varepsilon_{\text{b}}t_{\text{cd}}t_{\text{de}}-8t_{\text{ac}}t_{\text{ae}}t_{\text{bc}}t_{\text{bd}}t_{\text{de}}-4t_{\text{ac}}^{2}\varepsilon_{\text{b}}t_{\text{de}}^{2}\\\
-4t_{\text{ae}}^{2}\varepsilon_{\text{b}}t_{\text{cd}}^{2}+8t_{\text{ae}}^{2}t_{\text{bc}}t_{\text{bd}}t_{\text{cd}}-4t_{\text{ae}}^{2}t_{\text{bc}}^{2}\varepsilon_{\text{d}}\\}+\cos^{2}\left(\frac{ka}{4}\right)^{2}-16\varepsilon_{\text{a}}t_{\text{bc}}^{2}t_{\text{de}}^{2}+32t_{\text{ab}}t_{\text{ac}}t_{\text{bc}}t_{\text{de}}^{2}-32t_{\text{ab}}t_{\text{ae}}t_{\text{bc}}t_{\text{cd}}t_{\text{de}}-16t_{\text{ab}}^{2}\varepsilon_{\text{c}}t_{\text{de}}^{2}+32t_{\text{ab}}^{2}t_{\text{cd}}t_{\text{ce}}t_{\text{de}}\\\
-16t_{\text{ab}}^{2}t_{\text{cd}}^{2}\varepsilon_{\text{e}})+\varepsilon_{\text{a}}\varepsilon_{\text{c}}\varepsilon_{\text{e}}(t_{\text{bd}}^{2}-\varepsilon_{\text{b}}\varepsilon_{\text{d}})+\varepsilon_{\text{a}}\varepsilon_{\text{b}}t_{\text{ce}}^{2}\varepsilon_{\text{d}}-\varepsilon_{\text{a}}t_{\text{bd}}^{2}t_{\text{ce}}^{2}-2t_{\text{ac}}t_{\text{ae}}\varepsilon_{\text{b}}t_{\text{ce}}\varepsilon_{\text{d}}+2t_{\text{ac}}t_{\text{ae}}t_{\text{bd}}^{2}t_{\text{ce}}+t_{\text{ac}}^{2}\varepsilon_{\text{e}}(\varepsilon_{\text{b}}\varepsilon_{\text{d}}-t_{\text{bd}}^{2})+t_{\text{ae}}^{2}\varepsilon_{\text{c}}(\varepsilon_{\text{b}}\varepsilon_{\text{d}}-t_{\text{bd}}^{2}).$
(18)

## References

  * (1) K. S. Novoselov, et. al., Nature 490, 192-200 (2012).
  * (2) M. C. Lemme, et. al., Nat. Commun. 13, 1392 (2022).
  * (3) C. Liu, et. al., Nat. Nanotechnol. 15, 545–557 (2020).
  * (4) M. Xu, et. al., Chem. Rev. 113, 3766 (2013).
  * (5) A. Quellmalz, et. al., Nat. Commun. 12, 917 (2021).
  * (6) Y. Y. U. Illarionov, et. al., Nat. Commun. 11, 3385 (2020).
  * (7) V. Passi, et. al., Adv. Mater. Interfaces 6, 1801285 (2019).
  * (8) M. Fuechsle, et. al., Nat. Nanotechnol. 7, 242–246 (2012).
  * (9) S. Acharjee, et. al., Chaos 33, 013136 (2023).
  * (10) A. A. Balandin, et. al., Nano Lett. 8(3), 902-907 (2008).
  * (11) C. Lee, et. al., Science 321(5887), 385-388 (2008).
  * (12) N. Zhang, et. al., 2D Materials 5(4), 045004 (2018).
  * (13) M. Sang, et. al., Nanomaterials (Basel). 9(3), 374. (2019).
  * (14) E. Pop, V. Varshney, and A. Roy, MRS Bulletin, 37(12), 1273-1281 (2012).
  * (15) D. G. Papageorgiou, I. A. Kinloch and R. J. Young, Prog. in Mat. Sci. 90, 75-127 (2017).
  * (16) A. H. C. Neto, Rev. Mod. Phys. 81, 109 (2009).
  * (17) E. McCann and M. Koshino, Rep. Prog. Phys. 76 056503 (2013).
  * (18) T. Ando, NPG Asia Mater 1, 17–21 (2009).
  * (19) M. Akhtar, et. al., npj 2D Mater. Appl. 1, 5 (2017).
  * (20) P. Vogt, et. al., Phys. Rev. Lett. 108(15), 155501 (2012).
  * (21) F. F. Zhu, et. al., Nat. Mater. 14(10), 1020 (2015).
  * (22) J. Ji, et. al., Nat. Commun. 7, 13352 (2016).
  * (23) I. K. Drozdov, et. al., Nat. Phys. 10, 664 (2014).
  * (24) M. Davila, et. al., New. J. Phys. 16(9), 095002 (2014).
  * (25) A. Splendiani, et. al., Nano. Lett. 10, 1271 (2010).
  * (26) A. J. Mannix, et. al., Science 350, 1513-1516 (2015).
  * (27) Z. Zhang, et. al., Sci. Adv. 5, eaax0246 (2019).
  * (28) B. Feng, et. al., Phys. Rev. B 94, 041408 (2016).
  * (29) Z. Zhang, et. al., Adv. Funct. Mater. 27, 1605059 (2017).
  * (30) S. M. Mozvashi, et al., Sci. Rep. 11, 7547 (2021).
  * (31) j. Yua, et. al., RSC Adv., 7, 8654-8660 (2017).
  * (32) X. Wang et. al., Mater. Res. Express 8, 065003 (2021).
  * (33) V. Shukla, et. al., J. Phys. Chem. C 121, 26869–26876 (2017).
  * (34) C. S. Huang, et. al., J. Phys. Chem. C 122(26), 14665-14670 (2018).
  * (35) C. Hou, C., et al. Nano Res. 15, 2537–2544 (2022).
  * (36) J. Shen, et. al., Appl. Sur. Sci. 504, 144412 (2020).
  * (37) B. Peng, et. al., Mater. Res. Lett. 5(6), 399 (2017).
  * (38) P. T. T. Le, T. C. Phong and M. Yarmohammadi, Phys. Chem. Chem. Phys. 21, 21790 (2019).
  * (39) M. M. Nobahari, RSC Adv. 12, 648 (2022).
  * (40) F. Opoku, P. P. Govender, J. Phys. Chem. A 124(11), 2288-2300 (2020).
  * (41) R. Tsu, Superlattice to Nanoelectronics (Second Edition), (2011).
  * (42) J. Xia, et. al., Nat. Nanotechnol. 4, 177 (2019).
  * (43) M. Hassanzadazar and K. Hadidi, Mater. Res. Express 7, 065006 (2020).
  * (44) G. L. Yu, et. al., Proc. Natl. Acad. Sci. 110, 3282-3286 (2013).
  * (45) J. K. Li, et. al., Appl. Phys. Lett. 122, 063102 (2023).
  * (46) A. Inhofer, Phys. Rev. Applied 9, 024022 (2018).
  * (47) M. V. Cheremisin, Physica E: Low. Dimens. Syst. Nano. Struct. 69, 153-158 (2015).
  * (48) S. Dutta, Quantum transport: Atom to Transistor, Cambridge University Press. U. K. (2005).
  * (49) A. H. Pourasl, et. al. IEEE Sensors 19, 10 (2019).
  * (50) D. L. John, L. C. Castro, and D. L. Pulfrey, J. Appl. Phys. 96, 5180 (2004).
  * (51) H. Kolavada, et. al., Electrochim. Acta 439, 141589 (2023).


[en-US]ord=raise

# Realizing a tunable honeycomb lattice in ABBA-stacked twisted double bilayer
WSe2

Haining Pan  Department of Physics, Cornell University, Ithaca, NY, USA  Eun-
Ah Kim  Department of Physics, Cornell University, Ithaca, NY, USA  Radcliffe
Institute for Advanced Studies, Cambridge, MA, USA  Department of Physics,
Harvard University, Cambridge, MA, USA  Department of Physics, Ewha Womans
University, Seoul, South Korea  Chao-Ming Jian  Department of Physics, Cornell
University, Ithaca, NY, USA

###### Abstract

The ideal honeycomb lattice, featuring sublattice and SU(2) spin rotation
symmetries, is a fundamental model for investigating quantum matters with
topology and correlations. With the rise of the moiré-based design of model
systems, realizing a tunable and symmetric honeycomb lattice system with a
narrow bandwidth can open access to new phases and insights. We propose the
ABBA-stacked twisted double bilayer WSe2 as a realistic and tunable platform
for reaching this goal. Adjusting the twist angle allows the bandwidth and the
ratio between hopping parameters of different ranges to be tuned. Moreover,
the system’s small bandwidth and spin rotation symmetry enable effective
control of the electronic structure through an in-plane magnetic field. We
construct an extended Hubbard model for the system to demonstrate this
tunability and explore possible ordered phases using the Hartree-Fock
approximation. We find that at a hole filling of $\nu=2$ (two holes per moiré
unit cell), an in-plane magnetic field of a few Tesla can “dope” the system
from a semimetal to a metal. Interactions then drive an instability towards a
canted antiferromagnetic insulator ground state. Additionally, we observe a
competing insulating phase with sublattice charge polarization. Finally, we
discuss the experimental signatures of these novel insulating phases.

##  I Introduction

The ideal honeycomb lattice with both sublattice and SU(2) spin rotation
symmetries is an important platform for studying many-body physics with
topology and electronic correlations. Graphene has been widely studied as a
material realization of this lattice. Yet, the system’s wide bandwidth limits
the exploration of the full phase space of honeycomb lattice systems.
Realizing a symmetric honeycomb lattice with a narrow bandwidth and enhanced
tunability would significantly expand the accessible phase space, enabling the
exploration of higher doping levels and stronger correlations. Among other
possibilities, this advancement would facilitate experimental investigations
of various intriguing correlated phases proposed earlier, such as charge and
spin density waves [1, 2, 3], chiral superconductivity [4, 5], and topological
Mott insulators [6]. Over the past few years, moiré superlattices in
transition metal dichalcogenides (TMD) multilayer structures have emerged as a
highly-tunable platform for quantum simulation of two-dimensional lattice
models with narrow bandwidth [7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19,
20, 21]. Triangular superlattices [10, 7, 8, 9, 12, 13, 14, 11] and honeycomb
superlattices with spin-orbit coupling (SOC) or asymmetric sublattices [14,
19, 21, 22, 23] have been realized. In this study, we aim to find a twisted
TMD multilayer structure that realizes an ideal honeycomb lattice with both
sublattice and SU(2) spin rotation symmetries.

Theoretically, a strategy to create an ideal emergent honeycomb lattice with
sublattice and SU(2) spin symmetries in TMD multilayer structures follows
three criteria as first pointed out in Refs. [24, 25, 26]: (1) The
$\Gamma$-valley of each TMD layer, where SOC is negligible, appears at the
valence band edge (rather than the $K$/$K^{\prime}$-valleys), forming an
emergent lattice for the doped holes [24, 25, 26, 27]. (2) A symmetry relates
the independent high-symmetry-stacking MX and XM sites [see Fig. 1(b)], which
form the two sublattices of the emergent honeycomb lattice; (3) The energy at
the MX and XM sites is higher than the MM site (another high-symmetry stacking
location) [28]. However, previous papers only proposed these strategies in TMD
bilayers. The correct energetics (in criteria (1) and (3)) in these proposals
are yet to be confirmed and realized by experiments. For example, Refs. 29, 30
have experimentally identified the $\Gamma$-valley physics in twisted bilayer
WSe2 but found it to appear in energies much below the valence band edge,
which is located at the $K$/$K^{\prime}$-valleys.

Here, we pursue realizing an ideal emergent honeycomb lattice in ABBA-stacked
twisted double bilayer WSe2 (tdbWSe2), i.e., two AB- and BA-stacked WSe2
bilayers with a small relative twist angle [see Fig. 1(a)]. As experimentally
shown in a closely related system, the ABAB-stacked tdbWSe2, the increase in
the number of layers enhances the interlayer hybridization driving the
$\Gamma$-valley to the valence band edge [27, 31]. However, the ABAB-stacked
tdbWSe2 does not have the symmetry needed for criterion (2) and only realizes
a $\Gamma$-valley-based emergent triangular lattice. In this paper, we show
that the ABBA-stacked tdbWSe2 can pass the three criteria and realize a
tunable ideal honeycomb lattice with both sublattice and SU(2) spin rotation
symmetries.

The remainder of this paper is organized as follows. In Sec. II, we provide
the microscopic analysis that the ABBA-stacked tdbWSe2 can pass the three
criteria above and, hence, leads to an emergent honeycomb lattice with both
sublattice and SU(2) spin rotation symmetries. We construct a continuum model
for this tdbWSe2 system, enabling estimations of the hopping parameters of the
emergent honeycomb lattice. In Sec. III, we demonstrate the intriguing
tunability of this system by studying orders induced by an in-plane magnetic
field at the hole filling $\nu=2$, namely, the half-filling of the emergent
honeycomb lattice. The phase diagram is obtained using the Hartree–Fock
analysis of the extended Hubbard model on the emergent honeycomb lattice.
Finally, in Sec. IV, we summarize our results and experimental signatures of
the predicted orders.

##  II Modeling the emergent honeycomb lattice

The ABBA-stacked tdbWSe2 can satisfy the three criteria above for the
following reasons. First, this structure’s increased number of layers leads to
stronger interlayer tunnelings compared to a single bilayer, pushing the
$\Gamma$-valley to the top of the valence bands [32, 33]. Second, inherited
from the $D_{3h}$ point group symmetry of each WSe2 layer, the ABBA-stacked
tdbWSe2 manifests a three-fold rotation $C_{3z}$ around the $z$-axis, and a
two-fold rotation $C_{2y}$ around the in-plane $y$-axis, where the $y$-axis is
chosen to align with the $xy$-plane projection of a bond between the tungsten
atom and the selenium atom in the monolayer WSe2. These crystal symmetries
ensure a sublattice-symmetric honeycomb lattice composed of the MX and XM
sites. Third, from the ab initio tight-binding model in Appendix A, we find
the correct energy hierarchy amongst the MX, XM, and MM sites to ensure an
emergent honeycomb lattice.

To estimate the hopping parameters in the emergent honeycomb lattice based on
the three criteria, we need to obtain the dispersion of the two topmost moiré
valence bands. Adapting the approach in Ref. [31, 25], we begin with a four-
layer continuum Hamiltonian for our twisted ABBA-stacked system that preserves
$C_{2y}$, $C_{3z}$, and SU(2) spin rotation symmetries (due to negligible SOC
near the $\Gamma$-valley [34]):

$H_{\text{4L}}=-\frac{\hbar^{2}\bm{k}^{2}}{2m}+\begin{pmatrix}\Delta_{1}(\bm{r})&\Delta_{12}(\bm{r})&0&0\\\
\Delta_{12}^{\dagger}(\bm{r})&\Delta_{2}(\bm{r})&\Delta_{23}(\bm{r})&0\\\
0&\Delta_{23}^{\dagger}(\bm{r})&\Delta_{3}(\bm{r})&\Delta_{34}(\bm{r})\\\
0&0&\Delta_{34}^{\dagger}(\bm{r})&\Delta_{4}(\bm{r})\end{pmatrix}.$ (1)

$H_{\text{4L}}$ is identical for both spin species. The intralayer potentials
take the forms $\Delta_{1}(\bm{r})=\Delta_{4}(\bm{r})=V_{1}$ for the first and
fourth layer, and
$\Delta_{2}(\bm{r})=V_{2}^{(0)}+2V_{2}^{(1)}\sum\limits_{i=1,3,5}\cos(\bm{G}_{i}\cdot\bm{r}+\phi)$
and
$\Delta_{3}(\bm{r})=V_{2}^{(0)}+2V_{2}^{(1)}\sum\limits_{i=1,3,5}\cos(\bm{G}_{i}\cdot\bm{r}-\phi)$
for the second and third layer (the opposite sign of $\phi$ as a result of
$C_{2y}$ symmetry). The interlayer tunnelings take the forms
$\Delta_{12}(\bm{r})=\Delta_{34}(\bm{r})=V_{12}$ and
$\Delta_{23}(\bm{r})=V_{23}^{(0)}+2V_{23}^{(1)}\sum\limits_{i=1,3,5}\cos(\bm{G}_{i}\cdot\bm{r})$.
Here,
$\bm{G}_{i}=\frac{4\pi}{\sqrt{3}a_{M}}\left(\cos(\frac{(i-1)\pi}{3}),\sin(\frac{(i-1)\pi}{3})\right)$
for $i=1,2,...,6$ are the six first-shell moiré reciprocal lattice vectors
(see Fig. 1 (c)). The moiré lattice constant $a_{M}$ is determined by the
small twist angle $\theta$ through $a_{M}\approx a_{0}/\theta$ with the
lattice constant $a_{0}$ of monolayer WSe2 being 3.28Å [35]. At $\theta\sim
2^{\circ}$, $a_{M}\approx 10$ nm. The moiré reciprocal lattice vectors
$\bm{G}_{i}$ only show up in $\Delta_{2}=\Delta_{3}$ and $\Delta_{23}$ because
the small-angle twisting only appears between the second and third layer of
WSe2 in the twisted double bilayer structure.

We now estimate the parameters in the intralayer potentials and interlayer
tunnelings in Eq. (1). These potentials and tunnelings in the tdbWSe2
evaluated at the MM, and MX/XM sites should be fitted according to the band
structures of the _untwisted_ ABBA-stacked double bilayer WSe2 in the three
high-symmetry stacking configurations MM, and MX/XM [see Appendix A for more
details] [36]. The angle-independent parameters for the potentials and
tunnelings are estimated as follows:
$(V_{1},V_{2}^{(0)},V_{2}^{(1)})=(200,-159,-8)$ meV, $\phi=-0.17$, and
$(V_{12},V_{23}^{(0)},V_{23}^{(1)})=(184,356,-9)$ meV. We comment that these
parameter estimates are obtained based on assumptions on the interlayer
distances in the double bilayer WSe2 detailed in App. A. The assumed
interlayer distances produce a representative example of tdbWSe2 satisfying
the energetics criteria allowing us to estimate bandwidth and demonstrate the
tunability of the emergent honeycomb lattice, which are the main purposes of
our work. With the above parameters, we can solve Hamiltonian (1) for the
moiré valence bands of tdbWSe2. For example, Figure 1(d) shows the band
structure along the $\gamma$-$m$-$\kappa$-$\gamma$ path in the moiré Brillouin
zone (mBZ) at a twist angle of $\theta=2^{\circ}$. From the two topmost moiré
valence bands, especially the presence of a Dirac cone between them, we
confirm that ABBA-stacked tdbWSe2 with a hole filling less than 4 per moiré
unit cell captures an emergent honeycomb lattice with both sublattice and
SU(2) spin rotation symmetries. At a twist angle of $\theta=2^{\circ}$, the
bandwidth is around 10 meV (which is much smaller than the $\sim$2 eV
bandwidth of graphene [37]).

Figure 1:  (a) The schematic picture of the ABBA-stacked tdbWSe2. Within each
of the AB- and BA-stacked bilayer WSe2, the orientation of the two layers
differs by 180∘. The two bilayers are twisted so that the middle two layers
have a small twist angle $\theta$ between them. (b) The real-space moiré
pattern created by the twisted central two layers (blue and orange) with the
darker (lighter) for the metal (chalcogenide) atom. The three high-symmetry-
stacking sites in one moiré unit cell (black dashed lines) are MM (brown,
metal atoms on one layer aligned with metal atoms on the adjacent layer), XM
(blue, chalcogenide atom on one layer aligned with metal atoms on the adjacent
layer), and MX (grey, metal atom on one layer aligned with chalcogenide atoms
on the adjacent layer). Here, the alignment of atoms refers to the central two
layers in the double bilayer structure; (c) Moiré Brillouin zone with four
high-symmetry points $\gamma$, $m$, $\kappa$, and $\kappa^{\prime}$.
$G_{1,..,6}$ are the moiré reciprocal vectors; (d) Moiré band structure
obtained from the continuum model at twist angle $\theta=2^{\circ}$ along the
path $\gamma$-$m$-$\kappa$-$\gamma$ [blue dashed line shown in (c)]. The red
dashed line indicates Fermi energy at hole filling $\nu=2$ (two holes per
morié unit cell).

At hole filling $\nu=2$, i.e., two holes per moiré unit cell, the Fermi level
is precisely at the Dirac point protected by $C_{3z}$ and $C_{2y}$ symmetries.
This results in a semimetal (noninteracting) ground state, as shown in Fig.
1(d). The narrow bandwidth of the tdbWSe2-realized honeycomb lattice allows
for an intriguing way to tune the electronic structure: a realistic in-plane
magnetic field can “dope” the semimetal and create finite Fermi surfaces
located around both $\kappa$ and $\kappa^{\prime}$ points of the mBZ for both
spin species. Note that such tuning is impossible for the
$K$/$K^{\prime}$-valley-based TMD moiré system due to the strong SOC therein.
Without loss of generality, we assume the in-plane magnetic field $B$ is
applied along the $x$ direction, resulting in a Zeeman energy of
$h_{x}=\frac{1}{2}g\mu_{B}B$, where $\mu_{B}$ is the Bohr magneton
($5.79\times 10^{-2}$ meV/T,) and the $g$-factor is roughly 2 due to the
negligible SOC near the $\Gamma$-valley. For example, at a twist angle of
$\theta=2^{\circ}$, a magnetic field of 13 T, typically achievable
experimentally, creates a Zeeman splitting energy of 0.8 meV. This Zeeman
splitting leads to a finite Fermi surface of each spin species roughly halfway
in energy between the Dirac point (at the $\kappa$ point in the mBZ) and the
van Hove singularity (at the $m$ point), as shown in Fig. 2. In the presence
of the in-plane field, a magnetic particle-hole instability that opens a
charge gap is expected even at weak interactions because of the nesting
between the Fermi surfaces of the two spin species. This instability would
lead to a magnetic order that spontaneously breaks the remaining U(1) spin
rotation symmetry about the $x$-axis, namely, the in-plane field direction.

Conceptually, a similar in-plane-field-induced particle-hole instability
applies to graphene. The previous work [3] shows that the resulting state is a
canted Néel antiferromagnetic (AF) insulator, whose magnetizations projected
onto the plane perpendicular to the applied field are opposite in the two
sublattices of the honeycomb. However, realizing such a state in graphene
requires a very high in-plane field of $10^{2}\sim 10^{3}$ T [3]. As we show
below, a similar canted Néel AF insulator can be found in the tdbWSe2-realized
honeycomb lattice, whose bandwidth is engineered around the meV scale, at a
realistic in-plane field of a few Tesla.

Figure 2:  The two topmost moiré valence bands shown in Fig. 1(d) split by the
Zeeman energy shift due to an in-plane magnetic field of 13 T. The arrows
indicate the spin along the in-plane field direction. The red dashed line
indicates the Fermi energy at hole filling $\nu=2$, i.e., two holes per moiré
unit cell and half-filling for the emergent honeycomb lattice.

##  III Extended Hubbard model and phase diagram under in-plane field

Aside from the bandwidth distinction with graphene, the band structure of the
tdbWSe2-realized honeycomb lattice exhibits a more pronounced contribution
from the next-nearest-neighbor (NNN) hopping $t_{2}$ on the honeycomb lattice.
The effect of $t_{2}$ is manifested by the asymmetry (about the Dirac cone)
between the two topmost valence bands shown in Fig. 1 (d). In contrast to the
graphene-inspired model of [3] with only the nearest-neighbor (NN) hopping
$t_{1}$ and the on-site Hubbard interaction, we consider an extended Hubbard
model with longer-range hoppings and interactions pertaining to our
tdbWSe2-realized honeycomb lattice system. Our Hamiltonian consists of three
terms: the hopping $H_{\rm h}$, the Zeeman energy $H_{\rm Z}$ from the in-
plane magnetic field, and the extended Hubbard interaction $H_{\text{int}}$.
The term $H_{\rm
h}=-t_{1}\sum_{\expectationvalue{ij}}\sum_{\sigma}c_{i\sigma}^{\dagger}c_{j\sigma}-t_{2}\sum_{\expectationvalue{\expectationvalue{ij}}}\sum_{\sigma}c_{i\sigma}^{\dagger}c_{j\sigma}$
contains both NN hopping $t_{1}$ and NNN hoppings $t_{2}$. Here, $c_{i\sigma}$
is the fermion operator of the hole with spin $\sigma$ at the site $i$ of the
emergent honeycomb lattice. Both $t_{1}$ and $t_{2}$ are tunable by the twist
angle $\theta$. For instance, at $\theta=2^{\circ}$, we estimate $t_{1}=1.8$
meV and $t_{2}=0.2$ meV based on the moiré band structure obtained from the
continuum model Eq.(1). The ratio $t_{2}/t_{1}$ varies as we change the twist
angle: 0.018 for $\theta=1^{\circ}$, 0.11 for $\theta=2^{\circ}$, and 0.12 for
$\theta=3^{\circ}$. In the search for possible phases below, we treat
$t_{2}/t_{1}$ as a parameter to explore the effect of longer-range hopping.
The Zeeman energy $H_{\rm Z}$ is given by $H_{\rm
Z}=h_{x}\sum_{i}\left(c_{i\uparrow}^{\dagger}c_{i\uparrow}-c_{i\downarrow}^{\dagger}c_{i\downarrow}\right)$,
where the spin quantization axis conveniently is chosen along the $x$
direction. Under $H_{\rm Z}$, the spin symmetry is reduced to a U(1) spin
rotation in the $yz$ plane. We fix $h_{x}=0.4t_{1}$ in the following. For the
twist angle $\theta=2^{\circ}$, $h_{x}=0.4t_{1}$ (corresponding to a 13 T in-
plane magnetic field) sets the Fermi level roughly halfway between the Dirac
point at the $\kappa/\kappa^{\prime}$ point and the van Hove singularity at
the $m$ point of the mBZ.

The range of interactions in TMD-based superlattices typically extends beyond
on-site interactions. Here, we consider the extended Hubbard interaction

$H_{\text{int}}=U_{0}\sum_{i}c_{i\uparrow}^{\dagger}c_{i\uparrow}c_{i\downarrow}^{\dagger}c_{i\downarrow}+U_{1}\sum_{\langle
ij\rangle}\sum_{\sigma,\sigma^{\prime}=\uparrow,\downarrow}c_{i,\sigma}^{\dagger}c_{i,\sigma}c_{j,\sigma^{\prime}}^{\dagger}c_{j,\sigma^{\prime}},$
(2)

with the on-site interaction strength $U_{0}$, and the NN Hubbard interaction
of strength $U_{1}$. In the following study of the phase diagram, we treat
$U_{0}$ and $U_{1}$ as two free parameters to explore their effect on the
canted Néel AF insulator and other phases.

Figure 3:  (a) The color wheel shows the six spin orientations in the $yz$
plane. (b) The phase diagram as a function of $U_{0}$ and $t_{2}$. (c) The
phase diagram as a function of $U_{0}$ and $U_{1}$. Here, the size of the dot
indicates the number density of each site, while the color indicates the spin
orientation in the $yz$ plane as shown in the color wheel. The black dots in
the sublattice polarized state indicate zero spin polarization.

We start with understanding the effect of the NNN hopping $t_{2}$, and present
the ground-state phase diagram obtained using the Hartree-Fock method in Fig.
3(b). In this phase diagram, we consider a parameter range
$U_{0}\in[t_{1},9t_{1}]$ and $t_{2}\in[0,0.6t_{1}]$ while setting $U_{1}=0$.
Using the numerical Hartree-Fock method, we find three possible phases: the
canted Néel AF insulator, canted $\sqrt{3}\times\sqrt{3}$ spin density wave
(SDW), and a metal. At hole filling $\nu=2$, the Fermi surfaces of the two
spin species are always perfectly nested for $t_{2}\lesssim t_{1}/3$ at both
the $\kappa$ and $\kappa^{\prime}$ points of the mBZ. This perfect nesting
results in a Fermi surface instability towards a robust canted Néel AF
insulator, spontaneously breaking the U(1) spin rotation symmetry in the $yz$
plane as mentioned above. Beyond $t_{2}>t_{1}/3$, the system’s ground state is
metallic at small $U_{0}$. The metallic state is due to an extra Fermi surface
of the spin-down hole appearing around the $\gamma$ point when
$t_{2}>t_{1}/3$. Consequently, the Fermi surfaces of the two spin species are
no longer nested at $\nu=2$. As $U_{0}$ increases, the system can enter an
insulating SDW phase with an enlarged $\sqrt{3}\times\sqrt{3}$ unit cell. This
SDW phase exhibits a 120∘ AF order in the $yz$ plane on each sublattice of the
honeycomb. The chiralities of the two 120∘ AF orders are opposite. For the
tdbWSe2-realized honeycomb lattice, $t_{2}/t_{1}<1/3$ for twist angle
$\theta<4^{\circ}$. Hence, restricted to only on-site interactions, the canted
Néel AF insulator is the most relevant in-plane-field-induced phase for this
material at hole filling $\nu=2$.

In addition to the NNN hopping $t_{2}$, we also study the effect of the NN
Hubbard interaction $U_{1}$. We present the phase diagram obtained from the
Hartree-Fock calculation as a function of $U_{0}$ and $U_{1}$ as shown in Fig.
3(c). For this phase diagram, we set $t_{2}=0$ for simplicity. The canted Néel
AF insulator is robust at small $U_{1}$. As $U_{1}$ increases, the system is
driven to an insulating phase with the charge distribution polarized to one of
the sublattices, spontaneously breaking the system’s $C_{2y}$ symmetry. In
contrast to the canted Néel AF insulator, the $yz$-plane U(1) spin rotation
symmetry is preserved in this phase. The phase boundary between the AF
insulator and the sublattice-polarized insulator is roughly at
$U_{1}\sim\frac{1}{3}U_{0}+\frac{2}{3}h_{x}$. This is the exact phase boundary
in the strong interaction limit (where the hopping terms $H_{\rm h}$ are
neglected). In the tdbWSe2-realized honeycomb lattice, the screening of the
Coulomb interaction determines the relative strength between $U_{0}$ and
$U_{1}$. Stronger screening favors the canted Néel AF insulator, while weaker
screening favors the sublattice-polarized insulator. In Appendix B, we present
an analytical mean-field analysis in the weak-coupling limit $U_{0,1}\ll
t_{1}$ and find that the canted $\sqrt{3}\times\sqrt{3}$ SDW insulator can
also compete with the canted Néel AF insulator when $U_{1}/U_{0}$ increases.

##  IV Summary and Experimental Implication

In summary, we show that the ABBA-stacked tdbWSe2 can serve as a realistic and
tunable platform to simulate $\Gamma$-valley honeycomb lattice with both
sublattice and SU(2) spin rotation symmetries. We develop a continuum model to
estimate the relevant hopping parameters of this honeycomb lattice model,
including the NN and the NNN hoppings $t_{1,2}$. The small bandwidth, found to
be at the meV scale, and the ratio $t_{2}/t_{1}$ are both tunable via the
twist angle $\theta$. We show that the small bandwidth enables the effective
control of the electronic structure via an in-plane magnetic field. At a hole
filling $\nu=2$, an in-plane field of a few Teslas can dope this system from a
Dirac semimetal to a metal having finite Fermi surfaces with instabilities. To
demonstrate this tunability and understand the resulting instabilities, we
construct an extended Hubbard model for this honeycomb lattice and perform a
numerical Hartree-Fock calculation for the ground-state phase diagram at
$\nu=2$. We find a robust canted Néel AF insulating phase when the on-site
Hubbard interaction $U_{0}$ dominates. This phase is conceptually similar to
the magnetic-field-induced canted AF insulator studied in the context of
graphene by Ref. [3]. However, our finding resides in a more realistic
parameter regime. Moreover, a competing sublattice-polarized insulator is
found in our extended Hubbard model when the NN Hubbard interaction $U_{1}$
increases beyond a threshold.

The results presented in this work have direct experimental implications.
First, for ABBA-stacked tdbWSe2 at hole filling $\nu=2$, an in-plane magnetic
field is predicted to induce an insulating ground state. Consequently, at a
finite temperature, we expect a large magnetoresistance as a function of the
in-plane field around the Dirac semimetal state (at zero field). There are two
possible insulating states: the canted Néel AF insulator and the sublattice-
polarized insulator. The former can be detected through its spin texture
perpendicular to the in-plane field direction, potentially using spin-
polarized scanning tunneling microscopy [38]. The finite-temperature phase
transition into the canted Néel AF order is a Berezinskii–Kosterlitz–Thouless
transition associated with the remaining U(1) spin rotation symmetry under the
in-plane field. The sublattice-polarized insulator has no non-trivial spin
texture (other than the magnetization along the in-plane field direction).
Instead, the spontaneous charge polarization in the two sublattices of the
emergent honeycomb lattice leads to a finite electric polarization in the
tdbWSe2 sample along the $z$ direction. Consequently, a signature of this
sublattice-polarized insulating phase is the hysteresis of the $z$-axis
electric polarization as a function of the $z$-direcction electric field. The
finite-temperature phase transition into the sublattice-polarized phase is in
the Ising universality class. Switching between the two insulating phases can
be implemented by tuning the ratio between $U_{0}$ and $U_{1}$, which is
achievable via adjusting the twist angle and/or sample’s distance from the
gates.

Acknowledgements. We thank Kin Fai Mak, Jie Shan, and Liguo Ma for stimulating
discussions. H.P. acknowledges the helpful discussions with Fengcheng Wu, Ming
Xie, and Dan Mao. H.P. and E-A.K. acknowledge the funding support National
Science Foundation (Platform for the Accelerated Realization, Analysis, and
Discovery of Interface Materials (PARADIM)) under Cooperative Agreement No.
DMR-1539918. C.-M.J. is supported by a faculty startup grant at Cornell
University.

## References

  * Honerkamp [2008] C. Honerkamp, Density Waves and Cooper Pairing on the Honeycomb Lattice, Phys. Rev. Lett. 100, 146404 (2008).
  * Aleiner  _et al._ [2007] I. L. Aleiner, D. E. Kharzeev, and A. M. Tsvelik, Spontaneous symmetry breaking in graphene subjected to an in-plane magnetic field, Phys. Rev. B 76, 195415 (2007).
  * Bercx  _et al._ [2009] M. Bercx, T. C. Lang, and F. F. Assaad, Magnetic field induced semimetal-to-canted-antiferromagnet transition on the honeycomb lattice, Phys. Rev. B 80, 045412 (2009).
  * Nandkishore  _et al._ [2012] R. Nandkishore, L. S. Levitov, and A. V. Chubukov, Chiral superconductivity from repulsive interactions in doped graphene, Nature Physics 8, 158 (2012).
  * Nandkishore  _et al._ [2014] R. Nandkishore, R. Thomale, and A. V. Chubukov, Superconductivity from weak repulsion in hexagonal lattice systems, Phys. Rev. B 89, 144501 (2014).
  * Raghu  _et al._ [2008] S. Raghu, X.-L. Qi, C. Honerkamp, and S.-C. Zhang, Topological Mott Insulators, Phys. Rev. Lett. 100, 156401 (2008).
  * Regan  _et al._ [2020] E. C. Regan, D. Wang, C. Jin, M. I. Bakti Utama, B. Gao, X. Wei, S. Zhao, W. Zhao, Z. Zhang, K. Yumigeta, M. Blei, J. D. Carlström, K. Watanabe, T. Taniguchi, S. Tongay, M. Crommie, A. Zettl, and F. Wang, Mott and generalized Wigner crystal states in WSe2/WS2 moiré superlattices, Nature 579, 359 (2020).
  * Tang  _et al._ [2020] Y. Tang, L. Li, T. Li, Y. Xu, S. Liu, K. Barmak, K. Watanabe, T. Taniguchi, A. H. MacDonald, J. Shan, and K. F. Mak, Simulation of Hubbard model physics in WSe2/WS2 moiré superlattices, Nature 579, 353 (2020).
  * Xu  _et al._ [2020] Y. Xu, S. Liu, D. A. Rhodes, K. Watanabe, T. Taniguchi, J. Hone, V. Elser, K. F. Mak, and J. Shan, Correlated insulating states at fractional fillings of moiré superlattices, Nature 587, 214 (2020).
  * Wang  _et al._ [2020] L. Wang, E.-M. Shih, A. Ghiotto, L. Xian, D. A. Rhodes, C. Tan, M. Claassen, D. M. Kennes, Y. Bai, B. Kim, K. Watanabe, T. Taniguchi, X. Zhu, J. Hone, A. Rubio, A. N. Pasupathy, and C. R. Dean, Correlated electronic phases in twisted bilayer transition metal dichalcogenides, Nature Materials 19, 861 (2020).
  * Ghiotto  _et al._ [2021] A. Ghiotto, E.-M. Shih, G. S. S. G. Pereira, D. A. Rhodes, B. Kim, J. Zang, A. J. Millis, K. Watanabe, T. Taniguchi, J. C. Hone, L. Wang, C. R. Dean, and A. N. Pasupathy, Quantum criticality in twisted transition metal dichalcogenides, Nature 597, 345 (2021).
  * Zhou  _et al._ [2021] Y. Zhou, J. Sung, E. Brutschea, I. Esterlis, Y. Wang, G. Scuri, R. J. Gelly, H. Heo, T. Taniguchi, K. Watanabe, G. Zaránd, M. D. Lukin, P. Kim, E. Demler, and H. Park, Bilayer Wigner crystals in a transition metal dichalcogenide heterostructure, Nature 595, 48 (2021).
  * Li  _et al._ [2021a] H. Li, S. Li, E. C. Regan, D. Wang, W. Zhao, S. Kahn, K. Yumigeta, M. Blei, T. Taniguchi, K. Watanabe, S. Tongay, A. Zettl, M. F. Crommie, and F. Wang, Imaging two-dimensional generalized Wigner crystals, Nature 597, 650 (2021a).
  * Li  _et al._ [2021b] T. Li, S. Jiang, B. Shen, Y. Zhang, L. Li, Z. Tao, T. Devakul, K. Watanabe, T. Taniguchi, L. Fu, J. Shan, and K. F. Mak, Quantum anomalous Hall effect from intertwined moiré bands, Nature 600, 641 (2021b).
  * Gu  _et al._ [2021] J. Gu, L. Ma, S. Liu, K. Watanabe, T. Taniguchi, J. C. Hone, J. Shan, and K. F. Mak, Dipolar excitonic insulator in a moire lattice, arXiv:2108.06588 (2021).
  * Zhao  _et al._ [2022] W. Zhao, K. Kang, L. Li, C. Tschirhart, E. Redekop, K. Watanabe, T. Taniguchi, A. Young, J. Shan, and K. F. Mak, Realization of the Haldane Chern insulator in a moir\’e lattice (2022).
  * Xie  _et al._ [2022a] Y.-M. Xie, C.-P. Zhang, J.-X. Hu, K. F. Mak, and K. T. Law, Valley-Polarized Quantum Anomalous Hall State in Moiré ${\mathrm{MoTe}}_{2}/{\mathrm{WSe}}_{2}$ Heterobilayers, Phys. Rev. Lett. 128, 026402 (2022a).
  * Xie  _et al._ [2022b] Y.-M. Xie, C.-P. Zhang, and K. T. Law, Topological $p_{x}+ip_{y}$ inter-valley coherent state in Moir\’e MoTe2/WSe2 heterobilayers (2022b).
  * Tao  _et al._ [2022] Z. Tao, B. Shen, S. Jiang, T. Li, L. Li, L. Ma, W. Zhao, J. Hu, K. Pistunova, K. Watanabe, T. Taniguchi, T. F. Heinz, K. F. Mak, and J. Shan, Valley-coherent quantum anomalous Hall state in AB-stacked MoTe2/WSe2 bilayers (2022).
  * Liu  _et al._ [2021] E. Liu, T. Taniguchi, K. Watanabe, N. M. Gabor, Y.-T. Cui, and C. H. Lui, Excitonic and Valley-Polarization Signatures of Fractional Correlated Electronic Phases in a ${\mathrm{WSe}}_{2}/{\mathrm{WS}}_{2}$ Moiré Superlattice, Phys. Rev. Lett. 127, 037402 (2021).
  * Foutty  _et al._ [2023] B. A. Foutty, C. R. Kometter, T. Devakul, A. P. Reddy, K. Watanabe, T. Taniguchi, L. Fu, and B. E. Feldman, Mapping twist-tuned multi-band topology in bilayer WSe2 (2023).
  * Zhang  _et al._ [2021a] Y. Zhang, T. Devakul, and L. Fu, Spin-textured Chern bands in AB-stacked transition metal dichalcogenide bilayers, PNAS 118, 10.1073/pnas.2112673118 (2021a).
  * Wu  _et al._ [2019] F. Wu, T. Lovorn, E. Tutuc, I. Martin, and A. H. MacDonald, Topological insulators in twisted transition metal dichalcogenide homobilayers, Phys. Rev. Lett. 122, 086402 (2019).
  * Angeli and MacDonald [2021] M. Angeli and A. H. MacDonald, $\Gamma$ valley transition metal dichalcogenide moiré bands, Proceedings of the National Academy of Sciences 118, e2021826118 (2021).
  * Zhang  _et al._ [2021b] Y. Zhang, T. Liu, and L. Fu, Electronic structures, charge transfer, and charge order in twisted transition metal dichalcogenide bilayers, Phys. Rev. B 103, 155142 (2021b).
  * Xian  _et al._ [2021] L. Xian, M. Claassen, D. Kiese, M. M. Scherer, S. Trebst, D. M. Kennes, and A. Rubio, Realization of nearly dispersionless bands with strong orbital anisotropy from destructive interference in twisted bilayer MoS2, Nat Commun 12, 5644 (2021).
  * Movva  _et al._ [2018] H. C. P. Movva, T. Lovorn, B. Fallahazad, S. Larentis, K. Kim, T. Taniguchi, K. Watanabe, S. K. Banerjee, A. H. MacDonald, and E. Tutuc, Tunable $\mathrm{\Gamma}-K$ Valley Populations in Hole-Doped Trilayer ${\mathrm{WSe}}_{2}$, Phys. Rev. Lett. 120, 107703 (2018).
  * Pan  _et al._ [2020] H. Pan, F. Wu, and S. Das Sarma, Band topology, Hubbard model, Heisenberg model, and Dzyaloshinskii-Moriya interaction in twisted bilayer ${\mathrm{WSe}}_{2}$, Phys. Rev. Research 2, 033087 (2020).
  * Gatti  _et al._ [2022] G. Gatti, J. Issing, L. Rademaker, F. Margot, T. A. de Jong, S. J. van der Molen, J. Teyssier, T. K. Kim, M. D. Watson, C. Cacho, P. Dudin, J. Avila, K. C. Edwards, P. Paruch, N. Ubrig, I. Gutiérrez-Lezama, A. Morpurgo, A. Tamai, and F. Baumberger, Observation of flat $\Gamma$ moir\’e bands in twisted bilayer WSe2 (2022).
  * Pei  _et al._ [2022] D. Pei, B. Wang, Z. Zhou, Z. He, L. An, S. He, C. Chen, Y. Li, L. Wei, A. Liang, J. Avila, P. Dudin, V. Kandyba, A. Giampietri, M. Cattelan, A. Barinov, Z. Liu, J. Liu, H. Weng, N. Wang, J. Xue, and Y. Chen, Observation of $\Gamma$ -Valley Moiré Bands and Emergent Hexagonal Lattice in Twisted Transition Metal Dichalcogenides, Phys. Rev. X 12, 021065 (2022).
  * Foutty  _et al._ [2022] B. A. Foutty, J. Yu, T. Devakul, C. R. Kometter, Y. Zhang, K. Watanabe, T. Taniguchi, L. Fu, and B. E. Feldman, Tunable spin and valley excitations of correlated insulators in $\Gamma$-valley moir\’e bands, arXiv.2206.10631 10.48550/arXiv.2206.10631 (2022).
  * Mak  _et al._ [2010] K. F. Mak, C. Lee, J. Hone, J. Shan, and T. F. Heinz, Atomically Thin ${\mathrm{MoS}}_{2}$: A New Direct-Gap Semiconductor, Phys. Rev. Lett. 105, 136805 (2010).
  * Xu  _et al._ [2014] X. Xu, W. Yao, D. Xiao, and T. F. Heinz, Spin and pseudospins in layered transition metal dichalcogenides, Nature Phys 10, 343 (2014).
  * Fang  _et al._ [2015] S. Fang, R. Kuate Defo, S. N. Shirodkar, S. Lieu, G. A. Tritsaris, and E. Kaxiras, Ab initio tight-binding Hamiltonian for transition metal dichalcogenides, Phys. Rev. B 92, 205108 (2015).
  * He  _et al._ [2014] J. He, K. Hummer, and C. Franchini, Stacking effects on the electronic and optical properties of bilayer transition metal dichalcogenides ${\mathrm{MoS}}_{2}$, ${\mathrm{MoSe}}_{2}$, ${\mathrm{WS}}_{2}$, and ${\mathrm{WSe}}_{2}$, Phys. Rev. B 89, 075409 (2014).
  * Jung  _et al._ [2014] J. Jung, A. Raoux, Z. Qiao, and A. H. MacDonald, Ab initio theory of moir\’e superlattice bands in layered two-dimensional materials, Phys. Rev. B 89, 205414 (2014).
  * Jung and MacDonald [2013] J. Jung and A. H. MacDonald, Tight-binding model for graphene $\pi$-bands from maximally localized Wannier functions, Phys. Rev. B 87, 195450 (2013).
  * Pietzsch  _et al._ [2001] O. Pietzsch, A. Kubetzka, M. Bode, and R. Wiesendanger, Observation of Magnetic Hysteresis at the Nanometer Scale by Spin-Polarized Scanning Tunneling Spectroscopy, Science 292, 2053 (2001).

##  Appendix A Band structure calculation and parameter estimate

To show that criteria (1) and (3) introduced in the main text can be satisfied
by the tdbWSe2, the first task is to obtain its band structure. However, it is
not easy to derive the band structure of tdbWSe2 directly at the atomic level
because the emergent moiré superlattice contains up to ten thousand atoms even
within just one moiré unit cell. One way to approximate the band structure of
the twisted system is to use a continuum model with parameters fitted
according to the untwisted system in the high-symmetry stacking configurations
(MM, MX, and XM) that appear in the moiré unit cell of the twisted system
[36]. We will provide the details later in this section. Therefore, in this
section, we will first discuss the method to obtain the band structure of
untwisted systems.

We calculate the band structure of the untwisted system using the ab initio
tight-binding Hamiltonian [34]. The logic is to first construct a tight-
binding model for the monolayer by considering five $d$ orbitals of one
tungsten atom and six $p$ orbitals of the two selenium atoms in the monolayer
unit cell, and including spin-orbit coupling over the entire monolayer
Brillouin zone (BZ) (i.e., not just the small mBZ). It is called ab initio
because hopping parameters are evaluated from the density functional theory.
After constructing the tight-binding model for the monolayer, we can model the
interlayer tunneling between the adjacent layers by considering the hopping
between the two selenium atoms on the adjacent layers. Combining these
considerations, we can construct the Hamiltonian for the (untwisted) double
bilayer WSe2. Finally, we diagonalize the matrix of size 88 by 88 ($11$
orbitals $\times$ 2 spins $\times$ 4 layers) at each momentum to obtain the
band structure $E_{r}(k)$, where $r=\\{\text{MM},\text{XM},\text{MX}\\}$
represents the high-symmetry stacking configuration, and $k$ is the momentum
in the entire monolayer BZ. For the full recipe and the numerical details, we
refer to Ref. 34.

Our goal is to show that tdbWSe2 can satisfy criteria (1) and (3) so that a
honeycomb lattice with both sublattice and SU(2) spin rotation symmetries
emerges at the moiré scale. That is to say, we require
$E_{r}(\Gamma)>E_{r}(K)$ for $r=$MM and MX/XM, and
$E_{\text{MX}}(\Gamma)=E_{\text{XM}}(\Gamma)>E_{\text{MM}}(\Gamma)$. Obtaining
these energies from the model above requires us to make a choice for the
interlayer distances in the tdbWSe2 structure, which are unknown parameters.
Our main goal is to demonstrate the possibility and the tunability of the
emergent honeycomb lattice from tdbWSe2 and estimate the relevant energy
scale. We will not attempt to calculate the relaxed interlayer distances.
Instead, we will make assumptions on them based on the following reasoning. In
light of the $C_{2y}$ symmetry which flips the layers, the distance $d_{12}$
between the first and second layer should be the same as the distance $d_{34}$
between the third and fourth layer. For the distance between the central two
layers $d_{23}$, because the MM stacking corresponds to the case where all
tungsten atoms (also for selenium atoms) on both layers are stacked with
perfect alignment, respectively, the repulsion from the interlayer interaction
should be the strongest. Therefore, the value of $d_{23}$ in the MM stacking
configuration should be larger than $d_{12}$ and $d_{23}$ in the MX/XM
stacking configuration. With these considerations, we perturb the interlayer
distances $d_{12}$, $d_{23}$ for the MM stacking, and $d_{23}$ for the MX/XM
stacking, around the interlayer distance from the bulk WSe2 data, $\sim$ 6.5 Å
[35], and find that the system can satisfy the energetics hierarchy in
criteria (1) and (3) in a large parameter space around $d_{12}$=6 Å. As an
example that passes criteria (1) and (3), we choose the following interlayer
distances $(d_{12},d_{23},d_{34})$=(6,6.25,6) Å for MM stacking and
$(d_{12},d_{23},d_{34})$=(6,5.8,6)Å for MX and XM stacking for the estimation
of the parameters in the continuum model. Note that this choice has smaller
layer distances than the bulk data. In principle, one can consider adding
pressure to reduce the interlayer distance if needed.

Now, we present the valance band structure for the untwisted MM stacking in
Fig. 4(a), and untwisted MX/XM stacking in Fig. 4(b). The corresponding
valence band edge energies at the $\Gamma$ and $K$ valley for these two
stackings are listed in Table 1, which validate the two criteria.

To model the twisted moiré system from the untwisted system, we can represent
the energy at the high-symmetry-stacking sites in one moiré unit cell using
the energies in an untwisted system with the same type of stacking
configuration. Then we interpolate them using up to the first harmonics of the
moiré periodicity [36]. Therefore, by extracting the 4 topmost valence bands
at the $\Gamma$-valley for the untwisted systems in the MM and MX/XM stacking
configurations, we can estimate the unknown parameters in the moiré continuum
Hamiltonian Eq. (1) in the twisted system. To do that, we first will drop the
kinetic term in Eq. (1), as it is near $\Gamma$ valley, and then diagonalize
the potential term with $r=$MM and $r=$MX/XM, respectively, which gives four
discrete energies for each site. The fitting procedure is to adjust these four
energies at $r=$MM or $r=$MX/XM to match the energies obtained from the four
valence bands in the _ab initio_ tight-binding Hamiltonian at $\Gamma$ valley,
i.e., the first two columns in Table 1.

Figure 4: The valence band of tdbWSe2 in the monolayer Brillouin zone for (a) MM stacking configuration and (b) MX/XM stacking configuration Table 1: First four valence band (from the top to bottom) energies (eV) at $\Gamma$ and $K$-valley for MM and MX/XM stackings  Band | $E_{\text{MM}}(\Gamma)$ | $E_{\text{MX/XM}}(\Gamma)$ | $E_{\text{MM}}(K)$ | $E_{\text{MX/XM}}(K)$
---|---|---|---|---
0 | 0.341112 | 0.410080 | 0.334569 | 0.323034
1 | 0.244373 | 0.242830 | 0.286263 | 0.289225
2 | -0.050321 | 0.042969 | -0.182268 | -0.179475
3 | -0.553378 | -0.564056 | -0.232956 | -0.232795

##  Appendix B Phase competition in weak interaction limit

In this appendix, we present the analytical mean-field study of the Fermi
surface instabilities in our extended Hubbard model (with the in-plane
magnetic field) in the weak-coupling limit $U_{0,1}\ll t_{1}$ at hole filling
$\nu=2$. We find that, at zero temperature, dominant instabilities include
both the canted Néel AF order and the canted $\sqrt{3}\times\sqrt{3}$ SDW
order. Their competition is controlled by the ratio $U_{1}/U_{0}$. We also
discuss the relationship between our weak-coupling analysis and the Hartree-
Fock phase diagrams in Fig. 3 obtained in the stronger coupling regime
$U_{0}\gtrsim t_{1}$.

We start with the noninteracting part of the Hamiltonian
$H^{0}=H_{\text{h}}+H_{\text{Z}}$ containing the hopping term and the Zeeman
energy. At hole filling $\nu=2$, for small $t_{2}$, $h_{x}$ compared to
$t_{1}$, the Fermi surfaces are centered at the $\kappa$ and $\kappa^{\prime}$
points of the mBZ. We expand the dispersion around the
$\kappa$/$\kappa^{\prime}$ point:

$\begin{split}H^{0}_{\kappa,\sigma}&=\sum_{\bm{k}}\vec{c}_{\kappa,\sigma}^{\,\dagger}(\bm{k})\left(h_{x}\sigma\Lambda^{0}+v_{F}\vec{\Lambda}\cdot\bm{k}\right)\vec{c}_{\kappa,\sigma}(\bm{k}),\\\
H^{0}_{\kappa^{\prime},\sigma}&=\sum_{\bm{k}}\vec{c}_{\kappa^{\prime},\sigma}^{\,\dagger}(\bm{k})\left(h_{x}\sigma\Lambda^{0}-v_{F}\vec{\Lambda}^{*}\cdot\bm{k}\right)\vec{c}_{\kappa^{\prime},\sigma}(\bm{k}).\end{split}$
(3)

Here, ${\bm{k}}$ is the momentum measured from the $\kappa$/$\kappa^{\prime}$
point. The dispersion is expanded to the linear order of ${\bm{k}}$, which
suffices to study the Fermi surface instabilities. The Fermi velocity
$v_{F}=\sqrt{3}/{2}\absolutevalue{t_{1}}a_{M}$, which is $\sim 2.2\times
10^{4}$m/s at a twist angle of $\theta=2^{\circ}$.
$\vec{c}_{\tau,\sigma}=\left[c_{\tau,A,\sigma},c_{\tau,B,\sigma}\right]^{\intercal}$
is the fermionic spinor in the sublattice space at valley
$\tau=\kappa/\kappa^{\prime}$ with spin $\sigma=\uparrow,\downarrow$ along the
in-plane field direction $\hat{x}$. $\vec{\Lambda}$ comprises the Pauli-X and
Y matrices acting on the sublattice space. Note that $t_{2}$’s contribution
only starts to appear in the order-${\bm{k}}^{2}$ terms of the dispersion.
This contribution is independent of valley, spin, and sublattice. Hence, a
small $t_{2}$ is unimportant for the Fermi surface instability in the weak-
coupling limit.

At non-zero $h_{x}$, the Fermi surfaces of the spin-up fermions are perfectly
nested with those of the spin-down fermions. Consequently, particle-hole
instabilities towards magnetic ground states with spontaneous hybridization
between the two spin species are expected. The magnetic orders in these ground
states spontaneously break the U(1) spin rotation symmetry in the $yz$ plane
(perpendicular to the in-plane field direction). In the limit
$h_{x}\rightarrow 0$, the instabilities towards such magnetically-ordered
ground states vanish along with the Fermi surfaces in the noninteracting band
structure. Hence, these magnetic orders, to be discussed in detail below, are
dubbed in-plane-field-induced magnetic orders.

At hole filling $\nu=2$ (and with finite $h_{x}>0$), we begin our analysis of
possible instabilities by writing down the noninteracting bands crossing the
Fermi level:

$\begin{split}H_{\tau}=\sum_{\bm{k}}E_{+,\downarrow}(\bm{k})c_{\tau,+,\downarrow}^{\dagger}(\bm{k})c_{\tau,+,\downarrow}(\bm{k})\\\
+E_{-,\uparrow}(\bm{k})c_{\tau,-,\downarrow}^{\dagger}(\bm{k})c_{\tau,-,\downarrow}(\bm{k}),\end{split}$
(4)

for the two valleys $\tau=\kappa,\kappa^{\prime}$ and the two spin species
$\sigma=\uparrow,\downarrow$. The energies $E_{\pm,\sigma}(\bm{k})$ are given
by $E_{\pm,\sigma}(\bm{k})=h_{x}\sigma\pm v_{F}\absolutevalue{\bm{k}}$ with
$\pm$ playing the role of a band index. This band index $\pm$ is locked with
the spin index $\sigma$ for the bands that cross the Fermi level. The fermion
operators $c_{\tau,\pm,\sigma}$ are defined as

$\begin{pmatrix}c_{\kappa,+,\sigma}^{\dagger}(\bm{k})\\\
c_{\kappa,-,\sigma}^{\dagger}(\bm{k})\end{pmatrix}=\frac{1}{\sqrt{2}}\begin{pmatrix}1&e^{i\theta(\bm{k})}\\\
1&-e^{i\theta(\bm{k})}\\\
\end{pmatrix}\begin{pmatrix}c_{\kappa,A,\sigma}^{\dagger}(\bm{k})\\\
c_{\kappa,B,\sigma}^{\dagger}(\bm{k})\end{pmatrix},$ (5)

and

$\begin{pmatrix}c_{\kappa^{\prime},+,\sigma}^{\dagger}(\bm{k})\\\
c_{\kappa^{\prime},-,\sigma}^{\dagger}(\bm{k})\end{pmatrix}=\frac{1}{\sqrt{2}}\begin{pmatrix}-e^{i\theta(\bm{k})}&1\\\
e^{i\theta(\bm{k})}&1\\\
\end{pmatrix}\begin{pmatrix}c_{\kappa^{\prime},A,\sigma}^{\dagger}(\bm{k})\\\
c_{\kappa^{\prime},B,\sigma}^{\dagger}(\bm{k})\end{pmatrix},$ (6)

where $\theta(\bm{k})$ is the azimuth angle of $\bm{k}$. The spin-up (spin-
down) fermions have two particle-like (hole-like) Fermi surface pockets, one
around the $\kappa$ point and the other around the $\kappa^{\prime}$ point.
Hence, possible particle-hole instabilities (at weak couplings) can occur
through either intervalley or intravalley hybridizations of the two spin
species.

First, we consider the intervalley hybridization associated with the mean-
field order parameter
$\expectationvalue{c_{\tau,-,\uparrow}^{\dagger}(\bm{k})c_{\bar{\tau},+,\downarrow}(\bm{k})}$.
Here, $\bar{\tau}$ represents the opposite valley of $\tau$. The mean-field
decomposition of the extended Hubbard interaction corresponding to the
intervalley hybridization reads

$\begin{split}H_{\text{int}}^{\text{inter}}&=-\sum_{\begin{subarray}{c}\bm{k},\bm{k}^{\prime},\tau,\sigma\end{subarray}}\left(\frac{U_{0}}{2N}\cos(\theta(\bm{k})-\theta(\bm{k}^{\prime}))+\frac{3U_{1}}{2N}\right)\expectationvalue{c_{\bar{\tau},\bar{\sigma}}^{\dagger}(\bm{k}^{\prime})c_{{\tau},{\sigma}}(\bm{k}^{\prime})}c_{{\tau},{\sigma}}^{\dagger}(\bm{k})c_{\bar{\tau},\bar{\sigma}}(\bm{k})\\\
&+\sum_{\begin{subarray}{c}\bm{k},\bm{k}^{\prime},\tau\end{subarray}}\left(\frac{U_{0}}{2N}\cos(\theta(\bm{k})-\theta(\bm{k}^{\prime}))+\frac{3U_{1}}{2N}\right)\expectationvalue{c_{{\tau},\uparrow}^{\dagger}(\bm{k}^{\prime})c_{\bar{\tau},{\downarrow}}(\bm{k}^{\prime})}\expectationvalue{c_{\bar{\tau},\downarrow}^{\dagger}(\bm{k})c_{{\tau},{\uparrow}}(\bm{k})}.\end{split}$
(7)

Here, $\bar{\sigma}$ represents the opposite spin of $\sigma$. $N$ is the
system size. The band index is suppressed in the shorthands
$c_{{\tau},{\uparrow}}(\bm{k})=c_{{\tau},-,{\uparrow}}(\bm{k})$ and
$c_{{\tau},{\downarrow}}(\bm{k})=c_{{\tau},+,{\downarrow}}(\bm{k})$ due to the
locking of the spin and the band indices for the bands that cross the Fermi
level. In the derivation of $H_{\text{int}}^{\text{inter}}$, the terms
irrelevant to the intervalley hybridization of opposite spin species and the
terms of linear or higher order in $k$ are dropped. Note that the first
(second) term in the quadratic terms in the first line of Eq. (7) only couples
to the $p$-wave ($s$-wave) channel of the order parameter
$\expectationvalue{c_{\tau,-,\uparrow}^{\dagger}(\bm{k})c_{\bar{\tau},+,\downarrow}(\bm{k})}$.
We can consider the $s$-wave channel (with orbital angular momentum
$\absolutevalue{l}=0$) and the $p$-wave channel (with orbital angular momentum
$\absolutevalue{l}=1$) separately. At the mean-field level, the saddle point
equations in the $s$-wave channel only depend on $U_{1}$, while those in the
$p$-wave channel only depend on $U_{0}$.

For the intravalley hybridization associated with the mean-field order
parameter
$\expectationvalue{c_{\tau,-,\uparrow}^{\dagger}({\bm{k}})c_{\tau,+,\downarrow}({\bm{k}})}$,
the mean-field decomposition of interaction terms reads

$\begin{split}H_{\text{int}}^{\text{intra}}&=-\frac{U_{0}}{2N}\sum_{\begin{subarray}{c}\bm{k},\bm{k}^{\prime},\tau,\sigma\end{subarray}}\left(\expectationvalue{c_{{\tau},\bar{\sigma}}^{\dagger}(\bm{k}^{\prime})c_{{\tau},{\sigma}}(\bm{k}^{\prime})}-\expectationvalue{c_{\bar{\tau},\bar{\sigma}}^{\dagger}(\bm{k}^{\prime})c_{\bar{\tau},{\sigma}}(\bm{k}^{\prime})}\right)c_{{\tau},{\sigma}}^{\dagger}(\bm{k})c_{{\tau},\bar{\sigma}}(\bm{k})\\\
&-\frac{3U_{1}}{2N}\sum_{\begin{subarray}{c}\bm{k},\bm{k}^{\prime},\tau,\sigma\end{subarray}}\cos(\theta(\bm{k})-\theta(\bm{k}^{\prime}))\expectationvalue{c_{{\tau},\bar{\sigma}}^{\dagger}(\bm{k}^{\prime})c_{{\tau},{\sigma}}(\bm{k}^{\prime})}c_{\tau,\sigma}^{\dagger}(\bm{k})c_{{\tau},\bar{\sigma}}(\bm{k})\\\
&+\frac{U_{0}}{2N}\sum_{\begin{subarray}{c}\bm{k},\bm{k}^{\prime}\end{subarray}}\left(\expectationvalue{c_{\kappa,\uparrow}^{\dagger}(\bm{k}^{\prime})c_{\kappa,\downarrow}(\bm{k}^{\prime})}-\expectationvalue{c_{\kappa^{\prime},\uparrow}^{\dagger}(\bm{k}^{\prime})c_{\kappa^{\prime},\downarrow}(\bm{k}^{\prime})}\right)\left(\expectationvalue{c_{\kappa,\downarrow}^{\dagger}(\bm{k})c_{\kappa,\uparrow}(\bm{k})}-\expectationvalue{c_{\kappa^{\prime},\downarrow}^{\dagger}(\bm{k})c_{\kappa^{\prime},\uparrow}(\bm{k})}\right)\\\
&+\frac{3U_{1}}{2N}\sum_{\begin{subarray}{c}\bm{k},\bm{k}^{\prime},\tau\end{subarray}}\cos(\theta(\bm{k})-\theta(\bm{k}^{\prime}))\expectationvalue{c_{{\tau},{\uparrow}}^{\dagger}(\bm{k}^{\prime})c_{{\tau},{\downarrow}}(\bm{k}^{\prime})}\expectationvalue{c_{{\tau},{\downarrow}}^{\dagger}(\bm{k})c_{{\tau},{\uparrow}}(\bm{k})}.\end{split}$
(8)

Similar to the intervalley case, based on the form of
$H_{\text{int}}^{\text{intra}}$, we only need to consider the $s$-wave and the
$p$-wave channel of the order parameter
$\expectationvalue{c_{\tau,-,\uparrow}^{\dagger}({\bm{k}})c_{\tau,+,\downarrow}({\bm{k}})}$.
Also, it is natural to expect that the order parameters in the two valleys
$\tau=\kappa,\kappa^{\prime}$ share the same orbital angular momentum. With
this expectation, we find that the saddle point equations in the $s$-wave
channel only depend on $U_{0}$ while those in the $p$-wave channel only depend
on $U_{1}$.

At zero temperature, we find the two most dominant mean-field solutions with
the largest gaps opened around the Fermi surface. The corresponding orders are
then expected to have the highest mean-field transition temperatures. One of
these mean-field solutions is given by the intervalley hybridization in the
$s$-wave channel with
$\expectationvalue{c_{\kappa,-,\uparrow}^{\dagger}c_{\kappa^{\prime},+,\downarrow}}=\expectationvalue{c_{\kappa^{\prime},-,\uparrow}^{\dagger}c_{\kappa,+,\downarrow}}$,
which gives rise to the canted $\sqrt{3}\times\sqrt{3}$ SDW order, as
illustrated in Fig. 3(b). In this phase, the gap opened at the Fermi surface
is given by $2h_{x}\exp(-\frac{2}{3U_{1}\rho})$, where the density of states
$\rho$ can be approximated to a constant $\frac{4h_{x}}{3\pi t_{1}^{2}}$
within the energy interval between two Dirac points (at energies $-h_{x}$ and
$h_{x}$ relative to the Fermi level). The other dominant mean-field solution
is given by the intravalley hybridization in the $s$-wave channel with
$\expectationvalue{c_{\kappa,-,\uparrow}^{\dagger}c_{\kappa,+,\downarrow}}=-\expectationvalue{c_{\kappa^{\prime},-,\uparrow}^{\dagger}c_{\kappa^{\prime},+,\downarrow}}$,
which gives rise to the canted Néel AF order, as shown in Fig. 3 (b) and (c).
In this canted Néel AF phase, the gap size is given by
$2h_{x}\exp(-\frac{1}{U_{0}\rho})$.

By comparing the gaps, we expect the canted Néel AF phase to be the
energetically favorable ground state for $U_{1}<\frac{2}{3}U_{0}$ in the weak
coupling limit. When $U_{1}$ increases beyond $\frac{2}{3}U_{0}$, the ground
state goes through a first-order phase transition into the canted
$\sqrt{3}\times\sqrt{3}$ SDW phase. We comment that it is challenging to
extend our numerical Hartree-Fock calculation, which includes all the
interaction terms, to the weak-coupling limit because of the exponentially
small gap sizes. From the Hartree-Fock phase diagram with $U_{0}>t_{1}$ as
shown in Fig. 3 (c), the sublattice-polarized phase is more energetically
favorable than both the canted $\sqrt{3}\times\sqrt{3}$ SDW insulator and the
canted Néel AF insulator when $U_{1}\gtrsim\frac{1}{3}U_{0}+\frac{2}{3}h_{x}$.
However, there is no Fermi surface instability towards the sublattice-
polarized insulator in the weak coupling limit. It is interesting to
investigate the competition among the canted Néel AF insulator, the
$\sqrt{3}\times\sqrt{3}$ SDW insulator, and the sublattice-polarized insulator
in the regime with $U_{0}<t_{1}$, which we will leave for future study.