# Second-order and real Chern topological insulator in twisted bilayer
$\alpha$-graphyne

Bin-Bin Liu  School of Physics, Beihang University, Beijing 100191, China  Xu-
Tao Zeng  School of Physics, Beihang University, Beijing 100191, China  Cong
Chen  School of Physics, Beihang University, Beijing 100191, China  Department
of Physics, The University of Hong Kong, Hong Kong, China  Ziyu Chen
<EMAIL_ADDRESS>School of Physics, Beihang University, Beijing 100191,
China  Xian-Lei Sheng<EMAIL_ADDRESS>School of Physics, Beihang
University, Beijing 100191, China  Peng Huanwu Collaborative Center for
Research and Education, Beihang University, Beijing 100191, China

###### Abstract

The study of higher-order and real topological states as well as the material
realization have become a research forefront of topological condensed matter
physics in recent years. Twisted bilayer graphene (tbG) is proved to have
higher-order and real topology. However whether this conclusion can be
extended to other two-dimensional twisted bilayer carbon materials and the
mechanism behind it lack explorations. In this paper, we identify the twisted
bilayer $\alpha$-graphyne (tbGPY) at large twisting angle as a real Chern
insulator (also known as Stiefel-Whitney insulator) and a second-order
topological insulator. Our first-principles calculations suggest that the
tbGPY at 21.78∘ is stable at 100 K with a larger bulk gap than the tbG. The
non-trivial topological indicators, including the real Chern number and a
fractional charge, and the localized in-gap corner states are demonstrated
from first-principles and tight-binding calculations. Moreover, with
$\mathcal{C}_{6z}$ symmetry, we prove the equivalence between the two
indicators, and explain the existence of the corner states. To decipher the
real and higher-order topology inherited from the Moiré heterostructure, we
construct an effective four band tight-binding model capturing the topology
and dispersion of the tbGPY at large twisting angle. A topological phase
transition to a trivial insulator is demonstrated by breaking the
$\mathcal{C}_{2y}$ symmetry of the effective model, which gives insights on
the trivialization of the tbGPY as reducing the twisting angle to 9.43∘
suggested by our first-principles calculations.

##  I INTRODUCTION

The discovery of the topological insulator (TI) has stimulated a vibrant
research field in condensed matter physics Hasan and Kane (2010); Qi and Zhang
(2011); Shen (2012); Bansil  _et al._ (2016). A TI in $d$ dimensions has an
insulating bulk, but features topologically protected gapless states on its
$(d-1)$-dimensional (D) boundaries. Recently, the concept was extended to a
novel class of topological phase-the higher-order TI Zhang  _et al._ (2013);
Benalcazar  _et al._ (2017a); Langbehn  _et al._ (2017); Song  _et al._
(2017); Benalcazar  _et al._ (2017b); Schindler  _et al._ (2018a). An $n$th-
order TI has topological gapless states at its $(d-n)$D boundaries, but is
gapped otherwise. For instance, a second-order topological insulator (SOTI) in
2D hosts topological gapless states at its 0D corners between its edges that
are gapped. Previous works first reveal the higher-order TIs in 3D materials
Bradlyn  _et al._ (2017); Schindler  _et al._ (2018b, a); Yue  _et al._
(2019); Wang  _et al._ (2019); Zhang  _et al._ (2019); Vergniory  _et al._
(2019); Tang  _et al._ (2019a, b); Xu  _et al._ (2019). Then in 2D there are a
few material candidates being proposed, such as the graphdiyne family Sheng
_et al._ (2019); Liu  _et al._ (2019); Lee  _et al._ (2020); Chen  _et al._
(2021), twisted bilayer graphene (tbG) Park  _et al._ (2019); Liu  _et al._
(2021), black phosphorene Ezawa (2018); Hitomi  _et al._ (2021), Bi/EuO Chen
_et al._ (2020), monolayer group-V Radha and Lambrecht (2020); Huang and Liu
(2021), and group-IV materials Qian  _et al._ (2021); Pan  _et al._ (2022). It
is still a big challenge to find more realistic 2D SOTI materials.

With the spacetime inversion symmetry, the wavefunction over the Brillouin
zone (BZ) is real instead of complex. The so-called real topology is
characterized by the real Berry bundles over the BZ and can be indicated by
the Real Chern number (RCN)Zhao and Lu (2017) or the second Stiefel-Whitney
number Nakahara (2003). Twisting one layer of the AA-stacking periodic lattice
structures with translational symmetry, one can obtain a larger super cell
which forms the so-called Moiré pattern. In the field of twistronics, the tbG
is a pioneer material that has been investigated in many aspects including
topology. The tbG Moiré system is shown to be topological Song  _et al._
(2019) and features a nontrivial RCN or the second Stiefel-Whitney number Ahn
_et al._ (2019); Park  _et al._ (2019). An insulator with nontrivial RCN (or
second Stiefel-Whitney number) like tbG can be referred as a real Chern
insulator (RCI) Wang  _et al._ (2020) (or Stiefel-Whitney insulator Lee  _et
al._ (2020)). Besides, the tbG system at relatively large Moiré twisting
angle, say 21.78∘, is further proved to be a SOTI by explicitly demonstrating
the existence of the corner states Park  _et al._ (2019).

By inserting the acetylenic linkage into graphene lattice, many carbon
allotropes can be constructedBaughman  _et al._ (1987); Malko  _et al._
(2012). Except for graphene, other 2D carbon allotropes with twisted bilayer
structures are lacking of exploration, such as
($\alpha,~{}\beta,~{}\gamma$)-graphyne Malko  _et al._ (2012) and graphdiyne
Li  _et al._ (2010); Chen  _et al._ (2017); Jia  _et al._ (2017); Gao  _et
al._ (2019). These materials with hexagonal lattices all have the same space
group symmetry as graphene, but their properties are not identical to graphene
and vary in their own way due to the presence of the acetylenic linkageJana
_et al._ (2019). So their twisted bilayer structures do not necessarily behave
like the tbG.

Among these graphynes, the $\alpha$-graphyne is the most similar one to
graphene. As studied in Refs. Malko  _et al._ (2012); Kim and Choi (2012);
Bandyopadhyay and Jana (2020); Jana  _et al._ (2021), $\alpha$-graphyne has
eight carbon atoms per unit cell [see Fig. 1(a)], which can be regarded as
inserting an acetylenic linkage between every two atoms in graphene. The
symmetry group for $\alpha$-graphyne (and graphene) is the No. 191 space group
with point group D6h ($p$6m symmetry). The $\alpha$-graphyne shares some
properties with graphene. For example, the band structure of the
$\alpha$-graphyne features a similar linear gapless cone at the K in the BZ as
shown in Fig. 1(c). Therefore, it is natural to ask if the twisted bilayer
$\alpha$-graphyne (tbGPY for short) can inherit the second-order and real
topology from the Moiré heterostructure like the tbG at large twisting angle.

Figure 1: (a) Crystal structure and unit cell of the $\alpha$-graphyne with a
zigzag edge. (b) shows the high symmetry points in the first Brillouin zone
used in our calculations. (c) and (d) show the bulk and edge band structures
for the monolayer $\alpha$-graphyne.

In this paper, we address the issue by exploring the electronic and
topological properties of the tbGPY. Based on tight-binding and first-
principles calculations, we identified the tbGPY at large twisting angle as a
RCI and a SOTI, which has a larger bulk gap than the tbG. We confirmed the
bulk topology from both the nontrivial RCN $\nu_{R}$ Zhao and Lu (2017); Chen
_et al._ (2021) and a fractional charge $Q^{(2)}_{corner}$ Benalcazar  _et
al._ (2019) that indicates a filling anomaly. With the $\mathcal{C}_{6z}$
symmetry, we strictly proved that the nontrivial RCN equals an e/2 fractional
charge $Q^{(2)}_{corner}$, which in turn suggests the existence of corner
states on a hexagonal lattice. We further calculated the energy level of a
hexagonal nano-disk for the tbGPY at 21.78∘ and found six in-gap localized
corner states, which exemplifies the bulk-boundary correspondence of a RCI
with $\mathcal{C}_{6z}$ symmetry. The robustness of the corner states against
symmetry breaking disorders was demonstrated. We also performed an ab initio
molecular dynamics calculation which suggests that the tbGPY at 21.78∘ is
stable at 100 K.

Moreover, we built an effective tight-binding model to decipher the SOTI and
RCI state inherited from the Moiré heterostructure of tbGPY. The model
features a RCI as well as a SOTI with non-trivial fractional charge, and can
be transited to a trivial insulator by breaking $\mathcal{C}_{2y}$ symmetry.
To extend, the model can be applied to other twisted bilayer materials with a
Dirac cone at the K point in the monolayer, such as the tbG. Different from
the tbG, we observed that the RCI state in tbGPY becomes trivial as the
twisting angle becomes smaller. This trivialization can be simulated by our
effective model through the topological phase transition.

##  II Lattice and band structures

We begin by introducing the general geometry for the twisted bilayer Moiré
system with hexagonal lattice. To obtain the Moiré pattern, one may stack two
periodic lattice planes together (i.e. AA-stacking), and then twist one layer
with respect to the other with respect of certain commensurate angle. For the
hexagonal lattice unit cell, the commensurate angle can be formulated as
$\theta_{i}=\arccos(\frac{3i^{2}+3i+0.5}{3i^{2}+3i+1})$ (where $i$ is an
integer above zero) Lopes dos Santos  _et al._ (2007). The commensurate
condition does not depend on the inversion center, and is chosen as the middle
of the hexagon lattice [Fig. 1 (a) ] in this work. The Moiré pattern and the
associated lattice vectors only depend on the twisting angle or the index $i$.
In the following, we use the $i$ to label the commensurate angle.

Figure 2: The tbGPY ($i=1$) material and its properties. (a) Crystal structure
with lattice vectors and a specified armchair edge. (b) and (d) show the
gapped bulk and edge spectra with a gap around 13 meV for the bulk and 65 meV
for the edge. (c) The Wilson loop spectrum around $\theta=\pi$, which indicats
the non-trivial RCN ($\nu_{R}=1$). (e) The energy levels of the hexagonal
nano-disk in (f) near the fermi-level. (f) The real space distribution of the
six in-gap states marked as green dots in (e), suggesting the localized corner
states.

The lattice structure of the twisted bilayer $\alpha$-graphyne at the
commensurate angle 21.78∘ with $i=1$ is shown in Fig. 2(a), which is the main
focus of the paper. To test the dynamical stability, we performed a molecular
dynamics simulation for the tbGPY and find that it is dynamically stable under
at least 100 K. The simulation result is summarized in Appendix A. Its Moiré
supercell contains 112 carbon atoms, seven times as large as the AA-stacking
$\alpha$-graphyne. We obtained the layer distance of the tbGPY with a first-
principles relaxation calculation, which is $d_{0}=$3.444 Å, in agreement with
previous studiesLeenaerts  _et al._ (2013). See Appendix B.1 for calculation
details. Different from the monolayer $\alpha$-graphyne, the tbGPY lies in the
No. 177 space group with D6 point group symmetry. The generators for this
associated point group can be chosen as $\mathcal{C}_{6z}$ and
$\mathcal{C}_{2y}$.

