# Tunable chiral anomalies and coherent transport on a honeycomb lattice

Vasil A. Saroka,${}^{1,2,3\ast}\dagger$ Fanmiao Kong,4 Charles A. Downing,5
Renebeth B. Payod,6 Felix R. Fischer,7,8,9,10 Xiankai Sun,11
Lapo Bogani,4

1Center for Quantum Spintronics, Department of Physics,
Norwegian University of and Technology
NO-7491, Trondheim, Norway
2 Institute for Nuclear Problems, Belarusian State University,
Bobruiskaya 11, 220030 Minsk, Belarus
3 TBpack Ltd., 27 Old Gloucester Street,
London, WC1N 3AX, United Kingdom
4 Department of Materials, University of Oxford,
16 Parks Road, OX1 3PH Oxford, United Kingdom
5 Department of Physics and Astronomy, University of Exeter,
Exeter EX4 4QL, United Kingdom
6 Institute of Mathematical Sciences and Physics, University of the
Philippines,
Los Baños, Laguna 4031, Philippines
7 Department of Chemistry, University of California,
Berkeley, CA 94720, USA
8 Materials Sciences Division, Lawrence Berkeley National Laboratory,
Berkeley, CA 94720, USA
9 Kavli Energy NanoSciences Institute at the University of California Berkeley
and the Lawrence Berkeley National Laboratory,
Berkeley, California 94720, USA
10 Bakar Institute of Digital Materials for the Planet,
Division of Computing, Data Science, and Society, University of California,
Berkeley, CA 94720, USA.
11 Department of Electronic Engineering, The Chinese University of Hong Kong,
Shatin, New Territories, Hong Kong SAR, China

∗To whom correspondence should be addressed; E-mail<EMAIL_ADDRESS>
$\dagger$Present address: Department of Physics, University of Rome Tor
Vergata and INFN,
Via della Ricerca Scientifica 1, 00133 Roma, Italy

> The search for energy efficient materials is urged not only by the needs of
> modern electronics but also by emerging applications in neuromorphic
> computing and artificial intelligence. Currently, there exist two mechanisms
> for achieving dissipationless transport: superconductivity and the quantum
> Hall effect. Here we reveal that dissipationless transport is theoretically
> achievable on a honeycomb lattice by rational design of chiral anomalies.
> Breaking the usual assumption of commensurability leads to electronic modes
> exhibiting chiral anomalies capable of dissipationless coherent transport in
> the material bulk, rather than on the edge. The dissipationless transport
> undergoes a pitchfork bifurcation marking the onset of a novel cubic
> dispersion phase. These results open a new research avenue for the design of
> energy efficient information processing and higher-order dispersion
> materials.

## Introduction

Dissipationless charge transport is possible through superconductivity (?) and
the quantum Hall effect (?). Both phenomena have been reported for honeycomb
lattice materials. Mott-insulating and superconducting states have recently
been found in a moiré superlattices of graphene bilayers at the magic angle of
$1.1^{\circ}~{}$ (?, ?). The half-integer quantum Hall effect has been
reported for graphene (?). The former boils down to flat band and it carries a
lot of similarities with strongly correlated high temperature superconductors,
while the latter relies on gapless edge mode generation via time-reversal
symmetry breaking with magnetic field. These phenomena experimentally work at
low temperatures. Superconducting Cooper pairs are fragile not only with
respect to temperature ($>30$ K), but also magnetic field ($>0.2-20$ T),
magnetic impurities (?) and non-magnetic disorder (?). The Hall conductance
quantization is impressively stable against disorder, but its observation
still requires low temperatures ($\sim 4$ K) or very high magnetic fields
($\sim 30$ T) (?).

One alternative way to achieve dissipationless charge transport is suggested
by topological quantum matter. It has been realized that Haldane’s model of
the anomalous quantum Hall effect, when combined with spin-orbit interactions,
allows ballistic transport via helical states, that makes the quantum spin
Hall (QSH) effect possible even at zero magnetic field (?), see Fig.1A. This
problem has been also approached from a different route by making use of the
internal valley degree of freedom produced by the honeycomb lattice (?). The
reciprocal space of the honeycomb lattice features two non-equivalent Dirac
points conventionally denoted as K and K′. These regions of the reciprocal
space with linear energy bands centered at K and K′ points are called valleys.
The zigzag edges of the honeycomb lattice host the edge states that are well-
separated in the reciprocal space. However, these edge states do not
transverse the bulk energy gap. This can be seen by applying a staggered
potential that opens the band gap between the edge states. Nevertheless, one
can make the edge states to cross the bulk energy gap by creating a kink in
the staggered sublattice potential (?). By creating two kinks via back gate
potential transversing the ribbon, a valley Hall regime can be achieved (?).
In this regime, for any given edge state and valley, a counter propagating
state is found only at the opposite edge of the ribbon. Hence, in the absence
of the short range scatters causing intervalley mixing, the backscattering is
suppressed. The longitudinal electrical current along the ribbon causes a
transverse valley current – the quantum valley Hall (QVH) effect, see Fig. 1A.
The property of the edge modes which locks their direction of propagation to
the edge, spin or valley is referred to as chirality or helicity. Engineering
gapless chiral edge modes without external magnetic field is an important
milestone for the adaptation of quantum Hall dissipationless transport to low-
powered electronic circuits.

With respect to the transport properties and energy bands projected on the
transport direction, solids can be described as either gapped or gapless,
including the recent topological materials (?). Gapped materials like the
semiconducting Si and GaAs and the insulating diamond and SiO2, have a common
trait where their projected conduction and valence bands exhibit extrema at
the band edges. Around the band edge, an adequate expansion of the energy band
function requires at least second order terms, therefore the energy band
dispersion is quadratic for gapped structures. On the other hand, gapless
materials such as the typical metals and the interfaces of topological
insulators have projected bands that cross the Fermi level so that they
possess zero-energy states (ZES) pinned at the Fermi energy. In the ZES
vicinity, the energy band function has at least a first order expansion
resulting to a linear energy band dispersion. Thus, we can classify solids
into two major classes having either quadratic or linear dispersions of the
projected bands. Non-conventional cases aside from the ones described above
include the strongly correlated Mott insulators like the NiO and MnO which
exhibit dispersionless flat bands (?, ?), the Group V semimetals such as
bismuth with quadratic dispersion (?), and the Dirac (?) and Weyl (?)
semimetals that exhibit a linear dispersion relation in the vicinity of the
Fermi level. Materials with cubic or higher order projected dispersions have
not yet been reported.

In this paper, we search for metallic nanostructures favorable for the
dissipaitonless quantum transport and find an example of a completely flat
band that is an eigenmode of a chiral symmetry operator and that can be
transformed into a valley chiral charge conserving gapless mode exhibiting
cubic dispersion. This transformation to a gapless chiral mode is driven by an
external electric field allowing for continuous switching between strongly
correlated insulator-like and topological insulator-like dissipationless
transport regimes. In strike difference to topological insulators, however,
the reported dissipationless transport is realized by bulk rather than edge
states, see Fig. 1A. The onset of a new material phase with cubic dispersion
relation is manifested by a transitory state of the valley chiral mode that
demonstrates distinctive signatures of pitchfork bifurcation. These phenomena
offer unique possibility to investigate a broad variety of quantum transport
regimes in unconventional topological materials, such as a half-bearded
graphene nanoribbon (hbGNR). We unveil that within the graphene nanoribbons
(GNRs) incommensurate with the bulk graphene the said ribbon is an iconic
representative of a unique class exhibiting dissipationless transport and
cubic dispersion phase. Next, we introduce the notion of incomeensurability
from the point of graph theory and then we design the features outlined above.

## Main results

##### Incommensurability and strong pairing.

The hbGNR is incommensurate with the bulk graphene. This means the hbGNR unit
cell cannot be “tiled” with an integer number of graphene unit cells. This
statement can be formulated in terms of graph theory starting from a
celebrated Su-Schrieffer-Heeger (SSH) model. A periodic SSH chain presented in
Fig. 1B consists of a bipartite unit cell containing two atoms that can be
denoted as A (white) and B (gray). The chain is characterized by two
interatomic distances $d_{1}$ and $d_{2}$ and the corresponding hopping
integrals that are functions of $d$’s. The ratio between $d$’s determines the
topological phase of this system. The trivial and non-trivial topological
phases are characterized by 1D ${\bf Z}_{2}$ invariant that takes two distinct
values $0$ (trivial) and $1$ (non-trivial) (?, ?). When
$t_{1}(d_{1})>t_{2}(d_{2})$, the infinite SSH chain is a trivial insulator
with $\nu=0$, but for $t_{1}(d_{1})<t_{2}(d_{2})$ it is a topological
insulator with $\nu=1$. The bulk-boundary correspondence implies that a
boundary ZES forms at the interface between regions with different $\nu$’s. It
is easy to verify that, there are no ZESs, when a strong paring occurs between
intracell atoms ($t_{1}>t_{2}$) and the SSH chain is truncated between the
bulk unit cells. However, as shown in Fig. 1C, a pair of ZESs localized at the
edges of the chain arises when $\nu=1$ and $t_{1}<t_{2}$. However, this
consideration, that is valid for an even number of atoms in the SSH chain,
fails for an odd number of atoms. As seen from Fig. 1D, independently of the
$\nu$, the odd SSH chain features a sole ubiquitous ZES that localizes either
on one edge or another. Here the system with two edges hosts a single ZES.
Both even and odd SSH chains featuring ZESs can be interpreted as
incommensurable structures since a pair of strongly bonded atoms chosen as a
SSH model unit cell cannot fully tile these structures. At the same time, the
even SSH chain, that lacks ZESs, can be grouped into an integer number of
strongly bonded pairs representing a structure that is commensurable with the
periodic SSH unit cell.

##### Perfect matchings and topological phases.

An intuitive notion of commensurable and incommensurable SSH chains allow for
a precise description in terms of graph theory matchings (?). Consider a
finite SSH chain as an undirected regular multigraph $G(V,D)$, with $V$ and
$D$ being vertixes and edges, where parallel edges between the vertexes are
used to represent strong pairing between the atoms, while single edges are
used to represent inter pair connections (see Fig. 1E). The graph is regular
if each vertex has the same the number of connected edges – the vertex degree.
Hence, we assign an alternating pattern of single and parallel edges so that
to preserve the fixed vertex degree for all vertexes. This is, however,
impossible at the end vertixes of the graph. Therefore, we introduce ghost
edges that are the triangular features in Fig. 1E. A ghost edge attaches to a
sole vertex. It contributes into the vertex degree, but does not contribute
into the adjacency matrix of the graph. Thus, an even SSH chain corresponds to
a graph with each vertex end being supplemented with a ghost edge as presented
in Fig. 1E so that the regular vertex degree is $3$, i.e. the vertex degree is
$n+1$, where $n$ is the number of nearest neighbors in the interior of the
system. If both end vertexes of the even SSH chain graph are supplemented with
$2$ ghost edges, as shown in Fig. 1E for the incommensurate case, the regular
vertex degree $=3$ can be preserved only if the non-adjacent parallel edges
change their pattern. The two patterns of the parallel edges are, in fact,
matchings (?) of the simple underline graph shown in Fig. 1E. The commensurate
even SSH chain has a parallel edge distribution following a perfect matching
of the simple underline graph and it corresponds to a trivial insulator of
periodic SSH chain. In contrast, in the incommensurate even SSH chain, the
parallel edges are distributed according to a maximal matching of the
underline simple graph (but neither to a maximum nor to perfect matching) and
this corresponds to a 1D ${\bf Z}_{2}$ topological insulator of periodic SSH
chain. In other words, a perfect matching in the graph is an indicator of the
trivial topological phase. For any finite simple graph with even number of
vertixes the sufficient and necessary condition for the existence of a perfect
matching is given by Tutte’s theorem (?). This can be applied to identify if a
topological phase transition is possible in any complex graphs, which includes
those beyond condensed matter systems such as neural networks in deep learning
and neuromorphic computing (?). The given graph-theoretic scheme naturally
covers the case of the odd SSH chain. Since graphs with odd number of vertexes
lack perfect matchings, the odd SSH chains are intrinsically incommensurate
or, for the brevity of notation, simply incommensurate. The given theoretical
formulation bridges 1D ${\bf Z}_{2}$ topological insulators with topologically
frustrated polyaromatic hydrocarbons (?). It sheds light on a recent tailoring
of topological order to obtain metallic structures (?, ?, ?).

##### Nullity theorem.

Based on the above consideration we can suggest an analog of the bulk-boundary
correspondence in counting the number of in-gap states in non-trivial
topological phases. For simple bipartite hexagonal graphs, the number of zero-
energy eigenvalues of the graph adjacency matrix is given by the nullity
theorem (?): $\eta=\alpha-\beta$, where $\alpha$ and $\beta$ are the maximum
numbers of pairwise nonadjacent vertices and edges of $G$, respectively. While
the nullity theorem is valid for even SSH chain in the commensurate case, it
fails to count the number of zero-energy eigenvalues for the even SSH chain in
the incommesurate case presented in Fig. 1E. On the other hand, the number of
unsaturated vertices for any maximal matching (i.e. it covers maximum and
perfect matchings too) is observed to be a better criterion: the eigenspectrum
of the adjacency matrix of a multigraph contains in-gap zero or quasi-zero-
energy eigenvalues if the pairwise non-adjacent parallel edges follow a non-
perfect matching pattern. The number of such in-gap eigenvalues is equal to
the number of unsaturated vertexes in the multigraph. The generality of this
principle is demonstrated for two multigraphs in Fig. 1F hosting the two
different matching patterns – perfect and maximal ones – which can be found in
graphene and phosphorene quantum dots (QDs). This result is in agreement with
sophisticated multi-parameter calculations for monolayer phosphorene QDs (?).

##### Incommensurate graphene nanoribbons.

Now we apply the above principles for analysis of one-dimensional honeycomb
lattice structures – GNRs. An infinite regular multigraphs with the vertex
degree $=4$ can be constructed for zigzag and bearded GNRs. The only way to
get a multigraph with a perfect matching pattern for zigzag GNR (ZGNR) is to
introduce a parallel double edge within the unit cell by attaching a single
ghost edge to the outmost vertex of ZGNR, as shown in Fig. 1G, and then to
continue the pattern of alternating single and parallel double edges until the
opposite boundary vertex is reached for which a single ghost edge is
introduced. A similar pattern of perfect matching can be obtained for a
bearded GNR. However, one needs to attach two ghost edges to the outmost
vertexes on both sides of the graph for a bearded GNR as depicted in Fig. 1G.
Since perfect matchings are constructed for both graphs, the in-gap states are
not guaranteed in both zigzag and bearded GNRs so that both structures can be
classified as commensurate ones. To ensure stable in-gap states, a non-trivial
topology of the system is needed. As shown in Fig. 1G, the combination of
zigzag and bearded edges breaking the inversion symmetry in a hbGNR makes
perfect matching impossible. Hence, the hbGNR is intrinsically incommensurate.
In other words, a topological in-gap mode is expected by considering only the
graph topology.

##### Flat zero-energy modes.

As shown in Fig. 2A, for the hbGNR with $N_{r}=51$ atoms in the unit cell, a
perfectly flat zero-energy mode (ZEM) is observed at the Fermi level (see also
supplementary text S1 and S2 and fig. S1). The flat character is favorable for
strong many-body effects. It shall drive this system into a strongly
correlated regime since any weak Coulomb interaction is large compared to the
width of the perfectly flat band. This ZEM also possesses topological
stability such that it does not depend on the details of the physical model
and it cannot be removed from the gap unless the ribbon structure is modified.
In a graph theoretical sense presented above, the ZEM is topological. However,
being obtained from pure graphene, it originates neither from the Chern nor 2D
${\bf Z}_{2}$ topological invariant. The valley Chern number is not its
reason, too, since the staggered sublattice potential is absent. In fact,
incommensurate hbGNRs exhibit distinct behaviour even when driven into the
said topological phases (see supplementary text S3). For all the three flavors
of topological insulators, the bulk boundary correspondence should be
interpreted with care when dealing with incommensurate systems. In particular,
the standard phase diagrams need accounting for the chemical potential (see
figs. S2-S4 and their discussion in supplementary text S3). Nevertheless, the
wave function of hbGNR ZEM does exhibit several peculiar similarities to the
states that are characterized by the topological invariants mentioned above.

