# Network model for periodically strained graphene

Christophe De Beule  Department of Physics and Materials Science, University
of Luxembourg, L-1511 Luxembourg, Luxembourg  Department of Physics and
Astronomy, University of Pennsylvania, Philadelphia PA 19104  Võ Tiến Phong
Department of Physics and Astronomy, University of Pennsylvania, Philadelphia
PA 19104  E. J. Mele  Department of Physics and Astronomy, University of
Pennsylvania, Philadelphia PA 19104

###### Abstract

The long-wavelength physics of monolayer graphene in the presence of periodic
strain fields has a natural chiral scattering network description. When the
strain field varies slowly compared to the graphene lattice and the effective
magnetic length of the induced valley pseudomagnetic field, the low-energy
physics can be understood in terms of valley-polarized percolating domain-wall
modes. Inspired by a recent experiment, we consider a strain field with
threefold rotation and mirror symmetries but without twofold rotation
symmetry, resulting in a system with the connectivity of the oriented kagome
network. Scattering processes in this network are captured by a symmetry-
constrained phenomenological $S$ matrix. We analyze the phase diagram of the
kagome network, and show that the bulk physics of the strained graphene can be
qualitatively captured by the network when the $S$ matrix gains an appropriate
energy dependence. We also discuss the limitations of this approach to
properly account for boundary physics.

##  I Introduction

Superlattice engineering in low-dimensional materials by periodic modulations
that vary slowly compared to the microscopic lattice has been a highly
successful tool for creating artificial crystals that host new emergent
phenomena. One route is to impose periodic potentials by patterning
electrostatic gates Forsythe  _et al._ (2018) or by inducing periodic strain
fields with a substrate Mao  _et al._ (2020) or nano-pillars Jiang  _et al._
(2017). Another platform are moiré materials Andrei  _et al._ (2021), which
are layered materials subjected to rotational mismatch or heterostrain,
inducing a spatial modulation of the interlayer coupling. The archetypal
example being twisted bilayer graphene Lopes dos Santos  _et al._ (2007) where
the moiré structure induces isolated nearly-flat bands at the magic angle that
give rise to a rich phenomenology of correlated physics Cao  _et al._ (2018a,
b). Similar spectrally isolated and flattened minibands can also emerge in
monolayer graphene subjected to a periodic strain field Mao  _et al._ (2020);
Milovanović  _et al._ (2020); Manesco  _et al._ (2021); Manesco and Lado
(2021); Giambastiani  _et al._ (2022); Phong and Mele (2022).

Figure 1: (a) Pseudomagnetic field [Eq. (1) for $\phi=0$] superimposed on the
graphene lattice for $N=10$. (b) Nodal structure of (a). The PMF is positive
(negative) in the orange (green) regions and vanishes on the thick black
loops. Each loop supports a chiral mode which can scatter to neighboring
loops. Some processes are indicated by the dashed arrows. Here, dashed
hexagons highlight the connectivity of the nodal structure with the corners
acting as scattering centers (dots). (c) Reciprocal lattice for $N=4$, where
the large black hexagon is the original BZ of graphene, and the small orange
hexagons give the SBZ in the extended zone scheme. (d) Effective piecewise
constant PMF (kagome tiling).

In graphene, a strain field couples as a vector potential to the long-
wavelength Dirac excitations with an opposite sign for the two valleys Suzuura
and Ando (2002); Mañes (2007); Vozmediano  _et al._ (2010); Guinea  _et al._
(2010); Levy  _et al._ (2010). Hence, a periodic strain field in graphene
results in a periodic pseudomagnetic field (PMF), as shown in Fig. 1(a).
Wherever the PMF changes sign, which must occur in cases of zero spatially-
averaged flux, there are real-space valley Hall transitions. For example, the
PMF shown in Fig. 1(a) vanishes on disjoint closed loops illustrated in Fig.
1(b). Furthermore, for a PMF that varies slowly on the scale of the graphene
lattice and the effective magnetic length, it is meaningful to consider the
local Landau levels of regions separated by nodal lines. If the chemical
potential lies between the $m$th and $(m+1)$th local Landau level, the net
local valley Chern number is given by $\pm(m+1)/2$ since the lowest Landau
level of graphene only contributes Chern number 1/2 per valley Novoselov  _et
al._ (2005). Therefore, the local valley Chern number changes by an integer
amount when the PMF changes sign. It immediately follows that the nodal lines
of the PMF form domain walls that host chiral modes. Due to time-reversal
symmetry, these chiral modes propagate in opposite directions for different
valleys. However, intervalley backscattering is suppressed when the PMF varies
slowly on the scale of the graphene lattice. Furthermore, the bound chiral
modes can become dispersive by tunneling to neighboring loops. In our example,
the corresponding scattering centers form an effective kagome lattice, as
shown in Fig. 1(b). The single-valley low-energy physics of periodically
strained graphene can thus be understood in terms of percolating chiral modes
with the connectivity of an oriented scattering network. Moreover, for long-
wavelength physics, the smooth PMF can be replaced by an effective piecewise
constant field [see Fig. 1(d)]. This piecewise constant field has the same
symmetry and connectivity as the smooth PMF.

The above observations illustrate a more general paradigm in which the
superlattice modulates the local band topology to give to rise to a Chern
mosaic that hosts a chiral scattering network. Similar topological scattering
networks occur in moiré materials such as minimally twisted bilayer graphene
subjected to a perpendicular electric field San-Jose and Prada (2013); Efimkin
and MacDonald (2018); Fleischmann  _et al._ (2020); Tsim  _et al._ (2020);
Chou  _et al._ (2020); De Beule  _et al._ (2020); De Beule  _et al._ (2021a);
Vakhtel  _et al._ (2022) (triangular chiral network), as well as certain large
twist-angle graphene bilayers Pal  _et al._ (2019) (helical honeycomb
network), and double-aligned graphene-hexagonal boron nitride moirés Moulsdale
_et al._ (2022) (chiral kagome network).

In this work, we analyze the scattering network appropriate for periodically
strained graphene with zero total flux per valley. We find that the network
model captures many attributes of the spectral structure obtained from the
valley-projected Hamiltonian.

This paper is organized as follows. In Section II, we introduce our model for
the periodic pseudomagnetic field in graphene and calculate the resulting band
structure with a valley-projected continuum theory. We then argue that the
low-energy physics can be understood in terms of an effective chiral kagome
network. The phenomenological network model is developed in Section III where
we discuss its phase diagram and corroborate our results by finite-width
ribbon calculations. Finally, in Section IV, we demonstrate that the network
model can qualitatively reproduce the bulk low-energy physics of graphene
subjected to a periodic strain field and highlight its limitations in
describing boundary physics.

##  II Periodically strained graphene

To investigate monolayer graphene subjected to periodic strain, we consider a
pseudomagnetic field (PMF) with zero net flux that is commensurate with the
graphene lattice. Specifically, we consider the case with threefold rotation
symmetry ($\mathcal{C}_{3z}$) about the origin [see Fig. 1(a)] and mirror
symmetry across the $y$ axis ($\mathcal{M}_{x}$), but with broken twofold
rotation symmetry ($\mathcal{C}_{2z}$). Together with translations, these
symmetries form the wallpaper group 14 ($p3m1$) with point group $C_{3v}$. The
strain field thus reduces the graphene point group from $C_{6v}$ to $C_{3v}$,
which is the case for graphene on top of NbSe2 Mao  _et al._ (2020).
Alternatively, one can engineer this strain field using an artificial
substrate with a spatially-varying periodic height profile Jiang  _et al._
(2017); Phong and Mele (2022). In the first-star approximation, a PMF with
zero net flux that satisfies these symmetry constraints can be written as

$\bm{B}_{\nu}(\bm{r})=\bm{e}_{z}\nu
B_{0}\sum_{i=1}^{3}\cos\left(\bm{G}_{i}\cdot\bm{r}+\phi\right),$ (1)

where we take $B_{0}>0$ for concreteness, and which is shown in Fig. 1(a).
Note that Eq. (1) conserves time-reversal symmetry ($\mathcal{T}$) since it
has an opposite sign in valleys $K_{+}$ and $K_{-}$. We further note that
$\mathcal{C}_{2z}$ is conserved for $\phi=\pi/2$ and maximally broken for
$\phi=0$. Since we are interested in the latter, we set $\phi=0$ in the
remainder of this work. Here, $\bm{G}_{i}$ ($i=1,2,3$) are the three primitive
superlattice reciprocal lattice vectors related by $\mathcal{C}_{3z}$ symmetry
[see Fig. 1(c)] with $\bm{G}_{i}\cdot\bm{L}_{j}=2\pi\delta_{ij}$ ($i,j=1,2$)
and where

$\bm{L}_{1}=L(\sqrt{3}/2,1/2),\qquad\bm{L}_{2}=L(-1,0),$ (2)