We calculated the energy dispersion of the tbGPY with $i=1$ based on density
functional theory. The result is shown in the dashed lines in Fig. 2(b), which
is calculated along the line connecting the high symmetry points specified in
Fig. 1 (b). We found a sizable gap of around 12.8 meV. To compare, we also
calculate the energy spectrum for the twisted bilayer graphene with layer
distance 3.35 Å using the same calculation method, and find a gap of around
1.4 meV at K. See Appendix B.1 for explanations of the method. The opening of
the band gap for the tbGPY can be understood from the $U_{v}(1)$ valley
symmetry breaking for the bulk cones.

To proceed, we constructed a symmetry-based tight binding model with the
method introduced by Slater and KosterSlater and Koster (1954). The details
for the method are shown in Appendix B.2. Following the method, we constructed
a Slater Koster tight-binding (SKTB) model for tbGPY with one $p_{z}$ orbital
per site. The fitted energy dispersion for the tbGPY using the SKTB model is
shown in Fig. 2(b) with solid lines, which captures the main feature of the
energy dispersion of tbGPY in the low energy range.

The parameters for the SKTB model are illustrated as follows. The SKTB model
includes all the couplings that have bond lengths less than the range of $L=5$
Å. The nearest neighbor distance in the monolayer $\alpha-$graphyne is
$a_{0}=1.232$ Å, and the layer distance of the tbGPY is $d_{0}=$3.444 Å. The
nearest neighbor p-p $\pi$ coupling is fitted as $V_{pp\pi}^{0}=-4.45$ eV
while the nearest p-p $\sigma$ coupling is $V_{pp\pi}^{0}=0.32$ eV. Other
orbital couplings are determined relatively to those nearest ones by their
relative distance to $a_{0}$ and $d_{0}$ [See Eq. (B3) and Eq. (B4)]. The
decay length of the transfer integral is chosen as $\delta_{0}$=0.32 Å. The
onsite energy for each orbitals is fitted as -0.156 eV. In the following, we
investigate the topological properties of the material.

##  III Real Chern number

For a spinless system, under the spacetime inversion symmetry
$\mathcal{P}\mathcal{T}$, the wavefunction over the BZ is real instead of
complex. In such case, the topology is characterized by the real Berry bundles
over the BZ and can be termed as the real topology. The topological indicator
for the real topology is the the RCNZhao and Lu (2017), or the second Stiefel-
Whitney number Nakahara (2003). The RCN is the defining topological invariant
for the RCI as well as other exotic topological phase such as the second-order
nodal-line semimetalWang  _et al._ (2020); Chen  _et al._ (2022). In 2D the
defining spacetime inversion symmetry $\mathcal{P}\mathcal{T}$ can be replaced
with the $\mathcal{C}_{2z}\mathcal{T}$ symmetry. Originated from the $D_{6}$
point group, the tbGPY has the $\mathcal{C}_{2z}$ symmetry, while also
possessing the time reversal symmetry $\mathcal{T}$. So it has the combined
symmetry $\mathcal{C}_{2z}\mathcal{T}$ required for the real topological
classification.

In practice, there are two ways to compute the RCN Chen  _et al._ (2021). One
is the well-known Wilson loop method: calculating the Wilson loop along a
chosen direction e.g., $\mathbf{k}_{1}$, with fixed the $k_{2}\in[-\pi,\pi)$.
This results in N (the number of valence bands of the system) curves in the
$k_{2}-\theta$ diagram, which represents the evolution of the Wannier center.
The parity of the winding number or equivalently the RCN ($\nu_{R}$) can be
read off from the times of crossings ($\zeta$) of the Wilson loop spectrum
with the $\theta=\pi$ lineAhn  _et al._ (2018), i.e.,

$\displaystyle\nu_{R}=\zeta~{}\text{mod~{}2}.$ (1)

When not only the $\mathcal{PT}$ but also the $\mathcal{P}$ symmetry is
preserved, there is another more intuitive method to calculate the RCN, which
is by counting the parity eigenvalues of the valence bands at the four time
reversal invariant momenta (TRIM) points $\Gamma_{i}$. Practically, one may
follow

$\displaystyle(-1)^{\nu_{R}}=\prod_{i=1}^{4}(-1)^{\lfloor(n_{-}^{\Gamma_{i}}/2)\rfloor},$
(2)

to compute the RCN $\nu_{R}$, where the $n_{-}^{\Gamma_{i}}$ is the number of
the minus parity in the valence band at $\Gamma_{i}$.

We calculated the RCN for tbGPY with both methods. With the fist method, we
obtained the Wilson loop from the SKTB model of tbGPY as in Fig. 2(c), where a
single crossing of the spectrum with $\pi$ suggests the RCN to be nontrivial,
i.e., $\nu_{R}=1$. With the second parity counting method, we calculate the
parity eigenvalues of the $\mathcal{C}_{2z}$ at TRIM points from first-
principles calculations. Using Eq. (2), the RCN is also found to be
nontrivial. The parity counting method can be further simplified in our case.
As the tbGPY lattice possesses the $\mathcal{C}_{6z}$ symmetry, three of the
TRIM points (M points) are equivalent, so that we only need to calculate one M
point instead all three of them. Detailed information about the RCN
calculation is shown in Tab. D1 in Appendix D. Our results suggests that the
tbGPY system is a RCI at large commensurate angle.

The second parity counting method is more physically intuitive, because it
tells us information about band inversion. For example, a nontrivial RCN in
Tab. D1 corresponds to a double band inversion at the $\Gamma$ point.

A nontrivial RCN in a $\mathcal{P}\mathcal{T}$ invariant system could feature
a rich boundary correspondenceWang  _et al._ (2020). This means that despite
having a nontrivial RCN, the tbGPY is not necessarily a SOTI. In the
following, we explore the bulk-boundary correspondence for the tbGPY system
and prove that it is a SOTI with nontrivial fractional charge that features
localized corner states in a hexagonal lattice.

##  IV Bulk-boundary correspondence

In 2D, a TI has topologically protected gapless edge states, while for a SOTI
the edge state is gapped but the conducting corner charges in the edge gap
emerge. To proceed, we firstly calculate the edge spectrum for the tbGPY
($i=1$) and compare it with that of monolayer graphyne. In Fig. 2(d) and Fig.
1(d), we show the edge spectra for the tbGPY and the monolayer
$\alpha$-graphyne with the edges illustrated in Fig. 2(a) and Fig. 1(a).
Comparing the two, we find that instead of trivially stacking, the two bulk
Dirac cones in tbGPY are shifted to higher and lower energies with a gap in
between. Correspondingly, the edge states (bright yellow curves) that are
associated with the bulk Dirac cones are also shifted upward and downward,
resulting in an edge gap. The edge gap is found to be around 65 meV which is
larger than the bulk gap of 13 meV [Fig. 2(d)].

Then we explore the property of the second-order boundaries (corners) of tbGPY
on a hexagonal nano-disk geometry related by $\mathcal{C}_{6z}$ rotation
symmetry. Inside the edge gap, we find six localized corner states in the
nano-disk. Specifically, we calculate the energy level of a hexagonal nano-
disk of tbGPY with 33360 carbon atoms ($p_{z}$ orbitals), and find six states
inside the edge gap marked as green dots, as shown in Fig. 2(e). The observed
non-degeneracy of the energy levels of the six corner states is due to the
slight breaking of chiral symmetry in the tbGPY and the size effect of the
nano-disk. This is proved in our effective model analysis. We plot the charge
distribution of the six in-gap states and find that they are localized at six
corners of the hexagonal nano-disk in Fig. 2(f).

We now briefly introduce the concept of fractional charge to explain the
presence of the six corner states. The fractional charge is a bulk property
that indicates the filling anomaly of electrons and is defined under
crystalline symmetries. Under the $\mathcal{C}_{2z}$, a fractional charge
termed $Q^{(2)}_{corner}$ can be definedBenalcazar  _et al._ (2019) (see
Appendix C for detail). Further considering $\mathcal{C}_{6z}$ symmetry as in
the tbGPY case, we can relate the $Q^{(2)}_{corner}$ with the RCN in a neat
form

$\displaystyle Q^{(2)}_{corner}=e\frac{\nu_{R}}{2}.$ (3)

According to Ref. Benalcazar  _et al._ (2019), the fractional charge is a
secondary indicator which is well defined with the vanishing of the
polarization indicator $\bf P^{(2)}$. Remarkably, under the protection of
$\mathcal{C}_{6z}$ symmetry, we find

$\displaystyle\bf P^{(2)}=\bf 0.$ (4)

This is illuminating since it indicates a correspondence between the RCI and a
SOTI with localized corner states. A nontrivial RCN is a bulk property, which
in the presence of $\mathcal{C}_{6z}$ is equivalent to a nontrivial well-
defined fractional charge $Q^{(2)}_{corner}$ that corresponds to secondary
boundary states related by $\mathcal{C}_{2z}$. Here we derive the
correspondence (Eq. 3 and Eq. 4) with $\mathcal{C}_{6z}$, therefore there can
exist six localized corner states in general. We further demonstrate that the
tbGPY ($i=1$) owns nontrivial fractional charges of $Q^{(2)}_{corner}=e/2$
through first-principles. Therefore, the hexagonal nano-disk of tbGPY hosts
six corner states, as shown in Fig. 2 (f) and in Fig. 3 (f) in Sec. V.2.

The six corner charges are robust against disorders along the edge, provided
that the bulk and edge energy gaps remain open. This is because that the
disorder might shift the in-gap corner charges by only an integer, but leaves
the fractional part intact. The detailed calculations for the tbGPY nano-disk
with disorder can be found in Fig. C1 in Appendix C.

We have demonstrated the bulk-boundary correspondence of the tbGPY ($i=1$) by
showing that it is not only a RCI, but also a SOTI with localized in-gap
corner states.

##  V Effective four-band model

To decipher the higher-order and real topology inherited from the complex
Moiré supercell of tbGPY with numerous atoms, we provide a simpler effective
model with only four orbitals that captures the non-trivial topological
properties.

###  V.1 Model construction

We build the effective TB model starting from the symmetries of tbGPY. By
effective, we mean that the effective model captures the main features of the
bulk and edge dispersions, possesses the nontrivial topology, i.e., the
nontrivial RCN and fractional charge, and characterizes the localized corner
states.