The first peculiarity of ZEM is its chirality. By chirality we mean the
asymmetric distribution of the electron density for the band state wave
functions. In Fig. 2A, both the inverse participation ratio (IPR) and the
electron density show the shift of localization of ZEM from the bearded to
zigzag edge as $k$ changes from $0$ to $\pi$ (see methods for the IPR
definition). Thus, ZEM possesses a well-defined chirality everywhere except at
the Dirac K point, $k=2\pi/3$, where the chirality is ill defined since the
wave function extends over the whole width of the ribbon. As can be seen from
the IPR color scheme of Fig. 2A, the localization on a zigzag edge is higher
than that on the bearded. The hbGNR ZEM density behaviour differs from that of
a ZGNR having a unit cell with $N=52$ atoms, where the wave function is highly
localized on both edges, as presented in Fig. 1B. The hbGNR ZEM localization
is similar to the quantum valley Hall regime of the ZGNR presented for
comparison in Fig. 1C. Thus, the ZEM of the hbGNR is chiral. This result is
further supported by the continuum ${\bf k}\cdot{\bf p}$-model. The ZEM is the
Dirac equation solution that is an eigenstate of the chiral symmetry operator
given by the Pauli matrix $\sigma_{z}$: either
$\Psi_{+}(x)=(1,0)^{\mathrm{T}}\sqrt{\kappa_{y}/\sinh(\kappa_{y}L)}\,\exp\left(-\kappa_{y}x\right)$
or
$\Psi_{-}(x)=(0,1)^{\mathrm{T}}\sqrt{\kappa_{y}/\sinh(\kappa_{y}L)}\,\exp\left(\kappa_{y}x\right)$,
where $\kappa_{y}$ is the electron wave vector relative to the Dirac point and
$L$ is the ribbon width (see supplementary text S4 and fig. S5). In contrast
to a celebrated fermion-$1/2$ soliton (?), the ZEM solution is unbounded to
its topological defect which can be seen as a 1D analog of a 2D vortex of the
Kekulé pattern (?). This peculiarity is acquired because ZEM is normalized in
a finite space determined by the ribbon width rather than in the whole space,
as required for the soliton. We stress that such solutions are perfectly fine
in any bounded systems, therefore, under certain conditions they might even
result in some skin-effects for astrophysical superconducting strings (?).

The second peculiarity of hbGNR ZEM is the sublattice (pseudo-spin)
polarization throughout the whole Brillouin zone as clearly seen in Fig. 2A.
These numerical results are perfectly described by the given continuum model
solutions and are further supported by the analytical tight-binding model (see
supplementary text S2). The continuum model considered from the non-linear
dynamics and index theory viewpoint also reveals an unconventional topological
nature of this flat sublattice polarized mode. From this point of view,
however, the mode is characterized by a zero winding number but unlike the
bulk modes it resides in the vicinity of the unstable saddle fixed point of
the Dirac equation vector field (see supplementary text S4.3 and figs. S6 and
S7).

##### From flat modes to chiral anomalies.

Although the flatness of the mode is favorable for strong correlations and for
macroscopic coherent effects, such as superconductivity or ferromagnetism, it
is detrimental for the conventional electron transport. Hence, we speculate
the possibility of a dispersive in-gap mode arising form the sublattice mixing
introduced via an external electric field (?). Such a system can be considered
as the long-sought after 1D metal stable against the Peierls and Jahn–Teller
distortions. The external in-plane field causes a dispersion of the ZEM around
K and K′ points, as shown in the electronic band structure of hbGNR in Fig. 2D
(see also supplementary text S4.4 and fig. S8). The Fermi level can always be
adjusted to these points by the electrostatic doping through the back gate
voltage. As seen in Fig. 2E, the opposite direction of the electric field
flips the signs of the in-gap mode slopes in both valleys. The slopes are also
observed to be equal in magnitude but opposite in signs in the K and K′
valleys. The reciprocal space configuration of the ZEM is favorable for a
single channel dissipationless ballistic transport. The ZESs with the opposite
group velocities are far enough from each other and from the time reversal
(TR) invariant points in the $k$-space. Hence, the elastic backscattering is
suppressed, and this suppression is protected by the TR symmetry in a way
similar to the quantum valley Hall regime and topological kink states (?) (see
further continuum model topological accounts and how the chirality and TR
symmetry coexist here in supplementary text S4.5 and figs. S9-S10). Also, we
encounter here a momentum-valley locking: the charge current is valley
polarized so that this mode works as a perfect valley filter (?). Contrary to
the quantum valley Hall regime, the charge transport here does not necessarily
stick to the edges. Likewise, this system has no electrostatic domain walls
required for the topological kink states (?).

Applying a sufficiently large in-plane electric field, the dispersion of the
in-gap mode fills the bulk energy gap, thereby connecting the valence and
conduction bands. This regime is reminiscent of the edge modes in the
topological insulators as well as of the chiral anomaly of the Weyl semimetals
(WS) (?). However, in contrast to the topological insulators edge states, the
conductive mode here can be bulk which could be advantageous for managing heat
gradients in electronic devices, while in contrast to chiral anomaly of WS, in
the given case magnetic field is not needed, see Fig. 1A, which opens the way
to compact devices. The electric field required for filling the hbGNR bulk
energy gap can be reduced by increasing the hbGNR width. Figure 2F shows that
the ZEM band width is a linear function of the applied in-plane field in a
wide energy range. This linear dependence is also analytically justified
within the continuum model (see supplementary text S4.4 and fig. S8).

##### Cubic dispersion and pitchfork bifurcation.

A non-trivial behaviour is observed for wide ribbons and high external fields,
such as the $N_{r}=51$ with an in-plane electric field
$e\varepsilon/t_{1}>0.005$ Å-1. Once dispersive in-gap mode hybridizes with
the bulk subbands, its dispersion at the Dirac points changes from a linear to
cubic-like. Upon this hybridization, the dispersive in-gap mode breaks down
into three parts: two partly flat remnant modes and a new chiral gapless mode.
The gapless mode has a cubic dispersion without extrema and it features an
inflection point centered at the Dirac K point as shown in Fig. 3A. Upon
further increase of the external field and admixture of additional bulk bands,
the gapless mode transforms from having no extrema to having two local
extrema: one maximum (to the left/right) and one minimum (to the right/left)
in K/K′ valley, as shown in Fig. 3B. A detailed ZEM evolution towards a
gapless mode with cubic-type dispersion is presented in Fig. 3C. Thus, around
the K point, a smooth transformation from a completely flat to a linear and
finally a cubic character is observed. Here, the external field works as
control parameter which causes a pitchfork bifurcation. The dispersion curve
behaves as $E(k)=ak^{3}+ck$, which has two extrema for $c<0$ ($a>0$) and none
of them for $c>0$ ($a>0$). Not only does this breeding of the ZES with
opposite group velocities take similarity after the tripling of the fixed
points in the non-linear dynamical systems but it also replicates such effects
as critical slowing down that is characteristic for merging of the vector
field fixed points characterized by the opposite values of the topological
index (see supplementary text 4.6 and fig. S11). This analogy, however, is
approximate since, in effect, the hbGNR in the external field maps to a non-
stationary dynamical system featuring turbulent vector field (see
supplementary text S4.5 and Movie S1). The cubic nature of the gapless mode
dispersion is robust and supported by fitting the dispersion curve with a
cubic functional (see supplementary text S5 and fig. S12). Such a dispersion
is quite uncommon even for plasma waves even though it can be easily obtained
from the Korteweg–De Vries equation by neglecting its non-linearity (see
supplementary text S6). The only similar example is the dispersion of fast-
magnetosonic whistlers in the Earth’s bow shock frame (?). Similarly, tripling
of the Fermi levels has been only hypothesized for the edge states in a
peculiar context of the fractional quantum Hall effect (?).

The chiral character of the gapless mode persists despite all the drastic
transformations of the mode. In particular, the chiral charge defined as the
difference between the left and right movers at the Fermi level in each valley
is well-defined and preserved before and after the bifurcation (see
supplementary text 4.5). The chiral character of the gapless mode is also seen
from the electron density distributions before and after the bifurcation
presented in Figs. 3A and 3B. The characteristic feature of the bifurcation is
the change of the wave function localization from the zigzag to bearded edge
for $0<k<2\pi/3$ and from the bearded to zigzag edge for $2\pi/3<k<\pi$ as
compared to Fig. 2A,D, and E. The reported above rich properties of hbGNR ZEM
cannot be reproduced on a square lattice (see supplementary text S7 and fig.
S13).

##### Experimental detection.

The non-trivial behaviour exhibiting pitchfork bifurcation have consequences
in the density of states and quantum transport measurements. Figures 3D and 3E
present evolution the density of states (DOS) and transmission coefficients of
a hbGNR in the increasing in-plane electrostatic field.

The DOS tracks the flat bands and the pitchfork bifurcation. As seen from Fig.
3D, the initial DOS maximum from the fully flat ZEM splits into two lower
amplitude peaks. These peaks move towards the bulk bands and they eventually
mix with these bulk bands but do not lose their magnitude which indicates a
maintenance of the flatness feature in the band structure. The DOS between the
two peaks forms a nearly constant background. Upon reaching a critical in-
plane electric field a peak in the DOS arises at $E/t_{1}=0$ atop of the
uniform in-gap background. This peak corresponds to a cubic dispersion,
$E(k)\sim k^{3}$, presented in Fig. 3A. Further increase of the external field
splits the bifurcation DOS peak into two peaks with an inward looking
asymmetry, thereby marking the onset of the pitchfork bifurcation in Fig. 3B.
The two new peaks move into opposite directions, leaving behind an increased
DOS background between them.

The single channel electron transmission through the hbGNR changes from zero
to unity and finally reaches a multi-channel non-coherent transport regime. At
a zero external field, the transmission coefficient is zero throughout the
bulk gap of the hbGNR as shown in Fig. 3E. The absence of the transmission is
an indication of a Mott insulator-like regime. As the external field gradually
increases, the system develops a widening region of the perfect single channel
transmission originating from the dispersive in-gap mode discussed above. Once
the external field increases to the critical value corresponding to the
pitchfork bifurcation, a new transmission plateau appears on top of the in-gap
transmission pedestal which is proceeded by a short period when
$\partial_{k}E(k)_{k=\mathbf{K}}=0$ corresponding to a critical slowing down
in the quantum transport. The plateau is twice higher than the pedestal
signifying the onset of the multi-channel regime. As schematically presented
in Fig. 3E, the width of the plateau is bounded from left and right by the two
bifurcation peaks in DOS (dashed vertical lines in both Figs. 3D and 3E). The
width of the plateau determines the region of the cubic dispersion where the
energy within the valley is non-single valued, as schematically presented by
the dashed folded curve in Figs. 3D and 3E.

##### Classes of candidate GNRs.

The synthesis of one-dimensional organic metals remains challenging due to the
common characteristic of a polyradical ground state and consequently high
chemical reactivity (?, ?). In effect, a free-standing hbGNR which supports a
gapless chiral mode may be diffcult to experimentally prepare due to the
intrinsic high reactivity of localized radical states lining the edges. Thus,
we apply our prediction to a synthetically more realistic system that
recreates the hbGNR electronic structure – an asymmetrically hydrogenated ZGNR
(?). In this case, one zigzag edge of the ZGNRs is lined by trigonal planar
C–H groups while the opposite edge features tetrahedral methylene (CH2)
groups. This asymmetric substitution excludes one carbon atom per unit cell
from the $\pi$-network giving rise to the characteristic $\pi$-system of
hbGNRs. A ZGNR with $N=52$ featuring the asymmetric hydrogenation along the
zigzag edges is shown in Fig. 3F. Within this system, all the fundamental
results including the pitchfork bifurcation and electron density edge-to-edge
transfer for ZEM are perfectly replicated using density functional theory
(DFT) (see methods). The only difference between the DFT and tight-binding
results is an increase of the critical electric field by one order of
magnitude in the DFT calculations, which can be attributed to the screening
caused by the electric polarization of the metallic structure (?). Despite
this slight discrepancy, it is evident that the hbGNR $\pi$-orbital network
can in principle be modelled within the structure of asymmetrically
hydrogenated ZGNR.

Applying our basic concepts other families of incommensurate GNRs capable of
dissipationless coherent transport can be revealed. Namely, combining zigzag
and cove edges or zigzag and gulf edges in a single ribbon results in the
incommensurate zigzag-cove-edged GNRs (zcGNRs) and zigzag-gulf-edged GNRs
(zgGNRs), respectively. Even though the development of a practical synthesis
is still under investigation for zcGNRs and zgGNRs, recent advances in the
bottom-up synthesis of GNRs have shown that both classes of incommensurate
GNRs are synthetically accessible. One crucial difference between zcGNRs and
zgGNRs is in their translation period $T$. While zcGNRs possess a period
$T=2a$, the unit cell in zgGNRs is characterized by $T=3a$, where $a$ is the
graphene lattice constant. By using the cutting line method on the graphene
ﬁrst Brillouin zone (BZ) together with the zone folding, it easy to show that
the difference in the size of the unit cell between the ribbons will give rise
to two valleys centered around projected K and K′ points. In the case of
zcGNRs these valleys are shifted away from the TR invariant points $\Gamma$
and $X$ as demonstrated in Fig. 4A. In contrast, in the case of zgGNRs both K
and K′ points project to the TR invariant point $\Gamma$ as depicted in Fig.
4A so that these ribbons must exhibit only a single valley in their band
structures. These distinct features are evident in Fig. 4B, where DFT band
structures are presented for wide $N>50$ zcGNRs and zgGNRs. Moreover, both
zcGNRs and zgGNRs demonstrate flat ZEMs along with the above described valley
structures. In the presence of external in-plane electric field, both zcGNRs
and zgGNRs become metallic as seen in Fig. 4B. It is however only the zcGNRs
that exhibit the chiral anomaly structure similar to hbGNRs and are capable of
the dissipationless coherent transport. In general, the electronic band
structure of hbGNR-based incommensurate ribbon exhibits only a single valley
at $\Gamma$ if the translation period of an incommensurate ribbon is $T=3an$,
where $n$ is a natural number. Hence, such ribbons are ill-suited candidates
for the realization of a dissipationless conductor. Alternate families of
incommensurate GNRs with larger translation periods $T\neq 3an$ are
schematically shown in Fig. 4C. Wide ribbons of these GNR families featuring
the depicted edge geometries have three common properties: (i) two valley band
structures, (ii) a sufficiently large bulk energy gap, and (iii) a flat ZEM.
When subjected to an external in-plane electric ﬁeld, wide ribbons exhibit
chiral anomalies at low values of the ﬁeld while a pitchfork bifurcation can
be seen at higher values of the electric field. We have to note, however, that
further increase of the superlattice period of such incommensurate GNRs
results in a reduction of the two valley separation which eventually makes
inter-valley scattering more probable and is not favorable for dissipationless
transport.

##### Metamaterials perspective.

An alternative experimental verification can be obtained on analogous
ultracold atom lattices (?), photonic (?, ?, ?) and polaritonic (?, ?, ?)
crystals, and microwave metamaterials (?, ?, ?). Such artificially designed
systems do not have limitations of the synthetic chemistry, but the main
obstacle may be the imitation of the external electric field effect. The
proposed phenomena, however, can be readily verified in photonic and
nanomechanical metamaterials. In photonic crystals the pseudo-electric fields
can be created by the side-to-side laser writing (?). The second extremely
promising experimental playground is the platform of “nanomechanical
graphene”—a honeycomb lattice of free-standing Si3N4 nanomechanical membranes
with a modified lattice structure (?). Such nanomechanical graphene can be
fabricated with high precision by using electron-beam lithography and
subsequent dry and wet etching processes. The geometries of the nanomechanical
membranes are uniquely determined by the relative positions of the etched
holes $(r_{1},r_{2})=(r_{0}-\delta_{b},r_{0}+\delta_{b})$, as shown in Fig. 1C
of Ref. (?). The influence of an external in-plane electrostatic field to the
graphene honeycomb lattice can be imitated by tuning the on-site potential
wihtin the unit cell of a nanomechanical graphene ribbon. This can be realized
during the lithographic patterning by setting $\delta_{b}=0$ and linearly
varying $r_{0}$ across the width of the nanomechanical ribbon. As for
experimental characterization, one can electrically actuate the elastic waves
of the nanomechanical membranes and measure their propagating behaviors by
using optical interferometry. The energy band diagrams can be obtained by
first recording the real-space distribution of elastic waves along the desired
direction and then performing Fourier transform to project the signal to the
momentum space.

##### Conclusions.

In summary, a general graph-theoretic approach predicting in-gap zero-energy
modes is proposed to engineer dissipationless tranport regimes on the
honeycomb lattice. The bulk dissipationless transport is achievable in such
structures merely with an external electrostatic field, which demonstrates the
possibility to realize transport features of topological insulators and Weyl
semimetals under much more facile conditions. The identified dispersion
transformation at around the intrinsic Fermi level, ranging from fully flat to
cubic-like, paves the way to access a broad variety of notable quantum
transport regimes in a single experimental device. These results also show the
possibility of designing electric field tunable cubic materials suitable for
non-linear electronics and photonics. The cubic materials with the inflection
point in their dispersion available at the Fermi level may be promising
candidates for replacing complex semiconductor superlattices in observation of
Bloch oscillators and for implementation of submillimeter wavelength emitters
(?). Further studies may include other types of crystallographic system such
as kagome and Lieb lattices (?) or graphyne-based structures (?). Since the
reported flat band is intimately related to a fermion-$1/2$ soliton, it may
exhibit a fractional fermion number, which also should be a subject of future
research.

## Methods

