# Transport and spectral properties of magic angle twisted bilayer graphene
junctions based on local orbital models

M. Alvarado  Departamento de Física Teórica de la Materia Condensada C-V,
Condensed Matter Physics Center (IFIMAC) and Instituto Nicolás Cabrera,
Universidad Autónoma de Madrid, E-28049 Madrid, Spain  A. Levy Yeyati
Departamento de Física Teórica de la Materia Condensada C-V, Condensed Matter
Physics Center (IFIMAC) and Instituto Nicolás Cabrera, Universidad Autónoma de
Madrid, E-28049 Madrid, Spain

(August 27, 2024)

###### Abstract

The electronic properties of junctions defined electrostatically on twisted
bilayer graphene can be addressed theoretically using lattice models. Recent
works have introduced minimal local orbital models to describe twisted bilayer
graphene at the magic angle (MATBLG) with different degrees of approximation
and accounting for fragile topology. In the present work we use a Green’s
function formalism to obtain the spectral and transport properties for MATBLG
junctions based on these models. We introduce different symmetry breaking
perturbations to simulate the effect of interactions and characterize the
topology of the bulk bands by analyzing the corresponding Wilson loops. We
then analyze the spectral properties for different types of edges and in the
case of a domain wall in the sublattice symmetry breaking parameter. We
further consider a three region junction where one could control independently
the central and lateral regions doping level. In the limit where the central
region is fixed at the charge neutrality point and the lateral ones are
heavily doped, the spectral and the two terminal transport properties can be
understood in terms of the hybridization of chiral states along the junctions.
These properties are found to be extremely sensitive to the orientation of the
junctions along the moiré lattice.

##  I Introduction

The ability to tune the electronic properties of Van der Waals
heterostructures of two-dimensional materials by changing their relative angle
has given rise to the new field of "twistronics" [1, 2]. Twisted bilayer
graphene is by now the most prominent example of this field, exhibiting a
large variety of interesting phenomena. Near the so-called magic angle the
flat bands close to charge neutrality arising from the large scale moiré
pattern become extremely sensitive to electronic correlations and
interactions, giving rise to a rich phase diagram including superconducting
and insulating phases as a function of the doping level [3, 4].

The flat bands also exhibit topological properties [5] which manifest for
example in the appearance of anomalous quantum Hall effect [6, 7, 8]. More
recently, the possibility to produce tunable junctions on this material by
local electrostatic control has been demonstrated [9] thus opening the way
towards its application in transport devices.

From the point of view of theory, large efforts have been devoted to
understand the origin of correlation induced phenomena in MATBLG [10, 11, 12,
13, 14, 15, 16, 17, 18, 19]. A basic debate has emerged regarding the
appropriate microscopic modeling. While the continuous model of Bistritzer and
MacDonald [20] is broadly accepted as an accurate starting single-particle
description, real space or lattice models could be more appropriate for the
inclusion of interactions or for the study of transport in MATBLG junctions.
Lattice models suffer, however, of the so-called Wannier obstructions [21, 22,
23, 24, 25, 26] which prevent the description of topological bands in terms of
local orbitals. While simple low energy two bands models were proposed which
avoid this problem by neglecting certain symmetries [27], so called "faithful"
tight-binding (TB) models including higher energy bands are able to account
properly for the fragile topology [28, 29, 30, 31, 32, 33, 34] of the low
energy bands [35, 24]. Other models have been proposed mixing real and
momentum space descriptions [36], or including both valleys in a two-band
description [37].

The aim of the present work is to analyze spectral and transport properties
using lattice models for MATBLG junctions. We would like to address non-
homogeneous situations like the ones depicted in Fig. 1 where either local
spectral properties could be studied using STM techniques or two terminal
transport measurements could be performed. Inspired by the recent experiments
of Ref. [9] we study transport through regions with varying doping levels
which could be tuned by gates as schematically shown in Fig. 1 a).
Additionally, our approach should allow us to predict the local spectral
properties at such junctions or at the edges of TBLG samples, which could be
accessed experimentally as demonstrated in recent STM studies on NbSe2/CrBr3
moiré superlattices [38]. An important aspect of MATBLG which guides our study
is that the large size of its moiré pattern ($>10$ nm) would allow transport
experiments on junctions along well defined directions on the moiré lattice.
As we show in this work, the junction’s orientation plays an important role on
its transport and spectral properties. Eventually our work should serve as a
basis for analyzing transport and spectral properties in situations including
superconducting phases.

Our starting point are the lattice models from Refs [27] and [35] for which we
construct the boundary Green functions (bGFs) using an efficient analytical
method for different type of edges. As it is known, in general there is no
bulk-boundary correspondence for bands with fragile topology [39]. Our method
thus provide a direct way to test the appearance of topological edge states
within these models.

Figure 1: Different experimental setups based on MATBLG analyzed in this work.
a) Two terminal transport through a barrier defined by a backgate. The length
of the barrier region can vary between $\lambda_{b}\approx 300-50$ nm. b) STM
spectroscopy of a Domain Wall due to the interaction of MATBLG with the
substrate which is incorporated in the model through the parameter
$\delta_{S}$ (see text for details). The large size of the moiré pattern
($>10$ nm) in MATBLG would allow to study the effect of the alignement of the
junction or the domain wall along well defined directions on the moiré
lattice.

On the other hand, it is known that the one-electron bands in MATBLG are
greatly distorted near the magic angle by the effect of either electron-
electron interactions and by the interaction with the substrate [11, 14, 40,
41, 13, 12]. To simulate their effect we introduce in the models symmetry
breaking terms describing the coupling to the hBN substrate and possible
charge density wave (CDW) modulations. This allows us to study the effect of
most characteristic non-magnetic perturbations in MATBLG.

Interactions also lead to phases with broken valley and spin symmetries which
manifest in the so-called cascade transitions as a function of the doping
level [42, 7]. In this work we do not aim to make predictions on the system
phase diagram. We rather assume that valley and spin degrees of freedom can be
independently populated and study the resulting spectral and transport
properties for a given configuration.

The rest of the manuscript is organized as follows: in Sec. II, we discuss the
effect of Wannier obstruction and fragile topology of the nearly flat bands of
MATBLG in lattice models. We introduce the models from Refs. [27] and [35],
their relevant symmetries and topology. We implement different non-magnetic
symmetry breaking perturbations that can open topological gaps at the Fermi
energy. Furthermore, using Wilson loops we obtain the band’s Chern number and
study the competition among these perturbations in determining the bands
topology.

In Sec. III, we generalize the Green function formalism introduced in Ref.
[43] to 2D systems described by lattice models where we only assume finite
range hopping amplitudes. To compute the boundary Green functions for these
models we use an efficient procedure based on the Faddev-Leverrier algorithm
[44], to be described in a forthcoming publication [45]. In Sec. IV, we study
the spectral properties at edges of MATBLG starting with the spectral density
at a zig-zag boundary showing edge states at the CNP and at the moiré gap
between the nearly flat bands and the excited ones. Moreover we show the
appearance of chiral edge states at the CNP due to sublattice symmetry
breaking perturbation in the model of Ref. [35] unlike the gaped spectrum for
the perturbed model of Ref. [27]. Then we compute the spectral properties at a
domain wall for the sublattice symmetry breaking perturbation inducing the
appearance of two topological states with the same chirality at opposite
minivalleys. In Sec. V, we study spectral and transport properties in a three
region device where we can modulate the doping levels in those regions,
analyzing the influence of topological states at the boundaries of the central
region.

We finally offer, in Sec. VI, some conclusions summarizing the main results
and give an outlook of possible future work taking advantage of the methods
developed here. Technical details like the explicit description of the tight-
binding models or the recursive GF method to obtain the transport properties
in inhomogeneous devices are included in four appendices.

##  II Lattice models for MATBLG

Close to the magic angle a lattice model based on the $\pi$-orbitals at the
carbon atoms requires more than 10.000 states for each moiré unit cell. The
idea of minimal lattice models for MATBLG is to construct a basis of well-
localized Wannier orbitals purely consisting of flat band states. While this
idea is extremely attractive, simple two-bands lattice models fail to describe
the non trivial topology of the flat bands, which is linked to the fact that
they arise from two unperturbed Dirac cones coming from different graphene
layers but the same graphene valley, thus carrying the same helicity. The
topological character of the bands obstructs the construction of minimal
2-bands models based on Wannier orbitals with the correct topology, although
they can reproduce the bands dispersion [27]. Adding higher energy trivial
bands is necessary to consider all the emergent symmetries of the problem in
the hypothesis of fragile topology [35]. In spite of these limitations we
found it instructive to study first the case of the simplest two band model of
Ref. [27] and then to extend our study to the six band model of Ref. [35] that
accounts for fragile topology. We shall refer to these two models as two-band-
one-valley (2B1V) and six-band-one-valley (6B1V) models respectively.

###  II.1 2B1V model

Due to the fact that the flat bands are isolated from the rest of high energy
bands, this 2-band minimal model pursues to capture the electronic properties
of these bands in a basis of well-localized Wannier orbitals.

The Wannier orbitals in the flat bands have to be centered around the
$A=A_{1}B_{2}$ and $B=B_{1}A_{2}$ sites (see Fig. 2 a) in the emergent
honeycomb moiré lattice, where the sub-index refers to different pristine
graphene layers in TBG. This description reproduces the observed centering of
charge around $AA$ sites in the moiré lattice.

Figure 2: Real space representation and corresponding BZs for the 2B1V and
6B1V models. Panel a) shows the honeycomb lattice and the non-monotonic
variation with distance of the hopping integrals ($t_{1}$, $t_{2}$, $t_{4}$,
$t_{5}$) from a lattice site $B$ for the 2B1V model. b) Triangular (grey
shaded) and rectangular (white) moiré BZ showing the high symmetry points in
each one. The overline in the high symmetry points indicates that they belong
to the moiré BZ. Panels c) and d) correspond to the 6B1V model indicating the
armchair (panel c) and zigzag (panel d) boundaries in real space. The
$p$-orbitals are centered at the unit cell site red dots $(A,B)$ forming a
triangular lattice and the black dots corresponds to the $s$-orbitals placed
in a kagome lattice with three sites per unit cell. The shaded red region in
c) and d) corresponds to the doubling of the triangular primitive cell with
orthogonal lattice vectors producing the folding of the original triangular BZ
into a rectangular one as indicated in panel b).

From the localized Wannier description a tight binding Hamiltonian is derived
by calculating the hopping integrals between these well localized Wannier
orbitals, which behave as $(p_{x},p_{y})$-orbitals residing each one in one
site of the honeycomb lattice. To describe quantitatively the TBG band
dispersion it is necessary to take into account hopping elements between an
increasing number of neighbours due to rather extended character of the
Wannier orbitals representing the flat-bands. To assure protected Dirac nodes
in the flat bands it is necessary to add an additional $\mathbb{Z}_{2}$
symmetry corresponding to a twofold rotation with a change of valley.

This model has been intensively used [46, 47, 48, 49] due to its simplicity
but, although it exhibits naturally particle/hole symmetry, due to the Wannier
obstruction it lacks the characteristic symmetries of the system and leads to
moiré minivalleys with vanishing net chirality.