are the superlattice lattice vectors with $L=Na,$ $N>1$ a positive integer,
and $a$ the microscopic lattice constant of graphene [see Fig. 1(a)]. Because
the superlattice has a larger periodicity than the original graphene lattice,
the graphene Brillouin zone (BZ) is folded onto a smaller superlattice BZ
(SBZ) that fits $N^{2}$ times into the graphene BZ. This is illustrated in
Fig. 1(c). For $N\in 3\mathbb{N}$, the $K_{+}$ and $K_{-}$ points are folded
on top of each other and the system is an insulator, otherwise the system is a
semimetal Phong and Mele (2022). However, for $N\gg 1$, this distinction
becomes moot due to an emergent conservation of valley charge. In this case,
the PMF varies slowly on the scale of the graphene lattice such that the
valleys are effectively decoupled, and one can use an effective valley-
projected continuum theory.

###  II.1 Valley-projected theory

The valley-projected long-wavelength Hamiltonian in the presence of a
pseudovector potential is given by

$H_{\nu}=\frac{\hbar v_{F}}{V}\int
d^{2}\bm{r}\,\psi_{\nu}^{\dagger}\left[\left(-i\nabla+\frac{e\bm{A}_{\nu}}{\hbar}\right)\cdot\left(\nu\sigma_{x},\sigma_{y}\right)\right]\psi_{\nu},$
(3)

where $\psi_{\nu}(\bm{r})=(\psi_{\nu A}(\bm{r}),\psi_{\nu B}(\bm{r}))^{t}$ are
the field operators for sublattice $A$ and $B$, and $\sigma_{x,y}$ are Pauli
matrices that act in sublattice space. In dimensionless units, the only free
parameter is given by a ratio of length scales $L/l_{0}$ with
$l_{0}=\sqrt{\hbar/eB_{0}}$ the effective magnetic length at the maxima of the
PMF.

The symmetries of the single-valley theory are generated by
$\\{\mathcal{C}_{3z},\mathcal{M}_{x}\mathcal{T}\\}$ yielding the magnetic
point group $3m^{\prime}$. Note that twofold rotation $\mathcal{C}_{2z}$ is
not conserved since it sends
$\bm{B}_{\nu}(\bm{r})\mapsto-\bm{B}_{\nu}(\bm{r})$. In addition, there is a
chiral symmetry due to the absence of terms proportional to $\sigma_{0}$ or
$\sigma_{z}$ in Eq. (3). This is a model-dependent symmetry since such terms
are allowed and can be induced by, e.g., a strain-induced pseudo-electrostatic
potential or a constant sublattice-staggering term. However, in this work, we
will not consider these terms as they do not qualitatively change our results.
Since we consider the case where the net flux of the PMF vanishes, the
Hamiltonian (3) can be easily diagonalized in momentum space. We refer to App.
A for more details on the continuum model and its symmetries.

Figure 2: Band structure for valley $K_{+}$ (blue) and $K_{-}$ (red) of the
valley-projected theory along high-symmetry lines of the SBZ [green path in
Fig. 1(c)] for (a) $L/l_{0}=5$ and (b) $L/l_{0}=8$, given in units of $\hbar
v_{F}k_{0}$ where $k_{0}=|\bar{K}_{\nu}|=4\pi/3L$.

The resulting band structure is shown in Fig. 2 along high-symmetry lines of
the SBZ for two realistic values of $L/l_{0}\gg 1$, which is our regime of
interest. Owing to the chiral symmetry, the bands are symmetric about zero
energy and we label them as
$E_{n,\nu}(\bm{k})=\operatorname{sgn}(n)E_{|n|}(\nu\bm{k})$ with $n$ a nonzero
integer. Moreover, one can show that in the presence of chiral symmetry, the
bands near charge neutrality for a given valley always degenerate in pairs at
zero energy Phong and Mele (2022). However, this crossing can be avoided by
symmetry-allowed terms, but such terms are small and do not affect our results
qualitatively. Since $\mathcal{C}_{2z}$ is broken and the valley projection
implicitly breaks $\mathcal{T}$, the bands can have a nonzero Chern number. We
find that the two bands near charge neutrality carry a net Chern number Phong
and Mele (2022),

$\mathcal{C}_{+1}+\mathcal{C}_{-1}=\nu,$ (4)

while all remote low-energy bands that are shown in Fig. 2 are Chern trivial.
We note that this forms an obstruction to an effective lattice model for the
single-valley physics. The low-energy band manifold can thus be thought of as
a single Landau level that is half filled at charge neutrality. This situation
bears resemblance to the system of monolayer graphene in the presence of a
constant, real magnetic field. In that case, the $n=0$ Landau levels are also
half-filled at charge neutrality. Invoking chiral symmetry, one can argue that
the Hall conductivity must be the same in magnitude but opposite in sign in
gaps related by $E\mapsto-E$. This leads to the non-trivial prediction that
the gaps right above and below the $n=0$ Landau level must have Hall
conductivity equal to $\pm\frac{1}{2}\frac{e^{2}}{h}$ per valley and spin
Gusynin and Sharapov (2005); Peres  _et al._ (2006); Zhang  _et al._ (2005).
This half-integer quantum Hall effect stems from the $\pi$-quantized Berry
phase of the Dirac cones, and is related to the parity anomaly Haldane (1988);
Schakel (1991); Fialkovsky and Vassilevich (2012); Lapa (2019). Anomaly
cancellation is obtained when both valleys are simultaneously considered.

We also show the density $\rho_{A/B}(\bm{r})$ of the lowest conduction band
$(n=1)$ in Fig. 3. Here, we omit the valley index since the density is the
same for both valleys due to time-reversal symmetry. We observe that the
charge density on sublattice $B$ is localized with support on a triangular
lattice corresponding to the maxima of the PMF, while on sublattice $A$ it is
extended with support on the PMF minima which together form a honeycomb
lattice. Note also that $\rho_{A}$ breaks $\mathcal{C}_{2z}$ symmetry. The
separation of the density of the two sublattices can be understood in analogy
to the zeroth Landau level of graphene in a constant magnetic field, whose
wave function only has support on one of the two sublattices, depending on the
sign of the magnetic field and the valley Castro Neto  _et al._ (2009). For a
PMF, the valleys experience an opposite field and thus the support is
identical for both valleys.

Figure 3: Density of the lowest conduction band ($n=1$) for $L/l_{0}=5$
[energy bands shown in Fig. 2(a)] for sublattice $A$ ($\rho_{A}$, left) and
$B$ ($\rho_{B}$, right). The color scale is relative to the maximum which is a
factor $3$ larger for sublattice $B$.

###  II.2 Chiral interface modes

Next, we consider the nodal lines of the PMF, which are given here by disjoint
closed loops that encircle the maxima of the PMF, as shown in Fig. 1(b). In
the limit $L\gg l_{0}$, we can define a local valley Chern number inside and
outside of these loops. For example, if the chemical potential lies in between
the $m$th and $(m+1)$th local Landau levels, the net local valley Chern number
is given by $\operatorname{sgn}(B_{\nu})(m+1)/2$. The valley Chern number thus
changes by an integer when the PMF changes sign. This results in periodic
array of domain walls given by the nodal lines of the PMF, that host well-
defined chiral modes in the limit $L\gg l_{0}$. Moreover, the chiral modes can
percolate by scattering to neighboring loops, as illustrated in Fig. 1(b),
where the scattering centers form an effective kagome lattice. The low-energy
physics can thus be understood in terms of chiral modes bound to the nodal
lines of the PMF that percolate through the system with the connectivity of
the oriented kagome network.

Figure 4: (a) Spectrum of a pseudomagnetic domain wall [see Eq. (6)] for
valley $K_{+}$ (blue) and $K_{-}$ (red) where $\omega_{t}=\sqrt{2}v_{F}/l_{t}$
with $l_{t}=\sqrt{\hbar/e|B_{t}|}$. The column on the right shows the Landau
level and Chern number sequence in the hexagons (orange) and triangles (green)
of the kagome tiling [see Fig. 1(d)] as well as the number of chiral domain-
wall modes (arrows). (b) Normalized single-particle density for the state
indicated by an arrow in (a).

Since we consider long-wavelength physics, we can further replace the smooth
PMF by an effective piecewise constant PMF with the same symmetry and
connectivity. This yields a kagome tiling, as shown in Fig. 1(d). The values
of the constant field in the hexagons and triangles of the kagome tiling are
determined by the condition that the net flux should vanish,

$\Phi_{h}+2\Phi_{t}=0,$ (5)

which yields $B_{t}=-3B_{h}$. In terms of the smooth field, we find
$B_{h}=4B_{0}/\pi^{2}$. The corresponding sequence of Chern numbers and
interface modes for the kagome tiling are shown in Fig. 4(a). To demonstrate
the existence of these chiral modes explicitly, we performed a simple model
calculation with a piecewise constant pseudomagnetic field $\bm{B}_{\nu}=\nu
B(x)\bm{e}_{z}$ where

$B(x)=B_{t}\Theta(-x)+B_{h}\Theta(x),$ (6)

which is solved by the _ansatz_ $\psi(\bm{r})=e^{iky}\phi(x)$. A normalizable
solution can then be found in each region and the spectrum is obtained by
matching the spinor wave function on the boundaries. In Fig. 4(a), we show the
spectrum versus the transverse momentum $k$. Note that the net number of
chiral modes matches the sequence of Chern numbers. We show the single-
particle density for a given state in Fig. 4(b). As before, $\rho_{A}$
($\rho_{B}$) is mostly localized in the region with negative (positive) field.