To achieve the requirement, we find the symmetries for constructing the
effective model can be $\mathcal{C}_{2z}\mathcal{T}$ and crystalline
$\mathcal{C}_{6z}$ and $\mathcal{C}_{2y}$. Due to the presence of the combined
$\mathcal{C}_{2z}\mathcal{T}$ that defines the real topology, there exists a
Wannier obstruction for building a two band tight binding model Ahn  _et al._
(2019); Zou  _et al._ (2018). Following the strategy in Refs. Zou  _et al._
(2018); Song  _et al._ (2019), we build a four band tight binding (4BTB) model
in a honeycomb lattice, as shown in Fig. 3(a). The four orbitals of our 4BTB
model forms two layers with each layer two sublattices in a unit cell. Here,
we do not distinguish the orbital type, i.e., whether they are $s$ or $p_{z}$
orbitals. Without loss of generality, one may assume four $s$ orbitals with
each orbital per site. Incorporating the symmetries into the orbitals, we can
immediately specify the hoppings for the model as shown in Fig. 3(a).

Figure 3:  The topological 4BTB model and its properties. (a) Lattice,
orbitals and specified hoppings with a zigzag edge. (b) and (d) show the
gapped bulk and edge band structures. (c) The Wilson loop spectrum indicating
the nontrivial RCN, $\nu_{R}=1$. (e) The energy levels of the hexagonal nano-
disk in (f). (f) The charge distribution of the six in-gap states (marked as
green dots) in (e), indicating the presence of the localized corner states.
The parameters are chosen as $t_{t}=t_{b}=t=1/3$ eV, $\lambda=0.01t$ and
$M=0.1t$ for calculation.

Now we explain the coupling bonds for the 4BTB model. The black and red arrows
are the nearest neighbor hoppings for the top and bottom layers, which
generates two degenerated copies of massless Dirac cones at K point as that of
graphene. The amplitudes of the hoppings are denoted as $t_{t}$ (top) and
$t_{b}$ (bottom). The pink bond (denoted as $M$) connects each of the nearest
top and bottom orbitals as shown in the side view in Fig. 3(a). This bond
functions as the energy shifting term which shifts the two Dirac cones upward
and downward, which guarantees a gap for the edge spectrum. To have a bulk
band gap around K, we introduce the alternative hoppings $\lambda$ that
connect the top and bottom layer orbitals with solid orange and dashed dark
blue arrows. The alternating pattern of the $\lambda$ is that the sign of the
bond changes between positive and negative when rotating 60∘. These terms
effectively introduce a gap in the bulk, corresponding to the Moire twisting
in tbGPY. The 4BTB model in momentum space then reads

$\displaystyle
H^{\text{4BTB}}\left(\mathbf{k}\right)=\begin{bmatrix}G_{t}&M\sigma_{0}+iS\\\
M\sigma_{0}-iS&G_{b}\end{bmatrix}$ (5)

with

$\displaystyle G_{j}$
$\displaystyle=t_{j}\sum_{i=1}^{3}\left[\sigma_{x}\cos\left(\bm{\delta}_{i}\cdot\mathbf{k}\right)-\sigma_{y}\sin\left(\bm{\delta}_{i}\cdot\mathbf{k}\right)\right]~{}(j=t,b),$
$\displaystyle S$
$\displaystyle=2\lambda\sigma_{z}\sum_{i=1}^{3}\sin\left(\mathbf{d}_{i}\cdot\mathbf{k}\right).$

The $\sigma_{0}$ is a 2-by-2 identity matrix, and the $\sigma_{x,y,z}$ are
Pauli matrices acting on the sublattice degrees of freedom. The Hamiltonian in
Eq. 5 distinguishes top and bottom layers. The
$\bm{\delta}_{1}=\frac{1}{3}\mathbf{a}_{1}+\frac{2}{3}\mathbf{a}_{2}$,
$\bm{\delta}_{2}=-\frac{2}{3}\mathbf{a}_{1}-\frac{1}{3}\mathbf{a}_{2}$, and
$\bm{\delta}_{3}=\frac{1}{3}\mathbf{a}_{1}-\frac{1}{3}\mathbf{a}_{2}$ are the
nearest hopping vectors within one layer, while
$\mathbf{d}_{1}=\mathbf{a}_{1}$, $\mathbf{d}_{2}=\mathbf{a}_{2}$, and
$\mathbf{d}_{3}=-\mathbf{a}_{1}-\mathbf{a}_{2}$ are the second-nearest hopping
vectors within one layer. Note that the hopping parameters of the 4BTB model
are all real numbers due to an emergent $\mathcal{T}$ symmetry. Therefore the
4BTB model recovers the same symmetries as tbGPY.

###  V.2 Nontrivial bulk-boundary correspondence

To examine the validity of the effective 4BTB model, we check the criterion
mentioned at the beginning of the Sec.V.1. Specifically, we set
$t_{t}=t_{b}=t=1/3$ eV, $\lambda=0.01t$ and $M=0.1t$ for the model. The
parameters are obtained by simulating the main features of the bulk band of
tbGPY around K.

The bulk property concerns band dispersion as well as the topological
indicators. The band structure has a small gap around K with a valley
splitting feature, in agreement with the material calculations of tbGPY, as
shown in Fig. 3(b) and Fig. 2(b). The RCN as well as the fractional charges
$Q^{(2)}_{corner}$ are all nontrivial and agree with the tbGPY ($i=1$) as
shown in the $i=1$ row of Tab. 1. We also calculate the RCN with the Wilson
loop technique. We calculate the Wannier center evolution along the $k_{2}$
which has two curves and a single crossing point at $\theta=\pi$, which means
the RCN is nontrivial. See Fig. 3(c). Therefore, the topological 4BTB model
features a RCI with nontrivial fractional charges.

Table 1:  The RCNs and fractional charges for tbGPY material and effective 4BTB model. The first two columns indicate the systems with different $i$. The remaining columns are the RCN $\nu_{R}$ and the fractional charge of $Q^{(2)}_{corner}$ calculated by Eq. (2) and (3). The first two rows are tbGPY with $i=1$ calculated from first-principles and simulated by the topological 4BTB model. The last two rows are the tbGPY with $i=3$ calculated from first-principles and simulated by the trivial 4BTB model (which differs from the topological 4BTB model in $t_{b}=-t_{t}=-1/3$ eV).  $i$ | System | $\nu_{R}$ | $Q^{(2)}_{corner}$
---|---|---|---
$i=1$ | tbGPY | 1 | e/2
Topological 4BTB | 1 | e/2
$i=3$ | tbGPY | 0 | 0
Trivial 4BTB | 0 | 0

Then we check the boundary correspondence of the topological 4BTB model. We
investigate the edge spectrum along the zigzag direction specified in Fig.
3(a), and the result is shown in Fig. 3(d). One observes that each edge
connects two bulk cones, and leaves an edge gap in between, which has the same
feature as that in Fig. 2(d). To investigate the second-order edge, we
construct a nano-disk containing 30000 orbitals for the 4BTB model. There are
six degenerated zero modes at the Fermi level labeled by green dots, as shown
in Fig. 3(e), which are distributed at the six corners of the nano-disk as
shown in Fig. 3(f). In the model calculation, we are able to eliminate the
size effect and hence the corner modes are nearly degenerated at 0 eV. See
Appendix E for an explanation of topological origin of the corner state with
domain wall theory.

Therefore, we have constructed an effective model with only four orbitals for
the tbGPY at large commensurate angle. To extend, the 4BTB model may also be
applied to study other twisted bilayer system at large commensurate angle,
such as tbG which has the same symmetry as tbGPY and has a Dirac cone at K in
the monolayer Park  _et al._ (2019).

###  V.3 Topological phase transition

After demonstrating the nontrivial topology of the 4BTB model, it is natural
to explore the other side of the model, namely a trivial insulator. We here
provide a simple and intuitive way to trivialize the effective 4BTB model. To
achieve so, we simply reverse the sign of $t_{b}$, following that
$t_{b}=-t_{t}=-t$. In this way, we break the $\mathcal{C}_{2y}$ symmetry of
the 4BTB model and achieve a band inversion between the valence and
conductance bands. Since the $\mathcal{C}_{6z}$ symmetry is still remained in
the model, we can use Eq. 2 and Eq. 3 to calculate the RCN and the fractional
charge. The topological indicators for the trivial 4BTB model are both 0 as
shown in Tab. 1. Slight breaking of the $\mathcal{C}_{2y}$ symmetry will not
spoil the topological corner states as long as the bulk and edge gaps are not
closedSheng  _et al._ (2019). However, here the $\mathcal{C}_{2y}$ symmetry
breaking is accompanied with a band inversion, which is the reason that the
RCN and the fractional charge become trivial. As a consequence, there is no
localized corner state in the trivial 4BTB model, in contrast to that of the
topological 4BTB model as shown in Figs. 3(e, f).

For tbGPY with smaller commensurate angel, say 9.4∘ with $i=3$, the RCN
becomes trivial (see Appendix D for some details and discussions). This
indicates a band inversion. The trivial 4BTB model features the same RCN and
fractional charge as the material with $i=3$ as presented in Tab. 1, and gives
an understanding of the trivialization process from the symmetry breaking and
band inversion prospects.

##  VI Conclusion and outlook

In conclusion, we proposed the tbGPY at large commensurate angle as a RCI and
a SOTI from first-principles and tight-binding model calculations. The tbGPY
has a sizable bulk gap of around 13 meV larger than the tbG and is stable at
100 K, which is favorable for the experimental detection of the corner states.
To decipher the topological property inherited from the Moiré structure, we
constructed a simple but effective 4BTB model which captures the major
topological features of the tbGPY. The model may also apply to other twisted
bilayer systems with a Dirac cone at K in the monolayer, such as the tbG. A
phase transition of the effective model by breaking the $\mathcal{C}_{2y}$
symmetry is also studied, which gives insights on the trivialization of the
RCI state in the tbGPY. Moreover, in the presence of $\mathcal{C}_{6z}$
symmetry, we prove that the nontrivial RCN $\nu_{R}$ equals an e/2 fractional
charge $Q^{(2)}_{corner}$, which in turn explains the six localized corner
states in the tbGPY hexagonal nano-disk.