As it was pointed out by Ref. [35] this model does not include essential
symmetries of the continuum model like $C_{2}\mathcal{T}$, where $\mathcal{T}$
accounts for time reversal and $C_{2}$ is a twofold rotation about the z-axis,
and consequently, as mention before, has to introduce extra symmetries (i.e.
it needs to fine tune certain hoppings to zero) to assure the presence of
Dirac nodes. Another example of this failure in capturing the essential
symmetries of the continuum model is that a perpendicular electric field
$E_{z}$ opens a gap in this model [35].

Parameter | Meaning | Ratio to $t_{1}$
---|---|---
$t_{1}$ |  $(A-B)$ $1^{st}$-nn hopping | 1
$t_{2}^{\prime}$ |  $(A-A)$ $2^{st}$-nn hopping | 0.208
$t_{4}$ |  $(A-B)$ $2^{nd}$-nn hopping | 0.167
$t_{5}$ |  $(A-B)$ $3^{rd}$-nn hopping | 0.333
Table 1: Hopping amplitudes for the 2B1V model of Ref.[27] written in units of
the dominant energy scale $t_{1}=1.2$ meV assuming $t_{2}\approx
it_{2}^{\prime}$.

Following Ref. [46] we truncate the hopping terms in the 2B1V model by
considering just the four largest ones. As we can see in Table 1, due to the
particular form of the Wannier orbitals the hopping elements do not vary
monotonically with distance (i.e. for instance $t_{5}>t_{4}$). Notice also
that we neglect the $(A-A)\;1^{st}$-nn hopping $t_{3}$ which is negligible
compared to $t_{2}$ [27]. The relative size of the hopping elements considered
is illustrated in Fig. 2 a).

Within this approximation the bulk Hamiltonian within the 2B1V model is given
by

$\displaystyle\hat{\mathcal{H}}^{(2)}_{\pm}(\mathbf{k})=\begin{pmatrix}\pm
t_{2}^{\prime}g(\mathbf{k})&&h(\mathbf{k})\\\ h(\mathbf{k})^{*}&&\pm
t_{2}^{\prime}g(\mathbf{k})\end{pmatrix},$ $\displaystyle
h(\mathbf{k})=t_{1}f_{1}(\mathbf{k})+t_{4}f_{4}(\mathbf{k})+t_{5}f_{5}(\mathbf{k}),$
(1)

where the subindex $\xi=\pm$ indicate the valley, $t_{2}\approx
it_{2}^{\prime}$ and

$\displaystyle g(\mathbf{k})$ $\displaystyle=$ $\displaystyle
2\sin(\mathbf{k}(2\mathbf{L_{2}}-\mathbf{L_{1}}))+2\sin(\mathbf{k}(2\mathbf{L_{1}}-\mathbf{L_{2}}))$
$\displaystyle-2\sin(\mathbf{k}(\mathbf{L_{1}}+\mathbf{L_{2}})),$
$\displaystyle f_{1}(\mathbf{k})$ $\displaystyle=$ $\displaystyle
1+e^{i\mathbf{k}\mathbf{L_{1}}}+e^{i\mathbf{k}\mathbf{L_{2}}},$ $\displaystyle
f_{4}(\mathbf{k})$ $\displaystyle=$ $\displaystyle
e^{i\mathbf{k}(\mathbf{L_{1}}+\mathbf{L_{2}})}+2\cos(\mathbf{k}(\mathbf{L_{1}}-\mathbf{L_{2}})),$
$\displaystyle f_{5}(\mathbf{k})$ $\displaystyle=$ $\displaystyle
e^{i2\mathbf{k}\mathbf{L_{1}}}+e^{i2\mathbf{k}\mathbf{L_{2}}}+e^{-i\mathbf{k}\mathbf{L_{1}}}+e^{-i\mathbf{k}\mathbf{L_{2}}}$
(2)
$\displaystyle+e^{i\mathbf{k}(2\mathbf{L_{1}}-\mathbf{L_{2}})}+e^{i\mathbf{k}(2\mathbf{L_{2}}-\mathbf{L_{1}})},$

where the lattice vectors are $\mathbf{L_{1,2}}=L_{m}(\sqrt{3}/2,\pm 1/2)$ and
$L_{m}=a/(2\sin{\theta/2})$ is the moiré lattice vector (given in terms of the
graphene lattice parameter $a$ and the twist angle $\theta$). The list of the
parameters used are given in Table 1.

Within this model a gap may be opened by breaking the sublattice symmetry
which protects the Dirac nodes with a perturbation of the form
$\hat{\mathcal{H}}^{\prime}_{\pm}(\mathbf{k})=\hat{\mathcal{H}}(\mathbf{k})_{\pm}+\delta_{S}\sigma_{z}$.
This could be associated either to the effect of the substrate or arising from
interactions treated perturbatively in a mean field approximation. In this
work we shall consider $\delta_{S}$ as a phenomenological parameter.

Figure 3: Evolution of the nearly-flat bands in the hexagonal BZ for the 2B1V
(upper pannels) and 6B1V (lower pannels) models due to the sublattice symmetry
breaking parameter $\delta_{S}$, which takes the values 0, 1, 10 meV for each
column respectively. Color code: grey bands indicate trivial topological
character while colored ones have topological nature. Straight and dashed-
dotted lines correspond to different global valleys. Red (blue) is associated
to bands with Chern number +1 (-1) on the other hand, magenta means fragile
topology bands as in panel d).

The flat bands in the triangular moiré Brillouin zone (BZ) for the 2B1V model
are shown in the upper upper panels of Fig. 3. Full and dashed lines
correspond to the two graphene valleys and panels a), b) and c) illustrate the
effect of the $\delta_{S}$ parameter opening a gap at the Dirac nodes, which
in this model have opposite helicity and thus are trivial in nature.

###  II.2 6B1V model

Ref. [35] discusses how to build “faithful" TB models avoiding the topological
obstruction by adding different sets of trivial bands to the flat ones. These
models should satisfy the requirements imposed by the relevant symmetries of
the system plus the non-trivial topology of the flat bands; and also that the
complementary set of trivial bands describe quantitatively the higher energy
bands in MATBLG. For that purpose Ref. [35] obtains an atomic limit for the
sum of the topological flat bands plus the determined set of trivial bands
with the constraints of having two isolated bands with the appropriate
momentum-space symmetry representations. This results in a basis of $p_{z}$
and $p_{\pm}$-orbitals in a triangular lattice $(\tau,p_{z})$ and
$(\tau,p_{\pm})$ respectively, and three $s$-orbitals in a kagome lattice
$(\kappa,s)$, using a notation that indicates the lattice (triangular, $\tau$,
or kagome, $\kappa$) and the character of the orbitals ($s$ or $p$).

The model includes two $(\tau,p_{\pm})-$orbitals sitting on a triangular
lattice based on the charge centers in MATBLG at the flat bands placed at the
A-A bonds in a triangular moiré lattice to reproduce the Dirac nodes of the
flat bands, see Fig. 2 c) and d). To capture the behaviour of the flat bands
close to the $\Gamma$ point, however, it is necessary to hybridize these
$(\tau,p_{\pm})$-orbitals with the rest of the orbitals needed to fulfill the
constraints associated to symmetry.

It should be noticed that recent works question the idea of fragile topology
in MATBLG [29, 50]. These works suggest that particle/hole symmetry has to be
considered although it is an approximate symmetry inherited from pristine
graphene close at charge neutrality. In the models of Ref. [35] fragile
topology relies exclusively in the $C_{2}\mathcal{T}$ symmetry and lattice
translations which protects the Dirac nodes. According to Refs.[29, 50], with
the approximate particle/hole symmetry the system topology becomes stable
instead of fragile. However, no lattice model can reproduce all these
symmetries. Still these tight-binding models have been used to predict the
correct correlated insulating behaviour [17, 18, 51, 16, 49].

Parameter | Meaning | Ratio to $t_{\kappa}$
---|---|---
$t_{p_{z}}$ |  $(\tau,p_{z})$ nn hopping | 0.17
$t_{p_{\pm}}$ |  $(\tau,p_{\pm})$ nn intra-orbital hopping | -0.03
$t^{+}_{p_{\pm}p_{\pm}}$ |  $(\tau,p_{\pm})$ nn inter-orbital hopping | -0.065
$t^{-}_{p_{\pm}p_{\pm}}$ |  $(\tau,p_{\pm})$ nn inter-orbital hopping | -0.055
$t_{\kappa}$ |  $(\kappa,s)$ nn hopping | 1
$t^{\prime}_{\kappa}$ |  $(\kappa,s)$ nnn hopping | 0.25
$t^{+}_{p_{\pm}p_{z}}$ |  $(\tau,p_{\pm})\times(\tau,p_{z})$ nn hopping | 0.095
$t^{-}_{p_{\pm}p_{z}}$ |  $(\tau,p_{\pm})\times(\tau,p_{z})$ nn hopping | 0.085
$t^{+}_{\kappa p_{\pm}}$ |  $(\kappa,s)\times(\tau,p_{\pm})$ nn hopping | 0.6
$t^{-}_{\kappa p_{\pm}}$ |  $(\kappa,s)\times(\tau,p_{\pm})$ nn hopping | 0.2
$\delta_{p_{z}}$ |  $(\tau,p_{z})$ chemical potential | -0.2593
$\delta_{p_{\pm}}$ |  $(\tau,p_{\pm})$ chemical potential | -0.3628
$\delta_{\kappa}$ |  $(\kappa,s)$ chemical potential | 0.20
Table 2: 6B1V model parameters written in units of the dominant energy scale
$t_{\kappa}=27$ meV

Although in Ref. [35] 5, 6 and 10 bands models have been proposed, in the
present work we find it appropriate to concentrate in the 6 bands case. While
we are looking for the minimal possible models which accounts for the system
fragile topology the reason for this choice is that unlike the 5 bands model,
the 6 bands one is formulated in terms of a TB Hamiltonian with usual nearest-
neighbours hopping terms instead of quasi-orbitals, which makes it more
amenable for applying Green function techniques in real space.

Details on the 6B1V model as defined on a minimal triangular unit cell can be
found in Ref. [35]. Here we extend it to situations where there could be cell
doubling (cd) due to interactions producing charge transfer between
neighboring sites [40, 10, 15, 14, 49].

Figure 4: Evolution of the nearly-flat bands in the rectangular BZ for one
valley in the 6B1V model due to the translational symmetry breaking potential
$\delta_{CDW}=$ 0, 0.3, 0.7, 1 meV for each panel a), b), c), d),
respectively, and interaction with the substrate $\delta_{S}=0.1$ meV. Color
code: grey indicates trivial topological character while colored ones have
topological nature. Red (blue) is associated to bands with Chern number +1
(-1).