### Tight-binding calculations

In the tight-binding calculations, we follow an idealized picture where a
single-parameter tight-binding (TB) model is used. All numerical calculations
are performed in a $p_{z}$-orbital approximation with TBpack Mathematica
package (?). Since even for the nearest neighbour hopping $t_{1}$ a wide range
of phenomenological values is available (?), we make results independent of
this varying quantity by defining all energies in terms of $t_{1}$.

A homogeneous in-plane electric field applied across the ribbon width is
modeled as on-site potentials (?): $U=-e{\bf\varepsilon}\cdot{\bf r}/t_{1}$,
where $e$ is the elementary charge, ${\bf\varepsilon}$ is the strength of the
field and ${\bf r}$ is the position of the lattice site. The hopping integral,
$t_{1}$, is set to be independent of the external field, since the largest
field applied to the ribbon ($e\varepsilon/t_{1}=0.1$ Å-1) is much less than
the atomic field $1/a_{0}\approx 0.7$ Å-1, where $a_{0}=1.42$ Å is the
graphene lattice constant. Here we omit the electrostrictive deformation and
screening effects. While the electrostrictive deformation in graphene
nanoribbons is small ($2.6\%$) and nanoribbon structures can be considered as
rigid ones (?), the screening for gapless metalic ribbons is significant.
Hence, the external field shall be interpreted as an effective field
experience by the on-site electron after screening. While solving the
eigenproblem for the matrix Hamiltonian, we keep tracking the eigenenergy
branches by sorting the eigenvalues with respect to maximum overlap of the
eigenvectors at the neighboring $k$-points. The sorting starts from the center
of the Brillouin zone, $k=0$, where ZEM is flat and easily identifiable.

To visualize the electron wave function localization in the energy band
diagram we use a color scheme based on the inverse participation ratio (IPR)
(?, ?, ?, ?): $\mathrm{IPR}=\sum_{{\bf r}}\left|\Psi({\bf r})\right|^{4}$,
where $\int_{V}\left|\Psi({\bf r})\right|^{2}\,d{\bf r}=1$. The IPR
characterizes the volume of space, where the absolute value of the wave
function is essentially non-zero and as such can be imagined as the size of
the eigenstate. For a normalized eigenvector
$C_{j}(k)=\left(C_{j1}(k),C_{j2}(k),\ldots,C_{jN}(k)\right)^{\mathrm{T}}$,
which corresponds to the $j$-th eigenstate $E_{j}(k)$ of an $N\times N$ TB
Hamiltonian, the IPR is given by

$\mathrm{IPR}_{j}(k)=\sum_{i=1}^{N}\left|C_{ji}(k)\right|^{4}\,.$ (1)

This quantity is equal to one for a perfectly localized state in which an
electron density occupies only one atomic site. However, the IPR it is much
less than unity for a fully delocalized state, and it further approaches zero
in the thermodynamic limit ($N\to\infty$) as the size of the system increases.

### Quantum transport

The coherent quantum transport through a hbGNR is modeled within the non-
equilibrium Green’s function formalism. We partition the full system
Hamiltonian as in Ref. (?):

$H=\left(\matrix{H_{L}&H_{LC}&0\cr H_{CL}&H_{C}&H_{CR}\cr
0&H_{RC}&H_{R}}\right),$ (2)

where $H_{C}$ is a square matrix with size equal to the number of atoms in the
central scattering region of the device, while $H_{L,R}$ and
$H_{LC},H_{CL},H_{CR},H_{RC}$ are semi-infinite square and rectangular
matrices describing the semi-infinite electrodes and their coupling to the
central scattering region, respectively. Here, the scattering region is a
single hbGNR unit cell, while the leads are represented by semi-infinite
hbGNRs [cf. with Ref. (?)]. The Green’s function defined as

$(H-IE)\mathcal{G}(E)=1\,,$ (3)

where $I$ is the identity matrix, can be partitioned as

$\mathcal{G}=\left(\matrix{\mathcal{G}_{L}&\mathcal{G}_{LC}&\mathcal{G}_{LR}\cr\mathcal{G}_{CL}&\mathcal{G}_{C}&\mathcal{G}_{CR}\cr\mathcal{G}_{RL}&\mathcal{G}_{RC}&\mathcal{G}_{R}}\right)\,.$
(4)

In terms of partition blocks, the central region Green’s function is

$\mathcal{G}_{C}=\left[EI-H_{C}-\Sigma_{L}-\Sigma_{R}\right]^{-1}\,,$ (5)

where the self-energies $\Sigma_{X}$ ($X=L,R$) are given by

$\Sigma_{X}(E)=H_{CX}g_{X}(E)H_{XC}$ (6)

with $g_{X}(E)=\left(EI-H_{X}\right)^{-1}$ being Green’s functions of the
semi-infinite leads. We calculate the semi-infinite lead Green’s functions via
the surface Green’s functions by defining them in self-consistent equations:

$g_{\mathrm{surf}}(E)=\left((E+i\eta)I-H_{0}-V^{\dagger}g_{\mathrm{surf}}(E)V\right)^{-1}\,,$
(7)

where $H_{0}$ and $V$ are diagonal and first off-diagonal blocks of the lead
region, respectively, and $\eta$ is a small real and positive parameter that
sets $g_{\mathrm{surf}}(E)$ as a retarded Green’s function.

Finally, we define the coupling matrices

$\Gamma_{X}=i\left(\Sigma_{X}-\Sigma^{\dagger}_{X}\right)=-2\mathrm{Im}(\Sigma_{X})\,,$
(8)

so that the transmission coefficient can be obtained as

$T(E)=\mathrm{Tr}\left[\Gamma_{L}(E)\mathcal{G}_{C}(E)\Gamma_{R}(E)\mathcal{G}_{C}^{\dagger}(E)\right]\,.$
(9)

Equation (9) can be reduced to

$T(E)=\mathrm{Tr}\left[\Gamma_{L,11}(E)\mathcal{G}_{C}(E)\Gamma_{R,NN}(E)\mathcal{G}_{C}^{\dagger}(E)\right]\,,$
(10)

if the full Hamiltonian tri-block diagonal structure is taken into account
(?). This approach speeds up the calculations of the transmission coefficient.
In Eq. (10), $\mathcal{G}_{C}$ must be replaced with a block
$\mathcal{G}_{C,1N}$ when the scattering region consists of $N$ unit cells. In
practice, the Hamiltonian containing two unit cells for each lead and one unit
cell as a scattering region was constructed by TBpack (?) to extract $H_{0}$,
$V$ and $H_{C}$ as initial parameters needed in the calculations.

We noticed that the convergence of the self-consistent procedure is guaranteed
only if a small real $\eta>0$ is added not only to the initial approximation
$g_{\mathrm{surf},0}(E)=\left((E+i\eta)I-H_{0}\right)^{-1}$ but also to the
repeatedly applied function (7). Also, a narrow peak can be obtained at
$E/t_{1}=0$ in the calculations if large values of $\eta=10^{-3}$ are used.
However, this peak disappears for small $\eta=10^{-5}$, and as such it is
deemed to be a numerical artefact rather than a physical feature. We have
verified that the given algorithm reproduces the results reported in Ref. (?).

### Density functional theory calculations

The spin-restricted DFT calculation was carried out in Siesta (?) with SISL
(?) as postprocessing tool. We employed Perdew–Burke–Ernzerhof (PBE)
generalized gradient approximation (GGA) exchange and correlation functional
(?), and Monkhorst-Pack k-mesh of $21\times 1\times 1$, and energy cut-off of
400 Ry. The structure was relaxed until the force was below $0.01$ eV/Å.
Electron charge density was visualized in XCrySDen (?).

## References

  * 1. H. K. Onnes, KNAW, Proc. 13, 1274 (1911).
  * 2. K. V. Klitzing, G. Dorda, M. Pepper, Phys. Rev. Lett. 45, 494 (1980).
  * 3. Y. Cao, et al., Nature 556, 43 (2018).
  * 4. Y. Cao, et al., Nature 556, 80 (2018).
  * 5. K. S. Novoselov, et al., Nature 438, 197 (2005).
  * 6. A. S. Davydov, Phys. Rep. 190, 191 (1990).
  * 7. V. F. Gantmakher, V. T. Dolgopolov, Physics-Uspekhi 53, 1 (2010).
  * 8. K. S. Novoselov, et al., Science 315, 1379 (2008).
  * 9. X.-L. L. Qi, S.-C. C. Zhang, Rev. Mod. Phys. 83, 1057 (2011).
  * 10. Y. Ren, Z. Qiao, Q. Niu, Reports Prog. Phys. 79, 066501 (2016).
  * 11. Z. Wang, S. Cheng, X. Liu, H. Jiang, Nanotechnology 32, 402001 (2021).
  * 12. Z.-F. Liu, Q.-P. Wu, X.-B. Xiao, Surfaces and Interfaces 25, 101300 (2021).
  * 13. N. W. Ashcroft, D. N. Mermin, Solid State Physics (Saunders College Publishing, Philadelphia, 1976), first edn.
  * 14. C. Rödl, F. Fuchs, J. Furthmüller, F. Bechstedt, Phys. Rev. B 79, 235114 (2009).
  * 15. J. Gebhardt, C. Elsässer, J. Phys. Condens. Matter 35, 205901 (2023).
  * 16. J.-P. Issi, Aust. J. Phys. 32, 585 (1979).
  * 17. K. S. Novoselov, et al., Science 306, 666 (2004).
  * 18. B. Q. Lv, et al., Phys. Rev. X 5, 031013 (2015).
  * 19. J. K. Asbóth, L. Oroszlány, A. Pályi, A Short Course on Topological Insulators, vol. 919 of Lecture Notes in Physics (Springer, Cham, 2016), first edn.
  * 20. J. N. Fuchs, F. Piéchon, Phys. Rev. B 104, 235428 (2021).
  * 21. L. Lovász, M. Plummer, Matching Theory (North-Holland, New York, 1986).
  * 22. W. T. Tutte, J. London Math. Soc. s1-22, 107 (1947).
  * 23. C. D. Schuman, et al., Nat. Comput. Sci. 2, 10 (2022).
  * 24. S. Mishra, et al., Nat. Nanotechnol. 15, 22 (2020).
  * 25. B. Cirera, et al., Nat. Nanotechnol. 15, 437 (2020).
  * 26. D. J. Rizzo, et al., Science 369, 1597 (2020).
  * 27. R. D. McCurdy, et al., Engineering robust metallic zero-mode states in olympicene graphene nanoribbons (2023). arXiv:2302.08446.
  * 28. S. Fajtlowicz, P. E. John, H. Sachs, Croat. Chem. Acta 78, 195 (2005).
  * 29. V. A. Saroka, I. Lukyanchuk, M. E. Portnoi, H. Abdelsalam, Phys. Rev. B 96, 085436 (2017).
  * 30. R. Jackiw, C. Rebbi, Phys. Rev. D 13, 3398 (1976).
  * 31. C. Y. Hou, C. Chamon, C. Mudry, Phys. Rev. Lett. 98, 186809 (2007).
  * 32. E. Witten, Nuclear Physics B 249, 557 (1985).
  * 33. V. A. Saroka, K. G. Batrakov, Physics, Chem. Appl. Nanostructures, V. E. Borisenko, S. V. Gaponenko, V. S. Gurin, C. H. Kam, eds. (World Scientific, Singapore, 2015), pp. 240–243.
  * 34. A. Rycerz, J. Tworzydło, C. W. J. Beenakker, Nat. Phys. 3, 172 (2007).
  * 35. B. Yan, C. Felser, Annu. Rev. Condens. Matter Phys. 87, 337 (2017).
  * 36. D. Krauss-Varban, N. Omidi, J. Geophys. Res. 96, 17715 (1991).
  * 37. C. de C. Chamon, X. G. Wen, Physical Review B 49, 8227 (1994).
  * 38. K. Kusakabe, M. Maruyama, Phys. Rev. B 67, 092406 (2003).
  * 39. P. Gava, M. Lazzeri, A. M. Saitta, F. Mauri, Physical Review B 79, 165431 (2009).
  * 40. G. Jotzu, et al., Nature 515, 237 (2014).
  * 41. Y. Plotnik, et al., Nat. Mater. 13, 57 (2013).
  * 42. T. Ozawa, et al., Rev. Mod. Phys. 91, 015006 (2019).
  * 43. S. Xia, et al., Physical Review Letters 131, 013804 (2023).
  * 44. M. Milićević, et al., 2D Mater. 2, 034012 (2015).
  * 45. C. E. Whittaker, et al., Phys. Rev. B 99, 081402 (2019).
  * 46. C. E. Whittaker, et al., Nat. Photonics 15, 193 (2021).
  * 47. M. Bellec, U. Kuhl, G. Montambaux, F. Mortessagne, New J. Phys. 16, 113023 (2014).
  * 48. Y. N. Dautova, A. V. Shytov, I. R. Hooper, J. R. Sambles, A. P. Hibbins, Appl. Phys. Lett. 110, 261605 (2017).
  * 49. Y. N. Dautova, A. V. Shytov, I. R. Hooper, J. R. Sambles, A. P. Hibbins, Appl. Phys. Lett. 112, 191601 (2018).
  * 50. X. Liu, et al., Laser and Photonics Reviews 15, 2000563 (2021).
  * 51. X. Xi, J. Ma, S. Wan, C.-H. Dong, X. Sun, Sci. Adv. 7, eabe1398 (2021).
  * 52. L. Esaki, R. Tsu, IBM J. Res. Dev. 14, 61 (1970).
  * 53. L. Yan, P. Liljeroth, Adv. Phys. X 4, 1651672 (2019).
  * 54. J. Kang, Z. Wei, J. Li, ACS Appl. Mater. Interfaces 11, 2692 (2019).
  * 55. V. A. Saroka, TBpack (version 0.4.0 and higher) (2021).
  * 56. R. B. Payod, et al., Nat. Commun. 11, 82 (2020).
  * 57. C. Chang, et al., Carbon 44, 508 (2006).
  * 58. V. M. Martinez Alvarez, J. E. Barrios Vargas, M. Berdakin, L. E. F. Foa Torres, Eur. Phys. J. Spec. Top. 227, 1295 (2018).
  * 59. F. Wegner, Z. Phys. B 36, 209 (1980).
  * 60. D. Thouless, Phys. Rep. 13, 93 (1974).
  * 61. R. J. Bell, P. Dean, Discuss. Faraday Soc. 50, 55 (1970).
  * 62. F. Muñoz-Rojas, D. Jacob, J. Fernández-Rossier, J. J. Palacios, Phys. Rev. B 74, 195417 (2006).
  * 63. K. Wakabayashi, Phys. Rev. B 64, 125428 (2001).
  * 64. S. Compernolle, L. Chibotaru, A. Ceulemans, J. Chem. Phys. 119, 2854 (2003).
  * 65. J. M. Soler, et al., Journal of Physics: Condensed Matter 14, 2745 (2002).
  * 66. N. Papior, J. L. B., pfebrer, T. Frederiksen, S. S. Wuhl, zerothi/sisl: v0.11.0 (2021).
  * 67. J. P. Perdew, K. Burke, M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).
  * 68. A. Kokalj, Computational Materials Science 28, 155 (2003).
  * 69. V. A. Saroka, M. V. Shuba, M. E. Portnoi, Phys. Rev. B 95, 155438 (2017).
  * 70. D. Klein, Chem. Phys. Lett. 217, 261 (1994).
  * 71. K. Nakada, M. Fujita, G. Dresselhaus, M. S. Dresselhaus, Phys. Rev. B 54, 17954 (1996).
  * 72. V. Saroka, A. Pushkarchuk, S. Kuten, M. Portnoi, J. Saudi Chem. Soc. 22, 985 (2018).
  * 73. D. Cvetković, I. Gutman, N. Trinajstić, Theor. Chim. Acta 34, 129 (1974).
  * 74. C. T. White, J. Li, D. Gunlycke, J. W. Mintmire, Nano Lett. 7, 825 (2007).
  * 75. R. R. Hartmann, V. A. Saroka, M. E. Portnoi, J. Appl. Phys. 125, 151607 (2019).
  * 76. I. S. Gradshteyn, I. M. Ryzhik, Table of Integrals, Series, and Products (Academic Press, San Diego, 2000), 6th edn.
  * 77. F. D. M. Haldane, Phys. Rev. Lett. 61, 2015 (1988).
  * 78. C. L. Kane, E. J. Mele, Phys. Rev. Lett. 95, 226801 (2005).
  * 79. M. Z. Hasan, C. L. Kane, Rev. Mod. Phys. 82, 3045 (2010).
  * 80. C. L. Kane, E. J. Mele, Phys. Rev. Lett. 95, 146802 (2005).
  * 81. L. Fu, C. L. Kane, Phys. Rev. B 74, 195312 (2006).
  * 82. W. Yao, S. A. Yang, Q. Niu, Phys. Rev. Lett. 102, 096801 (2009).
  * 83. B. A. Benevig, T. L. Hughes, Topological Insulators and Topological Superconductors (Princeton University Press, Princeton, New Jersey, 2013), second edn.
  * 84. D. P. DiVincenzo, E. J. Mele, Phys. Rev. B 29, 1685 (1984).
  * 85. L. Brey, H. A. Fertig, Phys. Rev. B 73, 235411 (2006).
  * 86. A. Akhmerov, C. Beenakker, Phys. Rev. B 77, 085423 (2008).
  * 87. E. McCann, V. I. Falko, J. Phys. Condens. Matter 16, 2371 (2004).
  * 88. A. Niemi, G. Semenoff, Physics Reports 135, 99 (1986).
  * 89. R. Jackiw, Phys. Scr. T146, 014005 (2012).
  * 90. T. Yanagisawa, Symmetry 12, 373 (2020).
  * 91. S.-Q. Shen, Topological Insulators: Dirac Equation in Condensed Matters (Springer-Verlag, Berlin, Heidelberg, 2012), vol. 174 of Springer Series in Solid-State Sciences, p. 225.
  * 92. R. Jackiw, P. Rossi, Nucl. Phys. B 190, 681 (1981).
  * 93. S. Ryu, Y. Hatsugai, Phys. Rev. Lett. 89, 077002 (2002).
  * 94. J. Cayssol, J. N. Fuchs, J. Phys. Mater. 4, 034007 (2021).
  * 95. L. Fu, C. Kane, Phys. Rev. Lett. 100, 096407 (2008).
  * 96. C. Chamon, et al., Phys. Rev. Lett. 100, 110405 (2008).
  * 97. S. H. Strogatz, Nonlinear Dynamics and Chaos: With Applications to Physics, Biology, Chemistry and Engineering (Perseus Books, Reading, Massachusetts, 1994), first edn.
  * 98. M. Ezawa, New J. Phys. 16, 115004 (2014).
  * 99. M. M. Grujić, M. Ezawa, M. Ž. Tadić, F. M. Peeters, Phys. Rev. B 93, 245413 (2016).
  * 100. S. L. Adler, Phys. Rev. 177, 2426 (1969).
  * 101. J. S. Bell, R. Jackiw, Nuovo Cim. A 60, 47 (1969).
  * 102. H. Nielsen, M. Ninomiya, Phys. Lett. B 130, 389 (1983).