##  III Network model

In the previous section, we argued that the single-valley low-energy physics
of periodically strained graphene in the limit $L\gg l_{0}$ can be understood
in terms of a network of chiral modes. Moreover, at low energies, each link of
the network hosts a single chiral mode for a given valley. In this section, we
develop a phenomenological theory for the network, by assuming that all
microscopic details can be absorbed into an $S$ matrix that relates incoming
and outgoing amplitudes at each scattering center, as illustrated in Fig.
1(b). The $S$ matrix is then only constrained by unitarity and the valley-
preserving symmetries: $\mathcal{C}_{3z}$ and $\mathcal{M}_{x}\mathcal{T}$.

Figure 5: Unit cell of the chiral kagome network. Here, arrows represent
chiral modes propagating along oriented links and the three colored dots
correspond to the three subnodes. Figure 6: Network energy bands for valley
$K_{+}$ (blue) and $K_{-}$ (red) in the fundamental domain along high-symmetry
lines of the SBZ [green path in Fig. 1(c)] for constant scattering parameters
$(\varphi,P)$ as indicated on Fig. 7.

###  III.1 Scattering matrix

The unit cell of the chiral kagome network is shown in Fig. 5. It consists of
three scattering nodes which we refer to as _subnodes_. Each cell can be
labeled by a superlattice lattice vector
$\bm{R}=n_{1}\bm{L}_{1}+n_{2}\bm{L}_{2}$ with $n_{1,2}$ integers. Since the
$S$ matrix only relates local incoming and outgoing scattering amplitudes, the
total $S$ matrix of the unit cell is block diagonal in the subnodes:

$\mathcal{S}=\begin{pmatrix}S_{1}&0&0\\\ 0&S_{2}&0\\\ 0&0&S_{3}\end{pmatrix},$
(7)

such that $b_{\bm{R}}=\mathcal{S}a_{\bm{R}}$ where

$\displaystyle a$
$\displaystyle=\left(a_{11},a_{12},a_{21},a_{22},a_{31},a_{32}\right)^{t},$
(8) $\displaystyle b$
$\displaystyle=\left(b_{11},b_{12},b_{21},b_{22},b_{31},b_{32}\right)^{t},$
(9)

are the incoming and outgoing amplitudes of one cell, respectively, and which
are defined in Fig. 5. Threefold rotation symmetry $\mathcal{C}_{3z}$ yields
$S_{0}=S_{1}=S_{2}=S_{3}$ and $\mathcal{M}_{x}\mathcal{T}$ gives
$S_{0}=(S_{0})^{t}$. We thus have

$S_{0}=e^{iE_{0}}\begin{pmatrix}e^{i\varphi}\sqrt{P}&\sqrt{1-P}\\\
\sqrt{1-P}&-e^{-i\varphi}\sqrt{P}\end{pmatrix},$ (10)

where $0\leq P\leq 1$, and $\varphi$ and $E_{0}$ are phases. Here, $P$ is the
probability for a left turn and $\varphi$ is the relative phase difference of
the amplitudes for a left and right turn, as illustrated in Fig. 5. A detailed
derivation of the $S$-matrix is given in App. B.

###  III.2 Link phases

Neighboring scattering nodes in the network are connected by oriented links of
length $L/2$ along which chiral modes freely propagate. During propagation
between adjacent nodes, they acquire a phase $\lambda=\exp\left(-i2\pi
E/E_{L}\right)$ with $E_{L}=2hv_{F}/L$. Here, we have assumed that the modes
have linear dispersion with slope approximately given by $v_{F}$, see Fig.
4(a). Explicitly,

$\displaystyle a_{11\bm{R}}$ $\displaystyle=\lambda b_{21\bm{R}-\bm{L}_{3}},$
(11) $\displaystyle a_{12\bm{R}}$ $\displaystyle=\lambda b_{22\bm{R}},$ (12)
$\displaystyle a_{21\bm{R}}$ $\displaystyle=\lambda b_{31\bm{R}-\bm{L}_{2}},$
(13) $\displaystyle a_{22\bm{R}}$ $\displaystyle=\lambda b_{32\bm{R}},$ (14)
$\displaystyle a_{31\bm{R}}$ $\displaystyle=\lambda b_{11\bm{R}-\bm{L}_{1}},$
(15) $\displaystyle a_{32\bm{R}}$ $\displaystyle=\lambda b_{12\bm{R}}.$ (16)

Since the network has translation symmetry, we Fourier transform to momentum
space,

$b_{\tau j\bm{R}}=\frac{1}{\sqrt{N}}\sum_{\bm{k}}e^{i\bm{k}\cdot\bm{R}}b_{\tau
j\bm{k}},$ (17)

and similarly for incoming amplitudes. Note that we work in Bloch form. Here,
$\tau=1,2,3$ is the subnode index, and $j=1,2$ is similar to an orbital index.
We then obtain

$a_{\bm{k}}=\lambda\mathcal{M}_{\bm{k}}b_{\bm{k}},$ (18)

where

$\mathcal{M}_{\bm{k}}=\begin{pmatrix}0&0&e^{-ik_{3}}&0&0&0\\\ 0&0&0&1&0&0\\\
0&0&0&0&e^{-ik_{2}}&0\\\ 0&0&0&0&0&1\\\ e^{-ik_{1}}&0&0&0&0&0\\\
0&1&0&0&0&0\end{pmatrix},$ (19)

is the connectivity matrix with $k_{i}=\bm{k}\cdot\bm{L}_{i}$ ($i=1,2,3$) and
$\bm{L}_{3}=-(\bm{L}_{1}+\bm{L}_{2})$. Since the links always connect
different subnodes, $\mathcal{M}_{\bm{k}}$ is block antidiagonal.

###  III.3 Network bands

Combining $b_{\bm{k}}=\mathcal{S}a_{\bm{k}}$ with (18) yields

$\mathcal{U}_{\bm{k}}b_{\bm{k}}=e^{i2\pi E/E_{L}}b_{\bm{k}},$ (20)

where $\mathcal{U}_{\bm{k}}=\mathcal{S}\mathcal{M}_{\bm{k}}$. Solving for the
energy gives

$E_{n,s}(\bm{k})=\frac{E_{L}}{3}\left[n+s\,\frac{\arccos
f(\bm{k})}{2\pi}\right],$ (21)

with $s=\pm 1$ and integer $n$, and where we set $E_{0}=-\pi/6$ to obtain
symmetric bands. The result for the other valley is obtained by letting
$\bm{k}\rightarrow-\bm{k}$. Here, we have also defined

$f(\bm{k})=\sqrt{P}\left[\left(1-P\right)\sum_{i=1}^{3}\sin\left(k_{i}+\varphi\right)+P\sin
3\varphi\right].$ (22)

