# One-dimensional topological channels in heterostrained bilayer graphene

Nina C. Georgoulea  School of Physics, AMBER & CRANN Institute, Trinity
College Dublin, Dublin 2, Ireland  Nuala M. Caffrey  School of Physics,
University College Dublin, Dublin 4, Ireland  Centre for Quantum Engineering,
Science, and Technology, University College Dublin, Dublin 4, Ireland  Stephen
R. Power<EMAIL_ADDRESS>School of Physics, AMBER & CRANN Institute,
Trinity College Dublin, Dublin 2, Ireland  School of Physical Sciences, Dublin
City University, Glasnevin, Dublin 9, Ireland

###### Abstract

The domain walls between AB- and BA-stacked gapped bilayer graphene have
garnered intense interest as they host topologically-protected, valley-
polarised transport channels. The introduction of a twist angle between the
bilayers and the associated formation of a Moiré pattern has been the dominant
method used to study these topological channels, but heterostrain can also
give rise to similar stacking domains and interfaces. Here, we theoretically
study the electronic structure of a uniaxially heterostrained bilayer
graphene. We discuss the formation and evolution of interface-localized
channels in the one-dimensional Moiré pattern that emerges due to the
different stacking registries between the two layers. We find that a uniform
heterostrain is not sufficient to create one-dimensional topological channels
in biased bilayer graphene. Instead, using a simple model to account for the
in-plane atomic reconstruction driven by the changing stacking registry, we
show that the resulting expanded Bernal-stacked domains and sharper interfaces
are required for robust topological interfaces to emerge. These states are
highly localised in the AA- or SP-stacked interface regions and exhibit
differences in their layer and sublattice distribution depending on the
interface stacking. We conclude that heterostrain can be used as a mechanism
to tune the presence and distribution of topological channels in gapped
bilayer graphene systems, complementary to the field of twistronics.

††preprint: APS/123-QED

##  I Introduction

Graphene was the first 2D material to be experimentally isolated, and
monolayer graphene (MLG) remains the most well-known 2D material due to a
range of unique mechanical and electronic properties [1, 2, 3, 4]. MLG has
gapless, linear bands which form two Dirac cones, or valleys, near the Fermi
energy [5, 6]. The honeycomb lattice of MLG can be divided into two triangular
sublattices (A and B), which give rise to the pseudospin degree of freedom
underpinning exotic phenomena such as Klein tunneling [7]. Due to its high
elasticity and flexibility, it has also been suggested for potential
applications in flexible electronic devices [8], including touch screens [9]
and foldable organic light emitting diodes [10]. It is also the building block
for more complex multilayer structures [11], such as bilayer graphene (BLG),
which is formed from two stacked graphene layers coupled by weak van der Waals
interactions [12]. It can also be used in heterostructures, for example by
combining graphene with transition metal dichalcogenides [13].

BLG systems can be characterised by the relative stacking between their
layers. In the most favourable stacking, the carbon atoms from one sublattice
in each layer lie directly opposite the centre of a hexagon in the other
layer, whereas the atoms from the remaining sublattice in each layer lie
directly opposite each other [12, 14], forming so-called ‘dimers’. There are
two equivalent possibilities, AB and BA stacking, depending on which
sublattices from each layer form dimers, as shown in Fig. 1(a) and (b). The
least favourable stacking possibility is AA-stacking, where the two layers are
perfectly aligned and each atom sits opposite an atom from the same sublattice
in the opposite layer, as in Fig. 1(c) and (e). Between the low-energy AB and
high-energy AA extremes lie a range of stacking options which can be achieved
by varying the amount and directions by which one layer is shifted relative to
the other. The sublattice symmetric SP (or saddle point) stacking, shown in
Fig. 1(d) is of interest in this work as it occurs for shifts half way between
AB and BA stackings, and corresponds to a local energy maximum along this
direction.

The electronic properties of a BLG structure depend sensitively on the
stacking sequence, as shown in Fig. 1(f)–(h). For example, AB or BA-stacked
BLGs have parabolic electronic bands, whereas the bands of AA-stacked BLGs
remain linear, as in MLG [15, 1], but the two resultant cones are split in
energy. SP-stacked BLGs have also linear bands, but the resultant cones now
have relative shifts in both energy and momentum [16]. These differently
stacked systems also behave very differently under the application of an
interlayer bias. A band gap proportional to the bias is opened for AB/BA
stackings [17], while AA- and SP-stacked systems remain metallic [18].

The gaps opened by an interlayer bias in AB- and BA-stacked BLGs are equal in
magnitude, but are topologically non-equivalent [19, 20]. Interfaces between
AB and BA-stacked regions in a single system are therefore of interest as they
can host topologically-protected, valley-polarised transport channels. These
appear in the band structure of biased BLG systems as pairs of bands (in each
valley) which connect the otherwise gapped Dirac cones [19, 21, 22].
Interfaces between AB- and BA-stacked regions, and the associated topological
channels, can occur in a number of systems, with the simplest example being a
change in stacking registry caused by a grain boundary in one of the layers
[23, 24]. A similar effect occurs if the sign of the bias, instead of the
stacking registry, changes across the interface [25, 26, 27, 28]. AB/BA
interfaces also occur when a relative twist angle is introduced between the
two layers, to form twisted bilayer graphene (tBLG). This creates a Moiré
pattern where the stacking order has a periodic modulation, with alternating
domains of AA, AB/BA and SP stackings [29, 30]. The interfaces between AB and
BA domains form a triangular superlattice with AA-stacked vertices connected
by SP-stacked edges. The formation of localised, topological interface states
along the SP edges when a bias is applied is a topic of huge interest
recently, both theoretically [20, 31, 32, 33, 34] and experimentally [35, 36,
37, 38]. The localisation of electronic states in different domains can lead
to correlated insulating states and unconventional superconductivity for
specific twist angles [39, 40, 41].

Figure 1: Bilayer graphene with a 4% uniform heterostrain, applied in the AC-
direction to the bottom (orange) layer (L1). A one-dimensional Moiré pattern
emerges due to different stacking registries between the two layers, shown by
the zoomed regions in panels (a)–(e). The band structures for BLG systems,
with and without an interlayer bias, for uniform (f) AB-, (g) AA- and (h) SP-
stacking registries are shown in the bottom panels.