This unit-cell doubling leads to an orthogonal lattice with vectors
$L_{x}=\sqrt{3}L_{m}$ and $L_{y}=L_{m}$ (see figure 2) and the local fermion
operators are defined as $\hat{\Psi}=(\hat{\Psi}_{A}\;\hat{\Psi}_{B})^{T}$,
with
$\hat{\Psi}_{\mu}=(\hat{\tau}_{p_{z},\mu}\;\hat{\tau}_{p_{+},\mu}\;\hat{\tau}_{p_{-},\mu}\;\hat{\kappa}^{(1)}_{s,\mu}\;\hat{\kappa}^{(2)}_{s,\mu}\;\hat{\kappa}^{(3)}_{s,\mu}\;)^{T}$,
where $\mu\equiv A,B$ indicates the two sites within the orthogonal cell.
Within this basis the 6B1V Hamiltonian adopts the form

$\hat{\mathcal{H}}^{(6)}_{cd}(\textbf{k})=\begin{pmatrix}\hat{H}^{AA}&&\hat{H}^{AB}\\\
\hat{H}^{AB\dagger}&&\hat{H}^{AA}\end{pmatrix},$ (3)

where

$\displaystyle\hat{H}^{AA}=\begin{pmatrix}H_{p_{z}}^{AA}+\mu_{p_{z}}&&\hat{C}_{p_{\pm}p_{z}}^{AA}&&\hat{0}\\\
\hat{C}_{p_{\pm}p_{z}}^{AA\dagger}&&\hat{H}_{p_{\pm}}^{AA}+\mu_{p_{\pm}}\hat{\mathbb{I}}_{2}&&\hat{C}_{\kappa
p_{\pm}}^{AA}\\\ \hat{0}&&\hat{C}_{\kappa
p_{\pm}}^{AA\dagger}&&\hat{H}_{\kappa}^{AA}+\mu_{\kappa}\hat{\mathbb{I}}_{3}\end{pmatrix},$
$\displaystyle\hat{H}^{AB}=\begin{pmatrix}H_{p_{z}}^{AB}&&\hat{C}_{p_{\pm}p_{z},1}^{AB}&&\hat{0}\\\
\hat{C}_{p_{\pm}p_{z},2}^{AB}&&\hat{H}_{p_{\pm}}^{AB}&&\hat{C}_{\kappa
p_{\pm},1}^{AB}\\\ \hat{0}&&\hat{C}_{\kappa
p_{\pm},2}^{AB}&&\hat{H}_{\kappa}^{AB}\end{pmatrix}.$ (4)

The expressions for the matrices $H^{\mu\nu}_{\alpha}$ and
$C^{\mu\nu}_{\alpha\beta}$ in terms of the parameters in Table 2 are given in
Appendix A. The diagonal elements are set by
$\mu_{p_{z}}\equiv-6t_{p_{z}}+\delta_{p_{z}}$, $\mu_{p_{\pm}}\equiv
3t_{p_{\pm}}+\delta_{p_{\pm}}$ and
$\mu_{\kappa}\equiv-4(t_{\kappa}+t_{\kappa}^{\prime})+\delta_{\kappa}$.

In addition, we will consider the following symmetry breaking perturbations

$\hat{H}_{cd}^{S}=\delta_{S}\tau^{p_{\pm}}_{z}\sigma_{0}\quad,\quad\hat{H}_{cd}^{CDW}=\delta_{CDW}\hat{\mathbb{I}}_{6}\sigma_{z},$
(5)

where $\tau^{p_{\pm}}_{z}$ is a Pauli matrix acting in the
$(\tau,p_{\pm})_{\mu}$-orbitals with $\mu=(A,B)$, $\sigma_{0}$ is the identity
matrix and $\sigma_{\mu}$ are Pauli matrices which act in the $(A-B)$ doubled
unit cell sites space. The $\delta_{S}$ parameter is equivalent to a
sublattice staggered moiré potential as pointed out in Ref. [19].

The total Hamiltonian in the 6B1V model is then given by

$\hat{\mathcal{H}}^{(6)\prime}_{cd}(\textbf{k})=\hat{\mathcal{H}}^{(6)}_{cd}(\textbf{k})+\hat{H}_{cd}^{S}+\hat{H}_{cd}^{CDW}.$
(6)

Figure 5: Same as in Fig. 4 but with the translational symmetry breaking
potential $\delta_{CDW}=$ 0, 1, 4, 10.5 meV and strong coupling to the
substrate, $\delta_{S}=10$ meV. Panel d) shows a gap closing due to band
inversion between central Chern bands. Color code: grey bands shows trivial
topological character while colored ones have topological nature. Red (blue)
is associated to bands with Chern number +1 (-1), yellow in panel b)
corresponds to Chern number +2.

These perturbations can be either associated to the effect of the substrate
(i.e. hBN substrate induces sublattice symmetry breaking in the graphene
lattice placed on top) or to a CDW arising from interactions (i.e. breaking
translational symmetry) allowing us to explore a larger phase space [13, 12,
11, 40, 15, 14, 49].

On the lower panels of Fig. 3 we show the flat bands in the triangular BZ for
the 6B1V model for increasing values of the gap opening parameter
$\delta_{S}$. As can be observed, for $\delta_{S}=0$ the bands within the 2B1V
and the 6B1V models are quite similar despite the absent electron/hole
symmetry in 6B1V model. For sufficiently large $\delta_{S}$ their behavior
differ substantially. While in the 2B1V model the band minimum continues to
appear at $\bar{K}$ and $\bar{K}^{\prime}$, for the 6B1V model the minimum is
transferred into the $\bar{\Gamma}$ point. This behavior (which was also
pointed out in Ref. [40, 10]) is due to the fact that the states at the
$\bar{\Gamma}$ point in the 6B1V model are combination of the $(\tau,p_{z})$
and $(\kappa,s)$-orbitals, which are coupled to the perturbed $(\tau,p_{\pm})$
ones.

On the other hand, the models differ, as expected, on their topological
properties. The different colors in Fig. 3 indicate the bands’ Chern number as
extracted from the winding of the Berry curvature in the BZ of the
corresponding Wilson loops [52, 53, 54, 55]. When $\delta_{S}\neq 0$ the
positive and negative energy flat bands at the CNP in the 6B1V model are
characterized by a non-zero Chern number. A change in the global valley
induces a change in sign of the Chern number for each band due to time
reversal symmetry. The 2B1V bands remain trivial in all cases.

The effect of the translation symmetry breaking parameter $\delta_{CDW}$ and
the combined effect of $\delta_{CDW}$ and $\delta_{S}$ in the 6B1V model are
illustrated in Figs. 4 and 5. In Fig. 4 we show the low energy bands in the
orthogonal BZ for weak coupling with the substrate ($\delta_{S}=0.1$ meV) and
varying $\delta_{CDW}$. The former opens a gap at the CNP and the latter
splits each nearly-flat band as shown in panel b), redistributing the
topological charge of the original flat bands in a set of splitted trivial and
topological bands. The total topological charge at charge neutrality is still
$C=+1$. By effect of the competing perturbations a gap closing occurs in panel
c) due to a band inversion that produces in panel d) the trivialization of the
lower set of flat bands at the CNP.

On the other hand, in Fig. 5 we show the effect of translational symmetry
breaking in the case of strong coupling with the substrate ($\delta_{S}=10$
meV). As in Fig. 4 $\delta_{CDW}$ splits the flat bands and redistributes the
topological charge between the new set of bands. In this case the lower flat
bands are crossed by the set of trivial occupied bands. Although the Chern
number at charge neutrality in panel b) is $C=+1$, the redistribution of
topological charge is completely different due to the effect of the excited
bands which, individually considered, can carry some topological order. For
this reason we observe a flat band with $C=+2$. In the subsequent panels we
observe different topological transitions: first in panel c) a gap closing
between the lower topological bands that set $C=+1$ to the flat band close to
the charge neutrality and trivializes the rest of lowest bands and finally in
panel d) a band inversion occurs between the flat bands close to the CNP.

The competition between Chern insulating states and CDW in TBLG has been
recently analyzed in Refs. [13, 12, 18]. As we can see in Figs. 4 and 5,
starting from a Chern insulating phase the breaking of translational symmetry
can lead to a trivial phase, in agreement with these works.

##  III bGF method for 2D TB models

A central aim of the present work is to generalize the bGF technique develop
in [43] for 1D models to the case of 2D lattice models, particularly to the
ones describing bulk MATBLG by a $N\times N$ Hamiltonian in reciprocal space
$\hat{\mathcal{H}}(\textbf{k})$.

Figure 6: Open AC boundary LDOS for the 2B1V model showing trivial edge states
under the effect of the breaking of sublattice symmetry $\delta_{S}$. Left
column: spectral density for a right boundary. Right column: integrated LDOS
where straight (dot-dashed) line represents bulk (right boundary) LDOS. Panel
a): Weak coupling limit $\delta_{S}=0.1$ meV where a gap opens with a minimal
distorsion of the bands. It exhibits trivial edge states. Panel b): Strong
coupling limit $\delta_{S}=10$ meV where the bands are greatly distorted and
the trivial edge states disappear.

In order to compute the bGF we decompose the momenta into parallel and
perpendicular components $\textbf{k}=(k_{\parallel},k_{\perp})$ relative to
the boundary direction that we consider (in higher dimensional models the
parallel momentum component would be itself a vector
$\textbf{k}_{\parallel}$). The bulk Hamiltonian periodicity in both directions
is set by $(L_{\parallel},L_{\perp})$, such that
$\hat{\mathcal{H}}(\textbf{k}+2\pi\textbf{u}_{\perp}/L_{\perp})=\hat{\mathcal{H}}(\textbf{k})$,
where $\textbf{u}_{\perp}$ is the unitary vector in the perpendicular
direction. Thus, for a model with a triangular lattice like the 6B1V model we
have to double the primitive cell to obtain an orthogonal decomposition of the
momentum. Using this periodicity, the Hamiltonian can be expanded in a Fourier
series,
$\hat{\mathcal{H}}(\textbf{k})=\sum_{n}\hat{\mathcal{V}}_{n}(k_{\parallel})e^{ink_{\perp}L_{\perp}}$,
where $n$ is the number of neighbours and Hermiticity implies
$\hat{\mathcal{V}}_{-n}=\hat{\mathcal{V}}^{\dagger}_{n}$.

The advanced bulk GF is then defined as

$\hat{G}^{A}(\textbf{k},\omega)=\left[(\omega-i0^{+})\hat{\mathbb{I}}-\hat{\mathcal{H}}(\textbf{k})\right]^{-1},$
(7)

where the $N\times N$ matrix structure is indicated by the hat notation.
Fourier transforming along the perpendicular direction, the GF components are
given by

$\hat{G}^{A}_{jj^{\prime}}(k_{\parallel},\omega)=\frac{L_{\perp}}{2\pi}\int\limits_{-\pi/L_{\perp}}^{\pi/L_{\perp}}dk_{\perp}e^{i(j-j^{\prime})k_{\perp}L_{\perp}}\,\hat{G}^{A}(k_{\parallel},k_{\perp},\omega),$
(8)

where $j$ and $j^{\prime}$ are lattice site indices. By the identification
$z=e^{ik_{\perp}L_{\perp}}$, this integral is converted into a complex contour
integral,