We give an outlook on the experimental realization of the RCI and SOTI in
artificial systems with the 4BTB model. In the topological 4BTB model, all the
hopping amplitudes are real, either positive or negative, which preserves the
time-reversal symmetry $\mathcal{T}$. Therefore, the model can be realized in
a rich category of artificial systems, such as electrical circuit systemsImhof
_et al._ (2018), phonon latticesMa  _et al._ (2019); Xue  _et al._ (2020,
2019), photonic crystalsOzawa  _et al._ (2019); Mittal  _et al._ (2019),
mechanical systemsProdan and Prodan (2009), and cold atomsDalibard  _et al._
(2011). Moreover, a topological phase transition can be realized by simply
reversing the sign of the nearest neighbor hoppings in the bottom/top layer
$t_{b/t}$, which can be engineered to tune the existence of the corner states.
In addition, removing the $M$ term in the 4BTB model, the remaining spinless
model with pseudo-spin generated by the layer operator $\bf\tau$ behaves
effectively like the spinful model by Kane and MeleKane and Mele (2005). Here
the model can be realized with real hopping amplitudes and without spin degree
which is favorable for realizations in artificial systems.

## Acknowledgments

The authors thank Shengyuan A. Yang, Y. X. Zhao and Zhijun Wang for helpful
discussions. This work is supported by the National Natural Science Foundation
of China (Grants No. 12174018, No. 12074024, and No. 11774018).

##  Appendix A Molecular dynamics study of the tbGPY

To assess the dynamical stability of the tbGPY, we performed an ab initio
molecular dynamics (AIMD) simulation at 100K. With a time step of 1.5 fs, for
total 3 ps. We find the variation of the total energy is small, and the
structure of tbGPY and the Moiré pattern do not change too much after the
simulation, as shown in Fig. A1. Therefore, we can conclude that there are no
obvious distortions of the geometries for tbGPY heterostructure at 100K, which
suggests that it is dynamically stable at 100K.

Figure A1: AIMD simulation results for tbGPY ($i=1$). (a) Variation of the
total energy during the AIMD simulation at 100K. (b) Lattice structure at the
end of the simulation period.

##  Appendix B Methods

###  B.1 First-principles calculations

Our first-principles calculation is carried out based on the density-
functional theory implemented in the Vienna ab initio simulation package
(VASP) Kresse and Hafner (1994); Kresse and Furthmüller (1996). The projector
augmented wave (PAW) method Blöchl (1994) was used for treating the ionic
potentials. The generalized gradient approximation (GGA) with the Perdew-
Burke-Ernzerhof (PBE) Perdew  _et al._ (1996) realization was adopted for the
exchange-correlation functional. The plane-wave cutoff energy was set to 520
eV. Only the $\Gamma$ point is chosen for the BZ sampling in the self-
consistent calculations, considering the large unit cell of the Moiré
structure. A vacuum of around 20 Åis included in the simulation to reduce the
unwanted interaction between any graphyne bilayers. The crystal structure
optimization is stopped when the forces on the ions are less than 0.01 eV/Å.
The energy convergence criteria is set as 10-6 eV for the electronic self-
consistent calculations. The van der Waals corrections are carried out in the
relaxation calculations by the semiempirical density-functional theory (DFT)-
D3 methodGrimme  _et al._ (2010) with Becke-Johnson dampingGrimme  _et al._
(2011).

To obtain the layer distance of tbGPY, we performed several full-atom
structural relaxation with different initial interlayer distances ranging from
3.664 Åto 2.964 Åwith an interval of 0.1 Å. The free energy after relaxation
reaches a minimum of -938.3827 eV with the corresponding layer distance of
3.444 Å, as presented in the main text.

###  B.2 SKTB model

We study the edge spectra and corner states using the SKTB model, following
the Slater and Koster scheme Slater and Koster (1954); Moon and Koshino
(2013). We build the SKTB model with the $p_{z}$ orbital of the carbon atom in
tbGPY. The general Hamiltonian of the SKTB model is written as

$H=-\sum_{|\mathbf{r}_{i}-\mathbf{r}_{j}|<L}h(\mathbf{r}_{i}-\mathbf{r}_{j})|\mathbf{r}_{i}\rangle\langle\mathbf{r}_{j}|+\text{H.c.},$
(B1)

where $\mathbf{r}_{i}$ is the position of the lattice point i, and
$|\mathbf{r}_{i}\rangle$ represents the atomic state at site i. The $L$ is the
maximum length for the hopping considered. The
$h(\mathbf{r}_{i}-\mathbf{r}_{j})$ is the hopping amplitude (transfer
integral) between site i and j, which can be calculated using

$\displaystyle-h(\mathbf{d})=$ $\displaystyle
V_{pp\pi}\left[1-\left(\frac{\mathbf{d}\cdot\hat{\mathbf{e}}_{z}}{d}\right)^{2}\right]+V_{pp\sigma}\left(\frac{\mathbf{d}\cdot\hat{\mathbf{e}}_{z}}{d}\right)^{2},$
(B2)

with

$\displaystyle V_{pp\pi}=$ $\displaystyle
V_{pp\pi}^{0}\exp\left(-\frac{d-a_{0}}{\delta_{0}}\right),$ (B3)
$\displaystyle V_{pp\sigma}=$ $\displaystyle
V_{pp\sigma}^{0}\exp\left(-\frac{d-d_{0}}{\delta_{0}}\right),$ (B4)

where the $V_{pp\pi}^{0}$ is the nearest neighbor p-p $\pi$ coupling, and
$V_{pp\sigma}^{0}$ is the nearest neighbor p-p $\sigma$ coupling. The nearest
distance in the monolayer is $a_{0}$, and the layer distance is $d_{0}$. To
consider the interlayer effect, $L>d_{0}$ should be satisfied.

##  Appendix C Fractional charges and robustness of corner states against
symmetry-breaking disorders

The fractional charge quantization of corner states is due to the filling
anomaly between the number of electrons required to satisfy the charge
neutrality and the crystalline symmetry. The tbGPY preserves
$\mathcal{C}_{2z}$ and $\mathcal{C}_{6z}$ symmetries which correspond to
fractional charges of $Q^{(2)}_{corner}$ and $Q^{(6)}_{corner}$ defined as
Benalcazar  _et al._ (2019)

$\displaystyle Q^{(2)}_{corner}$
$\displaystyle=\frac{e}{4}\left[-n_{C_{2}}^{X}-n_{C_{2}}^{Y}+\left(n_{C_{2}}^{M}+n_{C_{2}}^{\Gamma}\right)\right]~{}\text{mod}~{}e,$
(C1) $\displaystyle Q^{(6)}_{corner}$
$\displaystyle=e\left(\frac{n_{C_{2}}^{M}-n_{C_{2}}^{\Gamma}}{4}+\frac{n_{C_{3}}^{K}-n_{C_{3}}^{\Gamma}}{6}\right)~{}\text{mod}~{}e,$
(C2)

where the $n_{C_{k}}^{{i}}$ denotes the number of eigenvalue 1 of
$\hat{C}_{k}$ (alone $\hat{z}$) at the point ${i}$ in the momentum space. In
the presence of $\mathcal{C}_{6z}$ symmetry, the $Q^{(2)}_{corner}$ admits a
more intuitive form shown in Eq. 3 with Eq. 4. As a result, the nontrivial RCN
which is a bulk property corresponds to the fractional charge
$Q^{(2)}_{corner}=e/2$ that features a unique boundary property. However, the
$Q^{(6)}_{corner}$ does not give a one-to-one correspondence to the RCN,

$\displaystyle Q^{(6)}_{corner}$
$\displaystyle=e\left(\frac{\text{sgn}\left(n_{C_{2}}^{M}-n_{C_{2}}^{\Gamma}\right)\nu_{R}}{2}+\frac{n_{C_{3}}^{K}-n_{C_{3}}^{\Gamma}}{6}\right)~{}\text{mod}~{}e.$
(C3)

We find $Q^{(6)}_{corner}=e/6$ for the nontrivial tbGPY.

Then we explain this robustness from the point of view of fractional charge.
In higher order topological insulator materials with disorders, if the
crystalline symmetry, e.g., the $\mathcal{C}_{2z}$ or $\mathcal{C}_{6z}$, is
slightly broken (which does not invert bands), it may cause a ground state
filling of the in-gap states. This could shift the corner charges by only an
integer, and therefore the fractional portion of the charge remains intact at
each corners. We demonstrated the above analysis of the robustness of corner
state by introducing disorders along the edge. Specifically, we removed
several atoms at the corners (outlined in Fig. C1. ) of the hexagonal nano-
disk. The disorder breaks the local $\mathcal{C}_{2z}$, $\mathcal{C}_{6z}$ and
the $\mathcal{C}_{2y}$ crystalline symmetries in the nano-disk, but does not
spoil the corner charges.

Figure C1: Robustness of corner states on a disordered hexagonal nano-disk of
tbGPY $i=1$. The disorders are introduced at the corners enclosed by circles.

##  Appendix D RCN calculations for tbGPY and tbG with different twisting
angles

Table D1:  The RCNs for tbGPY systems with i = 1–4 calculated from first principles. The first two columns are the commensurate twisting angles and their labels; the third column is the number of valence band for first-principles calculations. The forth and fifth columns are the number of minus parity at the $\Gamma$ and M counted from the valence bands. The last column suggests the RCNs calculated from the forth and fifth columns using Eq. (2).  i | twisted angle (∘) | Number of valence band | $n^{\Gamma}_{-}$ | $n^{M}_{-}$ | RCN
---|---|---|---|---|---
1 | 21.787 | 224 | 110 | 112 | 1
2 | 13.174 | 608 | 302 | 304 | 1
3 | 9.430 | 1184 | 592 | 592 | 0
4 | 7.341 | 1952 | 976 | 976 | 0
Table D2:  The RCNs for tbG systems with several $i=$ 1-10 calculated from DFT. The arrangement of table is the same as Tab. D1 i | twisted angle (∘) | Number of valence band | $n^{\Gamma}_{-}$ | $n^{M}_{-}$ | RCN
---|---|---|---|---|---
1 | 21.787 | 56 | 24 | 30 | 1
2 | 13.174 | 152 | 72 | 78 | 1
3 | 9.430 | 296 | 144 | 150 | 1
4 | 7.341 | 488 | 240 | 246 | 1
5 | 6.009 | 728 | 360 | 366 | 1
6 | 5.086 | 1016 | 504 | 510 | 1
7 | 4.409 | 1352 | 672 | 678 | 1
8 | 3.890 | 1736 | 864 | 870 | 1
9 | 3.481 | 2168 | 1080 | 1086 | 1
10 | 3.150 | 2648 | 1320 | 1326 | 1

Using Eq. 2 in the main text with DFT, we calculate the RCNs for tbGPY and tbG
which are presented in Tab. D1 and D2. We find a trivialization of RCN in the
tbGPY with $i=3$ and 4. For the tbG systems, however, we have not observed
such trivialization of RCN up to $i=10$.