Strain can also break the stacking registry and give rise to different
stacking domains and interfaces between them. A heterostrain, i.e. a
difference between the strains felt by each layer, introduces a mismatch
between the lattice constants in each layer. Local compression, stretching or
slipping of one layer can lead to the formation of strain solitons or more
complex topological defects [42, 43]. A pure uniaxial heterostrain, arising
from a uniaxial strain present in only one of the layers, leads to a one-
dimensional (1D) Moiré pattern [44, 45, 46, 47]. An example is shown in the
main panel of Fig. 1, where a 4% armchair-direction uniaxial tensile strain is
present in the bottom (orange) layer, with the top (blue) layer remaining
unstrained. The resulting lattice mismatch in this direction gives rise to an
alternating, repeating sequence of (from bottom to top in Fig. 1) AA, AB, SP,
and BA stackings. Such a situation can arise, for example, if the strain
applied to a bilayer system using a flexible substrate does not transfer
between the layer directly in contact with the substrate to the top layer [48,
49, 50]. In a recent work [46], we demonstrated that such a situation can
arise if the energetic penalty caused by the broken AB-stacking is less than
the energy cost required to strain the top layer by the same amount. For
applied heterostrains of approximately 1% or greater, we found that modulated
stacking case is energetically favourable, which suggests that such strains
can be used as a tool to tune the stacking of a BLG system.

In this work, we focus on the formation and evolution of interface states in
heterostrained BLGs. The 1D symmetry of the system gives rise to two types of
interface between AB and BA regions, namely those centred around the SP-
stacked and the AA-stacked regions, as shown in Fig. 1. In Sec. II, we
introduce the tight-binding methods use to calculate the electronic properties
of general BLG systems and a simple model to account for in-plane relaxations
and changes in domain size and interface sharpness. In Sec. III, we examine
the band structure and distribution of states in BLGs with a uniform
heterostrain, and compare the results to what is seen for uniformly-stacked
AA, AB and SP systems. Realistic values of bias and strain do not give clear
signatures of the expected topological interface channels, which we attribute
to the relatively wide interface regions and the relatively small AB and BA
domains that occur for a uniform heterostrain. Following this, in Sec. IV, we
use the simple model to mimic the effects of in-plane relaxation and the
resulting variances in interface width. This is motivated by previous studies
in twisted systems, which find that relaxation tends to minimise AA-stacked
domains, maximise AB/BA-stacked domains and sharpen SP-stacked interfaces [42,
51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 32, 61, 33, 62, 63]. These structural
deformations can affect the electronic band structure and, in particular, the
formation of flat bands at specific twist angles. In our heterostrained
systems, one layer is assumed to remain uniform due to its interaction with a
substrate. If strain is applied using a flexible substrate, the in-plane
relaxation should occur in the unstrained layer [64, 65, 66, 48]. However, if
strain is applied to the top layer using, for example, the tip of an atomic
force microscope [67], then we expect in-plane relaxation within this layer.
We therefore consider in-plane relaxations of both the strained and unstrained
layers independently in order to capture the geometries that emerge from both
processes. We find that in-plane relaxation, of either the strained or
unstrained layer, allows for more robust topological interfaces to emerge due
to larger, gapped AB/BA domains and sharper SP and AA interfaces. After
discussing the structural, energetic and band-structure considerations, we
examine the resulting topological interfaces in more detail in Sec. V.
Finally, we discuss our findings in the context of recent experimental works
and highlight the potential role of strain as a mechanism to tune the presence
and distribution of topological channels in BLG systems.

##  II Methods

###  II.1 Electronic structure

To model the electronic structure of a general multilayer graphene system with
a range of different stackings, strains, inter-layer biases and/or twists, a
tight-binding (TB) Hamiltonian of the form

$\hat{H}=\sum_{i}\epsilon_{i}\hat{c}^{\dagger}_{i}\hat{c}_{i}-t(\vec{d})\sum_{\langle
i,j\rangle}\big{(}\hat{c}^{\dagger}_{i}\hat{c}_{j}+\hat{c}^{\dagger}_{j}\hat{c}_{i}\big{)}$
(1)

can be used. Here $\hat{c}^{\dagger}_{i}$ and $\hat{c}$ are, respectively, the
creation and annihilation operators for an electron in the $p_{z}$ orbital at
site $i$. The effect of a simple interlayer bias $\Delta$ can be included by
setting the onsite parameters $\epsilon_{i}=\pm\frac{\Delta}{2}$, with the
choice of sign determined by the layer in which site $i$ is located. The
hopping parameter $t(\vec{d})$ depends on the distance vector
$\vec{d}=\vec{R}_{i}-\vec{R}_{j}$ between the two sites and is written in
terms of the appropriate Slater-Koster parameters [68]) as

$-t(\vec{d})=V_{pp\pi}(d)\Big{[}1-\Big{(}\dfrac{\vec{d}\cdot\vec{e}_{z}}{d}\Big{)}^{2}\Big{]}+V_{pp\sigma}(d)\Big{(}\dfrac{\vec{d}\cdot\vec{e}_{z}}{d}\Big{)}^{2}$
(2) $\displaystyle V_{pp\pi}(d)$
$\displaystyle=V^{0}_{pp\pi}e^{-\tfrac{d-a_{0}}{r_{0}}}$ (3) $\displaystyle
V_{pp\sigma}(d)$ $\displaystyle=V^{0}_{pp\sigma}e^{-\tfrac{d-d_{0}}{r_{0}}}.$
(4)

Here $V^{0}_{pp\pi}=-2.7$ eV and $V^{0}_{pp\sigma}=0.48$ eV are, respectively,
the nearest-neighbour in-plane and out-of-plane (dimer) transfer integrals for
unstrained graphene systems, $a_{0}=1.42$ Å and $d_{0}=3.35$ Å are the
corresponding unstrained in-plane and out-of-plane atomic separations,
$r_{0}=0.453$ Å is a decay length, and $\vec{e}_{z}$ is the unit vector in the
out-of-plane $z$ direction. We include all hopping terms whose corresponding
separation in the $xy$-plane is less than a cutoff of $r_{c}=1.82$ Å. Within a
single layer, this is equivalent to a nearest-neighbour TB model that can
account for local strains up to $\epsilon\sim 0.25$. This value is also
sufficient to account for the dimer ($\gamma_{1}$) and skew
($\gamma_{3},\gamma_{4}$) interlayer hoppings, up to similar levels of
deformation, in BLG systems. This type of TB model is suitable for systems
with inhomogeneous strains and stackings, as it can account for differing bond
lengths in both the in-plane and out-of-plane directions, and has previously
been applied in studies of twisted bilayers [69, 55, 33, 63].