$\hat{G}^{A}_{jj^{\prime}}(k_{\parallel},\omega)=\frac{1}{2\pi
i}\oint\limits_{|z|=1}\frac{dz}{z}z^{j-j^{\prime}}\hat{G}^{A}(k_{\parallel},z,\omega).$
(9)

Figure 7: Open ZZ boundary LDOS for the 6B1V model without symmetry breaking
perturbations. As in Fig. 6 the left column corresponds to the spectral
density for a right boundary while the right column shows the integrated LDOS
where straight (dot-dashed) line represents bulk LDOS $\langle\rho_{0}\rangle$
(right boundary LDOS $\langle\rho_{R}\rangle$). The upper and lower panels
show the same results but in different energy scale. In panel a) we can
observe trivial edge states inherit from pristine monolayer graphene appearing
within the nearly-flat bands. Panel b) illustrates the appearance of moiré
edge states in the gap between nearly-flat and the lower energy bands.

Further simplification can be obtained by introducing the roots
$z_{n}(k_{\parallel},\omega)$ of the secular polynomial in the complex-$z$
plane,

$\displaystyle P(k_{\parallel},z,\omega)$ $\displaystyle=$
$\displaystyle\mbox{det}\left[\omega\hat{\mathbb{I}}-\hat{\mathcal{H}}(k_{\parallel},z)\right]$
(10) $\displaystyle=$
$\displaystyle\frac{c_{m}}{z^{m}}\prod_{n=1}^{2m}\left[z-z_{n}(k_{\parallel},\omega)\right],$

where $m$ is the highest order of the polynomial and $c_{m}$ is the highest
order coefficient. In terms of these roots the contour integral in Eq. (9) can
be written as a sum over the residues of all roots inside the unit circle

$\hat{G}^{A}_{jj^{\prime}}(k_{\parallel},\omega)=\sideset{}{{}^{\prime}}{\sum}_{|z_{n}|<1}\frac{z_{n}^{p}\hat{A}(k_{\parallel},z_{n},\omega)}{c_{m}\prod_{l\neq
n}\left(z_{n}-z_{l}\right)},$ (11)

where $p=j-j^{\prime}+m-m^{\prime}-1$ and
$z^{-m^{\prime}}\hat{A}(k_{\parallel},z,\omega)$ is the adjugate matrix of
$\left[\omega\hat{\mathbb{I}}-\hat{\mathcal{H}}(k_{\parallel},z)\right]$ where
all the poles at zero were taken out of $\hat{A}$ as a common factor in
$z^{-m^{\prime}}$. Finally, $\sum^{{}^{\prime}}$ means that if $p=-1$, then we
include $z_{n}=0$ as a pole to take into account in the sum of residues (i.e.
in the non local GF components). When $p<-1$ higher order poles at zero appear
in the sum of residues. To simplify we can take advantage of the residue
theorem to avoid these poles and compute the integral as

$\hat{G}^{A}_{jj^{\prime}}(k_{\parallel},\omega)=-\sum_{|z_{n}|>1}\frac{z_{n}^{p}\hat{A}(k_{\parallel},z_{n},\omega)}{c_{m}\prod_{l\neq
n}\left(z_{n}-z_{l}\right)}.$ (12)

Figure 8: Same as Fig. 7 but under the effect of the coupling with the
substrate $\delta_{S}$ and broken translation symmetry $\delta_{CDW}$. Panel
a): Weak coupling with the substrate ($\delta_{S}=0.1$ meV) where interactions
open a gap without distorting the bands. It shows the appearance of
topological edge state between the flat bands at the $\bar{P}$ minivalley.
Panel b): Strong coupling with the substrate ($\delta_{S}=10$ meV) where
interactions distort the bands and the breaking of translational symmetry with
$\delta_{CDW}=1$ meV splits both flat bands. One can observe topological edge
states between the flat bands and another topological state with opposite
chirality between the splitted flat band and the lower energy bands due to
breaking of translation symmetry.

In [45] we describe an efficient method, based on Fadeev-LeVerrier algorithm
to construct the secular polynomial and the adjugate matrix which are needed
to evaluate the GF components in Eq. (11). On the other hand, in Appendix B we
give details on the method to obtain the bGF for the system with an open
boundary and to get the transport properties for junctions between different
TBLG regions.

##  IV Spectral properties at edges

The bGF $\hat{\mathcal{G}}_{L,R}$ of the semi-infinite system obtained as
described in Appendix B allows us to analyze the spectral properties at the
edges of MATBLG, encoded in the spectral densities
$\rho_{L,R}(k_{\parallel},\omega)=\mbox{Im}\hat{\mathcal{G}}_{L,R}(k_{\parallel},\omega)/\pi$
and the LDOS $\langle\rho_{L,R}(\omega)\rangle=\int
dk_{\parallel}/\Omega_{k_{\parallel}}\rho_{L,R}(k_{\parallel},\omega)$, where
$\Omega_{k_{\parallel}}=2\pi/L_{\parallel}$ accounts for the limits of
integration.

In Fig. 6 we illustrate the spectral properties at an armchair edge for the
2B1V model in the weak ($\delta_{S}=0.1$ meV) (a) and strong ($\delta_{S}=10$
meV) (b) coupling limits. In the AC boundary the $\bar{P}$ and $\bar{Q}$
minivalleys are projected onto the $\bar{\Gamma}$ point. As expected, no
chiral states appear when opening a gap in the 2B1V model. Only trivial states
can be observed in case (a). In the strong coupling limit the edge states just
merge with the continuum.

Figure 9: Topological domain wall with ZZ orientation in the 6B1V model. This
is produced by a change of sign of $|\delta_{S}|=0.1$ meV across the junction
and we also include translational symmetry breaking with $\delta_{CDW}=0.5$
meV. Panel a) shows the local spectral density at the domain wall where a
chiral state for each mini-valley with the same chirality can be observed.
Panel b) shows the integrated LDOS in log-scale along the domain wall where
the chiral states are localized around the center with a characteristic
penetration length of $\approx 200nm$. The vertical white dashed line
represents the domain wall frontier where the sign of $\delta_{S}$ changes.
The density of states are measured in [meV-1] units.

The edge spectral properties are radically different in the 6B1V model. Due to
the contribution of the set of higher energy bands in the 6B1V model we can
distinguish two types of edge states in Fig. 7: the ones inherited from the
pristine graphene appearing at the charge neutrality point (CPN) at each
minivalley ($\bar{Q}$ and $\bar{P}$) in panel (a) and the ones that emerge
entirely from the moiré structure at the gap between the flat bands and the
excited ones [47] (see panel b). One of these moiré edge states merges with
the flat bands at $\bar{\Gamma}$ and the other one is gaped and nearly
dispersionless. This gaped edge state also appears at the AC boundary and
might be related to higher order topology [56, 57, 58].

A gap opening at the CNP allows us to confirm the topological character of the
nearly flat bands. In Fig. 8 a) for a right ZZ boundary a trivial edge state
persists around the $\bar{Q}$ minivalley and there appears a chiral edge state
around $\bar{P}$ minivalley in the weak coupling limit with $\delta_{S}=0.1$
meV. In addition, the left boundary shows a topological edge state with
opposite chirality at $\bar{Q}$.

Figure 10: Two-terminal transport for a three region junction device with AC
orientation as a function of the height of the potential barrier $U_{b}$ in
the central region for a) 2B1V and b) 6B1V models. In both models we take
$\delta_{S}=1$ meV. Left column: momentum resolved transmission of the
junction. Right column: Solid red line represents the average transmission.
The length of the central region is $\lambda_{b}\approx 180$ nm and the Fermi
energy is kept at $E_{F}=-1$ meV on the lateral regions.

In panel b) we can analyze the effect of $\delta_{CDW}$ in the strong coupling
limit, $\delta_{S}=10$ meV. In the first place, there appears a chiral edge
state at charge neutrality between the new band minima around $\bar{\Gamma}$
and thus it loses the minivalley character which was present in panel a). As
we pointed out in Fig. 5 b), at CNP the occupied bands have $C=+1$. On the
other side, the gap between the lower flat band and the rest of lower energy
bands has $C=-1$. For that reason we can observe another chiral edge state
with opposite chirality with respect to the one at charge neutrality. Notice
that changing the global valley degree of freedom we obtain a similar picture
for each boundary where the topological edge states have the opposite
chirality and appears at the opposite minivalley (results not shown)

Figure 11: Finite size effects in the local spectroscopy for a three region
junction N-CNP-N with $\delta_{S}=0.1$ meV for both AC (left panels) and ZZ
(right panels) boundaries. Panel a) and e): LDOS across the central region of
the junction for AC and ZZ boundaries respectively. We can observe
anticrossings at low energies due to the hybridization of the edge states
located at the ends of the central region. The central panels show the
spectral density at both ends of the central region of each junction, upper
(lower) ones shows the left (right) end in arbitrary units. Panels d) and h)
show the LDOS asymmetry dimensionless factor $\Lambda_{LR}$ between the region
ends at the CNP as a function of the sublattice symmetry breaking potential
$\delta_{S}$. The length of the central region is $\lambda_{b}\approx 270$ nm
and its doping level is kept at the CNP, the lateral N regions are heavily
doped with $E_{F}=-140$ meV.

In Fig. 9 we show the spectral properties at a domain wall formed by an abrupt
change in the sign of the coupling with the substrate $\delta_{S}$ as
illustrated in Fig. 1 b). As can be observed in panel a) there appear two
chiral edge states at different minivalleys propagating in the same direction.
These correspond to states on the two sides of the domain wall which do not
annihilate as they “live" at different minivalleys. Notice that the chirality
of the valley polarized edge states depends on the sign of the sublattice
perturbation $\delta_{S}$ and on the general valley.

In panel b) we show the LDOS as a function of the distance from the domain
wall located at the center of the plot. It illustrates how the chiral edge
state inside the gap decays on a scale of 100 nm and the bulk gap is
reconstructed.

##  V Transport at N-N’-N junctions

In this section we analyze two terminal transport in different types of MATBLG
junctions. In Appendix. C we extend the GF method to take into account a non-
uniform potential profile in the central region, as in the situation depicted
in Fig. 1 a).

The two-terminal conductance properties which are predicted by the 2B1V and
6B1V models in the situation of Fig. 1 a) are illustrated in Fig. 10 as a
function of the height of the barrier, $U_{b}$, which defines the doping level
in the central region. The lateral regions are kept at $E_{F}=-1$ meV. For
both models we include a finite gap opening parameter of the same size
($\delta_{S}$ = 1 meV). As can be observed, in Fig. 10 there is a qualitative
and quantitative agreement between the two models, which reflects the fact
that the topological character of the bands in the 6B1V model does not play a
role in this case. In agreement with results of Ref. [46] for an abrupt
barrier the momentum resolved transmission exhibits oscillations due to finite
size effects. This spectral density oscillations are linked to the finite
length of the central region $\lambda_{B}$. We can also observe that the
normalized total transmission exhibits a maximum at the Van Hove singularity.
A lower transmission is obtained when the barrier induces injection of states
through the valence band.