## Acknowledgments

The authors thank K. Batrakov, E. Thingstad, J. Zheng, V. Demin, A.
Qaiumzadeh, A. Ferreira, and C. Beenakker for useful stimulating discussions
and J. Danon for providing computational facilities at NTNU.

##### Funding:

X.S. acknowledges funding support by the Research Grants Council of Hong Kong
(project no. 14209519). C.A.D. is supported by the Royal Society University
Research Fellowship (URF/ R1/ 201158). F.R.F. was partly supported by the
Heising-Simons Faculty Fellows Program at the UC Berkeley. V.A.S. was partly
supported by the Research Council of Norway Center of Excellence funding
scheme (project no. 262633, “QuSpin”) and EU HORIZON-MSCA-2021-PF-01 (project
no. 101065500, TeraExc).

##### Authors contributions:

V.A.S. and C.A.D. conceived the initial idea. V.A.S. developed graph
theoretical description and performed TB calculations. R.B.P. verified the
results. V.A.S and C.A.D. interpreted TB results within the continuum
effective mass model. X.S. proposed metamaterial platform and experimental
arrangement for testing predictions. F.R.F. proposed a zcGNR family for
investigation and realization by synthetic chemistry methods. F.K., V.A.S. and
L.B. performed DFT calculations and interpreted the results. V.A.S., R.B.P.,
F.K, X.S. and L.B. wrote the first draft of the manuscript. All authors
contributed to the further editing of the manuscript.

##### Competing interests:

The authors declare no competing interests.

##### Data and materials availability:

All data are available in the manuscript or the supplementary materials.
Additional data supporting the findings of this study are available from the
corresponding authors upon reasonable request. The code developed for this
study is available under MIT License from GitHub (?).

## Supplementary materials

Supplementary Text
Figs. S1 to S13
References (67-102)
Movie S1

Fig. 1: Zero-energy states and graph topology. (A) Dissipationless transport
schemes. Ivn: breaking of the inversion symmetry; E and B: electric and
magnetic fields. (B) Periodic SSH model. $t_{1,2}$ are hopping integrals,
$d_{1,2}$ are inter site distances. (C) Energy levels of the finite even SSH
chain as function of $t_{1}/t_{2}$ and the wave functions of $E=0$ states (red
curves). The light blue and white backgrounds mark regions of topological and
trivial phases, respectively. The negative and positive on-site weights of the
wave function are denoted as blue and red. For visibility the weights are
scaled by the factor shown to the left. The thin and thick gray lines
connecting the SSH chain sites visualize strengths of $t_{1}$ and $t_{2}$. (D)
Same as (C) but for the finite odd SSH chain. (E), The regular multigraphs for
a commensurate and incommensurate finite even SSH chain and a corresponding
simple graph. (F) Counting the in-gap states for two regular multigraphs and
matching their patterns corresponding to two chemical systems. Plot:
eigenvalues of adjacency matrices versus the normalized eigenvalue counting
number. (G) Regular multigraphs for various GNRs with translation invariance
along to the zigzag crystallographic direction. The dot dashed rectangle in
(B) and (G) marks the unit cells of periodic structures.  Fig. 2: Converting a
flat ZEM into a dispersive in-gap mode. (A) Energy bands and electron density
distribution over the unit cell of a ribbon at $3$ depicted $k$-points in the
K valley of the hbGNR $N_{r}=51$. Colour bars: the inverse participation ratio
(IPR) of each eigenstate. Circles: on-site electron densities scaled by a
factor shown to the left. (B) Same as (A), but for ZGNR $N=52$. (C) Same as
(B), but for the quantum valley Hall regime: $t_{1}=1$; $\Delta=0.2t_{1}$;
$U=0.3t_{1}$; $U_{l}=-U$; $U_{r}=U$, where $\Delta$ is the amplitude of the
on-site staggered potential, $U$ is the side back gate potential amplitude,
$U_{l,r}$ are the left and right side back gate potential that is applied to
the $12$ atoms on each side of the ribbon. (D) Same as (A), but for hbGNR
$N_{r}=51$ in the in-plane electric field $e\varepsilon/t_{1}=0.005$ Å-1. (E)
Same as (D), but for the opposite direction of the applied field and focused
at the zero energy of the band structure. (F) A zero-energy band width as a
function of the in-plane electrostatic field for hbGNR $N_{r}=11$, $21$, $31$,
$41$ and $51$. Scheme: the in-plane electric field applied across the ribbon.
Fig. 3: Chiral gapless mode pitchfork bifurcation. (A) The energy bands of
hbGNR $N_{r}=51$ in the electric field $e\varepsilon/t_{1}=0.012$ Å-1 together
with the electron density distribution over the ribbon unit cell for the
selected zero-energy mode energy levels at several $k$-points. (B) Same as (A)
but for $e\varepsilon/t_{1}=0.015$ Å-1. The circles representing on-site
electron densities are scaled by factors shown to the left. In (A) and (B)
dashed blue dispersion curves are the remnant partially flat bands and light
red curves are the chiral gapless modes. (C) hbGNR $N_{r}=51$ ZEM evolution in
the electric fields from $e\varepsilon/t_{1}=0.0$ to $0.015$ Å-1 with the step
$0.0015$ Å-1 marked with circled numbers $1$-$11$, respectively. (D), (E) The
density of states and the transmission coefficient in the electric fields in
(C). Red dots in (D) mark the peaks used for normalization of all the curves.
The spectra are shifted vertically for clarity. The natural horizontal drift
to the left is also retained for clarity. The wiggly dashed brown curve
represents the dispersion of the chiral gapless mode for
$e\varepsilon/t_{1}=0.015$ Å-1. The vertical dashed gray lines mark the
energies of the peaks in the density of states arising after the bifurcation
onset. (F) DFT energy bands and electron densities in a few $k$-points for
ZGNR $N=52$ asymmetrically hydrogenated as shown in the scheme above the plots
and subjected to $0$, $0.6$ and $1$ V/Å in-plane electric fields. The same
isovalue is used for imaging the charge densities. Fig. 4: Dissipationless
conductors. (A) The folding scheme explaining $3a$-periodicity of the single
valley appearance in intrinsically incommensurate GNRs. $a=2.46$ Å is the
graphene lattice constant. Black hexagons represent the first Brillouin zone
(BZ) of graphene. Light green lines are the cutting lines resulting from the
transverse momentum quantization in a GNR. BZs of GNRs with different
translation period $T$ are denoted by light blue curly brakets. (B) DFT energy
bands for zcGNR and zgGNR species with $25$ zigzag chains in the backbone ZGNR
and subject to the in-plane electric fields. (C) Schematic illustrations of
intrinsically incommensurate GNRs with the dissipationless transport
properties in the presence of an external in-plane electric field.

Supplementary Materials for

Tunable chiral anomalies and coherent transport on a honeycomb lattice

Vasil A. Saroka, Fanmiao Kong, Charles A. Downing, Renebeth B. Payod, Felix R.
Fischer, Xiankai Sun, Lapo Bogani

Corresponding author: Vasil Saroka<EMAIL_ADDRESS>
##### The PDF file includes:

  * Supplementary Text

  * Figs. S1 to S13

  * References

##### Other Supplementary Materials for this manuscript include the following:

  * Movie S1

## Supplementary Text

##  S1 Zero-energy mode and iso-spectral properties of hbGNRs and carbon
nanotubes

Let us take a closer look at the electronic properties of the half-bearded
graphene nanoribbons (hbGNRs). Figure S1A shows the energy band structure of
hbGNRs compared with that of an armchair single-wall carbon nanotube (aSWCNT).
An aSWCNT is a good reference since its energy band curves contour the full
spectrum of a bulk graphene. This feature is conveniently used to track
peculiar edge states that evade the bulk graphene spectrum (?). Indeed, for
the hbGNR a perfectly flat zero-energy band is seen in fig. S1A. This flat
band is fully outside the bulk graphene spectrum, which is a deep contrast to
the picture when zigzag and bearded graphene nanoribbons (GNRs) are compared
to aSWCNTs and partly flat bands only partially escape outside the bulk
graphene spectrum (?, ?, ?). A couple of other observations shall be mentioned
with respect to the energy bands of hbGNRs and aSWCNTs. The bulk energy bands
of the hbGNR are isospectral with the double degenerate energy bands of the
chosen single-wall carbon nanotube. This effect is of the same type as that
for $(1,4)$-divinyl benzene and $2$-phenyl butadiene and is rooted in the
spectral graph theory (?). The observed effect for hbGNRs differs from that of
both zigzag and bearded GNRs, thus making hbGNRs more similar to armchair
GNRs. (cf. Fig. 3 in Ref. (?), Fig. 2 in Ref. (?), and Fig. 6 in Ref. (?)). It
follows from fig. S1A that for an aSWCNT and a pair of hbGNRs there is a rule
of decomposition presented in fig. S1B. Namely, $N_{t}=2N_{r}+2$, where
$N_{t}(N_{r})$ is the number of atoms in the aSWCNT (hbGNR) unit cell. This
rule is distinct from the known rules for zigzag and armchair GNRs/carbon
nanotubes, $N_{t}=2N_{r}+4$, and bearded GNRs/aSWCNTs, $N_{t}=2N_{r}$ (?),
thereby making hbGNRs a special class of material as is expected from their
unique structural properties. Later we shall use this iso-spectral property to
find a correct description of the hbGNR zero-energy mode within the continuum
model.

##  S2 Analytical theory of hbGNRs

This section provides a detailed analytical description of the zero-energy
mode (ZEM) of hbGNR when an external electric field is absent. On the other
hand, it demonstrates the persistent character of the zero-energy solution in
the hbGNR when an external electrostatic field is applied in-plane across the
ribbon width.

###  S2.1 Pristine hbGNR

In this subsection, we shall treat the problem of hbGNR energy bands and wave
functions with the transfer matrix method from first principles without making
assumptions about the zero-energy states. This treatment is based on that for
zigzag GNR in Ref. (?). We start from the tight-binding Hamiltonian that can
be obtained from the Hamiltonian of a zigzag (Z) or bearded (B) GNR by its
extension to the odd size:

$\displaystyle H_{\mathrm{ZGNR}}$ $\displaystyle=$
$\displaystyle\left(\matrix{0&t_{1}q&0&0\cr t_{1}q&0&t_{1}&0\cr
0&t_{1}&0&t_{1}q\cr 0&0&t_{1}q&0}\right)\,,$ (S1) $\displaystyle
H_{\mathrm{BGNR}}$ $\displaystyle=$ $\displaystyle\left(\matrix{0&t_{1}&0&0\cr
t_{1}&0&t_{1}q&0\cr 0&t_{1}q&0&t_{1}\cr 0&0&t_{1}&0}\right)\,,$ (S2)
$\displaystyle H_{\mathrm{hbGNR}}$ $\displaystyle=$
$\displaystyle\left(\matrix{0&t_{1}q&0&0&0\cr t_{1}q&0&t_{1}&0&0\cr
0&t_{1}&0&t_{1}q&0\cr 0&0&t_{1}q&0&t_{1}\cr 0&0&0&t_{1}&0}\right)\,,$ (S3)

where $q=2\cos\left(k/2\right)$, $k\equiv ka$ with $a=2.46$ Å being graphene
lattice constant, and $t_{1}$ is the nearest-neighbor hopping integral. The
eigenproblem for hbGNR Hamiltonian (S3) can be presented in the transfer matix
form with the two transfer matrices $T_{1}$ and $T_{2}$:

$\left(\matrix{c_{N}\cr
c_{N+1}}\right)=T_{1}T^{(N_{r}-1)/2}\left(\matrix{c_{0}\cr c_{1}}\right)\,,$
(S4)

where $N_{r}$ is the number of atoms in the hbGNR unit cell, $T=T_{1}T_{2}$,
and

$\displaystyle T_{1}$ $\displaystyle=$
$\displaystyle\left(\matrix{0&1\cr-1/q&\alpha/q}\right)\,,$ (S5)
$\displaystyle T_{2}$ $\displaystyle=$
$\displaystyle\left(\matrix{0&1\cr-q&\alpha}\right)\,,$ (S6)

with $\alpha=E/t_{1}$ being the dimensionless energy. Since $N_{r}$ is odd for
hbGNR, the power $\left(N_{r}-1\right)/2$ is actually the parameter $w=N/2$
for ZGNR in Ref. (?). Hence, the power of the transfer matrix $T$ can be
evaluated by using Eq. 18 in Ref. (?) and later substituting
$w=\left(N_{r}-1\right)/2$. For $T^{w}$, we explicitly have

$T^{w}=\frac{(-1)^{w}}{q\sin\theta}\left(\matrix{-\left\\{\sin\left(w\theta\right)+q\sin\left[\left(w+1\right)\theta\right]\right\\}&\alpha\sin\left(w\theta\right)\cr-\frac{1+q^{2}-2q\cos\theta}{\alpha}\sin\left(w\theta\right)&\sin\left(w\theta\right)-q\sin\left[\left(w+1\right)\theta\right]}\right)\,,$
(S7)

where the new parameter $\theta$ reduces eigenvalues of $T$ to the exponential
form $\lambda_{1,2}=-e^{\mp i\theta}$, i.e. it results from
$(\alpha^{2}-q^{2}-1)/(2q)=-\cos\theta$ that corresponds to a direct band
numbering in terminology of Ref. (?). Symbolically, the full transfer matrix
$\mathbf{T}=T_{1}T^{w}$ is

$\displaystyle\mathbf{T}$ $\displaystyle=$
$\displaystyle\left(\matrix{0&1\cr-\frac{1}{q}&\frac{\alpha}{q}}\right)\left(\matrix{\left(T^{w}\right)_{11}&\left(T^{w}\right)_{12}\cr\left(T^{w}\right)_{21}&\left(T^{w}\right)_{22}}\right)\,,$
(S8) $\displaystyle=$
$\displaystyle\left(\matrix{\left(T^{w}\right)_{21}&\left(T^{w}\right)_{22}\cr\left(-\left(T^{w}\right)_{11}+\alpha\left(T^{w}\right)_{21}\right)/q&\left[-\left(T^{w}\right)_{21}+\alpha\left(T^{w}\right)_{22}\right]/q}\right)\,,$

whence the secular equation for the hard-wall boundary condition
$\mathbf{T}_{22}=0$ is

$\displaystyle\mathbf{T}_{22}$ $\displaystyle=$
$\displaystyle\frac{1}{q}\left(-\frac{\alpha(-1)^{w}}{q\sin\theta}\sin\left(w\theta\right)+\frac{\alpha(-1)^{w}}{q\sin\theta}\left\\{\sin\left(w\theta\right)-q\sin\left[\left(w+1\right)\theta\right]\right\\}\right)\,,$
(S9) $\displaystyle=$
$\displaystyle\frac{(-1)^{(w+1)}\alpha}{q\sin\theta}\sin\left[\left(w+1\right)\theta\right]=0\,,$