The main geometry considered in this work is a BLG sheet with a uniaxial
heterostrain of 1% applied along the armchair direction in one of the layers
(L1 – the ‘bottom’ layer), with the other layer (L2 – the ‘top’ layer)
unstrained. As discussed in the introduction, the layer which is actually
strained in experiment will depend on the method used to apply the strain. In
the absence of a Poisson contraction perpendicular to the applied strain, this
creates a regular one-dimensional Moiré pattern. This is shown for an
exaggerated strain of 4% in Fig. 1. As the stacking modulation does not affect
periodicity in the zigzag direction, the periodic unit cell required for
calculations of this system consists of a single chain of atoms from each
layer. We assume an initial AA-stacking, which remains unchanged at $y=0$ and
therefore the stacking must return to AA at the top of the unit cell in order
to maintain periodicity in the $y$-direction. This requires a total extension
which is a multiple of $y_{u}=3a_{0}$, the width of a four-atom graphene unit
cell. If L1 contains $N$ strained cells and has a total extension of $My_{u}$,
then L2 contains $N+M$ unstrained cells. This restricts calculations to
rational strain values $\varepsilon=\frac{M}{N}$ and requires larger unit
cells to simulate smaller strains. The unit cell for the 1% heterostrain case
discussed throughout this work corresponds to $M=1$, $N=100$ and contains 804
carbon atoms. The band structures presented in the following sections are
calculated using exact diagonalization of the corresponding Hamiltonian
matrices with the appropriate Bloch phases. As we are largely concerned with
the emergence of states confined to $x$-direction interfaces, we only show
bands as a function of $k_{x}$, with $k_{y}=0$.

###  II.2 Simple model of relaxation

In-plane relaxation has been shown to play a vital role in determining the
electronic properties of twisted Moiré systems [55, 58, 32, 59, 33, 62], with
the main effect of such relaxation an increase in the size of AB and BA
domains at the expense of AA and SP regions. In the presence of an interlayer
bias, this affects the localization of states along the sharper SP interfaces
between gapped regions of opposite mass. To consider how in-plane relaxations
could affect our heterostrained geometry, we employ a simple model that allows
the sharpness of the AA- and SP- interfaces, and hence the effective width of
AB/BA regions, to be continuously tuned. In our model, we assume that only one
layer of the system is allowed to relax. We first consider the case where
relaxation occurs within L1, the strained layer. A uniform $y$-direction
strain $\varepsilon=\frac{1}{N}$ applied to L1 produces an extension which is
linearly-dependent on the unstrained y-coordinate $y_{0,L1}$:

$\Delta y_{L1}\,(y_{0,L1})=\varepsilon\,y_{0,L1}\,.$ (5)

The linear extension and uniform strain, as a function of position, are shown
by the grey curves in Fig. 2(a) and (b) respectively. The resultant lattice
mismatch with L2 (the unstrained layer) gives rise to high-symmetry AA, AB, SP
and BA stackings when $\Delta y_{L1}=0$, $a_{0}$, $\tfrac{3}{2}a_{0}$ and
$2a_{0}$ respectively, which correspond to $y_{0,L1}^{AA_{1}}=0$,
$y_{0,L1}^{AB}=\frac{Ny_{u}}{3}$, $y_{0,L1}^{SP}=\frac{Ny_{u}}{2}$ and
$y_{0,L1}^{BA}=\frac{2Ny_{u}}{3}$. The choice of a rational strain values, as
discussed above, ensures that the stacking returns to AA again at $\Delta
y_{L1}=y_{u}=3a_{0}$ and $y_{0,L1}^{AA_{2}}=Ny_{u}$. (Note that, due to
periodicity, AA1 and AA2 correspond to the same point, but it is useful to
consider them separately in the context of a finite unit cell.) The positions
of these high-symmetry stacking locations are shown by the grey curve in Fig.
2(c) and the local stacking near them by the insets (a)–(e) in Fig. 1.

Figure 2: Comparison of the local (a) displacement, (b) strain and (c)
stacking order in heterostrained bilayer graphene as a function of position
along the strained direction ($y$). The grey curve shows the case of a uniform
uniaxial strain in layer L1, whereas the orange and purple curves correspond
to ‘relaxed’ structures with different values of the interface smoothness
parameter $\alpha$.

To simulate sharper interfaces, and broader gapped regions, we replace the
linear extension from Eq. (5) with a sum of sigmoid functions, $\Delta
y_{L1}^{\mathrm{rel.}}\,(y_{0,L1})=\sum_{i}S_{i}-\tfrac{S_{{AA}_{1}}}{2}$,
each of the form

$S_{i}(y_{0,L1})=\dfrac{\Delta
Y_{i}}{1+e^{-(y_{0,L1}-y_{0,L1}^{i})/(\alpha\Delta Y_{i})}}\,,$ (6)

with $i\in\\{AA_{1},SP,AA_{2}\\}$. These ‘S’-shaped functions are centred at
the AA and SP interfaces (see orange and purple curves in Fig. 2(a)), where
they locally increase the strain (panel (b)) and sharpen the interfaces (panel
(c)). This leads to much lower strains, and more uniform stackings, in the AB
and BA regions. The ‘smoothness’ or ‘sharpness’ of an interface is determined
by the parameter $\alpha$, which also sets the local maximum value of strain
at the interface. $\alpha$ relates the characteristic width of an interface to
the total stacking shift ($\Delta Y_{i}$) that occurs across it. We note that
the AA-interface has twice the interlayer shift of the SP-interface ($\Delta
Y_{AA}=2a_{0}$, $\Delta Y_{SP}=a_{0}$), but that $\alpha$ acts as a single
adjustable parameter to allow relaxation, with similar levels of strain, in
both interfaces. Low values of $\alpha$ give sharp interfaces, with more
uniform strain returning as $\alpha$ is increased. However, the model assumes
that the effects of each interface are independent and breaks down if $\alpha$
is increased to values where the interface regions begin to merge.

The above approach can also be used to consider relaxation occurring instead
in the initially unstrained layer (L2), as may be expected to occur if strain
is applied to the other layer via a flexible substrate. In this case, a
uniform tensile strain is first applied to L2 to match that in the strained
layer, before a non-uniform compressive strain is applied using the ‘sum-of-
sigmoids’ approximation above to return the average strain in the unstrained
layer to zero. The structural and electronic effects of both types of
relaxation are discussed in detail in Sec. IV.

##  III Electronic properties of heterostrained systems