The topological properties, however, do play a role in the case of N-N’-N
junctions where N’ is kept at the CNP and the side regions are heavily doped.
The outer regions in this junction exhibit smaller DOS than the flat bands,
differing from the metallic behaviour of heavily doped pristine graphene. The
spectral properties for this type of junctions are illustrated in Fig. 11. We
consider two different orientations of the junctions with respect to the moiré
lattice: AC (panels a-d) and ZZ (panels e-h). As illustrated by these results,
the orientation of the junctions leads to substantial differences even when
the gate potential profile is kept constant.

In the case of the armchair orientation we observe unexpected but mild left-
right (LR) asymmetry in the LDOS. This is due to the sublattice symmetry
breaking by $\delta_{S}$ which acts in the directional
$(\tau,p_{\pm})$-orbitals. Due to the large penetration depth of the chiral
edge states located at the ends of the central region we observe minigaps
associated to hybridization of these states. The spectral densities of the
boundaries show the hybridized opposite chirality end states. To analyze the
evolution of the LR asymmetry as a function of $\delta_{S}$ we define

$\Lambda_{LR}(\omega)=\frac{\langle\rho_{L}(\omega)\rangle-\langle\rho_{R}(\omega)\rangle}{\langle\rho_{L}(\omega)\rangle+\langle\rho_{R}(\omega)\rangle},$
(13)

which corresponds to a LDOS LR asymmetry dimensionless factor. Panel d) shows
that in the limit $\delta_{S}\rightarrow 0$ the AC junction does not show any
LR asymmetry because without perturbations both AC boundaries are equivalent.
Valley symmetry does not lift this LDOS asymmetry and only induces an
inversion of the spectral density in $k_{\parallel}$ (i.e. changing the
chirality and minivalley of the topological edge states). The step-like
increase of the LDOS at $\omega\approx-0.15$ meV is due to the top of valence
band. The bottom of the conduction band, placed at $\omega\approx 0.2$ meV.

On the other hand, ZZ junctions exhibit a much stronger LR asymmetry which
survives even in the limit $\delta_{S}\rightarrow 0$, see panel h). Within the
6B1V model this is due to the definition of the unit cell and the effect of
kagome orbitals at the boundaries, as can be seen in Fig. 2 d). Left (right)
boundaries are defined on triangular (kagome) sites and for that reason both
boundaries are naturally non-equivalent. This unperturbed limit for the LR
asymmetry relies in the fact the the charge centers are displaced from the
triangular lattice sites in the direction perpendicular to the ZZ edge by
effect of the kagome sites. Nevertheless, this is not exclusively a feature of
the 6B1V model but a general property of MATBG with Wannier orbitals
resembling a fidget spinner [27, 10, 59].

In addition, as for ZZ edges minivalleys are not coupled and opposite
chirality states are placed at different minivalleys, the hybridization
between chiral states is also asymmetric inducing anticrossings at different
places in the BZ. Panel f) shows an interplay between the chiral states at
$\bar{P}$ and the trivial ZZ edge state around $\bar{Q}$ at the left boundary,
which hybridizes with the chiral states at the right boundary at $\bar{Q}$,
see panel g).

Finally in Fig. 12 a) we illustrate the transport properties for a gate
voltage potential profile as in Fig. 11 but with a reduced length of the
central region to $\lambda_{b}\approx 60$ nm in order to obtain a larger
overlap between the chiral edge states at the ends of this narrow central
region. We observe in this case a small but finite transmission in the gap
region which can be associated with the chiral states hybridization. On the
other hand, for the strong coupling case the overlap diminishes because the
penetration length decays with $\delta_{S}$ and thus the transmission in the
gap region becomes negligible, see Fig. 12 b).

Figure 12: Two-terminal transport properties for a three region junction
N-CNP-N with ZZ orientation. Panel a): Weak coupling limit with the substrate
$\delta_{S}=0.1$ meV. Panel b): Strong coupling limit with the substrate
$\delta_{S}=10$ meV plus breaking of translational symmetry with
$\delta_{CDW}=1$ meV. Left column: spectral representation of the transmission
at the junction in log-scale. Right column: solid red line represents the
average transmission of the nearly-flat bands. The length of the central
region is $\lambda_{b}\approx 60$ nm and its Fermi energy is kept at the CNP,
while the normal regions N are heavily doped with $E_{F}=-140$ meV.

##  VI Conclusions and outlook

In this work we have analyzed electronic and transport properties of MATBLG
junctions using lattice models. We have considered both the 2B1V model of Ref.
[27] and the 6B1V model of Ref. [35], the last one accounting for fragile
topology. For that purpose we have implemented a boundary Green function
approach by extending the method developed in Ref. [43] to 2D lattices with
hopping elements between arbitrary distant neighbors.

Using this approach we have analyzed the spectral properties at different
types of edges showing the appearance of topological states for the 6B1V model
in the presence of symmetry breaking perturbations. We have also analyzed the
LDOS at a domain wall defined by the sublattice symmetry breaking potential
showing the coexistence of topological edge states at different minivalleys
with the same chirality. This analysis allows us to characterize the
penetration length of the chiral states.

Regarding the transport properties between MATLBG regions doped at the flat
bands through a barrier, we have found that these are not sensitive to
topology and thus both the 2B1V and the 6B1V models give qualitatively similar
results. In contrast, for a three region junction where the central one is
fixed at the CNP and the lateral ones are heavily doped, transport and LDOS
are determined by the hybridization of chiral states running in opposite
directions at the interfaces. As we have shown, an asymmetry in the LDOS
between the left and right interface appears which could allow one to
distinguish the orientation of the junction with respect to the moiré lattice.

As an outlook, we foresee several possible applications of the method
developed in the present work. On the one hand, an extension to a Bogoliubov
de Gennes formulation would allow to study Josephson and Andreev transport in
MATBLG junctions including superconducting correlations [60, 61, 62], in line
with recent experiments [9]. This type of calculations could help to identify
the symmetry of the order parameter in superconducting MATBLG which is at
present an open issue attracting great interest [48, 63, 64, 65, 66]. Similar
calculations have been implemented in the past by some of the authors for the
case of pristine graphene layers in proximity with unconventional and
topological superconductors [67]. Another extension of our work could be
oriented to analyze the possible higher-order topological insulating (HOTI)
states in MATBLG and the appearance of corner states associated to the gapped
moiré bound states [26, 30, 56, 58, 57].

###### Acknowledgements.

We acknowledge fruitful discussions with L. Brey, P. Burset, Ch. de Beule, P.
Recher and T. Stauber and thank all them for useful comments on this
manuscript. This project has been funded by the Spanish MICINN through Grant
No. FIS2017-84860-R; and by the María de Maeztu Programme for Units of
Excellence in n Research and Development Grant No. MDM-2014-0377.

##  Appendix A 6B1V Hamiltonian with unit cell doubling

As defined in the main text, for the 6B1V model we consider a doubled unit
cell leading to an orthogonal lattice with vectors $L_{x}=\sqrt{3}L_{m}$ and
$L_{y}=L_{m}$ and local fermion operators
$\hat{\Psi}=(\hat{\Psi}_{A}\;\hat{\Psi}_{B})^{T}$, with
$\hat{\Psi}_{\mu}=(\hat{\tau}_{p_{z},\mu}\;\hat{\tau}_{p_{+},\mu}\;\hat{\tau}_{p_{-},\mu}\;\hat{\kappa}^{(1)}_{s,\mu}\;\hat{\kappa}^{(2)}_{s,\mu}\;\hat{\kappa}^{(3)}_{s,\mu})^{T}$,
where $\mu\equiv A,B$ indicates the two sites within the doubled unit-cell.
The 6B1V Hamiltonian in this basis adopts the general form

$\displaystyle\hat{\mathcal{H}}^{(6)}_{cd}(\textbf{k})=\begin{pmatrix}\hat{H}^{AA}&&\hat{H}^{AB}\\\
\hat{H}^{AB\dagger}&&\hat{H}^{AA}\end{pmatrix},$ (14)

where

$\hat{H}^{AA}=\begin{pmatrix}H_{p_{z}}^{AA}+\mu_{p_{z}}&&\hat{C}_{p_{\pm}p_{z}}^{AA}&&\hat{0}\\\
\hat{C}_{p_{\pm}p_{z}}^{AA\dagger}&&\hat{H}_{p_{\pm}}^{AA}+\mu_{p_{\pm}}\hat{\mathbb{I}}_{2}&&\hat{C}_{\kappa
p_{\pm}}^{AA}\\\ \hat{0}&&\hat{C}_{\kappa
p_{\pm}}^{AA\dagger}&&\hat{H}_{\kappa}^{AA}+\mu_{\kappa}\hat{\mathbb{I}}_{3}\end{pmatrix}\quad,\quad\hat{H}^{AB}=\begin{pmatrix}H_{p_{z}}^{AB}&&\hat{C}_{p_{\pm}p_{z},1}^{AB}&&\hat{0}\\\
\hat{C}_{p_{\pm}p_{z},2}^{AB}&&\hat{H}_{p_{\pm}}^{AB}&&\hat{C}_{\kappa
p_{\pm},1}^{AB}\\\ \hat{0}&&\hat{C}_{\kappa
p_{\pm},2}^{AB}&&\hat{H}_{\kappa}^{AB}\end{pmatrix}.$ (15)

We now define all the contributions to the total Hamiltonian classified for
its orbital character. First the nearest-neighbor standard hopping for
$(\tau,p_{z})$-orbitals in a cd basis

$H_{p_{z}}^{AA}=t_{p_{z}}(\phi_{0\bar{1}}+\phi_{01})\quad,\quad
H_{p_{z}}^{AB}=t_{p_{z}}(1+\phi_{0\bar{1}}(1+\phi_{10})+\phi_{10}),$ (16)

where $\phi_{10}=e^{-ik_{x}L_{x}}$ and $\phi_{01}=e^{-ik_{y}L_{y}}$ with
$\phi_{\bar{ij}}=\phi_{ij}^{\dagger}$.

For the triangular lattice $(\tau,p_{\pm})-$orbitals we have the following
intra-orbital terms

$\hat{H}_{p_{\pm},0}^{AA}=t_{p_{\pm}}(\phi_{0\bar{1}}+\phi_{01})\hat{\mathbb{I}}_{2}\quad,\quad\hat{H}_{p_{\pm},0}^{AB}=t_{p_{\pm}}(1+\phi_{0\bar{1}}(1+\phi_{10})+\phi_{10})\hat{\mathbb{I}}_{2},$
(17)

and the inter-orbital terms