Now we discuss briefly the trivialization of the RCN in tbGPY with $i=3$.
Despite the similarities between graphene and $\alpha$-graphyne, the existence
of the acetylenic linkage makes some differences in the $\alpha$-graphyne. For
instance, the acetylenic linkage in the $\alpha$-graphyne inverts the energy
band, which in turn reverses the chirality between the graphene and the
$\alpha$-graphyne at K (or K’)Kim and Choi (2012); Jana  _et al._ (2019). The
trivialization in the RCN for tbGPY clearly indicates a double band inversion
as reducing the twisting angle, while for graphene, such band inversion does
not happen down to at least 3.15∘. This difference may be due to the existence
of the acetylenic linkage in the $\alpha$-graphyne.

##  Appendix E Higher-order topological origin from domain wall theory

The existence of corner states can also be understood in light of the domain
wall theory. Previous studies Wang and Moore (2019); Sheng  _et al._ (2019);
Chen  _et al._ (2020, 2021) suggest that a general SOTI system has the edge of
the form

$H_{edge}({k})=vk\sigma_{z}+m\sigma_{x}.$ (E1)

If the mass term is odd under a reflection symmetry which connects two edges,
then there should appear a domain wall state at the corner of the two edges.
In the 4BTB model the $\pi$-rotation about the y axis represented as
$\mathcal{C}_{2y}=\tau_{x}\sigma_{x}$ can serve as the reflection. Under
$\mathcal{C}_{2y}$, the velocity term reverses its sign with the mass term
unchanged. Therefore, it suggests that the mass term is odd under the
$\mathcal{C}_{2y}$, meaning that a domain wall state should appear at the
corner where the two edges connected by the $\mathcal{C}_{2y}$ cross. Further
considering the $\mathcal{C}_{3z}$ symmetry, there exist six corner states in
the hexagonal nano-disk.

## References

  * Hasan and Kane (2010) M. Z. Hasan and C. L. Kane, Rev. Mod. Phys. 82, 3045 (2010).
  * Qi and Zhang (2011) X.-L. Qi and S.-C. Zhang, Rev. Mod. Phys. 83, 1057 (2011).
  * Shen (2012) S.-Q. Shen,  _Topological Insulators_ (Springer, Berlin, 2012).
  * Bansil  _et al._ (2016) A. Bansil, H. Lin, and T. Das, Rev. Mod. Phys. 88, 021004 (2016).
  * Zhang  _et al._ (2013) F. Zhang, C. L. Kane, and E. J. Mele, Phys. Rev. Lett. 110, 046404 (2013).
  * Benalcazar  _et al._ (2017a) W. A. Benalcazar, B. A. Bernevig, and T. L. Hughes, Science 357, 61 (2017a).
  * Langbehn  _et al._ (2017) J. Langbehn, Y. Peng, L. Trifunovic, F. von Oppen, and P. W. Brouwer, Phys. Rev. Lett. 119, 246401 (2017).
  * Song  _et al._ (2017) Z. Song, Z. Fang, and C. Fang, Phys. Rev. Lett. 119, 246402 (2017).
  * Benalcazar  _et al._ (2017b) W. A. Benalcazar, B. A. Bernevig, and T. L. Hughes, Phys. Rev. B 96, 245115 (2017b).
  * Schindler  _et al._ (2018a) F. Schindler, A. M. Cook, M. G. Vergniory, Z. Wang, S. S. P. Parkin, B. A. Bernevig, and T. Neupert, Science Advances 4, eaat0346 (2018a).
  * Bradlyn  _et al._ (2017) B. Bradlyn, L. Elcoro, J. Cano, M. Vergniory, Z. Wang, C. Felser, M. I. Aroyo, and B. A. Bernevig, Nature 547, 298 (2017).
  * Schindler  _et al._ (2018b) F. Schindler, Z. Wang, M. G. Vergniory, A. M. Cook, A. Murani, S. Sengupta, A. Y. Kasumov, R. Deblock, S. Jeon, I. Drozdov, H. Bouchiat, S. Guéron, A. Yazdani, B. A. Bernevig, and T. Neupert, Nat. Phys. 14, 918 (2018b).
  * Yue  _et al._ (2019) C. Yue, Y. Xu, Z. Song, H. Weng, Y.-M. Lu, C. Fang, and X. Dai, Nature Physics 15, 577 (2019).
  * Wang  _et al._ (2019) Z. Wang, B. J. Wieder, J. Li, B. Yan, and B. A. Bernevig, Phys. Rev. Lett. 123, 186401 (2019).
  * Zhang  _et al._ (2019) T. Zhang, Y. Jiang, Z. Song, H. Huang, Y. He, Z. Fang, H. Weng, and C. Fang, Nature 566, 475 (2019).
  * Vergniory  _et al._ (2019) M. Vergniory, L. Elcoro, C. Felser, N. Regnault, B. A. Bernevig, and Z. Wang, Nature 566, 480 (2019).
  * Tang  _et al._ (2019a) F. Tang, H. C. Po, A. Vishwanath, and X. Wan, Nature Physics 15, 470 (2019a).
  * Tang  _et al._ (2019b) F. Tang, H. C. Po, A. Vishwanath, and X. Wan, Nature 566, 486 (2019b).
  * Xu  _et al._ (2019) Y. Xu, Z. Song, Z. Wang, H. Weng, and X. Dai, Phys. Rev. Lett. 122, 256402 (2019).
  * Sheng  _et al._ (2019) X.-L. Sheng, C. Chen, H. Liu, Z. Chen, Z.-M. Yu, Y. X. Zhao, and S. A. Yang, Phys. Rev. Lett. 123, 256402 (2019).
  * Liu  _et al._ (2019) B. Liu, G. Zhao, Z. Liu, and Z. F. Wang, Nano Letters 19, 6492 (2019).
  * Lee  _et al._ (2020) E. Lee, R. Kim, J. Ahn, and B.-J. Yang, npj Quantum Materials 5, 1 (2020).
  * Chen  _et al._ (2021) C. Chen, W. Wu, Z.-M. Yu, Z. Chen, Y. X. Zhao, X.-L. Sheng, and S. A. Yang, Phys. Rev. B 104, 085205 (2021).
  * Park  _et al._ (2019) M. J. Park, Y. Kim, G. Y. Cho, and S. B. Lee, Phys. Rev. Lett. 123, 216803 (2019).
  * Liu  _et al._ (2021) B. Liu, L. Xian, H. Mu, G. Zhao, Z. Liu, A. Rubio, and Z. F. Wang, Phys. Rev. Lett. 126, 066401 (2021).
  * Ezawa (2018) M. Ezawa, Phys. Rev. B 98, 045125 (2018).
  * Hitomi  _et al._ (2021) M. Hitomi, T. Kawakami, and M. Koshino, Phys. Rev. B 104, 125302 (2021).
  * Chen  _et al._ (2020) C. Chen, Z. Song, J.-Z. Zhao, Z. Chen, Z.-M. Yu, X.-L. Sheng, and S. A. Yang, Phys. Rev. Lett. 125, 056402 (2020).
  * Radha and Lambrecht (2020) S. K. Radha and W. R. L. Lambrecht, Phys. Rev. B 102, 115104 (2020).
  * Huang and Liu (2021) H. Huang and F. Liu, National Science Review (2021), https://doi.org/10.1093/nsr/nwab170.
  * Qian  _et al._ (2021) S. Qian, C.-C. Liu, and Y. Yao, Phys. Rev. B 104, 245427 (2021).
  * Pan  _et al._ (2022) M. Pan, D. Li, J. Fan, and H. Huang, npj Computational Materials 8, 1 (2022).
  * Zhao and Lu (2017) Y. X. Zhao and Y. Lu, Phys. Rev. Lett. 118, 056401 (2017).
  * Nakahara (2003) M. Nakahara,  _Geometry, Topology and Physics, 2nd ed._ (Institute of Physics, Bristol, 2003).
  * Song  _et al._ (2019) Z. Song, Z. Wang, W. Shi, G. Li, C. Fang, and B. A. Bernevig, Phys. Rev. Lett. 123, 036401 (2019).
  * Ahn  _et al._ (2019) J. Ahn, S. Park, and B.-J. Yang, Phys. Rev. X 9, 021013 (2019).
  * Wang  _et al._ (2020) K. Wang, J.-X. Dai, L. B. Shao, S. A. Yang, and Y. X. Zhao, Phys. Rev. Lett. 125, 126403 (2020).
  * Baughman  _et al._ (1987) R. H. Baughman, H. Eckhardt, and M. Kertesz, The Journal of Chemical Physics 87, 6687 (1987).
  * Malko  _et al._ (2012) D. Malko, C. Neiss, F. Viñes, and A. Görling, Phys. Rev. Lett. 108, 086804 (2012).
  * Li  _et al._ (2010) G. Li, Y. Li, H. Liu, Y. Guo, Y. Li, and D. Zhu, Chem. Commun. 46, 3256 (2010).
  * Chen  _et al._ (2017) Z. Chen, C. Molina-Jirón, S. Klyatskaya, F. Klappenberger, and M. Ruben, Annalen der Physik 529, 1700056 (2017).
  * Jia  _et al._ (2017) Z. Jia, Y. Li, Z. Zuo, H. Liu, C. Huang, and Y. Li, Accounts of Chemical Research 50, 2470 (2017).
  * Gao  _et al._ (2019) X. Gao, H. Liu, D. Wang, and J. Zhang, Chem. Soc. Rev. 48, 908 (2019).
  * Jana  _et al._ (2019) S. Jana, A. Bandyopadhyay, and D. Jana, Phys. Chem. Chem. Phys. 21, 13795 (2019).
  * Kim and Choi (2012) B. G. Kim and H. J. Choi, Phys. Rev. B 86, 115435 (2012).
  * Bandyopadhyay and Jana (2020) A. Bandyopadhyay and D. Jana, Reports on Progress in Physics 83, 056501 (2020).
  * Jana  _et al._ (2021) S. Jana, A. Bandyopadhyay, S. Datta, D. Bhattacharya, and D. Jana, Journal of Physics: Condensed Matter 34, 053001 (2021).
  * Benalcazar  _et al._ (2019) W. A. Benalcazar, T. Li, and T. L. Hughes, Phys. Rev. B 99, 245151 (2019).
  * Lopes dos Santos  _et al._ (2007) J. M. B. Lopes dos Santos, N. M. R. Peres, and A. H. Castro Neto, Phys. Rev. Lett. 99, 256802 (2007).
  * Leenaerts  _et al._ (2013) O. Leenaerts, B. Partoens, and F. M. Peeters, Applied Physics Letters 103, 013105 (2013).
  * Slater and Koster (1954) J. C. Slater and G. F. Koster, Phys. Rev. 94, 1498 (1954).
  * Chen  _et al._ (2022) C. Chen, X.-T. Zeng, Z. Chen, Y. X. Zhao, X.-L. Sheng, and S. A. Yang, Phys. Rev. Lett. 128, 026405 (2022).
  * Ahn  _et al._ (2018) J. Ahn, D. Kim, Y. Kim, and B.-J. Yang, Phys. Rev. Lett. 121, 106403 (2018).
  * Zou  _et al._ (2018) L. Zou, H. C. Po, A. Vishwanath, and T. Senthil, Phys. Rev. B 98, 085435 (2018).
  * Imhof  _et al._ (2018) S. Imhof, C. Berger, F. Bayer, J. Brehm, L. W. Molenkamp, T. Kiessling, F. Schindler, C. H. Lee, M. Greiter, T. Neupert, and R. Thomale, Nature Physics 14, 925 (2018).
  * Ma  _et al._ (2019) G. Ma, M. Xiao, and C. T. Chan, Nature Reviews Physics 1, 281 (2019).
  * Xue  _et al._ (2020) H. Xue, Y. Ge, H.-X. Sun, Q. Wang, D. Jia, Y.-J. Guan, S.-Q. Yuan, Y. Chong, and B. Zhang, Nature Communications 11, 2442 (2020).
  * Xue  _et al._ (2019) H. Xue, Y. Yang, F. Gao, Y. Chong, and B. Zhang, Nature materials 18, 108 (2019).
  * Ozawa  _et al._ (2019) T. Ozawa, H. M. Price, A. Amo, N. Goldman, M. Hafezi, L. Lu, M. C. Rechtsman, D. Schuster, J. Simon, O. Zilberberg, and I. Carusotto, Rev. Mod. Phys. 91, 015006 (2019).
  * Mittal  _et al._ (2019) S. Mittal, V. V. Orre, G. Zhu, M. A. Gorlach, A. Poddubny, and M. Hafezi, Nature Photonics 13, 692 (2019).
  * Prodan and Prodan (2009) E. Prodan and C. Prodan, Phys. Rev. Lett. 103, 248101 (2009).
  * Dalibard  _et al._ (2011) J. Dalibard, F. Gerbier, G. Juzeliūnas, and P. Öhberg, Rev. Mod. Phys. 83, 1523 (2011).
  * Kane and Mele (2005) C. L. Kane and E. J. Mele, Phys. Rev. Lett. 95, 1 (2005).
  * Kresse and Hafner (1994) G. Kresse and J. Hafner, Phys. Rev. B 49, 14251 (1994).
  * Kresse and Furthmüller (1996) G. Kresse and J. Furthmüller, Phys. Rev. B 54, 11169 (1996).
  * Blöchl (1994) P. E. Blöchl, Phys. Rev. B 50, 17953 (1994).
  * Perdew  _et al._ (1996) J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).
  * Grimme  _et al._ (2010) S. Grimme, J. Antony, S. Ehrlich, and H. Krieg, The Journal of Chemical Physics 132, 154104 (2010).
  * Grimme  _et al._ (2011) S. Grimme, S. Ehrlich, and L. Goerigk, Journal of Computational Chemistry 32, 1456 (2011).
  * Moon and Koshino (2013) P. Moon and M. Koshino, Phys. Rev. B 87, 205404 (2013).
  * Wang and Moore (2019) Y.-Q. Wang and J. E. Moore, Phys. Rev. B 99, 155102 (2019).