where we have used Eq. (S7). As seen from Eq. (S9), the secular equation is
satisfied if either (i) $\alpha=0$ or (ii) $\theta=\pi j/(w+1)$. The (i)
option yields the ZEM that is independent of $k$ and fully flat. The (ii)
option describes the bulk bands that are

$E=\pm
t_{1}\sqrt{4\cos^{2}\left(\frac{k}{2}\right)-4\cos\left(\frac{k}{2}\right)\cos\theta+1}\,.$
(S10)

One can notice that Eq. (S10) is identical to that for aSWCNTs (see Appendix C
of Ref. (?)) if

$\displaystyle\theta_{t,j}$ $\displaystyle=$ $\displaystyle\theta_{r,j}\,,$
$\displaystyle\frac{2\pi j}{N_{t}/2}$ $\displaystyle=$ $\displaystyle\frac{\pi
j}{(N_{r}-1)/2+1}\,,$ $\displaystyle N_{t}$ $\displaystyle=$ $\displaystyle
2N_{r}+2\,,$ (S11)

where we have used $w=(N_{r}-1)/2$. Equations (S10) and (S11) explain the
isospectral properties of hbGNRs and aSWCNTs presented in fig. S1A. The found
$\theta$’s allow us to obtain the eigenvectors of the hbGNR bulk bands. The
components of the eigenvectors can be found from the following recursive
relations [cf. with Eq. (S4)]:

$\displaystyle\left(\matrix{c_{j}\cr c_{j+1}}\right)$ $\displaystyle=$
$\displaystyle T_{1}T^{(j-1)/2}\left(\matrix{c_{0}\cr
c_{1}}\right)\,,\qquad\mathrm{odd}\,j$ $\displaystyle\left(\matrix{c_{j}\cr
c_{j+1}}\right)$ $\displaystyle=$ $\displaystyle T^{j/2}\left(\matrix{c_{0}\cr
c_{1}}\right)\,,\qquad\mathrm{even}\,j$ (S12)

where $c_{0}=0$ is required by the hard-wall boundary condition and $c_{1}$
can be found from the eigenvector normalization. The latter means that
$c_{j}=\left(T_{1}T^{(j-1)/2}\right)_{12}c_{1}$, when $j=2p+1$, and
$c_{j}=\left(T^{j/2}\right)_{12}c_{1}$, when $j=2p$, while $p=0,\ldots,w$ in
both cases. After reordering $c_{2p+1}\rightarrow c_{2p-1}$, the components
can be explicitly written as

$\displaystyle c^{(j)}_{2p-1}$ $\displaystyle=$ $\displaystyle
N_{j}\left\\{\sin\left[\left(p-1\right)\theta_{j}\right]-q\sin\left(p\theta_{j}\right)\right\\}\,,$
$\displaystyle c^{(j)}_{2p}$ $\displaystyle=$ $\displaystyle
N_{j}\alpha_{j,\pm}\sin\left(p\theta_{j}\right),\qquad p=1,\ldots,w+1\,,$
(S13)

where $\theta_{j}=\pi j/(w+1)$, $w=(N_{r}-1)/2$,
$\alpha_{j,\pm}=\pm\sqrt{q^{2}-2q\cos\theta_{j}+1}$ and

$\displaystyle N_{j}$ $\displaystyle=$
$\displaystyle\left(\sqrt{S^{(1)}_{\theta_{j},w}-2qS^{(2)}_{\theta_{j},w}+\left(2q^{2}-2q\cos\theta_{j}+1\right)S^{(3)}_{\theta_{j},w}}\right)^{-1}\,,$
(S14)

with $S^{(1,2,3)}_{\theta_{j},w}$ given by

$\displaystyle S^{(1)}_{\theta_{j},w}$ $\displaystyle=$
$\displaystyle\frac{1}{4}\left(1+2w-\frac{\sin\left[\left(2w+1\right)\theta_{j}\right]}{\sin\theta_{j}}\right)\,,$
(S15) $\displaystyle S^{(2)}_{\theta_{j},w}$ $\displaystyle=$
$\displaystyle\frac{1}{4}\left(2w\cos\theta_{j}+\frac{\sin
2\theta_{j}}{\sin\theta_{j}}-\frac{\sin\left[2\left(w+1\right)\theta_{j}\right]}{\sin\theta_{j}}\right)\,,$
(S16) $\displaystyle S^{(3)}_{\theta_{j},w}$ $\displaystyle=$
$\displaystyle\frac{1}{4}\left(3+2w-\frac{\sin\left[\left(2w+3\right)\theta_{j}\right]}{\sin\theta_{j}}\right)\,.$
(S17)

It should be noted that $\theta_{j}$ does not depend on $k$.

The ZEM case is more tricky, since we do not have $\theta_{j}$ to substitute
into the recursive relations given by Eq. (S12). In principle, one can notice
that the secular equation solution $\alpha=0$ formally defines $\theta$ from
$\pm\sqrt{q^{2}-2q\cos\theta+1}=0$ as
$\theta=\arccos\left((q^{2}+1)/2q\right)$. However,
$\left|(q^{2}+1)/2q\right|\geq 1$ for any $q\in\left(-\infty,+\infty\right)$,
making the obtained formal result meaningless. Similar to the treatment of
ZGNR edge states in Ref. (?), the found pathology can be cured by the analytic
continuation: $\theta\rightarrow i\beta$. This substitution necessarily
changes the trigonometric functions in Eqs. (S13) to the hyperbolic ones,
thereby enforcing the localization of eigenvector on a few or even one
component. Noticing also that $\alpha=0$ sets all even components of the
eigenvector to zero, we finally have

$\displaystyle c^{(j)}_{2p-1}$ $\displaystyle=$ $\displaystyle
N_{0}\left\\{\sinh\left[\left(p-1\right)\beta_{j}\right]-q\sin\left(p\beta_{j}\right)\right\\}\,$
$\displaystyle c^{(j)}_{2p}$ $\displaystyle=$ $\displaystyle 0,\qquad
p=1,\ldots,w+1\,,$ (S18)

where $\beta=\mathrm{arcosh}\left((q^{2}+1)/2q\right)$, and the normalization
constant is

$N_{0}=\left(\sqrt{\sum^{w+1}_{p=1}\left\\{\frac{\sinh\left[(p-1)\beta\right]}{\sinh\left(p\beta\right)}-q\right\\}^{2}\sinh^{2}\left(p\beta\right)}\right)^{-1}\,.$
(S19)

A simplified expression of Eq. (S19) with summations similar to the Eq. (S14)
can be obtained. Those expressions, however, are not practical from a
numerical point of view because of the need to operate with large numbers
resulting from the hyperbolic functions when $q\rightarrow 0$ and
$\beta\rightarrow\pm\infty$. In this case, the lost of precision leads to a
non-regular oscillatory behavior in the region close to $k=\pi$. In addition,
for an increasing width of the ribbon this numerically non-stable region
promptly expands from $k=\pi$ towards the Dirac point $k=2\pi/3$, thereby
making the primary region of interest inaccessible. The procedure bypassing
this obstacle is to keep the summation in Eq. (S19) and to set its terms to
zeros whenever the expression in $\left\\{\ldots\right\\}$ is small. Equations
(S18) explain the sublattice polarization of the flat ZEM wave function seen
in Fig. 2A of the main text.

###  S2.2 hbGNR with an in-plane external field

The purposes of this subsection are twofold. Firstly, to provide an insight
into the topological stability of the hbGNR ZEM. Here by topological stability
we mean an impossibility to remove this mode from the gap by application of an
external electrostatic potential applied in-plane of the GNR across its width.
Secondly, to give an approximate analytical description of the hbGNR ZEM in
the external potential. This analytical description implies obtaining
expressions for the ZEM electron wave function as function of the external
potential at the Dirac point as the main point of interest.

In order to achieve these goals, we slightly change the applied external
potential in such a way that it ensures the same increment, $\Delta$, of the
potential between the atomic sites. This approach is an accurate model of the
finite non-dimerized Su-Schrieffer-Heeger chain subjected to the external
electrostatic field. For GNRs, however, such a potential configuration may
seem more challenging since it requires a precise engineering of the potential
applied to the ribbon within its unit cell length scale. Nevertheless, this
configuration is not a hypothetical model, since the precise potential
engineering has been demonstrated on nanomechanical systems (?). It can be
easily checked numerically that the modified potential leads to the same
results as those presented in the main text of the manuscript, where the
potential is generated by a homogeneous electrostatic field. This
insensitivity to the specific details of the applied external potential
already indicates the topological stability of the hbGNR ZEM.

The Hamiltonian of a hbGNR in the model field is

$\displaystyle H_{\mathrm{hbGNR}}$ $\displaystyle=$
$\displaystyle\left(\matrix{\Delta&t_{1}q&0&0&\ldots&0\cr
t_{1}q&2\Delta&t_{1}&0&\ldots&0\cr
0&t_{1}&3\Delta&\ddots&\ddots&\vdots\cr\vdots&\ddots&\ddots&\ddots&t_{1}q&0\cr
0&\ldots&0&t_{1}q&(N_{r}-1)\Delta&t_{1}\cr
0&\ldots&0&0&t_{1}&N_{r}\Delta}\right)\,,$ (S20)

where $N_{r}$ is the number of atoms in the ribbon unit cell that is an odd
integer and $\Delta$ is the difference of the electrostatic potentials induced
on the neighboring sites by an external field. By specialising to the case of
the Dirac point, $k=2\pi/3$, we get $q=1$, so that the Hamiltonian

$\displaystyle H_{\mathrm{hbGNR}}$ $\displaystyle=$
$\displaystyle\left(\matrix{\Delta&t_{1}&0&0&\ldots&0\cr
t_{1}&2\Delta&t_{1}&0&\ldots&0\cr
0&t_{1}&3\Delta&\ddots&\ddots&\vdots\cr\vdots&\ddots&\ddots&\ddots&t_{1}&0\cr
0&\ldots&0&t_{1}&(N_{r}-1)\Delta&t_{1}\cr
0&\ldots&0&0&t_{1}&N_{r}\Delta}\right)\,$ (S21)

has a simpler diagonal band structure. By placing the coordinate system origin
in the center of the ribbon unit cell, the main diagonal of Eq. (S21) can be
made symmetric:

$\displaystyle H_{\mathrm{hbGNR}}$ $\displaystyle=$
$\displaystyle\left(\matrix{-w\Delta&t_{1}&0&0&\ldots&0\cr
t_{1}&-(w-1)\Delta&t_{1}&0&\ldots&0\cr
0&t_{1}&-(w-2)\Delta&\ddots&\ddots&\vdots\cr\vdots&\ddots&\ddots&\ddots&t_{1}&0\cr
0&\ldots&0&t_{1}&(w-1)\Delta&t_{1}\cr 0&\ldots&0&0&t_{1}&w\Delta}\right)\,,$
(S22)

where $w=\left(N_{r}-1\right)/2$. The eigenspace problem for the Hamiltonian
given by the Eq. (S22) can be re-scaled. Since any diagonalizable matrices $M$
and $\widetilde{M}$, such that $\widetilde{M}=M/s$, where $s\neq 0$, have the
same set of eigenvectors and are characterized by eigenvalues that are related
as $\widetilde{\lambda}=\lambda/s$, the said eigenspace problem can be
substituted by the eigenspace problem for the Hamiltonian

$\displaystyle\widetilde{H}_{\mathrm{hbGNR}}$
$\displaystyle=\left(\matrix{-2w&\tau&0&0&\ldots&0\cr\tau&-2(w-1)&\tau&0&\ldots&0\cr
0&\tau&-2(w-2)&\ddots&\ddots&\vdots\cr\vdots&\ddots&\ddots&\ddots&\tau&0\cr
0&\ldots&0&\tau&2(w-1)&\tau\cr 0&\ldots&0&0&\tau&2w}\right)\,,$ (S23)

where $\tau=t_{1}/(\Delta/2)$ is a single variable of this problem. The factor
of $2$ in the $\tau$ definition and on the main diagonal of Eq. (S23) is
introduced for the purpose of further analytical treatment. We are going to
show that $\widetilde{H}_{\mathrm{hbGNR}}$ always admits a zero eigenvalue as
valid solution with a well-defined eigenvector.

The set of simultaneous equations for the eigenspace problem in question is

$\displaystyle\left(\matrix{-2w&\tau&0&0&\ldots&0\cr\tau&-2(w-1)&\tau&0&\ldots&0\cr
0&\tau&-2(w-2)&\ddots&\ddots&\vdots\cr\vdots&\ddots&\ddots&\ddots&\tau&0\cr
0&\ldots&0&\tau&2(w-1)&\tau\cr
0&\ldots&0&0&\tau&2w}\right)\left(\matrix{c_{1}\cr c_{2}\cr c_{3}\cr\vdots\cr
c_{N_{r}-1}\cr c_{N_{r}}}\right)$ $\displaystyle=$ $\displaystyle
E\left(\matrix{c_{1}\cr c_{2}\cr c_{3}\cr\vdots\cr c_{N_{r}-1}\cr
c_{N_{r}}}\right)\,,$ (S24)

or equivalently

$\displaystyle\left(-2w-E\right)c_{1}+\tau c_{2}$ $\displaystyle=$
$\displaystyle 0\,,$ (S25) $\displaystyle\tau
c_{j-1}+\left[-2\left(w-j+1\right)-E\right]c_{j}+\tau c_{j+1}$
$\displaystyle=$ $\displaystyle 0\,,\qquad j=2\ldots N_{r}-1\,,$ (S26)
$\displaystyle\tau c_{N_{r}-1}+\left(2w-E\right)c_{N_{r}}$ $\displaystyle=$
$\displaystyle 0\,.$ (S27)

If $E=0$, then the set of equation for $j=2,\ldots,N_{r}-1$ in Eqs. (S26) is
automatically satisfied using a well-known Bessel function tri-term recurrence
relation (?): $zJ_{n-1}(z)+zJ_{n+1}(z)=2nJ_{n}(z)$. Namely, we shall choose
$c_{j}=J_{w-j+1}(\tau)$, where $j=1,2,\ldots,N_{r}$. Note that as $c_{j}$’s
index $j$ goes from $1$ to $N_{r}$, the diagonal index on the main diagonal of
Eq. (S24) changes from $-w$ to $w$, while the order of the Bessel function
solution of $c_{j}$’s changes from $w$ to $-w$. This choice of $c_{j}$
transforms each of the two remaining Eqs. (S25) and (S27) to the following
equation: $-(t/s)\,J_{(N_{r}+1)/2}(\tau)=0$. The found equation imposes
restriction on the validity of the suggested eigenvector. Ideally, one would
prefer an eigenvector that converts both Eqs. (S25) and (S27) into identities
for any value of $\tau$. The suggested eigenvector, however, is an accurate
solution only for such $\tau_{i}$’s that are zeros of $J_{(N_{r}+1)/2}(z)$:
$J_{(N_{r}+1)/2}(\tau_{i})=0$. Physically, this means that the proposed
eigenvector is valid for only a discrete set of values of the external field.
Owing to this limitation, the solution shall be referred to as quasi-exact.
Since the above mentioned re-scaling of the eigenproblem by $\Delta/2$ does
not affect the zero-energy, the eigenvalue is insensitive to the presence of
an external field. However, the corresponding eigenvector shows evidences that
exposing hbGNR to the external field eliminates the sublattice polarization of
the pristine ZEM. These observations are in agreement with the results in Fig.
2D of the main text and with considerations of the continuum ${\bf k}\cdot{\bf
p}$-model in the supplementary text section S4.

###  S2.3 Summary

In this section, we have shown that hbGNR ZEM always contains a state that is
pinned to the zero-energy at the Dirac point, which is unaffected by the
presence of an external field. This insensitivity to the external field
originates from the ZEM topological stability. The sublattice polarization of
the zero-energy state is sensitive to the external field and it can be tuned
by changing the applied field.

##  S3 Haldane and Kane-Mele models in commensurate and incommensurate ribbons

In this section, we analyse the classical topological edge states and their
transport properties in commensurate and incommensurate systems. We consider
two standard models: Haldane model of a Chern insulator (?) and Kane-Mele
model of a $\mathbf{Z}_{2}$ topological insulator (?). Also, our discussion is
supplemented with the case of topological kink states that arise from the
quantum valley Hall effect and are attributed to valley Chern insulators (?).
In each case, we also pay particular attention to the bulk-boundary
correspondence principle, the phase diagrams and an edge state dispersion
transformation as a function of the topological transition driving parameter.