$\displaystyle\hat{C}_{p_{\pm}p_{\pm}}^{AA}=\begin{pmatrix}0&&t^{+}_{p_{\pm}p_{\pm}}\phi_{0\bar{1}}+t^{-}_{p_{\pm}p_{\pm}}\phi_{01}\\\
t^{+}_{p_{\pm}p_{\pm}}\phi_{01}+t^{-}_{p_{\pm}p_{\pm}}\phi_{0\bar{1}}&&0\end{pmatrix},$
(18)
$\displaystyle\hat{C}_{p_{\pm}p_{\pm}}^{AB}=\begin{pmatrix}0&&t^{+}_{p_{\pm}p_{\pm}}(w+w^{\prime}\phi_{10})+t^{-}_{p_{\pm}p_{\pm}}(w^{\prime}\phi_{0\bar{1}}+w\phi_{0\bar{1}}\phi_{10})\\\
t^{+}_{p_{\pm}p_{\pm}}(w^{\prime}\phi_{0\bar{1}}\phi_{10}+w\phi_{0\bar{1}})+t^{-}_{p_{\pm}p_{\pm}}(w^{\prime}+w\phi_{10})&&0\end{pmatrix}$
,

where $w=e^{i2\pi/3}$ and $w^{\prime}$ is the complex conjugate of $w$. This
phase accounts the different orientation of $(\tau,p_{\pm})$-orbitals. We can
change valley with the transformation $w\rightarrow w^{\prime}$. The final
expression of the total Hamiltonian for this subspace has the form

$\hat{H}_{p_{\pm}}^{AA}=\hat{H}_{p_{\pm},0}^{AA}+\hat{C}_{p_{\pm}p_{\pm}}^{AA}\quad,\quad\hat{H}_{p_{\pm}}^{AB}=\hat{H}_{p_{\pm},0}^{AB}+\hat{C}_{p_{\pm}p_{\pm}}^{AB},$
(19)

which contains the fundamental symmetries that protect the Dirac nodes in the
triangular lattice. Next we consider the hopping terms concerning the kagome
lattice $(\kappa,s)$-orbitals which include first and second nearest neighbour
contributions within this subspace in the form

$\displaystyle\hat{H}_{\kappa}^{AA}=t_{\kappa}\begin{pmatrix}0&&1&&1\\\
1&&0&&1+\phi_{0\bar{1}}\\\
1&&1+\phi_{01}&&0\end{pmatrix}+t^{\prime}_{\kappa}\begin{pmatrix}0&&\phi_{01}&&\phi_{0\bar{1}}\\\
\phi_{0\bar{1}}&&0&&0\\\ \phi_{01}&&0&&0\end{pmatrix},$
$\displaystyle\hat{H}_{\kappa}^{AB}=t_{\kappa}\begin{pmatrix}0&&1&&\phi_{0\bar{1}}\\\
\phi_{0\bar{1}}\phi_{10}&&0&&0\\\
\phi_{10}&&0&&0\end{pmatrix}+t^{\prime}_{\kappa}\begin{pmatrix}0&&\phi_{0\bar{1}}&&1\\\
\phi_{10}&&0&&\phi_{0\bar{1}}(1+\phi_{10})\\\
\phi_{10}\phi_{0\bar{1}}&&1+\phi_{10}&&0\end{pmatrix}.$ (20)

On the other hand, the coupling terms between the $(\tau,p_{z})$ and the
$(\tau,p_{\pm})$ -orbitals are given by

$\displaystyle\hat{C}_{p_{\pm}p_{z}}^{AA}=i\begin{pmatrix}-t^{+}_{p_{\pm}p_{z}}\phi_{0\bar{1}}+t^{-}_{p_{\pm}p_{z}}\phi_{01}\\\
t^{+}_{p_{\pm}p_{z}}\phi_{01}-t^{-}_{p_{\pm}p_{z}}\phi_{0\bar{1}}\end{pmatrix}^{T}\quad,\quad\hat{C}_{p_{\pm}p_{z},1}^{AB}=i\begin{pmatrix}-t^{+}_{p_{\pm}p_{z}}(w+w^{\prime}\phi_{10})+t^{-}_{p_{\pm}p_{z}}(w^{\prime}\phi_{0\bar{1}}+w\phi_{0\bar{1}}\phi_{10})\\\
t^{+}_{p_{\pm}p_{z}}(w\phi_{0\bar{1}}+w^{\prime}\phi_{0\bar{1}}\phi_{10})-t^{-}_{p_{\pm}p_{z}}(w^{\prime}+w\phi_{10})\end{pmatrix}^{T},$
$\displaystyle\hat{C}_{p_{\pm}p_{z},2}^{AB}=i\begin{pmatrix}t^{+}_{p_{\pm}p_{z}}(w\phi_{0\bar{1}}+w^{\prime}\phi_{0\bar{1}}\phi_{10})-t^{-}_{p_{\pm}p_{z}}(w^{\prime}+w\phi_{10})\\\
-t^{+}_{p_{\pm}p_{z}}(w+w^{\prime}\phi_{10})+t^{-}_{p_{\pm}p_{z}}(w^{\prime}\phi_{0\bar{1}}+w\phi_{0\bar{1}}\phi_{10})\end{pmatrix},$
(21)

and finally the coupling orbital terms in the $(\tau,p_{\pm})\times(\kappa,s)$
sub-spaces

$\displaystyle\hat{C}_{\kappa p_{\pm}}^{AA}=\begin{pmatrix}0&&-t^{-}_{\kappa
p_{\pm}}w&&t^{+}_{\kappa p_{\pm}}w^{\prime}\\\ 0&&t^{+}_{\kappa
p_{\pm}}w^{\prime}&&-t^{-}_{\kappa
p_{\pm}}w\end{pmatrix}\quad,\quad\hat{C}_{\kappa
p_{\pm},1}^{AB}=\begin{pmatrix}t^{+}_{\kappa
p_{\pm}}\phi_{0\bar{1}}\phi_{10}-t^{-}_{\kappa
p_{\pm}}\phi_{10}&&t^{+}_{\kappa p_{\pm}}w\phi_{10}&&-t^{-}_{\kappa
p_{\pm}}w^{\prime}\phi_{0\bar{1}}\phi_{10}\\\ t^{+}_{\kappa
p_{\pm}}\phi_{10}-t^{-}_{\kappa
p_{\pm}}\phi_{0\bar{1}}\phi_{10}&&-t^{-}_{\kappa
p_{\pm}}w^{\prime}\phi_{10}&&t^{+}_{\kappa
p_{\pm}}w\phi_{0\bar{1}}\phi_{10}\end{pmatrix},$ $\displaystyle\hat{C}_{\kappa
p_{\pm},2}^{AB}=\begin{pmatrix}t^{+}_{\kappa p_{\pm}}-t^{-}_{\kappa
p_{\pm}}\phi_{0\bar{1}}&&t^{+}_{\kappa p_{\pm}}\phi_{0\bar{1}}-t^{-}_{\kappa
p_{\pm}}\\\ t^{+}_{\kappa p_{\pm}}w^{\prime}\phi_{0\bar{1}}&&-t^{-}_{\kappa
p_{\pm}}w\phi_{0\bar{1}}\\\ -t^{-}_{\kappa p_{\pm}}w&&t^{+}_{\kappa
p_{\pm}}w^{\prime}\end{pmatrix}.$ (22)

##  Appendix B boundary Green Function and Dyson Equations

To simplify notation, we omit the superscript ‘$A$’ in advanced GFs from now
on. Given the real-space components of the bulk GF in Eq. (11), we next extend
the method of Refs. [68, 69, 60] to derive the bGF characterizing a _semi-
infinite_ 2D dimensional system. To that effect, we add an impurity potential
line $\epsilon$ localized at the frontier region [70]. Taking the limit
$\epsilon\to\infty$, the infinite system is cut into disconnected semi-
infinite subsystems with $j\leq-1$ (left side, $L$) and $j\geq 1$ (right side,
$R$). Using the Dyson equation, the local GF components of the cut subsystem
follow as [69]

$\hat{\mathcal{G}}_{jj}=\hat{G}^{(0)}_{jj}-\hat{G}^{(0)}_{j0}\left[\hat{G}^{(0)}_{00}\right]^{-1}\hat{G}^{(0)}_{0j},$
(23)

where $\hat{G}^{(0)}$ are the bulk GF and $\hat{\mathcal{G}}$ are the semi-
infinite GF. The bGF for the left and right semi-infinite system,
respectively, are with Eq. (23) given by

$\hat{\mathcal{G}}_{L}(k_{\parallel},\omega)=\hat{\cal
G}_{\bar{1}\bar{1}}(k_{\parallel},\omega)\quad,\quad\hat{\mathcal{G}}_{R}(k_{\parallel},\omega)=\hat{\cal
G}_{11}(k_{\parallel},\omega).$ (24)

Using this boundary Green function we can compute the spectral properties of
open systems but also the spectral and transport properties for different type
of junctions defined on MATBLG.

For instance, to determine the local spectral density at a domain wall [71,
72] where the coupling to the substrate changes sign, we use the following
expression

$\hat{G}_{nn}=\hat{\mathcal{G}}_{nn}+\hat{\mathcal{G}}_{n1}\hat{\Sigma}^{\dagger}_{LR}\left[\hat{\mathbb{I}}-\hat{\mathcal{G}}_{L}\hat{\Sigma}_{LR}\hat{\mathcal{G}}_{R}\hat{\Sigma}^{\dagger}_{LR}\right]^{-1}\hat{\mathcal{G}}_{L}\hat{\Sigma}_{LR}\hat{\mathcal{G}}_{1n},$
(25)

where $\hat{\Sigma}^{\dagger}_{LR}$ is the coupling term between the different
regions that the define wall and requires the non-local GF for the right semi-
infinite system given by

$\displaystyle\hat{\mathcal{G}}_{n1}=\hat{G}^{(0)}_{n1}-\hat{\mathcal{G}}_{R}\hat{T}^{\dagger}\hat{G}^{(0)}_{n1}\left[\hat{G}^{(0)}_{00}\right]^{-1}\hat{G}^{(0)}_{01},$
$\displaystyle\hat{\mathcal{G}}_{1n}=\hat{G}^{(0)}_{1n}-\hat{G}^{(0)}_{10}\left[\hat{G}^{(0)}_{00}\right]^{-1}\hat{G}^{(0)}_{1n}\hat{T}\hat{\mathcal{G}}_{R},$
(26)

where $\hat{T}$ is the coupling term between consecutive sites of each semi-
infinite system and

$\displaystyle\hat{G}^{(0)}_{1n}=\hat{G}^{(0)}_{00}\left(\hat{T}\hat{\mathcal{G}}_{R}\right)^{n-1},$
$\displaystyle\hat{G}^{(0)}_{n1}=\left(\hat{\mathcal{G}}_{R}\hat{T}^{\dagger}\right)^{n-1}\hat{G}^{(0)}_{00}.$
(27)

For computing transport properties in junctions we start from the general
expression for the current between two leads using Keldysh Green functions in
[73, 69], which we adapt to the 1 valley MATBLG case

$\displaystyle I_{LR}(k_{\parallel})=\frac{2e}{h}\int d\omega\;$ tr
$\displaystyle\\{\hat{\Sigma}_{LR}\hat{G}_{RR}^{-+}\hat{\Sigma}_{RL}\hat{\cal
G}_{L}^{+-}$ (28)
$\displaystyle-\hat{\Sigma}_{LR}\hat{G}_{RR}^{+-}\hat{\Sigma}_{RL}\hat{\cal
G}_{L}^{-+}\\},$