We begin by considering the electronic structure of BLG with a 1% strain
applied to one of the layers, in the absence of either interlayer potentials
or relaxation effects. The 804-atom ($43$ nm wide) unit cell of this system is
four times wider than the schematic structure shown in Fig. 1, with a minimum
separation of approximately 14 nm between the centres of the AB and BA
domains. The electronic structure of this system along the $k_{x}$ direction
is shown in Fig. 3(a) and reveals a complex series of subbands due to the
large real-space unit cell and associated band-folding in reciprocal space.
The parabolic bands and semi-metallic, zero-bandgap behaviour expected for AB-
stacked BLG, shown schematically in Fig. 1(f), cannot be clearly distinguished
here. Instead, a large number of bands cross the Fermi energy near $E=0$ and
give the system a metallic character. This suggests that modulated stacking
introduces features from the band structures of uniformly AA (Fig. 1 (g)) and
SP-stacked (Fig. 1 (h)) BLGs. Unlike AB-stacking, these cases resemble the
conical band structure of MLG, but with two copies of the linear MLG cone
separated in either energy (AA), or energy and momentum (SP).

Figure 3: (a) $K$-valley band structure of a BLG system with 1% uniform
heterostrain. Panels (b), (c), (d) show the band structure projected onto 2.0
nm wide strips around the perfect AB-, AA- and SP-stacking points,
respectively. Most of the states are distributed across multiple regions, with
exceptions for low-energy, high-$k$ states which show confinement to regions
with particular stackings.

We now explore the relationship between the band structure of the
heterostrained structure and those of regularly-stacked systems in more
detail. Fig. 3(b–d) show the bands from Fig. 3(a) projected onto narrow 2.0 nm
strips centred around $y^{AB}$, $y^{AA}$, and $y^{SP}$, with the colour of
each point showing the weight of the associated heterostrained bilayer state
in these regions. These plots reveal the effective band structures in the AB,
AA and SP regions of the heterostrained BLG. Bands which are only strongly
coloured in one region indicate states that are largely localised in that
region. Similarly, a strong correspondence between the bands in each region
and their bulk counterparts in Fig. 1(f)–(h) would indicate that bands in a
region with a particular stacking resemble those in a uniform bilayer with the
same stacking. The majority of the heterostrained systems states, and
particularly the higher-order subbands, are not localised in this manner.
Instead they have weight in multiple regions, indicating states that are
distributed across large parts of the system and not confined to regions with
a particular stacking. These states reside in regions of $k_{x}$–$E$ common to
all three stackings, and can be viewed as hybridisations between states with
similar momentum and energy that occur in regions with different stackings.
The only significant exceptions to this trend is for states outside the region
of $k_{x}$–$E$ space spanned by the bands of AB-stacked BLG, e.g., larger $k$
values at low energies. These states are highly quenched in AB regions,
preventing the hybridization of states in different regions, and thus leading
to bands that are highly-localised in either AA or SP regions. Although they
have a similar spatial distribution to that expected for topological interface
states, we note that there is no gap in this system, and that these states
emerge instead from simple confinement effects.

Figure 4: Similar to Fig. 3, but with an additional interlayer bias of
$\Delta=-200$ meV between the layers. No clear band gap opens, even for states
in the AB-stacked region, despite a sizeable band gap opening for the
corresponding uniformly AB-stacked system in Fig. 1(f).

An interlayer bias opens a band gap only in BLG systems with certain
stackings, with the largest gap occurring for AB-stacking and no gap opened
for AA- or SP-stacked systems. This leads to a distribution of gapped and
conducting domains, and the formation of topological interface states, in
systems with modulated stacking, such as tBLG [20]. Similar behaviour should
emerge in heterostrained structures and, in particular, the AA and SP strips
separating AB and BA regions should host one-dimensional topologically
protected channels. Fig. 4 shows the band structure and projections of the
system from Fig. 3 when an interlayer bias of $-200$ meV is applied. We note
that here the negative potential is applied to the strained layer – a slightly
different result is obtained for the opposite case due to the strain-induced
layer asymmetry and the broken electron-hole symmetry resulting from longer-
ranged hopping terms in the Hamiltonian. This bias is sufficient to open a
clear gap for a uniform AB-stacked structure (Fig. 1(f)), but has little
qualitative effect on the gapless bands of AA- and SP-stacked systems (Fig.
1(g,h)). We note that, unlike other one-dimensional interface cases such as
grain boundaries or sharp bias flips [23, 24, 25, 26, 27, 28], the application
of an interlayer bias in Fig. 4 does not give rise to clearly defined
topological channels. In these other systems, an interlayer potential opens a
bulk band gap which is bridged only by pairs of chiral boundary modes with
opposite propagation directions in the K and K′ valleys [19]. For our system,
this would correspond to a total of four valley-protected topological modes in
each valley, two each for the interfaces along the AA and SP domains. Instead
of the clear emergence of a bulk gap and topological modes, we instead see
only minor changes compared to the unbiased system discussed in Fig. 3.

##  IV Role of interface smoothness

Figure 5: (a) Evolution of the band structure of the $\Delta=-200$ meV
heterostrained BLG system as the AB and SP interfaces are sharpened by
inhomogeneously relaxing the strained layer. The ‘pseudogap’ window, where
only four band crossings are present, is show in red. (b) Pseudogap magnitude
for different values of interlayer bias ($\Delta$) and interface smoothness
($\alpha$).

The absence of gap formation and topological modes in the previous section can
be explained by the continuously varying stacking order in the system
considered. This results in only very small regions with perfect AB or BA
stackings and prevents the formation of a band gap. However, the continuously
varying stacking order is an artefact of the homogeneous strain applied to L1,
and is not necessarily representative of realistic systems where non-uniform
strain profiles are adopted to minimise the total energy of the system [42,
51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 32, 61, 33, 62, 63]. We now consider
the simple relaxation model from Sec. II.2 to determine how in-plane
relaxation could affect the electronic properties of a biased, heterostrained
BLG system. This allows us to continuously adjust the interface smoothness and
redistribute the strain from the AB/BA-stacked regions to the SP and AA
interfaces. Fig. 5(a) shows the evolution of the band structure of the biased
system from the previous section as the interface smoothness $\alpha$ is
decreased. As the interfaces become sharper, the number of bands crossing the
Fermi energy decreases until only two pairs of crossings remain. This is the
expected number of band crossings for a BLG system with two topological
interfaces between AB and BA-type regions, i.e., with the formation of
counter-propagating valley-polarised interfaces states along both the AA and
SP strips in our system [19]. The energy window in which only these four bands
are present, which we refer to as the _pseudogap_ , is shown by the red shaded
areas for $\alpha=7$ and $\alpha=3$ in Fig. 5(a), and increases as the
interfaces become sharper. The pseudogap only opens below a critical value of
$\alpha=11$, as shown by the red dashed line, and is not present for unrelaxed
or smoother interfaces. This evolution of the pseudogap is shown explictly by
the orange curve in Fig. 5(b), together with corresponding results for smaller
and larger values of the interlayer bias. The critical smoothness required to
open a pseudogap and observe interface states varies with the magnitude of the
bias, with smaller biases requiring a sharper interface to open a pseudogap.
This is consistent with the different stacking profiles in these systems.
Sharper interfaces give rise to larger regions with nearly perfect AB or BA
stacking where a smaller bias is able to effectively quench bulk states.
Although precise AA or SP stackings still occur at the sharper interfaces,
they do not persist over wide enough areas for the gapless bulk-like bands
associated with uniform AA- or SP-stacked systems to emerge.