Before moving forward, let us first recollect the standard definition of a
topological insulator and define three main types of the topological
insulators in terms of their corresponding topological invariants. An integer
quantum Hall effect is a regime, where a material becomes a topological
insulator, that is, it is insulating in the bulk but conducting at the surface
due to the presence of the edge states. These states are explained by the
topological band theory and the bulk-boundary correspondence principle (?).
The bulk-boundary correspondence says that the interface between a topological
insulator and normal insulator produces a chiral gapless mode. The difference
between left and right moving gapless modes, is given by the difference
between the values of a topological invariant calculated for the bulk systems
to the left and right from the interface. The proper topological invariant for
the quantum Hall effect is the Chern number (?):
$\mathcal{C}_{n}=(2\pi)^{-1}\int_{\Omega}\nabla_{{\bf k}}\times A_{n}({\bf
k})\,d{\bf k}$, where $A_{n}({\bf k})=i\langle u_{n}({\bf
k})\left|\nabla_{{\bf k}}\right|u_{n}({\bf k})\rangle$ is the Berry connection
defined on the periodic parts of Bloch functions $\left|u_{n}({\bf
k})\right.\rangle$ and the integration is carried out over the whole Brillouin
zone $\Omega$. For the quantum spin Hall effect, the Chern invariant is zero.
The proper topological invariant is the 2D $\mathbf{Z}_{2}$ invariant (?, ?):
$\nu=(2\pi)^{-1}\left(\oint_{\partial\Omega_{1/2}}A({\bf k})\,d{\bf
k}-\int_{\Omega_{1/2}}\nabla_{{\bf k}}\times A({\bf k})\,d{\bf
k}\right)\mathrm{mod}\;2$, where $\Omega_{1/2}$ is the half of the Brillouin
zone connected to another half by the time reversal symmetry,
$\partial\Omega_{1/2}$ is the edge of $\Omega_{1/2}$ and $A({\bf
k})=\sum_{n}A_{n}({\bf k})$, where the sum runs over all occupied states. For
the quantum valley Hall effect the proper topological invariant is a
topological charge (?) or the valley Chern number (?): $\mathcal{C}_{n,{\bf
K}({\bf K}^{\prime})}=(2\pi)^{-1}\int_{\Omega_{{\bf K}({\bf
K}^{\prime})}}\nabla_{{\bf k}}\times A_{n}({\bf k})\,d{\bf k}$, where
$\Omega_{{\bf K}({\bf K}^{\prime})}$ stand for the K and K′ valley regions.
Thus, all gapless edge modes used as channels for dissipationless transport
can be associated with a distinct topological invariant. For the clarity of
the following discussion, we remind that term mode is used here as a synonym
of a band rather than that of a state. The mode has a dispersion. We refer to
the mode as the zero-energy one if it contains a zero-energy state which is a
state at the Fermi level. The flavor of the topological insulator is
determined by its topological invariant.

###  S3.1 Chern insulator

Let us now switch on Haldane’s $t_{2}$ term adiabatically and observe the
changes in zigzag (commensurate) and half-bearded (incommensurate) ribbons.
From the phase diagram of Haldane’s model presented in Fig. 2 of Ref. (?), we
expect the system to maintain two topological edge modes, since we use zero
staggered potential. The bulk system is topologically non-trivial and
characterized by Chern number $\mathcal{C}=1$. Therefore, when two interfaces
with vacuum are formed, each of two edges must carry a “one way” moving edge-
localized mode. Both modes cross the bulk gap and pass through the zero
energy. In agreement with the bulk-boundary correspondence, the difference
between the number of left and right moving modes for each edge is equal to
the difference of the Chern numbers across the interface:
$N_{L}-N_{R}=\Delta\mathcal{C}=\mathcal{C}-\mathcal{C}_{\mathrm{vac}}=1$,
where $\mathcal{C}_{\mathrm{vac}}=0$ is the Chern number of the vacuum.

####  S3.1.1 Commensurate ribbon

Figure S2 presents the energy bands and electron densities for the chosen
$k$-points for small, $t_{2}=0.01t_{1}$, and large, $t_{2}=0.1t_{1}$, values
of Haldane’s parameter for zigzag and half-bearded GNRs. In a zigzag GNR,
introducing a small next-to-nearest neighbor hopping parameter, $t_{2}$,
results in a splitting of the partially flat energy bands at the Dirac points
K and K′, while band degeneracy is preserved in the Kramer’s time reversal
(TR) invariant point $k=\pi$. In this latter point, the conduction band from K
(K′) valley connects to the valence band from K′(K) valley, forming a pair of
ZEMs filling the bulk energy gap. As seen from the distribution of the
electron densities in fig. S2A, the wave functions of these two ZEMs are
sharply localized at the opposite edges of the ribbon. This is exactly what is
expected from the bulk-boundary correspondence principle. We should also
notice that in $k$-space both ZEMs extend throughout the whole Brilloiun zone.
Moreover, the TR symmetry links $k=\pi-q$ state in each mode to $k=\pi+q$
states in another mode, where $q$ is a relative momentum measured form
$k=\pi$. In the vicinity of the TR invariant point $k=\pi$, the two modes have
linear dispersion with opposite group velocities. Since the Haldane’s term
breaks the TR symmetry of the Hamiltonian, the scattering amplitude for any
pair of TR symmetry related states is not vanishing; see Chapter 4 of Ref. (?)
and our further discussion of $\mathbf{Z}_{2}$ topological insulators. Thus,
an elastic backscattering is theoretically possible. However, the wave
functions of the TR symmetry related states, circled 2-5 and circled 3-4, are
spatially distant. This means that wider ribbons have wave functions of pair
states, which do not overlap. Hence, the elastic backscattering is suppressed
in such systems and dissipationless ballistic transport is realized. However,
if the ribbon is narrow enough, then there is no physical reason forbidding
the backscattering. Hence, dissipative transport shall be observed for narrow
ribbons. For the commensurate ribbons, the wave functions spatial separation
is a crucial factor that enables dissipantionless transport. This picture does
not change for larger values of $t_{2}$ as shown in fig. S2C.

####  S3.1.2 Incommensurate ribbon

In a half-bearded GNR, the same Haldane’s term makes the fully flat ZEM
dispersive, as shown in fig. S2B. In this case, the initial ZEM represented by
a straight horizontal line turns into a broken line due to two factors: (i)
the shifting of energy states at K and K′ point towards conduction and valence
bands, respectively, and (ii) the pinning of the energy states at the TR
invariant points $k=0$ and $\pi$. In contrast to the commensurate ZGNR, the
broken line in hbGNR does not link valence and conduction bands for small
values of $t_{2}$. We arrive at similar conclusions when using larger values
of $t_{2}$, as shown in fig. S2D. The inset of fig. S2D shows that for large
$t_{2}$ such as $t_{2}=0.1t_{1}$, the ZEM is still detached from the
conduction and valence bands by a small avoided crossing. This means that this
system can be turned into a true insulator as a result of doping by
positioning the Fermi level within this avoided crossing. We have checked that
the avoided crossing persists even for $t_{2}$ as large as $t_{1}$, though for
such extra large values the bands overlap at different $k$’s so that the
positioning of the Fermi level within avoided crossing gap is not possible. In
all given cases, the ZEM is subdivided by K and K′ points into two submodes,
where wave functions are localized at the opposite edges of the ribbon. The
bulk-boundary correspondence predicts two edge modes in a system with two
edges, therefore we have to associate it with the two edge localized submodes.
Another important observation is that there is a threshold value for $t_{2}$
to achieve an overlap between the ZEM and bulk bands. This implies that below
the threshold $t_{2}$ the standard phase diagram of Haldane’s model must have
an electrostatic doping level as one of its dimensions, when dealing with
incommensurate systems. Although such doping cannot change the Chern number
and, therefore, the topological class of the material, it significantly
affects the way a finite piece of topological material manifests itself in the
transport. For instance, when Haldane’s model hexagonal lattice is doped by
applying an out-of-plane electric field, then trimming this lattice into an
hbGNR can result in a system that is not conductive neither in the bulk nor at
the edges. In this case, despite being originating from a topological
insulator the system is physically equivalent to a trivial insulator.
Concurrently, choosing a different level of doping involves the edges of the
system into electronic transport again, which from the physical point of view
corresponds to a standard well-known presentation of a topological insulator.

In the vicinity of the two TR invariant points, $k=0$ and $\pi$, the
dispersion of the ZEM is linear. Contrary to the commensurate ZGNR, the states
of the ZEM in hbGNR do not have TR partners. The ZEM states at $k=\pi-q$ lack
the same energy partners located at $\pi+q$, as seen from the points circled
2-5 and circled 3-4 in Figs. S2B and S2D. Hence, the elastic backscattering is
by default forbidden. In this case, the edge localization of the ZEM states
wave functions and the width of the ribbon do not play any role in achieving
the ballistic transport. Therefore, we infer that dissipationless transport in
incommensurate ribbons can be realized via bulk modes too.

###  S3.2 Z2 topological insulator

Now that we have revealed the main features of Chern insulators with respect
to commensurate and incommensurate systems, we proceed to $\mathbf{Z}_{2}$
topological insulators. We adiabatically switch on the Kane-Mele spin-orbit
term via the $t_{2}$ parameter. To highlight the difference between spin up
and spin down states, we introduce a finite staggered on-site potential
$\Delta=2t_{2}$ while we set the Rashba spin-orbit term to zero. Figure S3
presents the energy bands and electron spin polarization densities for the
chosen $k$-points in zigzag and half-bearded GNRs. According to the phase
diagram of the Kane-Mele model presented in Fig. 1 of Ref. (?), the bulk
system is a 2D $\mathbf{Z}_{2}$ topological insulator characterized by $\nu=1$
as the spin-orbit term is switched on. For the $\mathbf{Z}_{2}$ topological
insulator, the bulk-boundary correspondence says that
$\Delta\nu=(N_{K}\mathrm{mod}\;2)$, where $\Delta\nu$ is the change of
$\mathbf{Z}_{2}$ index across an interface and $N_{K}$ is the number of the
Kramers pairs at the Fermi level for the edge states.

Since by default $\mathbf{Z}_{2}$ topological insulators deal with spinful
particles, we briefly recall what Kramers pairs are and how they participate
in elastic scattering. For half-integer spin particles, the TR symmetry
operator $\mathcal{T}$ acting on a Bloch state $\psi$ flips not only the
particle momentum but also its spin. Hence, Kramers pairs are formed by
$\psi_{k,\uparrow}\leftrightarrow\psi_{-k,\downarrow}$ and
$\psi_{k,\downarrow}\leftrightarrow\psi_{-k,\uparrow}$. Note that for one-
dimensional periodic structures the momentum flip around $k=0$ is the same as
the momentum flip around $k=\pi$; for the simplicity of the following
discussion we will always refer to $k\leftrightarrow-k$ flip even though fig.
S3 shows Brillouin zones that are more adapted for $\pi+q\leftrightarrow\pi-q$
notation, where $q$ is the momentum shift from the $k=\pi$ point. For
$\mathcal{T}$-symmetric Hamiltonians,
$[\mathcal{T},H(k)]=0\Rightarrow\mathcal{T}H(k)\mathcal{T}^{-1}=H(-k)$. Thus,
the scattering matrix element between the TR invariant Kramers pairs, i.e.
when $k=-k$, is vanishing:

$\displaystyle\langle\psi\left|H(k)\mathcal{T}\right|\psi\rangle$
$\displaystyle=$
$\displaystyle\langle\psi\left|\mathcal{T}H(k)\mathcal{T}^{-1}\mathcal{T}\right|\psi\rangle=\langle\psi\left|\mathcal{T}H(k)\right|\psi\rangle=$
(S28) $\displaystyle=$
$\displaystyle\langle\psi\left|H^{\dagger}(k)\mathcal{T}^{\mathrm{T}}\right|\psi\rangle\stackrel{{\scriptstyle\left\\{T^{\mathrm{T}}=-T\right\\}}}{{=}}-\langle\psi\left|H^{\dagger}(k)\mathcal{T}\right|\psi\rangle\stackrel{{\scriptstyle\left\\{H^{\dagger}=H\right\\}}}{{=}}$
$\displaystyle\stackrel{{\scriptstyle\left\\{H^{\dagger}=H\right\\}}}{{=}}$
$\displaystyle-\langle\psi\left|H(k)\mathcal{T}\right|\psi\rangle\,,$

where we have used the fact that $H$ is self-adjoint (Hermitian) operator and
$\mathcal{T}^{\mathrm{T}}=-\mathcal{T}$. The latter relation follows from the
definition of TR operator

$\mathcal{T}=UK\,,$ (S29)

where $U$ is a unitary matrix, i.e. $UU^{\dagger}=U^{\dagger}U=I$, and $K$ is
the complex conjugate. Since $U$ is antisymmetric, $-U=U^{\mathrm{T}}$, the
transposition of the TR operator yields

$\mathcal{T}^{\mathrm{T}}=(UK)^{\mathrm{T}}=U^{\mathrm{T}}K=-UK=-\mathcal{T}\,.$
(S30)

The antisymmetric character of $U$ comes from the condition imposed on
$\mathcal{T}$ by the half-integer spin:
$\mathcal{T}^{2}=-1=(UK)(UK)=UU^{\ast}KK=UU^{\ast}$. Accounting for the
unitarity of $U$, we then arrive at $-U^{\ast}=U^{-1}=U^{\dagger}$. Taking the
complex conjugate from the left- and right-hand sides of the latter equality
leads to the antisymmetric condition $-U=U^{\mathrm{T}}$ used in Eq. (S30).

The fundamental property of a particle with a half-integer spin is that
$\mathcal{T}^{2}=-1$. It originates from the fact that $\mathcal{T}$ flips the
spin of the particle thereby rotating it around some axis by angle $\pi$. This
rotation for a particle with the spin equal to $1/2$ can be realized by
choosing $U=i\sigma_{y}$, where

$\sigma_{y}=\left(\matrix{0&-i\cr i&0}\right)$

is Pauli matrix. This yields
$\mathcal{T}^{2}=(i\sigma_{y}K)(i\sigma_{y}K)=i\sigma_{y}(-i)\sigma^{\ast}_{y}KK=\sigma_{y}\sigma^{\ast}_{y}=-\sigma_{y}^{2}=-1$.
The proven relation in Eq. (S28) shows that the elastic backscattring is
forbidden between all Kramers pairs of a given pair of bands so that only a
spin-conserving scattering is allowed:
$\psi_{k,\uparrow}\leftrightarrow\psi_{-k,\uparrow}$ or
$\psi_{k,\downarrow}\leftrightarrow\psi_{-k,\downarrow}$. This is an important
corollary that we will use for analysis of the transport properties of
$\mathbf{Z}_{2}$ topological insulators with respect to the commensurate and
incommensurate GNRs.

####  S3.2.1 Commensurate ribbon

For a ZGNR, the effect of the spin-orbit term is similar to that in Haldane’s
model. For both spin polarized pairs of partially flat bands, we observe
splittings at the Dirac points K and K′, while their degeneracy is preserved
at the TR invariant point $k=\pi$. When the staggered potential is zero
($\Delta=0$), a fourfold degeneracy is found at $k=\pi$ point, since ZEM $1$,
ZEM $2$, ZEM $3$, and ZEM $4$ pass through zero energy in this point. As seen
from fig. S3A, a finite $\Delta=2t_{2}$ shifts the crossing between ZEM $1$
and ZEM $4$ from $k=\pi$ towards the K valley. The crossing between ZEM $2$
and ZEM $3$ is shifted towards K′ valley. As expected from the bulk boundary
correspondence, each edge hosts two modes with opposite spin polarizations.
Having opposite spin polarizations, ZEM $1$ and ZEM $3$ form a Kramers pair.
As seen from spin polarization densities in points circled 3 and circled 4 in
fig. S3A, the ZEM $1$ and ZEM $3$ are both localized at the same left edge.
They cross the Fermi level $E_{f}=0$ at the $k$ points related by the TR
symmetry. Thus, the number of Kramers pairs for the left edge localized ZEMs
is $N_{K}=1$. Then $(N_{K}\mathrm{mod}\;2)=1$, which is exactly equal to the
change of the $\mathbf{Z}_{2}$ index across the left ribbon edge:
$\Delta\nu=\nu-\nu_{\mathrm{vac}}=1-0=1$, where $\nu$ and $\nu_{\mathrm{vac}}$
are the $\mathbf{Z}_{2}$ index of the bulk Kane-Mele system and the vacuum,
respectively. Similar reasoning applies to the right edge of the ribbon, where
ZEM $2$ and ZEM $4$ are localized as seen from spin polarization densities in
points circled 2 and circled 5 in fig. S3A. Thus, four ZEMs are expected for
the system with two edges, which is indeed observed.

The four topological ZEMs have peculiar transport properties due to the
preserved TR symmetry. Equation (S28) states that the backscattering between
Kramers pairs is strictly forbidden as long as TR symmetry is respected by the
Hamiltonian. Thus, no scattering must be present between ZEM $1$ and ZEM $3$
as well as between ZEM $2$ and ZEM $4$. Equation (S28), however, does not
forbid spin-conserving elastic scattering between ZEM $1$ and ZEM $4$ as well
as ZEM $2$ and ZEM $3$. For $\Delta=0$, the same energy states in these pairs
of modes are symmetrically positioned at $k$ and $-k$. Thus, the momentum and
energy conservation laws in the elastic backscattering between these pairs of
modes is respected. The probability of this process, however, is given by the
overlap of the wave functions of the corresponding states. It is seen from
points circled 4 and circled 5 in fig. S3A that the wave functions in these
spin-conserving pairs are localized at the opposite edges of the ribbon, so
that the overlap is negligible. Hence, the elastic backscattering does not
occur. In summary, for commensurate ZGNR, the TR symmetry and wave function
spatial distancing are important in achieving a reliable dissipationless
ballistic transport.