# Living-off-the-Land Abuse Detection Using Natural Language Processing and
Supervised Learning

Ryan Stamp  School of Computer Science
University of Guelph
Guelph, ON, Canada
<EMAIL_ADDRESS>
###### Abstract

Living-off-the-Land is an evasion technique used by attackers where native
binaries are abused to achieve malicious intent. Since these binaries are
often legitimate system files, detecting such abuse is difficult and often
missed by modern anti-virus software. This paper proposes a novel abuse
detection algorithm using raw command strings. First, natural language
processing techniques such as regular expressions and one-hot encoding are
utilized for encoding the command strings as numerical token vectors. Next,
supervised learning techniques are employed to learn the malicious patterns in
the token vectors and ultimately predict the command’s label. Finally, the
model is evaluated using statistics from the training phase and in a virtual
environment to compare its effectiveness at detecting new commands to existing
anti-virus products such as Windows Defender.

###### Index Terms:

Living-off-the-land, Natural language processing, Regular expressions,
Supervised learning, Intrusion detection system, Machine learning, Endpoint
security.

##  I Introduction

As malware detection methods become more sophisticated, adversaries constantly
look for new approaches to remain undetected in target networks [1]. One such
method that has become popular in recent years is living-off-the-land (LOL).
This technique involves leveraging existing benign binaries (LOLBins) to
achieve the desired result. These binaries are typically legitimate system
files that are required for normal operating system execution; however, they
also possess the ability to perform particular procedures desired by
adversaries, such as downloading files from a remote source.

In a study by Barr-Smith et al. [2], a large range of malware samples were
analyzed (totaling nearly 32 million samples), and of them, 9.6% were found to
utilize native binaries. This statistic is non-negligible and demonstrates the
real threat LOL techniques possess.

However, detecting and/or preventing LOLBin abuse presents unique challenges
compared to other malware detection forms. The primary issue is that these
binaries are legitimate system files, so existing techniques such as signature
detection will not work here. Instead, more sophisticated methods must be
used.

In this paper, we propose a novel command string-based detection algorithm.
Our hypothesis states that certain tokens, or certain combinations of tokens,
will be used more frequently in malicious commands than in benign commands. By
employing natural language processing (NLP) techniques such as regular
expressions and one-hot encoding, we can parse and represent the tokens of a
command string as numerical vectors [3]. From here, we train a supervised
learning model to predict how malicious a token of a command string is likely
to be and, ultimately, how malicious the command itself is. Finally, we
evaluate the model to determine the effectiveness of the training set and its
ability to detect new threats not already detectable by existing anti-virus
products.

##  II Background and Related Works

Nowadays, cybersecurity is a hot topic for many researchers in the academic
and industrial sections. We can find a lot of security mechanisms and defenses
in attack detection and threat hunting in [4, 5, 6, 7]. Moreover, a
considerable number of studies have proposed other security solutions like
machine learning and blockchain to design and implement detection models [8,
9]. A few of these security solutions are the detection of outliers and
anomalous data and behavior in Internet of Things (IoT) environments, both of
which have been of interest to the academic and industrial research
communities, respectively [10, 11]. Studies have looked at the analysis of IoT
attacks and the countermeasures based on the devices connected to the
intelligent networks[12]. The detection of outliers in IoTs, Industrial
Internet of Things (IIoT), and similar environments is well known for its
notable performance using Artificial Intelligence (AI)-based methods,
including machine learning-based and deep learning-based approaches.

In order to identify anomalous behavior on IIoT sites, Hasan et al. [13] have
compared the performance of machine learning models by measuring precision,
accuracy, f1-score, and recall. For techniques like Artificial Neural Networks
(ANN), Random Forest (RF), and Decision Tree (DT), 99.4% accuracy was reported
by the system. Metrics such as f1-score and recall are better performed by the
RF model than other techniques. An ensemble technique based on the Average
One-Dependence Estimator (AODE) and RF model is presented by Jabbar, and his
colleague in [14]. Based on the Kyoto data, the proposed technique showed a
high accuracy of 90.51 percent and a low False Alarm Rate (FAR) of 0.14
percent. Furthermore, AL-Hawawreh et al. [15] have developed a deep learning-
based ensemble model for detecting anomalies in IIoT traffic. With the help of
TCP/IP data packets, they have used a feed-forward and AE architecture for
training and validating.

Research regarding LoLBin detection has mostly been conducted in recent years.
The main concern with LoLBin detection is that the adversary may not create
any malicious files for anti-virus software to detect [16][17]; nor can these
binaries simply be blocked since they provide necessary functionally to the
operating system. A common example is a certutil.exe binary native to Windows
operating systems. This tool is used for managing certificates but possesses
the ability to download files that adversaries can exploit. To overcome this
issue, researchers have turned to AI to detect the malicious use of otherwise
benign binaries.

In a study by Ongun et al. [16], the issues with current methods are outlined,
relying on command line parsing by regular expressions that yield high false-
positive rates. The study proposes an active learning approach where an AI
will select samples for an analyst to classify. The AI works under the same
hypothesis that we propose; that certain tokens in the command appear more
frequently for malicious commands. Although the results from the study are
promising, basing decisions on the tokens alone is very binary-specific as
different binaries will often use different tokens in their commands to
express the same functionality. Our research will improve on this limitation
by developing custom parsing rules to determine the functionality of specific
commands so that the same patterns can be recognized across different
binaries.

We also want to briefly mention that LOLBin detection also has use in the
mobile space. In a paper by Bellizzi et al. [18], Just-in-Time Memory
Forensics (JIT-MF) is leveraged to detect and fight against LOLBin abuse in
Android systems. However, although mobile forensics and security are arguably
just as important as forensics in desktop operating systems, this paper and
its proposed method will solely focus on desktop environments.

##  III Methodologies

This section will discuss the methodologies behind our LOLBin detection. It
will discuss the origin of the dataset, the pre-processing and tokenization
techniques, the building of the feature vectors and feature matrix, and
finally, the classifier model and its prediction interface.

###  III-A Dataset

The dataset consists of raw command strings for a variety of known LOLBins
labeled as either benign or malicious. Since there are likely hundreds of
LOLBins that malicious parties could abuse, we chose to focus on only a small
subset of the most popular. In total, our model supports 16 different LOLBins
listed as follows: bitsadmin, certutil, cmstp, csc, cscript, mmc, msiexec,
msxsl, reg, regsvcs, regsvr32, rundll32, schtasks, sqlps, wmic, and wscript.
Additionally, our methodologies were designed to make the process of adding
new LOLBins relatively trivial, and the such process will be discussed in the
later sections.

Benign commands for the LOLBins are rather plentiful. Our dataset they were
collected from our partner organization’s raw dataset of collected process
events from numerous client machines. Malicious samples, by comparison, are
harder to obtain. The primary struggle is a lack of readily available samples
that have been hand-labeled by experienced analysts. However, for our model,
we are less interested in learning specific commands and more interested in
learning specific patterns.

LOLBAS [19] is an excellent community-run project for documenting the various
LOLBins and providing examples of how they can be abused. Since our pre-
processing stages will extract the patterns from command strings, we can
collect the malicious samples for our dataset from these LOLBAS examples.

The DFIR Report [20] is another great resource for obtaining malicious
commands. This site hosts reports on real intrusions and lists the many
commands that were executed. Any commands involving LOLBins from this resource
were collected and added to the validation dataset, which is used later during
the evaluation of the model.

The final challenge with our dataset is the imbalance of malicious and benign
samples. Since the vast majority of commands encountered in the wild will be
benign, it follows that our dataset will also possess the same or similar
ratio. However, this is an issue as imbalanced datasets can heavily degrade
the effectiveness of the model [21].