The degree to which intralayer relaxation will occur can be simulated in a
number of ways, including molecular dynamics [51, 32, 54, 63] or ab-initio
simulations [61], and the Frenkel Kontorova model and related continuum
approaches [42, 55, 57, 33]. The specifics of the relaxation are also likely
to be affected by experimental conditions, such as substrate effects and the
manner in which strain is applied. While the specifics may vary, relaxation
occurs principally to minimise the energy costs introduced by straining the
graphene lattice and breaking uniform AB-stacking registry. Both contributions
can be calculated on a per-atom basis from ab-initio simulations of uniform
systems [46], allowing the energetics of larger, non-uniform systems to be
estimated by summing over local stacking and strain costs throughout the
system. Inspired by the possible experimental processes, we consider two
possible relaxation scenarios in Fig. 6: either (a) the heterostrained layer
or (b) the unstrained layer relaxes under the simple model in Sec. II.2. In
each case, the atomic positions in the other layer are unaffected by the
relaxation.

Figure 6: Energy costs (per supercell) associated with breaking AB-stacking
registry (red curves) and strain (blue curves) in heterostrained systems. The
panels show the effect of the simple relaxation model applied to either the
(a) strained or (b) unstrained layer to give sharper interfaces. The black
curves show the total energy from both terms, and the square symbols show the
corresponding values for uniformly strained (unrelaxed) systems.

Fig. 6 shows the energy costs (per supercell) due to breaking uniform AB-
stacking (red curves) and due to strain (blue curves) in heterostrained
systems with the relaxation model applied. The total energy cost is shown by
the black curve, with the corresponding energy values for a uniformly
heterostrained system (without relaxation) shown by the square symbols. The
stacking profiles and associated energy costs are identical if either the
strained (Fig. 6(a)) or unstrained layer (Fig. 6(b)) are relaxed, with the
energy costs decreasing linearly for sharper interfaces (smaller $\alpha$). In
contrast, the strain energy cost is always higher if we chose to relax the
initially unstrained layer. In this case, large portions of both layers
require strain values near $1\%$ in order to maintain the stacking order in
the larger AB/BA regions found for smaller $\alpha$. For both relaxation
scenarios, and for all values of $\alpha$, the strain energy cost is greater
than that of the uniformly-strained system (blue squares). This is trivial if
the unstrained layer is relaxed, as the cost of the additional strain
introduced by the relaxation of this layer must be added to the cost of the
uniformly strained layer. When the strained layer is relaxed, the strain is
redistributed throughout the layer such that the average strain in the layer
is conserved. As the energy cost of strain increases superlinearly [46], the
additional energy cost of higher strain near the interfaces is not compensated
by the reduced strain in the AB and BA domains.

Regardless of which layer is relaxed, the total energy curve is relatively
flat over a wide range of interface smoothness values, until it increases
sharply for $\alpha\lesssim 8$. This is due to the significant strain energy
cost for sharper interfaces, which becomes much larger than the energy savings
due to more preferable stacking configurations. There is a slight reduction in
the total energy, compared to the uniformly strained case (black square
symbol), if the strained (unstrained) layer is relaxed using $\alpha\gtrsim 8$
($10$). This is near the cut-off where a pseudogap, and associated topological
interface states, emerge for smaller values of bias, as shown in Fig. 5(b).
However, a reasonable energy window of order $\sim 100$ meV in which to
observe topological states should occur for a larger interlayer bias of
$\Delta=-400$ meV.

Before examining interface states in more detail, it is worth considering
additional factors which may affect strain redistribution in heterostrained
bilayer graphene. The relaxation model employed here is a simplification that
allowed us to examine, in a general manner, the role of interface sharpness in
such systems. Substrate effects are neglected, except under the simplifying
assumption that relaxation only occurs in one layer of the system. The
specific nature of the interaction between the substrate, or any tips or
contacts used for either the application of strain or other measurements, may
also effect the local distribution of strain in the system. Similarly, we have
not considered the role which out-of-plane deformations could play in reducing
the total energy of heterostrained systems. Indeed, a recent work shows that
such deformations can affect flat band formation and local chemical reactivity
in this type of system [47]. The local stacking energy will also be sensitive
to modulations of the interlayer separation introduced by such deformations
[56]. Finally, we note that for simplicity we have only considered a single
value of $\alpha$ within each structure so that the local strain at the AA and
SP interfaces is the same. However, these interfaces have different local
stacking energetics, relative shifts and widths, and so may be able to sustain
different levels of local strain. We do not expect these additional
considerations to dramatically alter the main trends in our results. The main
effect is likely to be that these more complicated relaxed geometries could
reduce the energy costs due to strain beyond the values predicted by our
simple model, leading to an increased preference for sharper interfaces. To
conclude this section, we re-emphasise that uniform heterostrain and bias
alone are not sufficient to create one-dimensional topological channels: an
accompanying relaxation which sharpens the interfaces and expands the Bernal-
stacked regions is also required.

##  V Interface states

Figure 7: (a) Band structure of a heterostrained BLG system with an interlayer
bias $\Delta=-400$ meV, where the interfaces are sharpened by a relaxation
with $\alpha=9$ in the strained layer. (b), (c), (d) show the band structure
projected onto AB, AA and SP regions of the system, as in Figs. 3 and 4. Clear
topological channels are formed at the $AA$ and $SP$ interfaces and are marked
in panels (c) and (d). The distributions of the these states around their
respective interfaces are shown in panels (e)–(g).