####  S3.2.2 Incommensurate ribbon

For the case of a hbGNR, the Kane-Mele spin-orbit term shifts the spin up and
down states at K (K′) points in the opposite directions while preserving their
degeneracy at the TR invariant Kramers points $k=0$ and $\pi$. This spin state
shifting breaks both spin down ZEM $1$ and spin up ZEM $2$ into two regions:
region I located between K and K′ points, and region II that is a complement
of the region I to the full Brillouin zone. In each region, the wave functions
of the spin up and down ZEM states are localized at the same ribbon edge as
shown by the spin polarization density in points circled 1, circled 3 and
circled 4 of fig. S3B. At the boundary of regions I and II the wave functions
are extended throughout the whole width of the ribbon. This picture differs
from the hbGNR case in the Haldane’s model, where the wave function preserves
the localized character while being distributed between the opposite ribbon
edges in equal fractions. In the given hbGNR case, the ZEM localized at the
edges of the ribbon do not connect the valence and conduction bands. As can be
seen in the inset of fig. S3B, there is an avoided crossing between ZEM $1$
and the valence bands. Despite considerable spin-orbit interaction that drives
the bulk system to the non-trivial $\mathbf{Z}_{2}$ topological regime with
$\nu=1$, the hbGNR being cut from the bulk system can be converted to a true
insulator by positioning the Fermi level within the avoided crossing.

Similar to the considered case in Haldane’s model, the bulk-boundary
correspondence can be formally satisfied if we subdivide both ZEM $1$ and ZEM
$2$ into submodes residing within the above mentioned regions I and II. Within
the region I, we have a Kramers pair formed by points circled 3 and circled 4
in fig. S3B. These points correspond to states with different spin
polarization, and they are positioned at TR symmetry related momenta $k$ and
$-k$. The wave functions of the corresponding states are localized at the
zigzag (left) edge of the ribbon. Hence, $N_{K}=1$, therefore
$(N_{K}\mathrm{mod}\;2)=1=\Delta\nu$. In the same way, we could consider a
pair of states in region II, thereby obtaining two additional submodes and a
total of four $\mathbf{Z}_{2}$ topological edge submodes for the system with
two edges. Also, we have checked that ZEMs do not touch any of the valence and
conduction bands even for as large spin-orbit coupling as $t_{2}=t_{1}$. The
dispersion of the ZEMs at the TR invariant points $k=0$ and $\pi$ is almost
linear. In contrast to the Haldane’s model, hbGNR ZEMs in Kane-Mele model have
TR partners thereby forming Kramers pairs. However, as we have shown in Eq.
(S28), the elastic backscattering is forbidden between the Kramers pairs as
long as the TR symmetry is respected. On the other hand, the considered ZEMs
lack spin-conserving partners for which the backscattering would be allowed.
Again, we come to the conclusion that the edge localization is not important
for achieving ballistic transport in this incommensurate GNR. Thus, the
dissipationless transport could be realized by bulk modes.

###  S3.3 Quantum valley Hall effect

Similar considerations can be given to the case of the quantum valley Hall
(QVH) regime and topological kink states in the commensurate zigzag and
incommesurate half-bearded GNRs. In analysing these systems, one should take
into account that TR symmetry operator switches valleys, and the suppression
of elastic backscattering comes from valley distancing in $k$-space rather
than from the vanishing of the scattering matrix elements between the Kramers
pairs. The case of commensurate GNR is presented in main text Fig. 2C, while
the case of incommensurate ribbon is shown in fig. S4. Two ingredients are
essential in achieving QVH regime in monolayer ribbons: (1) the breaking of
inversion symmetry by introducing a staggered sublattice potential leading to
the energy gap opening, and (2) creating two kinks in the staggered potential
by the side gate voltages. The former induces a non-zero Berry curvature
within each of the two valleys of graphene lattice, making the valley Chern
number a well-defined topological invariant. The side gate voltage pushes the
partially flat bands within the bulk gap, thereby creating channels for
ballistic transport. Notice that for the hbGNR there is only one band that is
localized on a single edge. Thus, we again arrive at the same conclusion above
that dissipationless transport in incommensurate systems can be realized via
bulk bands. The in-plane electric field applied to the hbGNR in the main text
is a single ingredient that achieves the same ballistic regime. Moreover, when
combined with an electrostatic doping, the in-plane electric field allows
tuning between the bulk and edge mode character of the dissipationless
transport.

###  S3.4 Summary

To conclude, we have seen that incommesurate systems possess peculiarities in
both electronic and transport properties. The behaviour of incommensurate
systems considerably deviates from that of commensurate systems for all
standard flavours of the topological insulators: Haldane’s, Kane-Mele’s and
Liu’s models representing standard Chern, $\mathbf{Z}_{2}$ and valley Chern
topological insulators. Not only can incommensurate systems provide
dissipationless tranport via edge mode, but also dissipationless transport via
bulk modes. In all the considered above cases, incommensurate systems can be
driven into true insulators by electrostatic doping, making the standard phase
diagrams less practical without accounting for doping dimension. For all
above-considered systems, the dispersion of ZEMs at the intrinsic Fermi level,
$E_{\mathrm{F}}=0$, is quasi-linear for all values of the primary model
parameter $t_{2}$ [Haldane’s and Kane-Mele’s models] or $U$ [Liu’s model of
the QVH effect in monolayer ribbons (?)]. Moreover, no cubic-type dispersion
is observed.

##  S4 Low-energy theory of hbGNR in continuum model

In condensed matter systems, the Dirac equation arises from the low-energy
limit of the Schrödinger equation around special degeneracy points of the
energy bands. These special points represent topological objects acting as the
Berry curvature sources in the $k$-space and resembling Dirac monopoles. The
graphene K and K′ points are classical examples of Berry curvature monopoles
referred to as Dirac points. The low-energy description of graphene has been
known for the long time (?). When supplemented by hard-wall boundary
conditions, it has also been used to describe the electronic properties of
graphene nanoribbons (?, ?) and finite length metallic single-wall carbon
nanotubes (?). In particular, the two best known GNR geometries with zigzag
and armchair edges have been considered and the derivation of partially flat
zero-energy bands featuring edge localized states has been reported (?).
However, a similar systematic description is currently missing for the hbGNRs
fully flat zero-energy band. Therefore, here we develop a relevant analytic
treatment of hbGNRs zero-energy bands as well as deepen understanding of their
fundamental behavior in an external in-plane electric field and supporting
claims of the main text with respect to the intimate relation of ZEM to the
quantum field theory, topology and chaos and non-linear dynamics.

###  S4.1 Bulk energy bands

Following the case of zigzag GNR in Ref. (?), we start from the Dirac equation
and amend the boundary condition in order to get a secular equation for the
hbGNR. For a GNR oriented along $y$-axis, the electron momentum along this
axis is a good quantum number. Hence, in K valley the Dirac equation for a
two-component spinor reads

$\displaystyle-i\partial_{x}\phi_{B}+i\kappa_{y}\phi_{B}$ $\displaystyle=$
$\displaystyle\varepsilon\phi_{A}\,,$ (S31)
$\displaystyle-i\partial_{x}\phi_{A}-i\kappa_{y}\phi_{A}$ $\displaystyle=$
$\displaystyle\varepsilon\phi_{B}\,,$ (S32)

where $\kappa_{y}$ is the electron momentum along the ribbon measured with
respect to the Dirac point. Expressing $\phi_{A}$ from Eq. (S31) as
$\phi_{A}=-\left(i/\varepsilon\right)\partial_{x}\phi_{B}+\left(i\kappa_{y}/\varepsilon\right)\phi_{B}$
and putting it into Eq. (S32), we arrive at a single second order equation:

$-\partial_{xx}\phi_{B}+\kappa_{y}^{2}\phi_{B}=\varepsilon^{2}\phi_{B}\,.$
(S33)

Equation (S33) has a general solution $\phi_{B}=C_{1}e^{zx}+C_{2}e^{-zx}$,
with $z=\sqrt{\kappa_{y}^{2}-\varepsilon^{2}}$ and $C_{1,2}$ as constants
determined from the boundary and normalization conditions. For a zigzag GNR,
the outermost atomic sites at the opposite edges belong to different
sublattices. The wave function must vanish on the next auxiliary atomic sites
(?), therefore the boundary condition for the ribbon of width $L$ is
$\phi_{B}(0)=\phi_{A}(L)=0$. In contrast, the outermost edge atoms of a hbGNR
belong to the same sublattice, say sublattice $A$. Hence, the wave function
must vanish on the next auxiliary atomic sites belonging to sublattice $B$.
Therefore, we set for the hbGNR

$\phi_{B}(0)=\phi_{B}(L)=0\,.$ (S34)

This condition yields the following set of equations:

$\displaystyle C_{1}e^{z\cdot 0}+C_{2}e^{-z\cdot 0}$ $\displaystyle=$
$\displaystyle 0\,,$ $\displaystyle C_{1}$ $\displaystyle=$ $\displaystyle-
C_{2}\,,$ (S35) $\displaystyle C_{1}e^{zL}+C_{2}e^{-zL}$ $\displaystyle=$
$\displaystyle 0\,,$ $\displaystyle 2C_{1}\sinh\left(zL\right)$
$\displaystyle=$ $\displaystyle 0\,,$ (S36)

Equation (S35) links the two initially independent constants $C_{1}$ and
$C_{2}$, leaving only one degree of freedom represented by either constant
$C_{1}$ or $C_{2}$. This remaining degree of freedom is fixed by normalization
condition. Equation (S36) leads to the secular equation
$\sinh\left(zL\right)=0$. This secular equation has only one solution $z=0$,
whence it follows that $\varepsilon=\pm\left|\kappa_{y}\right|$. It is obvious
that the chosen method of solving equation does not provide the correct
description of the ZEM. The method leads to two principle problems: (i) there
are two bands marked by ‘$+$’ and ‘$-$’ that normally attributed to the
conduction and valence bands, (ii) the obtained dispersion is linear. In
contrast to the derived result, the ZEM in question is a sole mode that
pertains to neither conduction nor valence band. Since the external electric
field is absent, the ZEM must also be independent of $\kappa_{y}$ rather than
a mode which exhibits linear dispersion. Thus, a different approach is needed
even for the most basic description of the mode not mentioning its spectacular
transformations in the external field from fully flat to the cubic dispersion.
The failure in finding a correct description of the ZEM, however, does not
mean that the initial method is not valid for the bulk bands of hbGNRs. By
substituting $z=i\kappa_{x}$, one can transform the above-mentioned secular
equation to $\sin\left(\kappa_{x}L\right)=0$ [cf. with Eq. (S9)], which leads
to quantization of the transverse electron momentum $\kappa_{x}$:
$\kappa_{x}=\pi n/L$, with $n$ being integer. Using this momentum
quantization, the valence and conduction energy bands are described by
$\varepsilon=\pm\sqrt{\kappa_{y}^{2}+\left(\pi n/L\right)^{2}}$, where $n=\pm
1,\pm 2,\ldots$. The option $n=0$ shall be excluded from the energy band
consideration for it leads to the already discussed
$\varepsilon=\pm\left|\kappa_{y}\right|$.

While the failure to achieve a valid description of the ZEM within the given
formalism can make us skeptical on the presented bulk states, we shall notice
the following. One more argument supporting the validity of the bulk bands
description of hbGNRs can be given based on their isospectral feature with
respect to the energy bands of the tubes. Indeed, imposing the periodic
boundary conditions on Eqs. (S31) and (S32), $\phi_{B}(0)=\phi_{B}(L)$ and
$\partial_{x}\phi_{B}(0)=\partial_{x}\phi_{B}(L)$, one can derive the
following secular equation for the circumferencial momentum of electrons,
$\kappa_{x}$: $e^{i\kappa_{x}C}=1$, where $C$ is the tube circumference. The
secular equation leads to $\kappa_{x}=2\pi n/C$, with $n$ being integer, which
results in the tube energy bands that are identical to hbGNR bulk energy bands
once $L=C/2$. The bulk bands derived from the low-energy description of hbGNR
are, indeed, isospectral with those of the corresponding armchair tube. Also,
the same structural relation takes place as in the tight-binding model between
the width and circumference of the isospectral ribbon and tube (see fig. S1B).

To find normalized spinor wave functions
$\Psi=\left(\phi_{A},\phi_{B}\right)^{\mathrm{T}}$ for the bulk bands of the
hbGNR, it is convenient to consider the ribbon to be positioned symmetrically
with respect to $x=0$. For this case, the boundary conditions are
$\phi_{B}(-L/2)=\phi_{B}(L/2)=0$ [cf. with Eqs. (S34)]:

$\displaystyle C_{1}e^{-zL/2}+C_{2}e^{zL/2}$ $\displaystyle=$ $\displaystyle
0\,,$ $\displaystyle C_{1}e^{zL/2}+C_{2}e^{-zL/2}$ $\displaystyle=$
$\displaystyle 0\,.$ (S37)

The set of Eqs. (S37) is a homogeneous system of linear equations with respect
to $C_{1}$ and $C_{2}$. Non-trivial solutions of such system can be found only
if the determinant of its matrix is zero. This condition defines the
quantization of $\kappa_{x}$, which holds the same form as before:
$\kappa_{x}=\pi n/L$, with $n$ being positive integer. From the second
equation in Eqs. (S37), we get $C_{2}=-C_{1}e^{i\kappa_{x}L}$ so that
$\phi_{B}=-2iC_{1}e^{i\kappa_{x}L/2}\sin\left[\kappa_{x}\left(L/2-x\right)\right]$
and $\phi_{A}=\pm
2C_{1}e^{i\kappa_{x}L/2}\cos\left[\kappa_{x}\left(L/2-x\right)+\varphi\right]$,
where ‘$+$’ and ‘$-$’ in $\phi_{A}$ are used for the conduction and valence
band, respectively, and
$\varphi=\arcsin\left(\kappa_{y}/\sqrt{\kappa_{y}^{2}+\kappa_{x}^{2}}\right)$.
Hence, the normalized spinor is

$\Psi=\left(\matrix{\phi_{A}\cr\phi_{B}}\right)=\frac{1}{\sqrt{L}}\left(\matrix{\pm\cos\left[\kappa_{x}\left(L/2-x\right)+\varphi\right]\cr-i\sin\left[\kappa_{x}\left(L/2-x\right)\right]}\right)\,,$
(S38)

where we have disregarded the global phase $e^{i\kappa_{x}L/2}$ obtained in
the derivation.

###  S4.2 Zero-energy mode

Since the seminal paper by Jackiw and Rebbi (?), the zero-energy modes (ZEMs)
of the Dirac fermions have been a subject of the intensive research in the
quantum field theory due to their relation to exotic concepts such as charge
fractionalization (?, ?) and Majorana fermions (?, ?) and superconducting
strings (?). Taking this into account, it is quite natural to make the
conjecture that the ZEM of hbGNR can be described in a similar approach (see
also Ref. (?)).

The Dirac equation itself is invariant with respect to the TR operation. The
hbGNR ZEM in the tight-binding model is also TR symmetric. Namely, the two
valleys around K and K′ points are related by the $\mathcal{T}$-operator. In a
single valley description of ZEM, the TR symmetry is locally broken while
globally it is preserved. Breaking of the TR symmetry implies the mode
chirality. The chiral mode can be described within the Dirac equation
formalism when the ZEM is an eigenspace mode of the particle conjugation or
chiral symmetry operator represented by the Pauli matrix $\sigma_{z}$ (see
Refs. (?, ?)). In this case, the two Eqs. (S31) and (S32) collapse into one.
Following this prescription, we can rewrite Eqs. (S31) and (S32) in a matrix
form:

$\left(-i\sigma_{x}\partial_{x}-\sigma_{y}\kappa_{y}\right)\Psi=\varepsilon\Psi\,,$
(S39)

where $\Psi=\left(\phi_{A},\phi_{B}\right)^{\mathrm{T}}$. By setting
$\varepsilon=0$ and multiplying the left and right hand sides by $\sigma_{x}$,
we arrive at

$\partial_{x}\Psi=-\sigma_{z}\kappa_{y}\Psi\,.$ (S40)