where the factor $2$ comes from the spin degree of freedom. Finally we expand
the non-equilibrium correlation functions in terms of advanced GF as using
Langreth rules

$\displaystyle\hat{G}^{+-/-+}_{RR}$ $\displaystyle=$
$\displaystyle(\hat{\mathbb{I}}+\hat{G}_{LR}^{\dagger}\hat{\Sigma}_{LR})\hat{\cal
G}^{+-/-+}_{R}(\hat{\mathbb{I}}+\hat{\Sigma}_{LR}^{\dagger}\hat{G}_{LR})+$
(29) $\displaystyle\hat{G}^{\dagger}_{RR}\hat{\Sigma}_{LR}^{\dagger}\hat{\cal
G}^{+-/-+}_{L}\hat{\Sigma}_{LR}\hat{G}_{RR},$

with

$\displaystyle\hat{\cal G}^{+-}_{j}=f_{j}(\omega-\mu_{j})(\hat{\cal
G}_{j}-\hat{\cal G}_{j}^{\dagger}),$ $\displaystyle\hat{\cal
G}^{+-}_{j}=\left(f_{j}(\omega-\mu_{j})-1\right)(\hat{\cal G}_{j}-\hat{\cal
G}_{j}^{\dagger}),$ (30)

$j=L/R$ and $f_{j}(\omega-\mu_{j})$ is the Fermi-Dirac distribution. Using
Dyson equation we can obtain the current in terms of the bGF

$\displaystyle\hat{G}_{LR}=\hat{\cal G}_{L}\hat{\Sigma}_{LR}\hat{G}_{RR},$
(31) $\displaystyle\hat{G}_{RR}=[\hat{\mathbb{I}}-\hat{\cal
G}_{R}\hat{\Sigma}_{LR}^{\dagger}\hat{\cal
G}_{L}\hat{\Sigma}_{LR}]^{-1}\hat{\cal G}_{R}$ .

After some algebra we get the Landauer formula

$I_{LR}(k_{\parallel})=\frac{2e}{h}\int
d\omega\;T(k_{\parallel},\omega)\left(f_{R}-f_{L}\right),$ (32)

where

$\displaystyle
T(k_{\parallel},\omega)=\textrm{tr}\left\\{\left(\hat{\mathbb{I}}+\hat{G}_{LR}^{\dagger}\hat{\Sigma}_{LR}\right)\left(\hat{\cal
G}_{R}-\hat{\cal
G}_{R}^{\dagger}\right)\left(\hat{\mathbb{I}}+\hat{\Sigma}_{LR}^{\dagger}\hat{G}_{LR}\right)\hat{\Sigma}_{LR}^{\dagger}\left(\hat{\cal
G}_{L}-\hat{\cal G}_{L}^{\dagger}\right)\right\\},$ $\displaystyle\langle
T(\omega)\rangle=\displaystyle\frac{1}{n_{k}}\sum\limits_{k_{\parallel}}T(k_{\parallel},\omega),$
(33)

where $T(k_{\parallel},\omega)$ defines a momentum resolved transmission
coefficient and $\langle T(\omega)\rangle$ is the normalized total
transmission function, which can be associated with a two-terminal zero
temperature conductance through $G=G_{0}\langle T(0)\rangle$ where
$G_{0}=2e^{2}/h$ is the conductance quantum.

##  Appendix C Recursive GF method

In cases where a non-uniform potential is applied the bGF can be calculated
recursively. We define the recursive the bGF at an dimensionless n-site as

$\displaystyle\left[\hat{\mathcal{G}}^{rc}_{R}(n)\right]^{-1}=\omega\hat{\mathbb{I}}-\hat{\mathcal{H}}_{0}(k_{\parallel})-\hat{\mathcal{H}}_{b}(n)-\Sigma_{R}(n),$
$\displaystyle\left[\hat{\mathcal{G}}^{rc}_{L}(n)\right]^{-1}=\omega\hat{\mathbb{I}}-\hat{\mathcal{H}}_{0}(k_{\parallel})-\hat{\mathcal{H}}_{b}(n)-\Sigma_{L}(n),$
(34)

where $\hat{\mathcal{H}}_{0}(k_{\parallel})$ is the local contribution defined
in a stripe and the recursive expression of the self energy takes the form

$\displaystyle\Sigma_{R}(n)$ $\displaystyle=$
$\displaystyle\hat{T}_{LR}\left[\hat{\mathcal{G}}^{rc}_{R}(n-1)\right]^{-1}\hat{T}^{\dagger}_{LR},$
$\displaystyle\Sigma_{L}(n)$ $\displaystyle=$
$\displaystyle\hat{T}^{\dagger}_{LR}\left[\hat{\mathcal{G}}^{rc}_{L}(n-1)\right]^{-1}\hat{T}_{LR}.$
(35)

As we can see the self-energy couples the n-site with the previous one and $n$
goes from $n=1$ to $n=N_{x}=L/2L_{m}$ where $L$ is the total length of the
system. The self energy at the first site $\Sigma_{L/R}(1)$ can be defined to
simulate the coupling to a doped normal lead with constant LDOS inducing a
characteristic broadening of the states through the junction.

Inspired by the experimental device described in Ref. [46] the potential
profile Hamiltonian for a three regions junction takes the form

$\displaystyle\hat{\mathcal{H}}_{b}(n)=\displaystyle\frac{U_{b}}{2}V_{x}(n)\hat{\mathbb{I}},$
$\displaystyle
V_{x}(n)=\tanh\left(\displaystyle\frac{n-(N_{x}-N_{b})-1}{d_{b}}\right)+1.$
(36)

For the 6B1V model and following the polynomial expansion in $z$ of the
Hamiltonian $\hat{\mathcal{H}}(z)$ as described in [45] takes the form

$\displaystyle\hat{\mathcal{H}}(z)$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{2m+1}\hat{H}_{i}z^{i-(m+1)}$ $\displaystyle=$
$\displaystyle\hat{H}_{1}z^{-m}+\dots+\hat{H}_{m+1}+\dots+\hat{H}_{2m+1}z^{m}.$

we can define the matrix that define the recursive method as

$\hat{\mathcal{H}}_{0}=\begin{pmatrix}\hat{H}_{2}&&\hat{H}_{3}\\\
\hat{H}_{1}&&\hat{H}_{2}\end{pmatrix}\quad,\quad\hat{T}_{LR}=\begin{pmatrix}\hat{0}&&\hat{H}_{1}\\\
\hat{0}&&\hat{0}\end{pmatrix}.$ (38)

We note that this expression is totally general and can be use for any TB
nearest neighbour Hamiltonian.

In the case of the 2B1V model, which goes beyond the nearest neighbors
approximation, we follow Ref. [46] to set the recursive equations from

$\displaystyle\hat{\mathcal{H}}_{0}=\begin{pmatrix}\hat{H}_{0}&&\hat{H}_{1}&&\hat{H}_{2}&&\hat{H}_{3}\\\
\hat{H}_{1}^{\dagger}&&\hat{H}_{0}&&\hat{H}_{1}&&\hat{H}_{2}\\\
\hat{H}_{2}^{\dagger}&&\hat{H}_{1}^{\dagger}&&\hat{H}_{0}&&\hat{H}_{1}\\\
\hat{H}_{3}^{\dagger}&&\hat{H}_{2}^{\dagger}&&\hat{H}_{1}^{\dagger}&&\hat{H}_{0}\end{pmatrix},$
$\displaystyle\hat{T}_{LR}=\begin{pmatrix}\hat{0}&&\hat{H}_{3}&&\hat{H}_{2}&&\hat{H}_{1}\\\
\hat{0}&&\hat{0}&&\hat{H}_{3}&&\hat{H}_{2}\\\
\hat{0}&&\hat{0}&&\hat{0}&&\hat{H}_{3}\\\
\hat{0}&&\hat{0}&&\hat{0}&&\hat{0}\end{pmatrix},$ (39)

where

$\displaystyle\hat{H}_{0}$ $\displaystyle=$
$\displaystyle\begin{pmatrix}2t_{2}^{\prime}\sin{k_{x}L_{x}}&&t_{1}e^{i\frac{k_{x}L_{x}}{2}}+t^{*}_{4}e^{-i\frac{k_{x}L_{x}}{2}}\\\
t_{1}e^{-i\frac{k_{x}L_{x}}{2}}+t_{4}e^{i\frac{k_{x}L_{x}}{2}}&&2t_{2}^{\prime}\sin{k_{x}L_{x}}\end{pmatrix},$
$\displaystyle\hat{H}_{1}$ $\displaystyle=$
$\displaystyle\begin{pmatrix}0&&t_{1}+t_{5}^{*}e^{ik_{x}L_{x}}\\\
t_{1}+t_{5}e^{-ik_{x}L_{x}}&&0\end{pmatrix},$ $\displaystyle\hat{H}_{2}$
$\displaystyle=$
$\displaystyle\begin{pmatrix}0&&t_{4}^{*}e^{i\frac{k_{x}L_{x}}{2}}+t_{5}^{*}e^{-i\frac{k_{x}L_{x}}{2}}\\\
t_{4}e^{-i\frac{k_{x}L_{x}}{2}}+t_{5}e^{i\frac{k_{x}L_{x}}{2}}&&0\end{pmatrix},$
$\displaystyle\hat{H}_{3}$ $\displaystyle=$
$\displaystyle\begin{pmatrix}-2t_{2}^{\prime}\sin{\frac{k_{x}L_{x}}{2}}&&t_{5}^{*}\\\
t_{5}&&-2t_{2}^{\prime}\sin{\frac{k_{x}L_{x}}{2}}\end{pmatrix}.$ (40)

However, we equivalently can express the relevant matrices for the recursive
method in terms of the polynomial expansion of the Hamiltonian.