Fig. 7(a) shows the band structure near the K point of a heterostrained BLG
system with the relaxation model applied to the strained layer, with
$\alpha=9$ and $\Delta=-400$ meV. From Figs. 5(b) and 6(a), these parameters
are consistent with a sizeable pseudogap and reasonable energetic costs. A
similar pseudogap could also be achieved for a lower bias, as shown for
$\Delta=-200$ meV in Fig. 5(a), but at the cost of a much higher energetic
cost due to larger local strains near a sharper interface. Conversely,
achieving a similar pseudogap with a smoother interface would require a larger
bias to be applied between the layers. The band structure in Fig. 7(a) shows a
pseudogap spanned by four dispersive bands. That these indeed correspond to
interface-localised channels is made clear by the band projections in panels
(b)–(d). The AB regions are gapped, whereas the AA- and SP- interfaces each
host a pair of boundary modes in the K valley. We note that the modes at a
particular interface propagate in only a single direction for this valley,
while the corresponding modes in the $K^{\prime}$ valley (not shown) propagate
in the opposite direction. We further note that the $K$ valley modes from the
two different interfaces have opposite propagation directions. This is due to
the opposite placements of the AB and BA regions relative to the two
interfaces. All these features are consistent with the formation of valley-
protected topological modes at the interfaces due to a change in the valley
Chern number between AB and BA regions in biased BLG [19].

Finally, we note that there are subtle differences between modes at AA- or SP-
interfaces. Firstly, the energy window for observing these states is slightly
different, due to the onset of higher-energy bands at different energies in
the AA and SP regions. We note that AA interface in our system is broader than
the SP interface, as discussed in Sec. II.2. This does not seem to correspond
to a wider interface state however. Fig. 7(e)–(h) show the projections of the
corresponding states in panels (c) and (d) onto each atom in a $20$ Å strip
around the relevant interface. A slightly broader spread is noted for the SP
interface modes in (f), (g) compared to the AA modes in panels (e), (h). There
is a more significant difference however in the layer (colour) and sublattice
(symbol) distribution of these modes. Each of the two AA modes (e, h) has
largest weight exactly at the interface on one of the two layers, with both
sublattices having very similar distributions in the dominant layer. The other
layer has a smaller, sublattice-split distribution with peaks to either side
of the interface. In contrast, the SP modes (f, g) appear to have a more
symmetric distribution between both layers and sublattices. We note that these
trends are likely to be sensitive to the strain levels at the interfaces, and
may not be present if, for example, relaxation results in different local
strains at each interface.

##  VI Conclusions

Twisted bilayer graphene systems and the associated Moiré superlattice have
come under intense scrutiny recently due to the emergence of flat electronic
bands, strongly correlated states, and unconventional superconductivity in
this system. However, twisting is just one way to generate a Moiré
superlattice. In this study we use tight-binding methods to study the
formation and evolution of interface states in heterostrained bilayer
graphene. Uniaxial heterostrain in this system gives rise to a one-dimensional
Moiré pattern with two interface types between AB and BA regions, namely those
with an SP-stacking and those with an AA-stacking. Although uniform
heterostrain creates domain boundaries between AB- and BA- regions, as in
twisted bilayer graphene (tBLG), they are too wide to show clear signatures of
topological interface channels. We demonstrate that in-plane relaxations,
which maximise AB/BA-stacked domains and sharpen AA- and SP-stacked
interfaces, are necessary for robust one-dimensional topologically protected
channels to emerge.

By comparing the electronic band structure of a biased heterostrained bilayer
to those of regularly-stacked systems, we identify bands that are highly-
localised in AA- or SP-stacked interfaces. These bands emerge in the pseudogap
which opens in the AB-stacked regions with the application of an interlayer
bias. We show that the size of the pseudogap and the emergence of interface
states depends on how sharp the interface is as well as the magnitude of
applied bias. We show that when the AB regions are gapped, the AA- and SP-
interfaces each host a pair of boundary modes in the K valley propagating
unidirectionally. The corresponding modes in the K′ valley propagate in the
opposite direction. This is consistent with the formation of valley-protected
topological modes at the interfaces due to a change in the valley Chern number
between AB and BA regions in biased BLG. Finally, unlike tBLG where only the
SP-like interfaces contribute to the network, the Moiré system generated by
the application of a uniaxial strain has two distinct interface states, namely
those at the AA- and SP-stacked interfaces. These interface states at the AA-
and SP-interfaces are not identical, with differences in the layer and
sublattice distribution of these modes. We expect the exact nature of the
distribution of these modes to depend sensitively on the strain levels at
these interfaces.

So far we have considered a pure uniaxial strain which creates a strictly one-
dimensional modulation of the stacking order. Allowing for a biaxial strain,
or a Poisson compression of the strained layer, produces a two-dimensional
Moiré pattern [70]. We are now investigating how this will modify the emergent
network of interface states. We expect that the channels associated with the
AA-stacked interface will no longer contribute to the network, as in tBLG,
while those associated with the SP-regions will remain. However, their
configuration is likely to be different to those in tBLG.

Heterostrained, untwisted bilayer graphene, as considered here, is a very
promising platform from which to tune the presence and distribution of
topological channels. This system could potentially exhibit phenomena akin to
those in tBLG while circumventing some of the limitations associated with
achieving precise twist angle control [67].