Equation (S40) has a standard form of the Jackiw-Rebbi ZEM equation:
$\partial_{x}\Psi=m(x)\sigma_{z}\Psi$, where $m(x)$ is the spatial mass-kink
profile (?, ?, ?). However, in contrast to the standard Jackiw-Rebbi equation,
see Refs. (?, ?, ?) and references therein, in Eq. (S40) the gap generating
mass does not depend on the spatial coordinates but rather on some external
parameter. In other words, we deal with the case
$\partial_{x}\Psi=m(\eta)\sigma_{z}\Psi$, where $m(\eta)$ is the mass term
profile in a parameter space. For solids, the reciprocal space is a natural
parameter space. When $m\equiv k$ in a 1D system, the reciprocal space origin
provides a sign change in the mass profile . Thus, in Eq. (S40) we deal with a
$k$-dependent mass. This mass is somewhat similar to the Haldane mass (?), but
it does not need inversion symmetry to have a sign change.

Equation (S40) reduces to a scalar equation if spinor $\Psi$ is an eigenvector
of $\sigma_{z}$: $\sigma_{z}V_{\lambda}=\lambda V_{\lambda}$, with
$\lambda_{+}=1,V_{+}=\left(1,0\right)^{\mathrm{T}}$ and
$\lambda_{-}=-1,V_{-}=\left(0,1\right)^{\mathrm{T}}$. Hence,
$\Psi_{\lambda}=V_{\lambda}\,\chi(x)$, where $\chi(x)$ is the scalar function
to be found from

$\partial_{x}\chi(x)=-\lambda\kappa_{y}\chi(x)\,.$ (S41)

The solution of this first order equation is
$\chi_{\pm}(x)=A\exp\left(-\lambda_{\pm}\kappa_{y}x\right)$, carrying only one
constant that can be fixed by the normalization condition. Thus, the spinor is

$\Psi_{\pm}(x)=V_{\pm}A\exp\left(-\lambda_{\pm}\kappa_{y}x\right)\,.$ (S42)

The spinor is exponentially decaying or growing for $\kappa_{y}\neq 0$ and it
is constant for $\kappa_{y}=0$, which aligns well with the behaviour of ZEM
known from the tight-binding model. Moreover, eigenvectors of $\sigma_{z}$
exhibit the two-component polarization which can be attributed to the
sublattice polarization of the ZEM. There are, however, substantial
difficulties in interpretation if we construct the ZEM as the original Jackiw-
Rebbi solution in Ref. (?). As follows from the tight-binding results, the
sublattice polarization of ZEM is preserved even when the localization of the
wave function changes form zigzag to bearded edge of the ribbon, i.e. when
$\kappa_{y}>0$ changes to $\kappa_{y}<0$ in K valley. However, to produce a
normalizible Jackiw-Rebbi solution, $\lambda$ should be chosen for
$\kappa_{y}>0$ and $\kappa_{y}<0$ so that to ensure that spinors
$\Psi_{\pm}(x\rightarrow\pm\infty)\rightarrow\mathbf{0}$:

$\displaystyle\kappa_{y}>0:\;\Psi_{+}$ $\displaystyle=$ $\displaystyle
AV_{+}\exp\left(-\kappa_{y}x\right)\,,$ $\displaystyle\kappa_{y}<0:\;\Psi_{-}$
$\displaystyle=$ $\displaystyle AV_{-}\exp\left(\kappa_{y}x\right)\,.$ (S43)

As one can see from Eqs. (S43), neither the sublattice polarization nor the
wave function localization can be properly described by thus constructed
solution. According to Eq. (S43) the sublattice polarization changes as
$\kappa_{y}$ changes sign. Also, the wave function is always localized at the
same interface $x=0$ and does not shift upon $\kappa_{y}$ variation.
Furthermore, the boundary condition
$\Psi_{\pm}(x\rightarrow\pm\infty)\rightarrow\mathbf{0}$ is incompatible with
that given by Eq. (S34).

The above reasoning is true only within the original Jackiw-Rebbi framework
from Ref. (?) that disregards the exponential solutions on the basis of their
non-normalizability. But this is true only for an infinite system. The
exponential solutions are normalizable on a finite region of space. Therefore,
if one accounts for the finite transverse size of the hbGNR, then the
counterpart of the Jackiw-Rebbi solution is a valid physical entity found in
real-world systems. For fully polarized modes provided by Eq. (S42) the
boundary conditions in Eq. (S34) are automatically satisfied. As shown in fig.
S5A, all outer atoms of hbGNR (blue) belong to the same sublattice $A$, while
the next lattice nodes (red) all belong to another sublattice $B$. The
hardwall boundary condition demands the wave function being zero on all red
sites that belong to $B$ sublattice in fig. S5A. This requirement is
automatically fulfilled once the spinor $\Psi$ is proportional to the
eigenvector of $\sigma_{z}$ with zero $B$ component:
$V_{+}=\left(1,0\right)^{\mathrm{T}}$. In this way, we correctly choose the
solution between the two spinors provided by Eq. (S42). It becomes obvious now
that the arbitrary constant from the general solution of the first order Eq.
(S41) corresponds only to the normalization degree of freedom. The proper
solution is

$\Psi_{+}(x)=\left(\matrix{1\cr
0}\right)C_{1}\exp\left(-\kappa_{y}x\right)\,.$ (S44)

It is convenient to normalize the solution given by Eq. (S44) assuming that
the ribbon is placed symmetrically with respect to $x=0$. Owing to the limited
width of the nanoribbon, the integration has a finite result. Hence, the
normalization constant is readily found:

$\displaystyle\int_{-L/2}^{L/2}\Psi_{+}^{\dagger}(x)\Psi_{+}(x)d\,x$
$\displaystyle=$ $\displaystyle 1\,,$
$\displaystyle\left|C_{1}\right|^{2}\int_{-L/2}^{L/2}\exp\left(-2\kappa_{y}x\right)d\,x$
$\displaystyle=$ $\displaystyle 1\,,$ $\displaystyle\left|C_{1}\right|$
$\displaystyle=$
$\displaystyle\sqrt{\frac{\kappa_{y}}{\sinh(\kappa_{y}L)}}\,.$ (S45)

Thus, the wave function is

$\Psi_{+}(x)=\left(\matrix{1\cr
0}\right)\sqrt{\frac{\kappa_{y}}{\sinh(\kappa_{y}L)}}\,\exp\left(-\kappa_{y}x\right)\,.$
(S46)

Equation (S46) describes the behaviour of the ZEM wave function familiar from
the tight-binding model in the main text Fig. 2A. As seen in fig. S5B, the
wave function for $\kappa_{y}>0$ is exponentially growing towards the left
edge which shows the high localization on the zigzag edge. A similar
localization is seen for the bearded edge as shown in fig. S5C, when
$\kappa_{y}<0$. Finally, the wave function is constant throughout the ribbon
interior when $\kappa_{y}=0$ as shown in fig. S5D. To summarize, the variation
of $\kappa_{y}$ around the Dirac point ($\kappa_{y}=0$) shifts the ZEM wave
function from one edge to another while preserving the sublattice polarization
and its chiral character. Only in the vicinity of the Dirac point the wave
function loses localization. In all of these cases, the wave function is still
the eigenfunction of the chiral symmetry operator $\mathcal{S}=\sigma_{z}$.
Thus, the asymmetric electron density distribution attributed to the ZEM
chirality in Fig. 2 is well-described in the continuum model and further
supported by a deeper argument based on the $\mathcal{S}$ operator
eigenfunctions. Note that once $B$ sublattice provides outer atoms of the
hbGNR, the bearded and zigzag edges swap their positions as seen in fig. S5A.
This is also accounted for by Eq. (S42). The $B$ sublattice analog of Eq.
(S46) carries not only a different eigenvector but also an additional minus
sign in the exponent:

$\Psi_{-}(x)=\left(\matrix{0\cr
1}\right)\sqrt{\frac{\kappa_{y}}{\sinh(\kappa_{y}L)}}\,\exp\left(\kappa_{y}x\right)\,,$
(S47)

that swaps the wave function localization from the left (right) to right
(left) side for $\kappa_{y}>0$ ($\kappa_{y}<0$). The swapping of the wave
function localization corresponds to the switching of zigzag and bearded
edges.

The solutions given by Eqs. (S46) and (S47) are equivalent to Eq. 9 describing
isolated Majorana bound states in Ref. (?) but they do not describe Majorana
fermions. For the $2\times 2$ Dirac equation $i\partial_{t}\Psi=H\Psi$, with
$H$ given by Eq. (S39), the particle-hole (charge conjugation) symmetry
operator is $\Xi=\sigma_{x}K$, where $K$ denotes complex conjugation. Acting
with this $\Xi$ operator on $\Psi_{\pm}$ in Eqs. (S46) and (S47), we see that
$\Xi\Psi_{\pm}\neq\Psi_{\pm}$, which means the found ZEM is not a Majorana
fermion.

It is known that chiral states appear whenever TR symmetry $\mathcal{T}$ is
broken. Usually, these chiral states require presence of an external magnetic
field or an intrinsic net zero magnetic flux like in the Haldane model.
However, the breaking of $\mathcal{T}$ occurs for the ZEM of a hbGNR due to
the peculiar hbGNR edges and boundary conditions. For the above $2\times 2$
Dirac equation $\mathcal{T}=i\sigma_{y}K$. Hence,
$\mathcal{T}\Psi_{+}(\kappa_{y})=-\Psi_{-}(-\kappa_{y})$ which means that
$\Psi_{+}$ mode does not have a TR partner. Although
$\mathcal{T}\Psi_{+}(\kappa_{y})$ is an eigenstate of $H(-\kappa_{y})\equiv
H_{\mathbf{K}^{\prime}}(\kappa_{y})$, with
$H_{\mathbf{K}^{\prime}}=-i\sigma_{x}\partial_{x}+\sigma_{y}\kappa_{y}$ being
Hamiltonian in K′ valley, this $\mathcal{T}\Psi_{+}(\kappa_{y})$ eigenstate
corresponds to imposing boundary conditions on another sublattice. A more
detailed picture is the following. In the K′ valley one finds ZEM solutions
$\Psi^{\mathbf{K}^{\prime}}_{+}=\left(\matrix{1\cr
0}\right)\sqrt{\frac{\kappa_{y}}{\sinh(\kappa_{y}L)}}\,\exp\left(\kappa_{y}x\right)$
and $\Psi^{\mathbf{K}^{\prime}}_{-}=\left(\matrix{0\cr
1}\right)\sqrt{\frac{\kappa_{y}}{\sinh(\kappa_{y}L)}}\,\exp\left(-\kappa_{y}x\right)$,
so that in general we have
$\mathcal{T}\Psi^{\mathbf{K}}_{-}=\Psi^{\mathbf{K}^{\prime}}_{+}$ and
$\mathcal{T}\Psi^{\mathbf{K}}_{+}=-\Psi^{\mathbf{K}^{\prime}}_{-}$, where
$\Psi^{\mathbf{K}}_{\pm}$ are given by Eqs. (S46) and (S47). In the hbGNR, the
TR partner of ZEM is absent due to the inversion of the sublattice
polarization by $\mathcal{T}$ which provides another argument for supporting
the chiral character of the ZEM.

The formally defined charge conjugation and TR symmetry operators are defined
with respect to the sublattice degree of freedom or pseudo-spin. In order to
see how this makes them different with respect to the real spin analogs it is
useful to consider symmetry classes as we do further.

At the level of continuum model, the system is described by fundamental
symmetries such that $\mathcal{T}^{2}=-\mathbf{1}_{2\times 2}$,
$\Xi^{2}=\mathbf{1}_{2\times 2}$, and $\mathcal{S}=\mathbf{1}_{2\times 2}$,
where $\mathbf{1}_{2\times 2}$ stands for a $2\times 2$ unit matrix. This
means the low-energy Hamiltonian of our one-dimensional system belongs to the
DIII Altland-Zirnbauer symmetry class (?). In contrast, the full Hamiltonian
from the tight-binding model belongs to BDI symmetry class, where it is
characterized by $\mathcal{T}^{2}=\mathbf{1}_{2\times 2}$,
$\Xi^{2}=\mathbf{1}_{2\times 2}$, and $\mathcal{S}=\mathbf{1}_{2\times 2}$.
This symmetry classification can be verified by taking the graphene
Hamiltonian

$H(\mathbf{k})=\left(\matrix{0&h(\mathbf{k})\cr
h^{\ast}(\mathbf{k})&0}\right)$

and noticing that $\mathrm{Re}[h(\mathbf{k})]=\mathrm{Re}[h(-\mathbf{k})]$ and
$\mathrm{Im}[h(\mathbf{k})]=-\mathrm{Im}[h(-\mathbf{k})]$. Then we clearly see
that $\mathcal{T}H(\mathbf{k})\mathcal{T}^{-1}=H(-\mathbf{k})$ for
$\mathcal{T}=\mathbf{1}_{2\times 2}K$, which is the definition of the TR
symmetry operator for a spinless particle since
$\mathcal{T}^{2}=\mathbf{1}_{2\times 2}$. The definition of the chiral
symmetry operator does not change for the full Hamiltonian,
$\mathcal{S}=\sigma_{z}$. However, the charge conjugation symmetry operator
becomes $\Xi=\mathcal{S}\mathcal{T}=\sigma_{z}K$ so that
$\Xi^{2}=\mathbf{1}_{2\times 2}$. Therefore, the BDI class is identified. In
the above-mentioned identification of the DIII symmetry class, the TR symmetry
operator acts on the sublattice pseudo-spin instead of the real spin of an
electron. This important difference between the continuum and tight-binding
models must be taken into account in interpretation of the $\Xi$\- and
$\mathcal{T}$-breaking results. For instance, the $\mathcal{T}$ acting on the
real spin is unbroken for the tight-binding model.

Having now solutions for both K and K′ valleys, we notice that if the
intervalley and inter-sublattice (i.e. inter-ribbon) mixings were achievable
simultaneously, one would be able to construct a Majorana pseudo-fermion (i.e.
characterized by a pseudo-spin in lieu of a real spin), since there would then
be a way of exploiting the solutions as follows:
$\Psi_{\mathrm{M}}=\Psi^{\mathbf{K}}_{+}+\Psi^{\mathbf{K}^{\prime}}_{-}=\Psi^{\mathbf{K}}_{+}+\Xi\Psi^{\mathbf{K}}_{+}=\left(\matrix{1\cr
1}\right)\sqrt{\frac{\kappa_{y}}{\sinh(\kappa_{y}L)}}\,\exp\left(-\kappa_{y}x\right)$,
where $\left(\matrix{1\cr 1}\right)$ is an eigenvector of $\Xi\sim\sigma_{x}$.
This relation between the hbGNR ZEM and a Majorana pseudo-fermion is similar
to the connection between the ZEM of a 2D Dirac fermion bounded to a real
space vortex and a Majorana fermion, cf. with Eq. (35) in Ref. (?). The hbGNR
ZEM also resembles the edge ZEM hosted in a finite piece of a square lattice
with a real space vortex, see Fig. 2e of Ref. (?). Since similar spatial
vortexes have been constructed for an infinite honeycomb lattice (?), the ZEM
analogs must also exist in finite structures on the honeycomb lattice. In some
sense hbGNR ZEM can be considered as a 1D analog of those modes to be found on
the honeycomb lattice.

A brief comment is needed on the exponential growth in Eq. (S46) when
$\kappa_{y}\rightarrow+\infty$ and $x>0$. Formally, $\kappa_{y}$ is not
bounded variable in the continuum model. This unbounded $\kappa_{y}$ is the
source of potential misinterpretations. In reality, the range of electron wave
vectors is limited by the Brillouin zone. Therefore, the assumption
$\kappa_{y}\rightarrow+\infty$ is not a valid physical case.

To conclude this subsection, the analytical description based on the Dirac
equation has revealed a range of distinctive features making the ZEM of hbGNR
very different from all previously reported ZEMs in condensed matter systems
and quantum field theory.

###  S4.3 Topological properties analysis: non-linear dynamics and index
theory perspective

A general intuition tells us that whenever one encounters an edge localized
state, it must be topological and characterized by some non-zero topological
invariant. It happens, however, that depending on a chosen topological
perspective, the situation may be quite the opposite. By noticing that the
Eqs. (S31) and (S32) depend only on one independent variable $x$, we can draw
an analogy with typical dynamic equations featuring time as an independent
variable. This allows one to make use of the whole apparatus of non-linear
dynamics and chaos theory (?) for the analysis of the given equations. In what
follows, we analyse both bulk and edge eigenstates of hbGNRs via phase
portraits and indexes of their phase space trajectories.

The phase portrait of a dynamic system is the graphical representation of a
vector field in phase space. The phase space is the space featuring the
coordinate and momentum as independent variables. Generally speaking, for any
second order ordinary differential equation, the phase space coordinates can
be defined as the sought function and its first derivative.

The principle difference between classical dynamic equations and Eqs. (S31)
and (S32) is that classical coordinates are always real while wave functions
admit complex values. However, the wave functions can always be chosen to be
real provided that the magnetic field is absent in the system. For the
equations in question this can be achieved by the following unitary transform
of the Hamiltonian:

$U=\left(\matrix{1&0\cr 0&-i}\right)$

and $\mathcal{H}(\kappa_{y})=UH(\kappa_{y})U^{\dagger}$, so that the new set
of equations for analysis is