## References

  * Carr  _et al._ [2017] S. Carr, D. Massatt, S. Fang, P. Cazeaux, M. Luskin, and E. Kaxiras, Physical Review B 95, 075420 (2017).
  * Ribeiro-Palau  _et al._ [2018] R. Ribeiro-Palau, C. Zhang, K. Watanabe, T. Taniguchi, J. Hone, and C. R. Dean, Science 361, 690 (2018).
  * 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).
  * Ma  _et al._ [2020] C. Ma, Q. Wang, S. Mills, X. Chen, B. Deng, S. Yuan, C. Li, K. Watanabe, T. Taniguchi, X. Du, _et al._ , Nano Letters 20, 6076 (2020).
  * Serlin  _et al._ [2020] M. Serlin, C. Tschirhart, H. Polshyn, Y. Zhang, J. Zhu, K. Watanabe, T. Taniguchi, L. Balents, and A. Young, Science 367, 900 (2020).
  * Wu  _et al._ [2021] S. Wu, Z. Zhang, K. Watanabe, T. Taniguchi, and E. Y. Andrei, Nature Materials , 1 (2021).
  * Stepanov  _et al._ [2020] P. Stepanov, M. Xie, T. Taniguchi, K. Watanabe, X. Lu, A. H. MacDonald, B. A. Bernevig, and D. K. Efetov, arXiv preprint arXiv:2012.15126 (2020).
  * Rodan-Legrain  _et al._ [2020] D. Rodan-Legrain, Y. Cao, J. M. Park, S. C. de la Barrera, M. T. Randeria, K. Watanabe, T. Taniguchi, and P. Jarillo-Herrero, arXiv preprint arXiv:2011.02500 (2020).
  * Guinea and Walet [2018] F. Guinea and N. R. Walet, Proceedings of the National Academy of Sciences 115, 13174 (2018).
  * Cea  _et al._ [2020] T. Cea, P. A. Pantaleón, and F. Guinea, Physical Review B 102, 155136 (2020).
  * Liu and Dai [2021] J. Liu and X. Dai, Physical Review B 103, 035427 (2021).
  * Wilhelm  _et al._ [2021] P. Wilhelm, T. C. Lang, and A. M. Läuchli, Physical Review B 103, 125406 (2021).
  * Bultinck  _et al._ [2020] N. Bultinck, S. Chatterjee, and M. P. Zaletel, Physical Review Letters 124, 166601 (2020).
  * Nuckolls  _et al._ [2020] K. P. Nuckolls, M. Oh, D. Wong, B. Lian, K. Watanabe, T. Taniguchi, B. A. Bernevig, and A. Yazdani, Nature , 1 (2020).
  * Stauber  _et al._ [2018] T. Stauber, T. Low, and G. Gómez-Santos, Physical Review B 98, 195414 (2018).
  * Kang and Vafek [2019] J. Kang and O. Vafek, Physical Review Letters 122, 246401 (2019).
  * Zhang  _et al._ [2020] M. Zhang, Y. Zhang, C. Lu, W.-Q. Chen, and F. Yang, Chinese Physics B 29, 127102 (2020).
  * Carr  _et al._ [2019] S. Carr, S. Fang, H. C. Po, A. Vishwanath, and E. Kaxiras, Physical Review Research 1, 033072 (2019).
  * Bistritzer and MacDonald [2011] R. Bistritzer and A. H. MacDonald, Proceedings of the National Academy of Sciences 108, 12233 (2011).
  * Thouless [1984] D. Thouless, Journal of Physics C: Solid State Physics 17, L325 (1984).
  * Soluyanov and Vanderbilt [2011] A. A. Soluyanov and D. Vanderbilt, Physical Review B 83, 035108 (2011).
  * Hejazi  _et al._ [2019] K. Hejazi, C. Liu, H. Shapourian, X. Chen, and L. Balents, Physical Review B 99, 035111 (2019).
  * Zou  _et al._ [2018] L. Zou, H. C. Po, A. Vishwanath, and T. Senthil, Physical Review B 98, 085435 (2018).
  * Po  _et al._ [2018a] H. C. Po, H. Watanabe, and A. Vishwanath, Physical Review Letters 121, 126402 (2018a).
  * Ahn  _et al._ [2019] J. Ahn, S. Park, and B.-J. Yang, Physical Review X 9, 021013 (2019).
  * Koshino  _et al._ [2018] M. Koshino, N. F. Yuan, T. Koretsune, M. Ochi, K. Kuroki, and L. Fu, Physical Review X 8, 031087 (2018).
  * Wang  _et al._ [2019] H.-X. Wang, G.-Y. Guo, and J.-H. Jiang, New Journal of Physics 21, 093029 (2019).
  * Bernevig  _et al._ [2020] B. A. Bernevig, Z. Song, N. Regnault, and B. Lian, arXiv preprint arXiv:2009.11301 (2020).
  * Wieder and Bernevig [2018] B. J. Wieder and B. A. Bernevig, arXiv preprint arXiv:1810.02373 (2018).
  * de Paz  _et al._ [2019] M. B. de Paz, M. G. Vergniory, D. Bercioux, A. García-Etxarri, and B. Bradlyn, Physical Review Research 1, 032005 (2019).
  * Kooi  _et al._ [2019] S. H. Kooi, G. Van Miert, and C. Ortix, Physical Review B 100, 115160 (2019).
  * Li  _et al._ [2020] Z. Li, H.-C. Chan, and Y. Xiang, Physical Review B 102, 245149 (2020).
  * Peri  _et al._ [2020] V. Peri, Z.-D. Song, M. Serra-Garcia, P. Engeler, R. Queiroz, X. Huang, W. Deng, Z. Liu, B. A. Bernevig, and S. D. Huber, Science 367, 797 (2020).
  * Po  _et al._ [2019] H. C. Po, L. Zou, T. Senthil, and A. Vishwanath, Physical Review B 99, 195455 (2019).
  * Hejazi  _et al._ [2021] K. Hejazi, X. Chen, and L. Balents, Physical Review Research 3, 013242 (2021).
  * Cao  _et al._ [2020] J. Cao, M. Wang, C.-C. Liu, and Y. Yao, arXiv preprint arXiv:2012.02575 (2020).
  * Kezilebieke  _et al._ [2020] S. Kezilebieke, M. N. Huda, V. Vaňo, M. Aapro, S. C. Ganguli, O. J. Silveira, S. Głodzik, A. S. Foster, T. Ojanen, and P. Liljeroth, Nature 588, 424 (2020).
  * Song  _et al._ [2020] Z.-D. Song, L. Elcoro, and B. A. Bernevig, Science 367, 794 (2020).
  * Pierce  _et al._ [2021] A. T. Pierce, Y. Xie, J. M. Park, E. Khalaf, S. H. Lee, Y. Cao, D. E. Parker, P. R. Forrester, S. Chen, K. Watanabe, _et al._ , arXiv preprint arXiv:2101.04123 (2021).
  * Choi  _et al._ [2021] Y. Choi, H. Kim, C. Lewandowski, Y. Peng, A. Thomson, R. Polski, Y. Zhang, K. Watanabe, T. Taniguchi, J. Alicea, _et al._ , arXiv preprint arXiv:2102.02209 (2021).
  * Zondiner  _et al._ [2020] U. Zondiner, A. Rozen, D. Rodan-Legrain, Y. Cao, R. Queiroz, T. Taniguchi, K. Watanabe, Y. Oreg, F. von Oppen, A. Stern, _et al._ , Nature 582, 203 (2020).
  * Alvarado  _et al._ [2020] M. Alvarado, A. Iks, A. Zazunov, R. Egger, and A. L. Yeyati, Physical Review B 101, 094511 (2020).
  * Householder [2013] A. S. Householder,  _The theory of matrices in numerical analysis_ (Courier Corporation, 2013).
  * [45] M. Alvarado and A. L. Yeyati, In preparation .
  * De Beule  _et al._ [2020] C. De Beule, P. G. Silvestrov, M.-H. Liu, P. Recher, _et al._ , Physical Review Research 2, 043151 (2020).
  * Fujimoto and Koshino [2020] M. Fujimoto and M. Koshino, arXiv preprint arXiv:2012.02537 (2020).
  * Wu  _et al._ [2020] X. Wu, W. Hanke, M. Fink, M. Klett, and R. Thomale, Physical Review B 101, 134517 (2020).
  * Chichinadze  _et al._ [2020] D. V. Chichinadze, L. Classen, and A. V. Chubukov, Physical Review B 102, 125120 (2020).
  * Song  _et al._ [2019] Z. Song, Z. Wang, W. Shi, G. Li, C. Fang, and B. A. Bernevig, Physical Review Letters 123, 036401 (2019).
  * Choi  _et al._ [2019] Y. Choi, J. Kemmer, Y. Peng, A. Thomson, H. Arora, R. Polski, Y. Zhang, H. Ren, J. Alicea, G. Refael, _et al._ , Nature Physics 15, 1174 (2019).
  * Alexandradinata  _et al._ [2014] A. Alexandradinata, X. Dai, and B. A. Bernevig, Physical Review B 89, 155114 (2014).
  * Alexandradinata  _et al._ [2016] A. Alexandradinata, Z. Wang, and B. A. Bernevig, Physical Review X 6, 021008 (2016).
  * Ahn  _et al._ [2018] J. Ahn, D. Kim, Y. Kim, and B.-J. Yang, Physical Review Letters 121, 106403 (2018).
  * Bradlyn  _et al._ [2019] B. Bradlyn, Z. Wang, J. Cano, and B. A. Bernevig, Physical Review B 99, 045140 (2019).
  * Park  _et al._ [2019] M. J. Park, Y. Kim, G. Y. Cho, and S. Lee, Physical review letters 123, 216803 (2019).
  * Park  _et al._ [2021] M. J. Park, S. Jeon, S. Lee, H. C. Park, and Y. Kim, Carbon 174, 260 (2021).
  * Liu  _et al._ [2021] B. Liu, L. Xian, H. Mu, G. Zhao, Z. Liu, A. Rubio, and Z. Wang, Physical Review Letters 126, 066401 (2021).
  * Po  _et al._ [2018b] H. C. Po, L. Zou, A. Vishwanath, and T. Senthil, Physical Review X 8, 031089 (2018b).
  * Burset  _et al._ [2008] P. Burset, A. L. Yeyati, and A. Martín-Rodero, Phys. Rev. B 77, 205425 (2008).
  * Burset  _et al._ [2009] P. Burset, W. Herrera, and A. Levy Yeyati, Phys. Rev. B 80, 041402 (2009).
  * Herrera  _et al._ [2010] W. J. Herrera, P. Burset, and A. L. Yeyati, Journal of Physics: Condensed Matter 22, 275304 (2010).
  * Guo  _et al._ [2018] H. Guo, X. Zhu, S. Feng, and R. T. Scalettar, Physical Review B 97, 235453 (2018).
  * Ray  _et al._ [2019] S. Ray, J. Jung, and T. Das, Physical Review B 99, 134515 (2019).
  * Christos  _et al._ [2020] M. Christos, S. Sachdev, and M. S. Scheurer, Proceedings of the National Academy of Sciences 117, 29543 (2020).
  * Peri  _et al._ [2021] V. Peri, Z.-D. Song, B. A. Bernevig, and S. D. Huber, Physical review letters 126, 027002 (2021).
  * Casas  _et al._ [2019] O. E. Casas, S. G. Páez, A. L. Yeyati, P. Burset, and W. J. Herrera, Physical Review B 99, 144502 (2019).
  * Arrachea  _et al._ [2009] L. Arrachea, G. S. Lozano, and A. Aligia, Physical Review B 80, 014425 (2009).
  * Zazunov  _et al._ [2016] A. Zazunov, R. Egger, and A. L. Yeyati, Physical Review B 94, 014502 (2016).
  * Pinon  _et al._ [2020] S. Pinon, V. Kaladzhyan, and C. Bena, Physical Review B 101, 115405 (2020).
  * Kwan  _et al._ [2020] Y. H. Kwan, G. Wagner, N. Chakraborty, S. H. Simon, and S. Parameswaran, arXiv preprint arXiv:2007.07903 (2020).
  * Padhi  _et al._ [2020] B. Padhi, A. Tiwari, T. Neupert, and S. Ryu, Physical Review Research 2, 033458 (2020).
  * Cuevas  _et al._ [1996] J. Cuevas, A. Martín-Rodero, and A. L. Yeyati, Physical Review B 54, 7366 (1996).