## References

  * Novoselov  _et al._ [2004] K. S. Novoselov, A. K. Geim, S. V. Morozov, D.-E. Jiang, Y. Zhang, S. V. Dubonos, I. V. Grigorieva, and A. A. Firsov, Science 306, 666 (2004).
  * Novoselov  _et al._ [2005] K. S. Novoselov, D. Jiang, F. Schedin, T. Booth, V. Khotkevich, S. Morozov, and A. K. Geim, Proceedings of the National Academy of Sciences 102, 10451 (2005).
  * Katsnelson [2007] M. I. Katsnelson, Materials today 10, 20 (2007).
  * Lee  _et al._ [2008] C. Lee, X. Wei, J. W. Kysar, and J. Hone, science 321, 385 (2008).
  * Wallace [1947] P. R. Wallace, Physical review 71, 622 (1947).
  * Castro Neto  _et al._ [2009] A. H. Castro Neto, F. Guinea, N. M. R. Peres, K. S. Novoselov, and A. K. Geim, Rev. Mod. Phys. 81, 109 (2009).
  * Katsnelson  _et al._ [2006] M. I. Katsnelson, K. S. Novoselov, and A. K. Geim, Nature physics 2, 620 (2006).
  * Novoselov  _et al._ [2012] K. S. Novoselov, V. Fal, L. Colombo, P. Gellert, M. Schwab, K. Kim, _et al._ , Nature 490, 192 (2012).
  * Bae  _et al._ [2010] S. Bae, H. Kim, Y. Lee, X. Xu, J.-S. Park, Y. Zheng, J. Balakrishnan, T. Lei, H. R. Kim, Y. I. Song, _et al._ , Nature nanotechnology 5, 574 (2010).
  * Han  _et al._ [2012] T.-H. Han, Y. Lee, M.-R. Choi, S.-H. Woo, S.-H. Bae, B. H. Hong, J.-H. Ahn, and T.-W. Lee, Nature Photonics 6, 105 (2012).
  * Geim and Grigorieva [2013] A. K. Geim and I. V. Grigorieva, Nature 499, 419 (2013).
  * McCann and Koshino [2013] E. McCann and M. Koshino, Reports on Progress in physics 76, 056503 (2013).
  * Azadmanjiri  _et al._ [2020] J. Azadmanjiri, V. K. Srivastava, P. Kumar, Z. Sofer, J. Min, and J. Gong, Applied Materials Today 19, 100600 (2020).
  * Rozhkov  _et al._ [2016] A. V. Rozhkov, A. Sboychakov, A. Rakhmanov, and F. Nori, Physics Reports 648, 1 (2016).
  * Ho  _et al._ [2006] J.-H. Ho, C. Lu, C. Hwang, C. Chang, and M.-F. Lin, Physical Review B 74, 085406 (2006).
  * Park  _et al._ [2015] C. Park, J. Ryou, S. Hong, B. G. Sumpter, G. Kim, and M. Yoon, Physical review letters 115, 015502 (2015).
  * Castro  _et al._ [2007] E. V. Castro, K. Novoselov, S. Morozov, N. Peres, J. L. Dos Santos, J. Nilsson, F. Guinea, A. Geim, and A. Castro Neto, Physical review letters 99, 216802 (2007).
  * Silva  _et al._ [2020] E. Silva, M. Santos, J. Skelton, T. Yang, T. Santos, S. Parker, and A. Walsh, Materials Today: Proceedings 20, 373 (2020).
  * Zhang  _et al._ [2013] F. Zhang, A. H. MacDonald, and E. J. Mele, Proceedings of the National Academy of Sciences 110, 10546 (2013).
  * San-Jose and Prada [2013] P. San-Jose and E. Prada, Physical Review B 88, 121408 (2013).
  * Vaezi  _et al._ [2013] A. Vaezi, Y. Liang, D. H. Ngai, L. Yang, and E.-A. Kim, Physical Review X 3, 021018 (2013).
  * Ju  _et al._ [2015] L. Ju, Z. Shi, N. Nair, Y. Lv, C. Jin, J. Velasco, C. Ojeda-Aristizabal, H. A. Bechtel, M. C. Martin, A. Zettl, _et al._ , Nature 520, 650 (2015).
  * Yasaei  _et al._ [2014] P. Yasaei, B. Kumar, R. Hantehzadeh, M. Kayyalha, A. Baskin, N. Repnin, C. Wang, R. F. Klie, Y. P. Chen, P. Král, _et al._ , Nature communications 5, 4911 (2014).
  * Jaskólski  _et al._ [2016] W. Jaskólski, M. Pelc, L. Chico, and A. Ayuela, Nanoscale 8, 6079 (2016).
  * Oostinga  _et al._ [2008] J. B. Oostinga, H. B. Heersche, X. Liu, A. F. Morpurgo, and L. M. Vandersypen, Nature materials 7, 151 (2008).
  * Zhang  _et al._ [2009] Y. Zhang, T.-T. Tang, C. Girit, Z. Hao, M. C. Martin, A. Zettl, M. F. Crommie, Y. R. Shen, and F. Wang, Nature 459, 820 (2009).
  * Li  _et al._ [2016] J. Li, K. Wang, K. J. McFaul, Z. Zern, Y. Ren, K. Watanabe, T. Taniguchi, Z. Qiao, and J. Zhu, Nature nanotechnology 11, 1060 (2016).
  * Rickhaus  _et al._ [2018] P. Rickhaus, J. Wallbank, S. Slizovskiy, R. Pisoni, H. Overweg, Y. Lee, M. Eich, M.-H. Liu, K. Watanabe, T. Taniguchi, _et al._ , Nano letters 18, 6725 (2018).
  * Bistritzer and MacDonald [2011] R. Bistritzer and A. H. MacDonald, Proceedings of the National Academy of Sciences 108, 12233 (2011).
  * Andrei and MacDonald [2020] E. Y. Andrei and A. H. MacDonald, Nature materials 19, 1265 (2020).
  * Efimkin and MacDonald [2018] D. K. Efimkin and A. H. MacDonald, Phys. Rev. B 98, 035404 (2018).
  * Fleischmann  _et al._ [2019] M. Fleischmann, R. Gupta, F. Wullschläger, S. Theil, D. Weckbecker, V. Meded, S. Sharma, B. Meyer, and S. Shallcross, Nano Letters 20, 971 (2019).
  * Tsim  _et al._ [2020] B. Tsim, N. N. Nam, and M. Koshino, Physical Review B 101, 125409 (2020).
  * De Beule  _et al._ [2021] C. De Beule, F. Dominguez, and P. Recher, Physical Review B 104, 195410 (2021).
  * Huang  _et al._ [2018] S. Huang, K. Kim, D. K. Efimkin, T. Lovorn, T. Taniguchi, K. Watanabe, A. H. MacDonald, E. Tutuc, and B. J. LeRoy, Phys. Rev. Lett. 121, 037702 (2018).
  * Kerelsky  _et al._ [2019] A. Kerelsky, L. J. McGilly, D. M. Kennes, L. Xian, M. Yankowitz, S. Chen, K. Watanabe, T. Taniguchi, J. Hone, C. Dean, _et al._ , Nature 572, 95 (2019).
  * Xu  _et al._ [2019] S. Xu, A. Berdyugin, P. Kumaravadivel, F. Guinea, R. Krishna Kumar, D. Bandurin, S. Morozov, W. Kuang, B. Tsim, S. Liu, _et al._ , Nature communications 10, 4008 (2019).
  * Verbakel  _et al._ [2021] J. Verbakel, Q. Yao, K. Sotthewes, and H. Zandvliet, Physical Review B 103, 165134 (2021).
  * Cao  _et al._ [2018a] Y. Cao, V. Fatemi, S. Fang, K. Watanabe, T. Taniguchi, E. Kaxiras, and P. Jarillo-Herrero, Nature 556, 43 (2018a).
  * Cao  _et al._ [2018b] Y. Cao, V. Fatemi, A. Demir, S. Fang, S. L. Tomarken, J. Y. Luo, J. D. Sanchez-Yamagishi, K. Watanabe, T. Taniguchi, E. Kaxiras, _et al._ , Nature 556, 80 (2018b).
  * Stepanov  _et al._ [2020] P. Stepanov, I. Das, X. Lu, A. Fahimniya, K. Watanabe, T. Taniguchi, F. H. Koppens, J. Lischner, L. Levitov, and D. K. Efetov, Nature 583, 375 (2020).
  * Popov  _et al._ [2011] A. M. Popov, I. V. Lebedeva, A. A. Knizhnik, Y. E. Lozovik, and B. V. Potapkin, Physical Review B 84, 045404 (2011).
  * Alden  _et al._ [2013] J. S. Alden, A. W. Tsen, P. Y. Huang, R. Hovden, L. Brown, J. Park, D. A. Muller, and P. L. McEuen, Proceedings of the National Academy of Sciences 110, 11256 (2013).
  * Choi  _et al._ [2010] S.-M. Choi, S.-H. Jhi, and Y.-W. Son, Nano letters 10, 3486 (2010).
  * Van der Donck  _et al._ [2016] M. Van der Donck, C. De Beule, B. Partoens, F. Peeters, and B. Van Duppen, 2D Materials 3, 035015 (2016).
  * Georgoulea  _et al._ [2022] N. C. Georgoulea, S. R. Power, and N. M. Caffrey, Journal of Physics: Condensed Matter 34, 475302 (2022).
  * Schleder  _et al._ [2023] G. R. Schleder, M. Pizzochero, and E. Kaxiras, “One-dimensional moiré physics and chemistry in heterostrained bilayer graphene,” (2023), arXiv:2306.09799 [cond-mat.mtrl-sci] .
  * Androulidakis  _et al._ [2020] C. Androulidakis, E. N. Koukaras, G. Paterakis, G. Trakakis, and C. Galiotis, Nature communications 11, 1595 (2020).
  * Frank  _et al._ [2012] O. Frank, M. Bousa, I. Riaz, R. Jalil, K. S. Novoselov, G. Tsoukleri, J. Parthenios, L. Kavan, K. Papagelis, and C. Galiotis, Nano Letters 12, 687 (2012).
  * Wang  _et al._ [2019] K. Wang, W. Ouyang, W. Cao, M. Ma, and Q. Zheng, Nanoscale 11, 2186 (2019).
  * Lin  _et al._ [2013] J. Lin, W. Fang, W. Zhou, A. R. Lupini, J. C. Idrobo, J. Kong, S. J. Pennycook, and S. T. Pantelides, Nano letters 13, 3262 (2013).
  * Van Wijk  _et al._ [2015] M. Van Wijk, A. Schuring, M. Katsnelson, and A. Fasolino, 2D Materials 2, 034010 (2015).
  * Dai  _et al._ [2016] S. Dai, Y. Xiang, and D. J. Srolovitz, Nano letters 16, 5923 (2016).
  * Jain  _et al._ [2016] S. K. Jain, V. Juričić, and G. T. Barkema, 2D Materials 4, 015018 (2016).
  * Nam and Koshino [2017] N. N. Nam and M. Koshino, Physical Review B 96, 075311 (2017).
  * Gargiulo and Yazyev [2017] F. Gargiulo and O. V. Yazyev, 2D Materials 5, 015019 (2017).
  * Carr  _et al._ [2018] S. Carr, D. Massatt, S. B. Torrisi, P. Cazeaux, M. Luskin, and E. Kaxiras, Physical Review B 98, 224102 (2018).
  * Lin  _et al._ [2018] X. Lin, D. Liu, and D. Tománek, Physical Review B 98, 195432 (2018).
  * Guinea and Walet [2019] F. Guinea and N. R. Walet, Physical Review B 99, 205134 (2019).
  * Yoo  _et al._ [2019] H. Yoo, R. Engelke, S. Carr, S. Fang, K. Zhang, P. Cazeaux, S. H. Sung, R. Hovden, A. W. Tsen, T. Taniguchi, _et al._ , Nature materials 18, 448 (2019).
  * Lucignano  _et al._ [2019] P. Lucignano, D. Alfè, V. Cataudella, D. Ninno, and G. Cantele, Physical Review B 99, 195419 (2019).
  * Nguyen  _et al._ [2021] V. H. Nguyen, D. Paszko, M. Lamparski, B. Van Troeye, V. Meunier, and J.-C. Charlier, 2D Materials 8, 035046 (2021).
  * Leconte  _et al._ [2022] N. Leconte, S. Javvaji, J. An, A. Samudrala, and J. Jung, Phys. Rev. B 106, 115410 (2022).
  * Yu  _et al._ [2008] T. Yu, Z. Ni, C. Du, Y. You, Y. Wang, and Z. Shen, The Journal of Physical Chemistry C 112, 12602 (2008).
  * Huang  _et al._ [2009] M. Huang, H. Yan, C. Chen, D. Song, T. F. Heinz, and J. Hone, Proceedings of the National Academy of Sciences 106, 7304 (2009).
  * Anagnostopoulos  _et al._ [2015] G. Anagnostopoulos, C. Androulidakis, E. N. Koukaras, G. Tsoukleri, I. Polyzos, J. Parthenios, K. Papagelis, and C. Galiotis, ACS applied materials & interfaces 7, 4216 (2015).
  * Kapfer  _et al._ [2022] M. Kapfer, B. S. Jessen, M. E. Eisele, M. Fu, D. R. Danielsen, T. P. Darlington, S. L. Moore, N. R. Finney, A. Marchese, V. Hsieh, _et al._ , arXiv preprint arXiv:2209.10696 (2022).
  * Slater and Koster [1954] J. C. Slater and G. F. Koster, Physical Review 94, 1498 (1954).
  * Trambly de Laissardière  _et al._ [2010] G. Trambly de Laissardière, D. Mayou, and L. Magaud, Nano letters 10, 804 (2010).
  * Parhizkar and Galitski [2022] A. Parhizkar and V. Galitski, Phys. Rev. Res. 4, L022027 (2022).