Note that shifting the origin to $\bar{K}_{\pm}$ is equivalent to sending
$\varphi\rightarrow\varphi\pm 2\pi/3$ in Eq. (22), respectively, such that all
unique cases are contained in $\varphi\in[-\pi/3,\pi/3[$.

When the scattering parameters $P$ and $\varphi$ are constant, the energy
enters only in the link phases and Eq. (20) is periodic in energy with period
$E_{L}$. However, the network energy bands (21), which are shown in Fig. 6,
have a smaller period $E_{L}/3$. This can be understood from the phase-
rotation symmetry Delplace (2020); Delplace  _et al._ (2017),

$D\mathcal{U}_{\bm{k}}D^{-1}=e^{-i2\pi/3}\mathcal{U}_{\bm{k}},$ (23)

where

$D=\textrm{diag}\left(1,e^{i2\pi/3},e^{i4\pi/3}\right)\otimes\mathds{1}_{2},$
(24)

such that $Db_{n,s}$ is an eigenstate with energy $E_{n+1,s}=E_{n,s}+E_{L}/3$.
Note that $D^{3}=1$, corresponding to one full period and $b_{n+3,s}=b_{n,s}$.
Hence, we only need to consider two bands, called the fundamental domain,
chosen here as $-E_{L}/6<E\leq E_{L}/6$. However, when $\mathcal{S}$ is energy
dependent, Eq. (20) becomes a nonlinear eigenvalue problem, which simplifies
here to a transcendental equation with $f\rightarrow f(E)$. Finally, we note
that (23) is a result of the cyclic structure of the network. It is
independent of symmetry constraints, and holds for any subnode block diagonal
$S$ matrix.

Figure 7: Phase diagram of the chiral kagome network. White (colored) regions
are trivial (Chern) gapped phases, where the Chern number $\mathcal{C}_{+}$ of
the band $E_{n,+}$ is shown. The curves correspond to gap closings at high-
symmetry points and the dots indicate the parameters for the bands shown in
Fig. 6.

###  III.4 Phase diagram

The phase diagram of the oriented kagome network with a single channel has
been discussed previously in the context of photonic systems Pasek and Chong
(2014). Here, we discuss the phase diagram from a different perspective and
give some new analytical results.

We start by identifying the gap closings between pairs of bands $E_{n,\pm}$
and between $E_{n,+}$ and $E_{n+1,-}$. In general, the gap closes for
$f(\bm{k})^{2}=1$ which occurs only at the high-symmetry points $\bar{\Gamma}$
and $\bar{K}_{\pm}$. This yields

$P_{\bar{\Gamma}}=\frac{1}{4\sin^{2}\varphi},\qquad
P_{\bar{K}_{\pm}}=\frac{1}{4\sin^{2}\left(\varphi\pm 2\pi/3\right)},$ (25)

which define gap-closing lines in the $(\varphi,P)$ plane, as shown in Fig. 7.

Now consider the two flatband limits of the chiral kagome network, namely
$P=0$ and $P=1$. In these limits, the bands are isolated and correspond to
classical loop configurations, such that the Chern number vanishes. This can
also be understood from a strong phase-rotation symmetry in the flatband
limit, which reduces the fundamental domain to a single band Delplace (2020);
Delplace  _et al._ (2017). Since the Chern number is invariant under unitary
transformations, it follows together with completeness that the sum of the
Chern numbers in each fundamental domain vanishes. At a strong-phase rotation
symmetric point, all bands are therefore Chern trivial. In the following, we
refer to the phases that are adiabatically connected to the $P=0$ and $P=1$
flatband limit as the $H$ and $T$ phase, respectively. Indeed, for $P=0$, the
flatbands correspond to hexagonal closed orbits in real space, while for
$P=1$, the network is localized in triangular orbits, as illustrated in Fig.
7. As such, the two flatband limits correspond to two distinct classical
limits. Further note that there are two $T$ phases, depending on which of the
two triangles of the kagome network the density is mostly localized, see App.
B. The $H$ and $T$ phase are connected by a metallic phase at two points in
the phase diagram, where two gap-closing lines cross, see Fig. 7. In the
context of scattering networks, this is called a percolation point. Hence, an
interface between the trivial $H$ and $T$ phases supports a chiral mode which
propagates in a snakelike fashion. This result seems to indicate that the
Chern trivial phases do posses a notion of relative topology, similar to a
massive Dirac electron subjected to a mass inversion on a domain wall.

Figure 8: Network ribbons of width $M=2$ with (a) triangle and (b) hexagon
edges. Figure 9: (a–c) Energy bands of a ribbon with triangle (a), hexagon
(b), and mixed (c) edges of width $M=12$ for a single valley in the $T$ phase
[$(\varphi,P)=(0,1/2)$]. Here, vertical gray lines indicate the projected
$\bar{K}_{\pm}$ points at $k=\mp 2\pi/3L$ and the relevant winding numbers
$(W_{I},W_{II})$ of the gaps are are indicated. (d–f) Same for the
$\mathcal{C}_{+}=+1$ phase [(c) in Fig. 7].

For the remaining two gapped phases, we numerically calculate the Chern number
and we find two Chern phases, see Fig. 7, where pairs of bands $E_{n,\pm}$
carry opposite Chern numbers $\mathcal{C}_{+}+\mathcal{C}_{-}=0$ with
$\mathcal{C}_{+}=\pm 1$, as expected from phase-rotation symmetry. However,
because of the unbounded spectrum of the network model which repeats
periodically in energy, the total Chern number of the ground state is an ill-
defined quantity. Indeed, there are two distinct gaps whose Chern number
oscillates between zero and $\pm 1$. By calculating the spectrum in a ribbon
geometry, we find that the boundary conditions pick out one particular
sequence. Hence, we need another invariant that is sensitive to the choice of
unit cell, namely a winding number Delplace  _et al._ (2017); Rudner  _et al._
(2013). The winding number is defined for each gap and a difference in winding
numbers equals the net Chern number of the intervening bands. One can
calculate the winding number explicitly by mapping the network to a Floquet
model De Beule  _et al._ (2021b). However, this mapping is not unique since it
depends on the choice of cell Delplace  _et al._ (2017). Consequently, certain
boundary configurations of the network have no counterpart in the Floquet
system once a choice has been made.

In addition to the Chern number, the chiral kagome network is thus also
characterized by a pair of winding numbers $(W_{I},W_{II})$, where $W_{I}$
corresponds to our choice of unit cell shown in Fig. 5, and $W_{II}$
corresponds to a cell whose subnodes all belong to a single hexagon.

###  III.5 Network ribbon

We now consider a ribbon periodic along the $x$ direction with period $2L$ and
finite along the $y$ direction with width $W=\sqrt{3}\,ML$ where
$M=1,2,\ldots$. Moreover, we consider two types of boundaries: triangle and
hexagon edges, illustrated in Fig. 8(a) and (b), respectively. These are the
simplest edge types that respect current conservation, i.e., at each
scattering node the total number of incoming modes should equal the total
number of outgoing modes. The energy bands for ribbons with triangle and
hexagon edges are shown in Figs. 9(a,b) and (d,e) for the trivial $T$ and
Chern phase, respectively. In the trivial phases, the existence of edge modes
depends on the type of edge. Indeed, in the flatband limit of the $T$ phase
($P=1$), a hexagon edge is decoupled from the bulk which is localized in
triangular orbits. Hence, the $T$ phase supports edge modes at a hexagon edge
while there are no edge modes at a triangle edge. The opposite situation then
holds for the $H$ phase. In the Chern phases, edge modes exist at any
boundary. However, as we mentioned previously, while differences of Chern
numbers are well-defined, the net Chern number is undetermined because of the
unbounded spectrum. Moreover, the net Chern number depends in principle on the
entire set of occupied bands. This includes bands that are not captured by the
network model, which is only valid in a finite energy window near charge
neutrality. It is thus a priori not clear which of the two gaps in the
fundamental domain hosts edge modes. This is resolved in the network model by
accounting for the winding number as discussed in the previous section. Here,
we use the bulk-boundary correspondence Rudner  _et al._ (2013) to obtain the
winding numbers from the finite-size calculation, as indicated on the figures.

Finally, we consider a ribbon with a hexagonal edge and a triangle edge on
opposite sides. In this case, each trivial phase hosts edge modes but only on
the edge that does not support a localized bulk mode. At the topological phase
transition, the edge mode moves to the opposite edge without changing
propagation direction. This is consistent with the fact that the hexagonal and
triangle edge naturally support edge modes of opposite chirality (see Fig. 8).
We find single-boundary edge modes in all gapped phases for ribbons with both
types of edges. The mixed-edge ribbon is thus characterized by three
invariants: the Chern number of one of the two bands in the fundamental domain
and a pair of winding numbers, one for each edge. The corresponding spectrum
is shown in Fig. 9(c,f), where we have indicated the boundary for each edge
mode. In a rectangular sample, it is not immediately clear what happens to the
single-boundary edge modes. We find that in a theory with local current
conservation, such a setup is only possible if the system is coupled to a
reservoir.

Figure 10: (a) Network energy bands for valley $K_{+}$ (blue) and $K_{-}$
(red) along high-symmetry lines of the SBZ [green path in Fig. 1(c)] for
$\varphi=0$ and $P(E)$ as shown in (b). The corresponding path in the
$(\varphi,P)$ plane is shown in Fig. 7 as the green dashed line. (c) Wilson
loop spectrum (solid) and the total Berry phase (dashed) for the two-band
manifold near charge neutrality of valley $K_{+}$.

##  IV Discussion

By inspection of the band structure of the valley-projected continuum theory,
shown in Fig. 2, it is clear that the $S$ matrix must gain an energy
dependence to make further progress. This is similar to a quantum Hall
transition, where the system goes through a percolation point at a critical
filling Chalker and Coddington (1988); Kramer  _et al._ (2005); Potter  _et
al._ (2020). In this case, the $S$ matrix traces out a path in $(\varphi,P)$
space as a function of energy and Eq. (20) becomes a nonlinear eigenvalue
problem which is solved self-consistently. Ideally, one obtains the energy-
dependence of the $S$ matrix from a microscopic theory. However, this is a
nontrivial task and instead we demonstrate that we can qualitatively reproduce
both the energy bands and their topology by taking an appropriate energy
dependence. To this end, we require that the resulting bands are approximately
symmetric and gapless at zero energy. Moreover, the net valley Chern number of
the two lowest bands should equal $\nu=\pm 1$. All these conditions are met in
10(a) by taking scattering parameters $\varphi=0$ and $P=P(E)$ which is shown
in Fig. 10(b).

At small energies relative to $E_{L}$, we have $P(E)=1/3+cE$. Hence, the
network goes through a percolation transition at zero energy, where $c$
determines the bandwidth of the low-energy bands. Moreover, we find that the
valley-projected two-band manifold at near zero energy carries a net unit
valley Chern number, while all other bands are Chern trivial. Here, we
numerically calculated the Chern number by opening a tiny gap at charge
neutrality with a small nonzero value for $\varphi$. We find that the sign of
$\varphi$ determines which of the two bands carries the Chern number, similar
to a sublattice-staggering term in the continuum theory Phong and Mele (2022).
We also calculated the winding of the Wilson loop spectrum of the low-energy
bands along a non-contractible loop in the SBZ, which remains well-defined for
$\varphi=0$. The Wilson loop spectrum is shown in Fig. 10(c). Intuitively, the
band topology is understood as follows. The bands at positive (negative)
energies correspond locally to the $T$ ($H$) trivial phase, as indicated on
Fig. 10(a). The unit Chern number of the two valley-projected bands near
charge neutrality corresponds to the difference in the nominal winding number
between the $H$ and $T$ gaps.

While the network model can qualitatively capture both the energy bands and
the topology of the bulk low-energy physics of the periodic PMF, it does not
reproduce the boundary physics of Ref. Phong and Mele, 2022. Indeed, we find
that the band topology in the network model implies that the trivial gaps
above and below the two-band manifold at charge neutrality carry different
winding numbers. Hence, only one of these gaps supports edge modes for a given
boundary, contrary to Ref. Phong and Mele, 2022. This is a consequence of the
coarse-grained approach that breaks down at the boundary. Indeed, in the
network model, the links correspond to structures on the scale of the strain
field, while all microscopics is smoothed over. Furthermore, although both the
network and atomistic models Phong and Mele (2022) support single-boundary
edge modes, there is no clear connection between these two cases. Indeed,
there is no well-defined bulk-boundary correspondence when both valleys are
taken into account because the total Chern number inevitably vanishes when
$\mathcal{T}$ symmetry is restored. As such, the fate of the single-valley
edge modes depends on intervalley scattering at the boundary.

In conclusion, we have demonstrated that the single-valley long-wavelength
physics of periodically strained graphene can be qualitatively understood with
an oriented network model. Here, the symmetries and connectivity of the
network are solely determined by the nodal lines of the corresponding
pseudomagnetic field which host chiral domain-wall modes. In this work, we
considered a commensurate strain field with $C_{3v}$ symmetry which is
relevant for graphene on NbSe2. In this case, the network corresponds to a
chiral kagome network constrained by the valley-preserving symmetries of the
pseudomagnetic field. Because the existence of the network only depends on the
nodal structure of the pseudomagnetic field, similar network physics holds for
any periodic pseudogmagnetic field in graphene with zero net flux.

###### Acknowledgements.

This research was funded in whole, or in part, by the Luxembourg National
Research Fund (FNR), grant reference 16515716. For the purpose of open access,
the author has applied a Creative Commons Attribution 4.0 International (CC BY
4.0) license to any Author Accepted Manuscript version arising from this
submission. Work by VTP and EJM is supported by the Department of Energy under
grant DE-FG02-84ER45118.

##  Appendix A Valley-projected theory

###  A.1 Continuum model

Here, we present more details on the continuum Hamiltonian given in Eq. (3).
For simplicity, we take the gauge $\bm{A}_{\nu}=\nu A(\bm{r})\bm{e}_{x}$ where

$A(\bm{r})=-B_{0}\sum_{i=1}^{3}\frac{\sin(\bm{G}_{i}\cdot\bm{r}+\phi)}{G_{iy}},$
(26)

with

$\bm{G}_{1}=\frac{4\pi}{\sqrt{3}L}\begin{pmatrix}0\\\
1\end{pmatrix},\qquad\bm{G}_{2,3}=\frac{4\pi}{\sqrt{3}L}\begin{pmatrix}\mp\sqrt{3}/2\\\
-1/2\end{pmatrix},$ (27)

such that $\bm{G}_{i}\cdot\bm{L}_{j}=2\pi\delta_{ij}$ ($i,j=1,2$). Up to a
gauge transformation, this gauge is equivalent to the gauge that is directly
related to the strain field Vozmediano  _et al._ (2010),

$\bm{A}_{\nu}(\bm{r})=\frac{\nu\hbar}{e}\frac{\sqrt{3}\beta}{2a}\begin{pmatrix}\varepsilon_{xx}-\varepsilon_{yy}\\\
-2\varepsilon_{xy}\end{pmatrix},$ (28)

where $\nu=\pm 1$ is the valley index, $\beta\sim 1$ a dimensionless material
constant, $a\approx 0.25$ nm the graphene lattice constant, and
$\varepsilon_{ij}$ the strain tensor. The corresponding pseudomagnetic field
(PMF) is given in Eq. (1) which has vanishing net flux,

$\frac{1}{V_{c}}\int_{cell}d^{2}\bm{r}\,\bm{B}_{\nu}\cdot\bm{e}_{z}=\nu
B_{\bm{G}=\bm{0}}=0,$ (29)

with $V_{c}=\sqrt{3}L^{2}/2$ the area of the unit cell.

Remember that $L=Na$ with integer $N\gg 1$ such that the PMF varies slowly on
the graphene lattice scale. In this case, the single-valley theory is a good
approximation. We diagonalize the continuum Hamiltonian in Eq. (3) by Fourier
transform:

$\psi_{\nu}(\bm{r})=\sum_{\bm{k}\in\textrm{SBZ}}\sum_{\bm{G}}e^{i\left(\bm{k}+\bm{G}\right)\cdot\bm{r}}c_{\nu\bm{k}+\bm{G}},$
(30)

where the sum over momentum has been written as a sum over the superlattice
Brillouin zone (SBZ) and a sum over reciprocal lattice vectors of the
superlattice. The valley-projected Hamiltonian given in Eq. (3) becomes

$\displaystyle H_{\nu}$ $\displaystyle=\hbar
v_{F}\sum_{\bm{k}\in\textrm{SBZ}}\sum_{\bm{G},\bm{G}^{\prime}}c_{\nu\bm{k}+\bm{G}}^{\dagger}\Big{[}$
$\displaystyle+\delta_{\bm{G},\bm{G}^{\prime}}\left(\bm{k}+\bm{G}-\bar{\bm{K}}_{\nu}\right)\cdot\left(\nu\sigma_{x},\sigma_{y}\right)$
(31)
$\displaystyle+\sum_{i=1}^{3}\frac{e^{-i\phi}\delta_{\bm{G}+\bm{G}_{i},\bm{G}^{\prime}}-e^{i\phi}\delta_{\bm{G}-\bm{G}_{i},\bm{G}^{\prime}}}{2il_{0}^{2}G_{iy}}\,\sigma_{x}\Big{]}c_{\nu\bm{k}+\bm{G}^{\prime}},$

where we have placed the momentum origin at
$\bar{\bm{K}}_{\nu}=\nu\left(4\pi/3L\right)\bm{e}_{x}$. Next, we truncate the
Bloch Hamiltonian by only taking $|\bm{G}|<k_{c}$ where $k_{c}$ is increased
until the energy bands under consideration converge. The bands are symmetric
about zero energy and we label them as
$E_{n,\nu}(\bm{k})=\operatorname{sgn}(n)E_{|n|}(\nu\bm{k})$ with $n$ a nonzero
integer. The corresponding Bloch wave function is written as
$\psi_{n\bm{k}}(\bm{r})=e^{i\bm{k}\cdot\bm{r}}u_{n\bm{k}}(\bm{r})/\sqrt{N}$
with

$u_{n\bm{k}}(\bm{r})=\frac{1}{\sqrt{V_{c}}}\sum_{\bm{G}}e^{i\bm{G}\cdot\bm{r}}\phi_{n\bm{k},\bm{G}},$
(32)

where $N$ is the number of cells and $\phi_{n\bm{k}}$ are eigenstates of the
Bloch Hamiltonian. The total density of sublattice $\sigma=A,B$ for band $n$
is

$\rho_{n\sigma}(\bm{r})=\frac{1}{N}\sum_{\bm{k}\in\textrm{SBZ}}u_{n\bm{k}}^{\dagger}(\bm{r})\left(\frac{\sigma_{0}\pm\sigma_{z}}{2}\right)u_{n\bm{k}}(\bm{r}).$
(33)

Throughout this paper, we work in dimensionless units by sending

$\displaystyle\bm{k}$ $\displaystyle\rightarrow k_{0}\bm{k},$ (34)
$\displaystyle E$ $\displaystyle\rightarrow\hbar v_{F}k_{0}E,$ (35)

with $k_{0}=4\pi/3L$. In these units, the Hamiltonian only contains one
dimensionless parameter

$\left(k_{0}l_{0}\right)^{-1}=\frac{3}{4\pi}\frac{L}{l_{0}},$ (36)

where (for $t\approx 3$ eV)

$\displaystyle\hbar v_{F}k_{0}=\frac{2\pi t}{\sqrt{3}N}$
$\displaystyle\approx\frac{11\,\textrm{eV}}{N},$ (37)
$\displaystyle\frac{L}{l_{0}}$
$\displaystyle\approx\frac{N}{100}\sqrt{\frac{B_{0}}{T}}.$ (38)

For example, for $N=60$ ($L\approx 15$ nm) and $B_{0}=100$ T, which are
realistic values for engineering a periodic strain field with a varying height
profile Phong and Mele (2022) of the order of $6$ Å, we obtain $L/l_{0}\approx
6$ and $\hbar v_{F}k_{0}\approx 180$ meV.

###  A.2 Symmetries

Here, we show how the (pseudo)vector potential is constrained by the
symmetries of the system. Note that the $\mathcal{C}_{3z}$ symmetry inherited
from the strain field is not the microscopic $\mathcal{C}_{3z}$ symmetry,
which has to be broken since otherwise the strain tensor is only given by its
trace (pure dilatation) and the pseudovector potential would be absent.
Indeed, one can think about the pseudovector potential in terms of a strain-
induced local shift of the Dirac point from the zone corner, to which it would
be fixed when microscopic $\mathcal{C}_{3z}$ were conserved. The remaining
symmetries that we consider are spinless time reversal ($\mathcal{T}$),
twofold rotation ($\mathcal{C}_{2z}$), and mirror symmetry across the $y$ axis
$(\mathcal{M}_{x})$. We take the following representation for the action of
the symmetries on the field operators:

$\displaystyle\mathcal{T}\psi_{\nu}(\bm{r})\mathcal{T}^{-1}$
$\displaystyle=\psi_{-\nu}(\bm{r})$ (39)
$\displaystyle\mathcal{C}_{2z}\psi_{\nu}(\bm{r})\mathcal{C}_{2z}^{-1}$
$\displaystyle=\sigma_{x}\psi_{-\nu}(-\bm{r})$ (40)
$\displaystyle\mathcal{M}_{x}\psi_{\nu}(x,y)\mathcal{M}_{x}^{-1}$
$\displaystyle=\psi_{-\nu}(-x,y)$ (41)
$\displaystyle\mathcal{C}_{3z}\psi_{\nu}(\bm{r})\mathcal{C}_{3z}^{-1}$
$\displaystyle=e^{i\nu(2\pi/3)\sigma_{z}}\psi_{\nu}(C_{3z}\bm{r}).$ (42)

Next, we investigate how these symmetries act on the (pesudo)vector potential.
Under time reversal, we find from $\mathcal{T}H\mathcal{T}^{-1}=H$ with
$H=H_{+}+H_{-}$,

$\mathcal{T}:\quad\bm{A}(\bm{r})\rightarrow-\bm{A}(\bm{r}),\quad\bm{A}_{\nu}(\bm{r})\rightarrow\bm{A}_{\nu}(\bm{r}),$
(43)

where we used $\mathcal{T}i\mathcal{T}^{-1}=-i$. We see that a PMF conserves
time-reversal symmetry as expected. On the other hand, under a twofold
rotation, we have

$\mathcal{C}_{2z}:\quad\bm{A}(\bm{r})\rightarrow-\bm{A}(-\bm{r}),\quad\bm{A}_{\nu}(\bm{r})\rightarrow\bm{A}_{\nu}(-\bm{r}),$
(44)

such that $\mathcal{C}_{2z}$ is conserved when
$\bm{A}(\bm{r})=-\bm{A}(-\bm{r})+\nabla f$ for a real magnetic field and
$\bm{A}_{\nu}(\bm{r})=\bm{A}_{\nu}(-\bm{r})+\nabla f$ for a PMF, where
$f(\bm{r})$ is a scalar function. Hence, the real (pseudo) magnetic field has
to be an even (odd) function of the position vector. The PMF that we consider
in this work [see Eq. (1)] generally does not conserve $\mathcal{C}_{2z}$
which is in fact broken by the substrate. Only the special case $\phi=\pi/2$
conserves $\mathcal{C}_{2z}$, yielding a gapless spectrum that can be mapped
to a triangular chiral network. Finally, the mirror $\mathcal{M}_{x}$ gives,

$\displaystyle\bm{A}(x,y)$
$\displaystyle\rightarrow\begin{bmatrix}-A_{x}(-x,y)\\\
+A_{y}(-x,y)\end{bmatrix},$ (45) $\displaystyle\bm{A}_{\nu}(x,y)$
$\displaystyle\rightarrow\begin{bmatrix}+A_{\nu x}(-x,y)\\\ -A_{\nu
y}(-x,y)\end{bmatrix},$ (46)

such that for a PMF, $\mathcal{M}_{x}$ is conserved if

$\displaystyle A_{\nu x}(x,y)$ $\displaystyle=+A_{\nu x}(-x,y)+\partial_{x}f,$
(47) $\displaystyle A_{\nu y}(x,y)$ $\displaystyle=-A_{\nu
y}(-x,y)+\partial_{y}f,$ (48)

and $\bm{B}_{\nu}(x,y)=\bm{B}_{\nu}(-x,y)$. Note that in general, the mirror
line is offset to $x=-2\phi/3$.

Finally, we note that the Hamiltonian in Eq. (3) also has a chiral symmetry,

$\mathcal{C}\psi_{\nu}(\bm{r})\mathcal{C}^{-1}=\psi_{\nu}^{\dagger}(\bm{r})\sigma_{z},\qquad\mathcal{C}i\mathcal{C}^{-1}=-i,$
(49)

such that $\mathcal{C}H_{\nu}\mathcal{C}^{-1}=H_{\nu}$. However, this is a
model-dependent symmetry due to the absence of terms proportional to
$\sigma_{0}$ or $\sigma_{z}$ in the Hamiltonian. Such terms are in principle
allowed by symmetry, e.g., a strain-induced pseudo-electrostatic potential or
a constant sublattice-staggering term, respectively. The latter is
proportional to the trace of the strain tensor Suzuura and Ando (2002); Low
_et al._ (2011). In the tight-binding framework, the pseudo-electrostatic
potential originates from a modulation of intra-sublattice hopping terms.

##  Appendix B Network model

###  B.1 Symmetry constraints

We demonstrate how the $S$ matrix is constrained by the valley-preserving
symmetries $\mathcal{C}_{3z}$ and $\mathcal{M}_{x}\mathcal{T}$.

Under the threefold rotation, the subnodes of the kagome network undergo a
cyclic permutation $1\rightarrow 2\rightarrow 3\rightarrow 1$ and therefore
$S_{0}=S_{1}\sim S_{2}\sim S_{3}$ where $\sim$ indicates that the equality
holds up to a diagonal unitary transformation. However, the resulting phases
can always be removed by choosing a suitable basis for the amplitudes, i.e.,
$\left(a_{1},a_{2},a_{3}\right)\rightarrow\left(a_{1},U_{2}a_{2},U_{3}a_{3}\right)$
and similar for outgoing amplitudes, where $U_{2}$ and $U_{3}$ are diagonal
$2\times 2$ unitary matrices.

The composite symmetry $\mathcal{M}_{x}\mathcal{T}$ is more tricky. First, we
consider the action of $\mathcal{M}_{x}$ and $\mathcal{T}$ separately. To this
end, we need to consider both valleys. For subnode $1$ (red dot in Fig. 5) we
have

$\begin{pmatrix}b_{1}\\\ b_{1}^{\prime}\end{pmatrix}=\begin{pmatrix}S_{1}&0\\\
0&S_{1}^{\prime}\end{pmatrix}\begin{pmatrix}a_{1}\\\
a_{1}^{\prime}\end{pmatrix},$ (50)

where the unprimed and primed quantities correspond to valleys $K_{+}$ and
$K_{-}$, respectively. Here, we have assumed that there is no intervalley
scattering. Under $\mathcal{M}_{x}$, the amplitudes transform as

$\displaystyle\mathcal{M}_{x}a_{1}\mathcal{M}_{x}^{-1}=a_{1}^{\prime},$ (51)
$\displaystyle\mathcal{M}_{x}b_{1}\mathcal{M}_{x}^{-1}=b_{1}^{\prime},$ (52)

which is clear from Fig. 5. We then obtain

$\begin{pmatrix}b_{1}^{\prime}\\\
b_{1}\end{pmatrix}=\mathcal{M}_{x}\begin{pmatrix}S_{1}&0\\\
0&S_{1}^{\prime}\end{pmatrix}\mathcal{M}_{x}^{-1}\begin{pmatrix}a_{1}^{\prime}\\\
a_{1}\end{pmatrix},$ (53)

or

$\mathcal{M}_{x}\begin{pmatrix}S_{1}&0\\\
0&S_{1}^{\prime}\end{pmatrix}\mathcal{M}_{x}^{-1}=\begin{pmatrix}S_{1}^{\prime}&0\\\
0&S_{1}\end{pmatrix}.$ (54)

When $\mathcal{M}_{x}$ is a symmetry, we have $S_{1}\sim S_{1}^{\prime}$. On
the other hand, time-reversal symmetry yields

$\displaystyle\mathcal{T}a_{1}\mathcal{T}^{-1}=\left(b_{1}^{\prime}\right)^{*},$
(55)
$\displaystyle\mathcal{T}b_{1}\mathcal{T}^{-1}=\left(a_{1}^{\prime}\right)^{*},$
(56)

or

$\mathcal{T}\begin{pmatrix}S_{1}&0\\\
0&S_{1}^{\prime}\end{pmatrix}\mathcal{T}^{-1}=\begin{pmatrix}S_{1}^{\prime}&0\\\
0&S_{1}\end{pmatrix}^{t},$ (57)

such that $S_{1}\sim\left(S_{1}^{\prime}\right)^{t}$ when $\mathcal{T}$ is
conserved. The combined symmetry $\mathcal{M}_{x}\mathcal{T}$ yields
$S_{1}\sim\left(S_{1}\right)^{t}$. Since the diagonal unitary only acts within
subnode $1$, it can removed independently of the phases that were removed
under the constraint given by $\mathcal{C}_{3z}$.

In summary, threefold rotation symmetry ($\mathcal{C}_{3z}$) and mirror
symmetry combined with time-reversal symmetry ($\mathcal{M}_{x}\mathcal{T}$)
constrains the $S$ matrix as follows:

$\mathcal{S}=\mathds{1}_{3}\otimes S_{0},\qquad S_{0}=\left(S_{0}\right)^{t},$
(58)

where an explicit form of $S_{0}$ is given in Eq. (10).

###  B.2 Classical loop configurations

Since there are two nonequivalent triangles in the kagome lattice, as
$\mathcal{C}_{2z}$ is broken, there are tree possible classical loop
configurations, depending on whether the loops are going along one of two
nonequivalent triangles or the hexagon. On average, for a single band, the
loops run only along one of the two triangles for $P\approx 1$. To this end,
we calculated the following “loop order parameters” for the network bands
$E_{n,s}(\bm{k})$:

$\displaystyle t_{1,s}$
$\displaystyle=\frac{3^{3}}{N}\sum_{\bm{k}}\left|b_{11}b_{21}b_{31}\right|^{2}_{s},$
(59) $\displaystyle t_{2,s}$
$\displaystyle=\frac{3^{3}}{N}\sum_{\bm{k}}\left|b_{12}b_{22}b_{32}\right|^{2}_{s},$
(60) $\displaystyle h_{s}$
$\displaystyle=\frac{6^{6}}{N}\sum_{\bm{k}}\left|b_{11}b_{12}b_{21}b_{22}b_{31}b_{32}\right|^{2}_{s},$
(61)

which are shown in Fig. 11. Here, $t_{1s}$ ($t_{2s}$) is close to unity if the
band corresponds to a classical loop configuration along the downward (upward)
pointing triangles in the kagome network and $h_{s}$ is close to unity if it
corresponds to loops around hexagons (see Fig. 7). We find that $h_{+}=h_{-}$
and $t_{1,s}=t_{2,-s}$. Note also that $t_{1,s}$ and $t_{2,s}$ do not vanish
entirely in the $H$ phase, because in that case, the amplitude is equal for
all links, and they attain the average value $3^{3}/6^{3}=1/8$.

Figure 11: Loop order parameters for the band $E_{n,+}$ shown in the
$(\varphi,P)$ scattering parameter space of the chiral kagome network, where
the white curves give the gap closing lines.

###  B.3 Network ribbons

The calculation for the network ribbons goes as follows: one first constructs
a large $S$ matrix that describes scattering within the unit cell of the
ribbon. Next, one constructs the matrix $\mathcal{M}(k)$ which contains the
Bloch phases. We then obtain

$\mathcal{S}\mathcal{M}(k)b_{k}=e^{i2\pi E/E_{L}}b_{k},$ (62)

where now $\mathcal{S}$ contains many copies of $S_{0}$ together with boundary
$S$ matrices for the top and bottom edge. The connectivity matrix
$\mathcal{M}(k)$ for the network ribbon is obtained in a similar way as for
the bulk system, see Eqs. (11)–(16).

###  B.4 Energy dependence of the $S$ matrix

When the $S$ matrix depends on the energy, Eq. (20) becomes a nonlinear
eigenvalue problem, which we solve self-consistently.

The energy dependence of the scattering parameters that we have used to obtain
the network band structure shown in Fig. 10(a) is explicitly given by

$P(E)=\left[\frac{1}{1+\left(3^{1/a}-1\right)e^{-bE}}\right]^{a},\qquad\varphi=0,$
(63)

with $a=1.9$ and $b=1.1$. These values are chosen to make the bands
approximately particle-hole symmetric.

###  B.5 Wilson loop

The Wilson loop is a generalization of the Berry phase for a number of
isolated bands. Since these bands are isolated, there is a finite energy gap
between these bands and other bands, while the gap between the constituent
bands can vanish. For example, the valley-projected two-band manifold at
charge neutrality in the main text consists of two isolated bands.

To define the Wilson loop, we first consider a closed non-contractible path
$\mathcal{C}$ in the SBZ. The Wilson loop along $\mathcal{C}$ is then defined
as

$\mathcal{W}(\mathcal{C})=\mathcal{P}\exp\left[-i\oint_{\mathcal{C}}d\bm{k}\cdot\mathcal{A}\right],$
(64)

which is a $M\times M$ unitary matrix with $M$ the number of isolated bands.
Here, $\mathcal{P}$ denotes path ordering and $\mathcal{A}$ is the Berry
connection matrix of the isolated manifold of bands. If we want to calculate
the Wilson loop numerically, we need to discretize the loop $\mathcal{C}$. For
a loop at $\bm{k}_{0}$, consider the discretization
$\\{\bm{k}_{1},\bm{k}_{2},\ldots,\bm{k}_{n}\\}$ where
$\bm{k}_{n}=\bm{k}_{0}+\bm{G}$ with $\bm{G}$ a reciprocal lattice vector. The
discretized Wilson loop along $\mathcal{C}$ is then defined as

$\mathcal{W}_{n}(\mathcal{C})=\Psi_{\bm{k}_{1}}^{\dagger}\Psi_{\bm{k}_{2}}\Psi_{\bm{k}_{2}}^{\dagger}\Psi_{\bm{k}_{3}}\cdots\Psi_{\bm{k}_{n-1}}^{\dagger}\Psi_{\bm{k}_{n}},$
(65)

where $\Psi_{\bm{k}}$ is a matrix of eigenstates of the $M$ isolated bands at
$\bm{k}$. In our case, we have

$\Psi_{\bm{k}}=\left(b_{1\bm{k}}\,b_{2\bm{k}}\,\cdots\,b_{M\bm{k}}\right).$
(66)

Note that since we work in Bloch form, we need to enforce
$\Psi_{\bm{k}_{n}}=\Psi_{\bm{k}_{1}}$. In the limit $n\rightarrow\infty$, Eq.
(65) becomes exact and
$\mathcal{W}(\mathcal{C})=\lim_{n\rightarrow\infty}\mathcal{W}_{n}(\mathcal{C})$.

We then investigate the band topology with the winding of the phase of the
Wilson loop eigenvalues (the so-called Wilson spectrum) as the loop
$\mathcal{C}$ is varied. For concreteness, if
$\bm{k}=k_{1}\bm{G}_{1}+k_{2}\bm{G}_{2}$ is a general momentum, we consider
the loops $\mathcal{C}(k_{1})$ for $0\leq k_{2}\leq 1$. We then plot the
Wilson spectrum as a function of $k_{1}$. This is shown in Fig. 10(c) for the
valley-projected two-band manifold near charge neutrality. We also show the
sum of the Wilson loop eigenvalues, which corresponds to the total Berry phase
of the two bands. The nontrivial winding of the total Berry phase indicates
that these bands carry a net Chern number of $\pm 1$. This is a well-defined
notion since the crossing point between the two bands is not protected by
symmetry.

## References

  * Forsythe  _et al._ (2018) C. Forsythe, X. Zhou, K. Watanabe, T. Taniguchi, A. Pasupathy, P. Moon, M. Koshino, P. Kim, and C. R. Dean, “Band structure engineering of 2d materials using patterned dielectric superlattices,” Nat. Nanotechnol. 13, 566–571 (2018).
  * Mao  _et al._ (2020) J. Mao, S. P. Milovanović, M. Anđelković, X. Lai, Y. Cao, K. Watanabe, T. Taniguchi, L. Covaci, F. M. Peeters, A. K. Geim, Y. Jiang, and E. Y. Andrei, “Evidence of flat bands and correlated states in buckled graphene superlattices,” Nature 584, 215–220 (2020).
  * Jiang  _et al._ (2017) Y. Jiang, J. Mao, J. Duan, X. Lai, K. Watanabe, T. Taniguchi, and E. Y. Andrei, “Visualizing strain-induced pseudomagnetic fields in graphene through an hbn magnifying glass,” Nano Lett. 17, 2839–2843 (2017).
  * Andrei  _et al._ (2021) E. Y. Andrei, D. K. Efetov, P. Jarillo-Herrero, A. H. MacDonald, K. F. Mak, T. Senthil, E. Tutuc, A. Yazdani, and A. F. Young, “The marvels of moiré materials,” Nat. Rev. Mater. 6, 201–206 (2021).
  * Lopes dos Santos  _et al._ (2007) J. M. B. Lopes dos Santos, N. M. R. Peres, and A. H. Castro Neto, “Graphene Bilayer with a Twist: Electronic Structure,” Phys. Rev. Lett. 99, 256802 (2007).
  * Cao  _et al._ (2018a) Y. Cao, V. Fatemi, A. Demir, S. Fang, S. L. Tomarken, J. Y. Luo, J. D. Sanchez-Yamagishi, K. Watanabe, T. Taniguchi, E. Kaxiras, R. C. Ashoori, and P. Jarillo-Herrero, “Correlated insulator behaviour at half-filling in magic-angle graphene superlattices,” Nature 556, 80–84 (2018a).
  * Cao  _et al._ (2018b) Y. Cao, V. Fatemi, S. Fang, K. Watanabe, T. Taniguchi, E. Kaxiras, and P. Jarillo-Herrero, “Unconventional superconductivity in magic-angle graphene superlattices,” Nature 556, 43–50 (2018b).
  * Milovanović  _et al._ (2020) S. P. Milovanović, M. Anđelković, L. Covaci, and F. M. Peeters, “Band flattening in buckled monolayer graphene,” Phys. Rev. B 102, 245427 (2020).
  * Manesco  _et al._ (2021) A. L. R. Manesco, J. L. Lado, E. V. S. Ribeiro, G. Weber, and D. Rodrigues Jr, “Correlations in the elastic landau level of spontaneously buckled graphene,” 2D Mater. 8, 015011 (2021).
  * Manesco and Lado (2021) A. L. R. Manesco and J. L. Lado, “Correlation-induced valley topology in buckled graphene superlattices,” 2D Mater. 8, 035057 (2021).
  * Giambastiani  _et al._ (2022) D. Giambastiani, F. Colangelo, A. Tredicucci, S. Roddaro, and A. Pitanti, “Electron localization in periodically strained graphene,” J. Appl. Phys. 131, 085103 (2022).
  * Phong and Mele (2022) V. T. Phong and E. J. Mele, “Boundary modes from periodic magnetic and pseudomagnetic fields in graphene,” Phys. Rev. Lett. 128, 176406 (2022).
  * Suzuura and Ando (2002) H. Suzuura and T. Ando, “Phonons and electron-phonon scattering in carbon nanotubes,” Phys. Rev. B 65, 235412 (2002).
  * Mañes (2007) J. L. Mañes, “Symmetry-based approach to electron-phonon interactions in graphene,” Phys. Rev. B 76, 045430 (2007).
  * Vozmediano  _et al._ (2010) M. Vozmediano, M. Katsnelson, and F. Guinea, “Gauge fields in graphene,” Phys. Rep. 496, 109–148 (2010).
  * Guinea  _et al._ (2010) F. Guinea, A. K. Geim, M. I. Katsnelson, and K. S. Novoselov, “Generating quantizing pseudomagnetic fields by bending graphene ribbons,” Phys. Rev. B 81, 035408 (2010).
  * Levy  _et al._ (2010) N. Levy, S. A. Burke, K. L. Meaker, M. Panlasigui, A. Zettl, F. Guinea, A. H. C. Neto, and M. F. Crommie, “Strain-Induced Pseudo–Magnetic Fields Greater Than 300 Tesla in Graphene Nanobubbles,” Science 329, 544–547 (2010).
  * Novoselov  _et al._ (2005) K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang, M. I. Katsnelson, I. V. Grigorieva, S. V. Dubonos, and A. A. Firsov, “Two-dimensional gas of massless dirac fermions in graphene,” Nature 438, 197–200 (2005).
  * San-Jose and Prada (2013) P. San-Jose and E. Prada, “Helical networks in twisted bilayer graphene under interlayer bias,” Phys. Rev. B 88, 121408(R) (2013).
  * Efimkin and MacDonald (2018) D. K. Efimkin and A. H. MacDonald, “Helical network model for twisted bilayer graphene,” Phys. Rev. B 98, 035404 (2018).
  * Fleischmann  _et al._ (2020) M. Fleischmann, R. Gupta, F. Wullschläger, S. Theil, D. Weckbecker, V. Meded, S. Sharma, B. Meyer, and S. Shallcross, “Perfect and Controllable Nesting in Minimally Twisted Bilayer Graphene,” Nano Lett. 20, 971–978 (2020).
  * Tsim  _et al._ (2020) B. Tsim, N. N. T. Nam, and M. Koshino, “Perfect one-dimensional chiral states in biased twisted bilayer graphene,” Phys. Rev. B 101, 125409 (2020).
  * Chou  _et al._ (2020) Y.-Z. Chou, F. Wu, and S. Das Sarma, “Hofstadter butterfly and Floquet topological insulators in minimally twisted bilayer graphene,” Phys. Rev. Res. 2, 033271 (2020).
  * De Beule  _et al._ (2020) C. De Beule, F. Dominguez, and P. Recher, “Aharonov-Bohm Oscillations in Minimally Twisted Bilayer Graphene,” Phys. Rev. Lett. 125, 096402 (2020).
  * De Beule  _et al._ (2021a) C. De Beule, F. Dominguez, and P. Recher, “Network model and four-terminal transport in minimally twisted bilayer graphene,” Phys. Rev. B 104, 195410 (2021a).
  * Vakhtel  _et al._ (2022) T. Vakhtel, D. O. Oriekhov, and C. W. J. Beenakker, “Bloch oscillations in the magnetoconductance of twisted bilayer graphene,” Phys. Rev. B 105, L241408 (2022).
  * Pal  _et al._ (2019) H. K. Pal, S. Spitz, and M. Kindermann, “Emergent Geometric Frustration and Flat Band in Moiré Bilayer Graphene,” Phys. Rev. Lett. 123, 186402 (2019).
  * Moulsdale  _et al._ (2022) C. Moulsdale, A. Knothe, and V. Fal’ko, “Kagome network of miniband-edge states in double-aligned graphene–hexagonal boron nitride structures,” Phys. Rev. B 105, L201112 (2022).
  * Gusynin and Sharapov (2005) V. P. Gusynin and S. G. Sharapov, “Unconventional integer quantum hall effect in graphene,” Phys. Rev. Lett. 95, 146801 (2005).
  * Peres  _et al._ (2006) N. M. R. Peres, F. Guinea, and A. H. Castro Neto, “Electronic properties of disordered two-dimensional carbon,” Phys. Rev. B 73, 125411 (2006).
  * Zhang  _et al._ (2005) Y. Zhang, Y.-W. Tan, H. L. Stormer, and P. Kim, “Experimental observation of the quantum Hall effect and Berry’s phase in graphene,” Nature 438, 201–204 (2005).
  * Haldane (1988) F. D. M. Haldane, “Model for a quantum hall effect without landau levels: Condensed-matter realization of the ”parity anomaly”,” Phys. Rev. Lett. 61, 2015–2018 (1988).
  * Schakel (1991) A. M. J. Schakel, “Relativistic quantum Hall effect,” Phys. Rev. D 43, 1428–1431 (1991).
  * Fialkovsky and Vassilevich (2012) I. Fialkovsky and D. Vassilevich, “Quantum field theory in graphene,” Int. J. Mod. Phys. A 27, 1260007 (2012).
  * Lapa (2019) M. F. Lapa, “Parity anomaly from the Hamiltonian point of view,” Phys. Rev. B 99, 235144 (2019).
  * Castro Neto  _et al._ (2009) A. H. Castro Neto, F. Guinea, N. M. R. Peres, K. S. Novoselov, and A. K. Geim, “The electronic properties of graphene,” Rev. Mod. Phys. 81, 109–162 (2009).
  * Delplace (2020) P. Delplace, “Topological chiral modes in random scattering networks,” SciPost Phys. 8, 081 (2020).
  * Delplace  _et al._ (2017) P. Delplace, M. Fruchart, and C. Tauber, “Phase rotation symmetry and the topology of oriented scattering networks,” Phys. Rev. B 95, 205413 (2017).
  * Pasek and Chong (2014) M. Pasek and Y. D. Chong, “Network models of photonic Floquet topological insulators,” Phys. Rev. B 89, 075113 (2014).
  * Rudner  _et al._ (2013) M. S. Rudner, N. H. Lindner, E. Berg, and M. Levin, “Anomalous Edge States and the Bulk-Edge Correspondence for Periodically Driven Two-Dimensional Systems,” Phys. Rev. X 3, 031005 (2013).
  * De Beule  _et al._ (2021b) C. De Beule, F. Dominguez, and P. Recher, “Effective Floquet model for minimally twisted bilayer graphene,” Phys. Rev. B 103, 195432 (2021b).
  * Chalker and Coddington (1988) J. T. Chalker and P. D. Coddington, “Percolation, quantum tunnelling and the integer Hall effect,” J. Phys. C Solid State Phys. 21, 2665–2679 (1988).
  * Kramer  _et al._ (2005) B. Kramer, T. Ohtsuki, and S. Kettemann, “Random network models and quantum phase transitions in two dimensions,” Phys. Rep. 417, 211–342 (2005).
  * Potter  _et al._ (2020) A. C. Potter, J. T. Chalker, and V. Gurarie, “Quantum Hall Network Models as Floquet Topological Insulators,” Phys. Rev. Lett. 125, 086601 (2020).
  * Low  _et al._ (2011) T. Low, F. Guinea, and M. I. Katsnelson, “Gaps tunable by electrostatic gates in strained graphene,” Phys. Rev. B 83, 195436 (2011).