Synthetic Minority Oversampling Technique (SMOTE) [21] is an algorithm for
balancing a dataset. In short, the algorithm generates new feature vectors by
creating line segments between minority class neighbors. These line segments
allow for the generation of new vectors to balance the dataset but without
introducing any new information that could alter the model’s training. In
practice, this technique proved very effective.

###  III-B Lexers and Tokenization

Since the model operates on raw command strings, custom-built lexers are first
employed to tokenize the input strings.

The lexers will first split the command string into tokens. Characters
separated by whitespace are considered different tokens unless a group of
whitespace-separated characters is wrapped in quotations. For example,
“`"Windows Defender"`” would be considered a single token. Further
normalization is applied for each token generated in this manner to decrease
the number of unique tokens the model will encounter. This normalization
involves removing the leading or trailing quotes, dashes, or slashes. The
tokens are also converted to lowercase to enable case insensitive matching
later on. Figure 1 demonstrates this process visually.

Figure 1: Tokens are normalized

To further decrease the token space, the lexers will next apply several
regular expressions to the tokens to detect particular patterns and replace
that token with a special token. For example, a url regex is employed to
detect if a token is an HTTP or HTTPS url and replaces it with the “$<$url$>$”
token [22, 23]. The file regex detects file names and replaces that token with
two special tokens, a generic “$<$file$>$” token and a specific extension
token such as “$<$ext_exe$>$”. Figure 2 demonstrates this process visually.
For a full list of custom tokens, see appendix -A.

Figure 2: Detected patterns are replaced with custom tokens

Do note that since the syntax for every binary can differ greatly, a unique
lexer is required for each unique binary. To accomplish this while reusing as
much logic as possible, a base class that implements independent binary
patterns, such as the URL and file regexes discussed above, was implemented
and can simply be instanced to fulfill the requirements for most binaries.
However, if necessary, the base class can be extended to more easily add
additional pattern recognition if a particular binary requires it. For
example, “rundll32” implements a subclass to detect javascript parameters.

###  III-C Feature Vectors and Feature Matrix

Building feature vectors from command tokens present a unique challenge.
Vectors for machine learning models must be numerical [24, 25]. However, our
tokenized strings are words. Thus, some encoding of words to floats must be
employed. Although many such techniques exist, the relatively simple one-hot
encoding technique was chosen.

Since there are theoretically a nearly infinite number of unique tokens, we
must determine which ones are the most important for learning benign versus
malicious behavior and only include those tokens in the one-hot encoding.
Additionally, we add a special “$<$rare$>$” token to act as a dust-bin
category for all unrecognized tokens.

To this extent, as with the lexers, we will need a unique feature builder for
each unique binary. The reason is that different binaries will have different
tokens that are important to determining the command’s label and therefore
need to be included in the one-hot encoding.

To build the corpus, the entire dataset for the binary is iterated through,
and each command is tokenized. Every token that appears in the malicious
commands is considered important and is included in the word corpus. Since
there are many more unique examples of benign commands compared to malicious
samples, only tokens that appear in at least three unique tokenized commands
are included in the corpus.

Figure 3 demonstrates visually how a sample tokenized command may be one-hot
encoded. Also, since building the word corpus can be expensive depending on
the dataset size for a binary, the word corpuses are stored alongside the
model artifacts when saving to disk so that they do not need to be regenerated
every time.

Figure 3: Example one-hot encoding for a collection of tokens

Now that we can encode our tokens numerically, we next need to build the
actual vectors. To achieve better model performance, we not only want to
encode what tokens are present in command but also the context of what other
tokens appear around them. To achieve this, we employ a novel token feature-
building algorithm. Note that this algorithm, as described below, builds a
feature vector for each token instead of one vector for the entire command.

Firstly, the feature corresponding to the token being examined is set to one.
Next, features corresponding to the tokens away from the target token are set
to 1/2. Two away are set assigned values of 1/3 and so on. The window size for
how many tokens to consider can be adjusted, but we used a window size of 2.
This process is also additive, so if the same token appears more than once in
the window around the target token, then this is reflected in the final
vector. Figure 4 demonstrates this algorithm visually.

Figure 4: Individual tokens are converted into feature vectors based on what
tokens are around them

###  III-D LOLBin Models

So far, we have needed to create unique lexers and feature builders for each
unique binary, and the final model is no different. Because the feature
vectors are different lengths due to our unique one-hot encoding per binary, a
separate model will need to be trained for each binary. At first, this may
seem problematic, but this approach boasts several advantages:

  1. 1.

If an existing LOLBin’s dataset needs to be updated, then the LOLBin’s model
can easily be retrained without having to touch the other models.

  2. 2.

Similarly, adding support for a new LOLBin can be done simply by adding the
new model and the performances of the other models will be completely
unaffected.

  3. 3.

Re-balancing the dataset using the SMOTE utility (see section III-A) can more
easily be done on a per binary basis instead of on the entire dataset.

The final training matrices consist of a token vector for every token in every
command in the dataset. The train-test split is done with a standard 80%
training and 20% testing. Note that some more obscure LOLBins have no benign
samples in the dataset. In this rare case, that binary model is simply trained
on the entire dataset with no splitting or testing.

Two different types of supervised learning models are utilized. SKlearn’s
multi-layer perceptron (MLP) is the preferred classifier; however, like all
neural network-based models, a sufficient amount of data is required to be
effective. Suppose the LOLBin in question does not meet the required
threshold, in our case, 500 samples. In that case, a Random Forest (RF)
classifier is used instead, which can achieve much better performance when
less data is available.

Additionally, since the model is trained on token vectors, predictions will be
in the form of token scores, not command scores. To obtain command scores, a
feature matrix is constructed by combining the token vectors for each token in
the command. Submitting the matrix to the model will produce a token score for
each row in the matrix. Aggregation techniques can now be employed to
interpret the results, such as min, max, and average pooling. In our testing,
the max value among the token scores proved to be the most reliable in
determining the command’s label. Figure 5 demonstrates this process visually.
Note the token scores used are only examples and may not reflect the true
token scores generated by the model.

Figure 5: Command score is calculated form multiple token scores

###  III-E Unimodel

Since we trained a separate model for each LOLBin, we have developed a
container class, dubbed the ”unimodal”, to encapsulate the models into a
single object. The unimodeladditionally provides methods for loading the model
artifacts from disk and a single interface for making predictions. This method
uses techniques such as regular expressions to extract the binary name from
the command string and determine which model to query. Figure 6 demonstrates
this process visually. The use of the unimodel proved very effective. It
solved any issues related to the organization that may arise from using
multiple models.

Figure 6: The unimodel encapsulates individual LOLBin models and provides
single interface for predictions

##  IV Evaluation

###  IV-A Unimodel Statistics

The statistics for each LOLBin model are listed in table I. Note that the
three LOLBins with no benign samples are excluded. It can be seen that the
models performed very well overall with near-perfect statistics.

TABLE I: Model evaluation in VM results | Classifier | Accuracy | Precision | Recall | F1 Score
---|---|---|---|---|---
bitsadmin | RF | 1.0 | 1.0 | 1.0 | 1.0
certutil | RF | 1.0 | 1.0 | 1.0 | 1.0
csc | MLP | 1.0 | 1.0 | 1.0 | 1.0
cscript | MLP | 1.0 | 1.0 | 1.0 | 1.0
mmc | RF | 1.0 | 1.0 | 1.0 | 1.0
msiexec | MLP | 0.9940 | 0.9884 | 1.0 | 0.9942
reg | MLP | 0.9987 | 1.0 | 0.9973 | 0.9986
regsvr32 | MLP | 0.9757 | 1.0 | 0.9451 | 0.9718
rundll32 | MLP | 0.9739 | 1.0 | 0.9480 | 0.9733
schtasks | RF | 0.9835 | 0.9836 | 0.9836 | 0.9836
sqlps | RF | 1.0 | 1.0 | 1.0 | 1.0
wmic | MLP | 0.9984 | 0.9968 | 1.0 | 0.9984
wscript | MLP | 1.0 | 1.0 | 1.0 | 1.0
AVERAGE | — | 0.9941 | 0.9976 | 0.9903 | 0.9938

These statistics also help to emphasize the advantages of modular models.
Since each model is trained specifically for one LOLBin, it can much more
easily obtain better accuracy than if it was trained for multiple LOLBins.
Separate stats additionally allows analysts to zero in on which models may be
underperforming. From here, their datasets can be updated to improve
performance; and as mentioned earlier, the models can be trained individually
without affecting the others.

###  IV-B Virtual Machine Simulation

Although statistics can give some idea of how well the model is performing,
its real value comes from its ability to detect new commands that are not
already detectable by existing solutions [26]. We also want to validate how
effective the model is at detecting new commands not present in the training
set, as most of the commands encountered in the wild will be new. For this
purpose, we utilize the validation set discussed in section III-A.

The model was evaluated against Windows Defender and our partner
organization’s own Endpoint Detection and Response (EDR) solution. The
evaluation consisted of running 83 malicious commands (spanning all LOLBins)
in a virtual machine (VM) with both Defender enabled and the EDR agent
installed. These 83 commands combine the validation dataset and select
commands from the training dataset.

The results are summarized in table II. As can be seen, Defender can currently
only detect a mere 8 of the malicious LOLBin commands. In contrast, our model
can detect all of them. Also, note that for confidential reasons, the EDR
results are not disclosed in this paper and are only shared internally with
the partner organization.

TABLE II: Model evaluation in VM results | Defender | Model
---|---|---
bitsadmin | 0/4 | 4/4
certutil | 2/8 | 8/8
cmstp | 0/2 | 2/2
csc | 0/2 | 2/2
cscript | 0/2 | 2/2
mmc | 0/1 | 1/1
msiexec | 0/5 | 5/5
msxsl | 0/2 | 2/2
reg | 0/19 | 19/19
regsvcs | 1/1 | 1/1
regsvr32 | 1/2 | 2/2
rundll32 | 2/8 | 8/8
schtasks | 1/10 | 10/10
sqlps | 0/1 | 1/1
wmic | 1/13 | 13/13
wscript | 0/3 | 3/3
TOTAL | 8/83 | 83/83

This extra level of evaluation provides much greater confidence in the model’s
effectiveness than utilizing the model statistics alone. Comparing against
existing solutions such as Windows Defender and the partner organization’s EDR
solution additionally helps to prove the model and algorithm’s worth and value
to detect new threats.

###  IV-C False Positives and the Base-Rate Fallacy

In this sub-section, we would like to bring attention to the base-rate fallacy
and how it affects the problem of detecting malicious commands.

Axelsson [27] excellently describes the base-rate fallacy using the following
example. Suppose a disease affects only 1/10000 people in a population, and a
test for that disease has a 1% false-positive rate. Then, if an individual
tests positive for the disease using the said test, there is still only a
1/100 chance that they actually have the disease. Intuitively, humans tend to
disregard the actual base rate and zero in on the low false-positive rate, so
the final results seem rather counter-intuitive. Hence the name “base-rate
_fallacy_ ”.

We now extend this example to the problem of malicious command detection.
Malicious commands make up a tiny percentage of all commands executed on a
system, and the rate at which commands are executed is very high. So, by the
base-rate fallacy, even if our model (or any other solution) has a very low
false positive rate, the resulting hits are still much more likely to be
benign than actually malicious.

As part of our evaluation, we ran our model on a live stream of command data
generated from numerous clients. The goal was to check for false positives and
adjust the models accordingly. However, the sheer amount of false positives
was much higher than expected, which we later determined was due to the base-
rate fallacy.

This result presents a problem for analysts who monitor alerts generated from
such systems. Since the flow of command data to the model is very high,
analysts will likely be flooded with false positive alerts that obscure true
positives and, in general, consume a lot of resources. Therefore, we recommend
that any solution to the malicious command detection problem (whether our
presented solution or another) not generate alerts outright, or at the very
least, also apply some form of filtering such as whitelisting.

This result also emphasizes the need for human analysts, or more specifically,
that this problem may never reach a point of complete automation. Although it
may be tempting to adjust the models more to further lower the false positive
rate, one would be doing so with the risk of increasing the false negative
rate. It is better to generate more false positives and have human analysts
available to sift through them than to lower the rate of incorrect alerts and
run the risk of missing true malicious activity.

##  V Extensibility

As stated previously, the unimodel and related procedures were designed to be
extendable. This feature is essential in the context of LOLBin detection, as
new LOLBins and techniques are being discovered regularly as adversaries
constantly adapt their methods to the changing landscape. If adding new
patterns and/or retraining the model was a lengthy process, it could act as a
bottleneck for the system. Additionally, making sure the unimodel stays
modular is important for upholding consistency by not affecting the
performance of one LOLBin’s detection by adjusting the model of another
unrelated LOLBin.

##  VI Conclusions and Future Work

This paper proposed a novel LOLBin command abuse detection algorithm using
both regular expressions and supervised learning. Model statistics demonstrate
high efficiency at detecting malicious patterns, and further evaluation in a
virtual environment shows the model’s ability to detect many commands not
already detectable by existing anti-virus software such as Windows Defender.
The base-rate fallacy was also briefly discussed, and its applications to the
LOLBin problem were analyzed.

In terms of future work, support for more LOLBins and/or more token patterns
could be added as needed. More sophisticated techniques, in general, may also
be required to handle more complex LOLBins such as cmd. Finally, other
techniques for reducing false positives could be explored, as was covered in a
previous section.

###  -A Custom Tokens

Common (used by all lexers):

_$ <$url$>$_ : http or https url
eg. `https://malicious[.]com`

_$ <$url$>$ $<$ext_*$>$_ : url that ends with a path to a file
eg. `https://malicious[.]com/file.exe`

_$ <$ads$>$ _: alternate data stream
eg. `benign.txt:malicious.exe`

_$ <share>$ $<$ext_*$>$_ : shared file
eg. `\\10.10.10.10\malicious\file.exe`

_$ <number>$_ : positive or negative integer
eg. `100`

_$ <decimal>$_ : positive or negative real number
eg. `10.5`

_$ <ip\\_addr>$_ : IP address
eg. `192.168.1.0`

_$ <guid>$_ : Globally Unique Identifier (GUID)
eg. `{e77b42d3-55a5-4b3e-9d08-d59047c2e4c8}`

_$ <benign\\_keyword>$_ : LOLBin specific keyword indicating benign behaviour

_$ <mal\\_keyword>$_ : LOLBin specific keyword indicating malicious behaviour

regsvr32:

_$ <script>$ *_ : local or remote script
eg. `i:file.sct`

rundll32:

_$ <javascript>$_ : javascript
eg. `javascript:*`

_$ <dll>$ * *_ : DLL file and entrypoint
eg. `file.dll,EntryPoint`

## References

  * [1] A. Yazdinejad, H. HaddadPajouh, A. Dehghantanha, R. M. Parizi, G. Srivastava, and M.-Y. Chen, “Cryptocurrency malware hunting: A deep recurrent neural network approach,” Applied Soft Computing, vol. 96, p. 106630, 2020.
  * [2] F. Barr-Smith, X. Ugarte-Pedrero, M. Grazian, R. Spolaor, and I. Martinovic, “Survivalism: Systematic analysis of windows malware living-off-the-land,” in Symposium on Security and Privacy (SP), pp. 1557–1574, IEEE, 2021.
  * [3] A. Yazdinejad, R. M. Parizi, A. Dehghantanha, and H. Karimipour, “Federated learning for drone authentication,” Ad Hoc Networks, vol. 120, p. 102574, 2021.
  * [4] A. Yazdinejad, A. Dehghantanha, R. M. Parizi, M. Hammoudeh, H. Karimipour, and G. Srivastava, “Block hunter: Federated learning for cyber threat hunting in blockchain-based iiot networks,” IEEE Transactions on Industrial Informatics, pp. 1–1, 2022.
  * [5] A. Corallo, M. Lazoi, M. Lezzi, and A. Luperto, “Cybersecurity awareness in the context of the industrial internet of things: A systematic literature review,” Computers in Industry, vol. 137, p. 103614, 2022.
  * [6] A. Yazdinejad, B. Zolfaghari, A. Azmoodeh, A. Dehghantanha, H. Karimipour, E. Fraser, A. G. Green, C. Russell, and E. Duncan, “A review on security of smart farming and precision agriculture: Security aspects, attacks, threats and countermeasures,” Applied Sciences, vol. 11, no. 16, p. 7518, 2021\.
  * [7] A. Yazdinejad, R. M. Parizi, A. Dehghantanha, Q. Zhang, and K.-K. R. Choo, “An energy-efficient sdn controller architecture for iot networks with blockchain-based security,” IEEE Transactions on Services Computing, vol. 13, no. 4, pp. 625–638, 2020.
  * [8] S. Nakhodchi, B. Zolfaghari, A. Yazdinejad, and A. Dehghantanha, “Steeleye: An application-layer attack detection and attribution model in industrial control systems using semi-deep learning,” in 2021 18th International Conference on Privacy, Security and Trust (PST), pp. 1–8, IEEE, 2021.
  * [9] A. Yazdinejad, R. M. Parizi, A. Dehghantanha, and K.-K. R. Choo, “Blockchain-enabled authentication handover with efficient privacy protection in sdn-based 5g networks,” IEEE Transactions on Network Science and Engineering, vol. 8, no. 2, pp. 1120–1132, 2019.
  * [10] X. Yan, Y. Xu, X. Xing, B. Cui, Z. Guo, and T. Guo, “Trustworthy network anomaly detection based on an adaptive learning rate and momentum in iiot,” IEEE Transactions on Industrial Informatics, vol. 16, no. 9, pp. 6182–6192, 2020.
  * [11] G. Han, J. Tu, L. Liu, M. Martínez-García, and Y. Peng, “Anomaly detection based on multidimensional data processing for protecting vital devices in 6g enabled massive iiot,” IEEE Internet of Things Journal (Early Access Article), pp. 1–1, 2021.
  * [12] A. Yazdinejad, R. M. Parizi, G. Srivastava, and A. Dehghantanha, “Making sense of blockchain for ai deepfakes technology,” in 2020 IEEE Globecom Workshops (GC Wkshps, pp. 1–6, IEEE, 2020.
  * [13] M. Hasan, M. M. Islam, M. I. I. Zarif, and M. Hashem, “Attack and anomaly detection in iot sensors in iot sites using machine learning approaches,” Internet of Things, vol. 7, p. 100059, 2019.
  * [14] M. Jabbar, R. Aluvalu, et al., “Rfaode: A novel ensemble intrusion detection system,” Procedia computer science, vol. 115, pp. 226–234, 2017\.
  * [15] A.-H. Muna, N. Moustafa, and E. Sitnikova, “Identification of malicious activities in industrial internet of things based on deep learning models,” Journal of information security and applications, vol. 41, pp. 1–11, 2018\.
  * [16] T. Ongun, J. W. Stokes, J. B. Or, K. Tian, F. Tajaddodianfar, J. Neil, C. Seifert, A. Oprea, and J. C. Platt, “Living-off-the-land command detection using active learning,” 2021.
  * [17] R. Tarek, S. Chaimae, and C. Habiba, “Runtime api signature for fileless malware detection,” in Advances in Intelligent Systems and Computing, vol 1129, 2020.
  * [18] J. Bellizzi1, M. Vella1, C. Colombo1, and J. Hernandez-Castro, “Responding to living-off-the-land tactics using just-in-time memory forensics (jit-mf) for android,” 2021.
  * [19] “Lolbas: Living off the land binaries, scripts and libraries.”
  * [20] “The dfir report.”
  * [21] N. V. Chawla, K. W. Bowyer, L. O. Hall, and W. P. Kegelmeyer, “Smote: Synthetic minority over-sampling technique,” in Journal of Artificial Intelligence Research 16, pp. 321–357, 2002.
  * [22] A. Yazdinejad, R. M. Parizi, A. Dehghantanha, and K.-K. R. Choo, “P4-to-blockchain: A secure blockchain-enabled packet parser for software defined networking,” Computers & Security, vol. 88, p. 101629, 2020.
  * [23] A. Yazdinejad, A. Bohlooli, and K. Jamshidi, “Efficient design and hardware implementation of the openflow v1. 3 switch on the virtex-6 fpga ml605,” The Journal of Supercomputing, vol. 74, no. 3, pp. 1299–1320, 2018.
  * [24] A. Yazdinejad, G. Srivastava, R. M. Parizi, A. Dehghantanha, K.-K. R. Choo, and M. Aledhari, “Decentralized authentication of distributed patients in hospital networks using blockchain,” IEEE journal of biomedical and health informatics, vol. 24, no. 8, pp. 2146–2156, 2020.
  * [25] R. M. Parizi, S. Homayoun, A. Yazdinejad, A. Dehghantanha, and K.-K. R. Choo, “Integrating privacy enhancing techniques into blockchains using sidechains,” in 2019 IEEE Canadian Conference of Electrical and Computer Engineering (CCECE), pp. 1–4, IEEE, 2019.
  * [26] A. Yazdinejad, R. M. Parizi, A. Dehghantanha, G. Srivastava, S. Mohan, and A. M. Rababah, “Cost optimization of secure routing with untrusted devices in software defined networking,” Journal of Parallel and distributed Computing, vol. 143, pp. 36–46, 2020.
  * [27] S. Axelsson, “The base-rate fallacy and the difficulty of intrusion detection,” in ACM Transactions on Information and System Security, Vol. 3, No. 3, pp. 186–205, IEEE, 2000.