# Exotic edge states of $C_{3}$ high-fold fermions in honeycomb lattices

L. Madail1,2, R. G. Dias2, J. Fernández-Rossier1111On permanent leave from
Departamento de Física Aplicada, Universidad de Alicante, Spain 1
International Iberian Nanotechnology Laboratory (INL), Av. Mestre José Veiga,
4715-330 Braga, Portugal
2 Department of Physics and I3N, University of Aveiro, 3810-193 Aveiro,
Portugal

###### Abstract

A generalization of the graphene honeycomb model to the case where each site
in the honeycomb lattice contains a $n-$fold degenerate set of eigenstates of
the $C_{3}$ symmetry has been recently proposed to describe several systems,
including triangulene crystals and photonic lattices. These generalized
honeycomb models are defined by $(n_{a},n_{b})$, the number $C_{3}$
eigenstates in the $a$ and $b$ sites of the unit cell, resulting in
$n_{a}+n_{b}$ bands. Thus, the $(1,1)$ case gives the coventional honeycomb
model that describes the two low-energy bands in graphene. Generalizations,
such as $(2,1)$, $(2,2)$ and $(3,3)$ display several non-trivial features,
such as coexisting graphene-like Dirac cones with flat-bands, both at zero and
finite-energy, as well as robust degeneracy points where a flat-band and a
parabolic band meet at the $\Gamma$-point. Here, we explore the edge states of
this class of crystals, using as reference triangulene crystals, and we find
several types of edge states absent in the conventional $(1,1)$ honeycomb
case, associated to the non-trivial features of the two-dimensional (2D) bands
of the high-fold case. First, we find dispersive edge states associated to the
finite-energy flat-bands, that occur both at the armchair and zigzag
termination. Second, in the case of non-centrosymmetric triangulene crystals
that lead to a $S=1$ Dirac band, we have a bonding-antibonding pair of
dispersive edge states, localized in the same edge so that their energy
splitting is reduced as their localization increases, opposite to the
conventional behavior of pairs of states localized in opposite edges. Third,
for the $(3,3)$ case, that hosts a gap separating a pair of flat conduction
and valence bands, we find non-dispersive edge states with $E=0$ in all edge
terminations.

##  I Introduction

It has long been known that the topology of the Fermi surface of materials has
a gigantic impact in their electronic properties. The discovery of graphene
promoted a way to rationalize our understanding of this interplay by analogy
with massless Dirac fermions [1]. This made it possible to understand a number
of graphene’s exotic properties such as electron-hole symmetry [2], the
unconventional Landau Level spectrum [1], Klein tunneling [3, 4], the finite
value of its minimal conductivity [5], the universal value of the optical
transmission coefficient [6], and the presence of zero-energy edge states in
zigzag boundaries [7].

This paradigm has been extended more recently to the case of three-dimensional
systems with diabolic points giving rise to Dirac cones, whose low-energy
quasiparticles are isomorphic to Weyl fermions [8]. The common theme is the
presence of degeneracy points in the band structure, that endows them with
finite Berry curvatures and non-trivial quantum geometry [9, 10].
Unsurprisingly, analogy with a pre-existing relativistic theory is not
necessary, and concepts such as new fermions [11] and high-fold fermions [12]
has been put forward in order to dub exotic quasiparticles with non-trivial
properties and without a relativistic analog.

It was recently proposed by one of us [13] a generalization of the graphene
paradigmatic model, where quantum particles move in a honeycomb lattice, to
the case where every site in the lattice hosts multiple orbitals that are
eigenstates of the $C_{3}$ symmetry operator. The corresponding eigenvalues
are given by $e^{i\nu\frac{2\pi}{3}}$,where $\nu=0,\pm 1$. Here we shall refer
to the quasiparticles in these systems as $C_{3}$ high-fold fermions, although
the analysis will be relevant to bosonic quasiparticles as well. The $C_{3}$
high-fold model provides an accurate description of the low-energy bands of
two similar systems, the junctions of a crystalline network of graphene
nanoribbons [14] and triangulene two-dimensional crystals [13], a honeycomb
lattice whose unit cell contains a pair of $[n]$-triangulenes – equilateral
triangle shaped graphene nanoislands with a lateral dimension of $n$ benzene
rings. Importantly, $[n]$-triangulenes have a sublattice imbalance of
$N_{A}-N_{B}=n-1$, leading to the existence of $n-1$ sub-lattice polarized
zero modes [15], that govern their low-energy properties and can be chosen as
eigenstates of the $C_{3}$ symmetry operator [16].

Using both density-functional based calculations and the standard tight-
binding (TB) model for graphene, it was found [13] that honeycomb lattices,
with a unit cell made of two triangulenes dimensions $n_{a},n_{b}$, denoted in
the following $[n_{a},n_{b}]$-triangulenes, have a set of $n_{a}+n_{b}-2$ low-
energy bands, providing a realization of generalized honeycomb model with
$(n_{a}-1,n_{b}-1)$. Hence, we found that the low-bands of
$[n_{a},n_{b}]$-triangulene crystals are described with the
$(n_{a}-1,n_{b}-1)$ generalized honeycomb lattice, with the constraint that
$\sum_{i=1,n}\nu_{i}=0$ for every triangulene. Specifically, the calculations
show that [13]:

  1. 1.

The two low-energy bands of the $[2,2]$-triangulene crystal were isomorphic to
graphene, with reduced bandwidth, in agreement with previous work [17].

  2. 2.

The three low-energy bands of the $[2,3]$-triangulene crystal are isomorphic
to the so-called $S=1$ Dirac model [18] and the Lieb lattice [19, 20], with a
single Dirac cone at the $\Gamma$-point, crossed at $E=0$ by a flat-band.

  3. 3.

The four low-energy bands of the $[3,3]$-triangulene crystal are isomorphic to
the $p_{x}-p_{y}$ honeycomb lattice [21, 22].

  4. 4.

The six low-energy bands of the $[4,4]$ are made of two pairs of Kagome-like
bands, separated by a gap, where both the conduction and valence band are
flat, in agreement with work by other authors [23, 24].

Most of theory work in triangulenes has focused on their interacting
properties associated to electron-electron interactions, for individual
molecules [15], dimers [25, 26, 27], monostrand chains [28, 29, 30, 31] and
two-dimensional crystals [32, 33, 34]. Very recently, the edge states of non-
interacting electrons in the $[4,4]$-triangulene crystal have been also
studied [23]. In this work we undertake the systematic study of edge states of
the four representative classes of crystals listed above and relate their
properties to the underlying isospin of the fermions in the low-energy bands
of the two-dimensional lattice. For that matter, we adopt the same strategy of
the seminal paper [35] on graphene edge states of Nakada et al. and we study
the band structure of one-dimensional triangulene ribbons. We do this at two
different levels. First, we use the same tight-binding model of the two-
dimensional lattice [13], for the complete atomistic model [see Fig. 2(a)].
Second, we use a reduced tight-binding model, that includes only the low-
energy states of the triangulenes and maps them into an effective lattice.
This comparison permits us to ensure that the edge states found in the
complete calculation are actually arising from the low-energy bands associated
to the zero-modes of the triangulenes.

We stress here that our results will be of interest not only for the case of
graphene, but also to any other system that can be described with a honeycomb
lattice where every site hosts several modes with $C_{3}$ symmetry. Thus, very
much like there are many artificial graphene systems, where honeycomb lattices
host one orbital per site, realized with adatoms [36], photons [37],
polaritons [38], surface electrons [39], electric circuits [40], our work
provides the generalization to analogous realizations of the multi-orbital
$C_{3}$ symmetric honeycomb lattice.

The rest of this work is organized as follows. In section II we introduce the
complete tight-binding model that describes triangulene crystals as well as
the reduced low-energy tight-binding models relevant for the $[2,2]$, $[2,3]$,
$[3,3]$ and $[4,4]$ crystals. In addition, we review the evanescent Bloch
states method to determine the dispersion of edge states and establish the
topological invariant used to classify them. In section III, IV, V and VI we
present the electronic structure for the four types of crystals, comparing the
results for the complete and the reduced models, focusing on the properties of
the low-energy states and on their topological regime. In section VII we
summarize our main results.

##  II Methods

###  II.1 Complete tight-binding model

Figure 1: Single particle properties of graphene triangular nanoislands. (a)
Example for $[2]$-triangulenes composed of (top) excess $A$ or (bottom) excess
$B$ sublattice sites. (b) $[n_{a},n_{b}]$-triangulene dimer. (c)
Representation of the zero modes of $n=2,3,4$ triangulenes, chosen as
eigenstates of $R_{2\pi/3}$ rotation operator with eigenvalues $e^{i\omega}$
and $\omega=0,\pm\frac{2\pi}{3}$. Circle size scales with the amplitude of the
eigenfunction while the color represent their phase.

The single particle properties of $[n_{a},n_{b}]$-triangulenes, with one
orbital per carbon atom on an unit cell dimer formed by two triangulenes of A
and B types, is described by the following bipartite tight-binding Hamiltonian

$H_{\text{TB}}=-t_{1}\sum_{<i,j>}a^{\dagger}_{i}b_{j}-t_{3}\sum_{\lll
i,j\ggg}a^{\dagger}_{i}b_{j},$ (1)

where $a^{\dagger}_{i}$ ($b^{\dagger}_{i}$) is the fermionic operator that
creates an electron at site $i$ belonging to the $A$($B$)-sublattice. The
first (second) sum is taken between first (third) neighbor sites in each
triangulene whose vicinity is defined in Fig. 1 and includes all intercell
hopping elements at the binding sites of the dimer. Following previous work
[13], we assume $t_{2}=0$, so that the model is bipartite ensuring thereby
local charge neutrality, in agreement with density functional based
calculations.

The single-particle spectrum of individual $[n]$-triangulenes described with
equation (1) has a set of electron-hole symmetric states with finite-energy,
well separated by a gap, and $n-1$ in-gap sublattice-polarized zero modes with
$E=0$. These modes can be chosen to be eigenvectors of the $C_{3}$ symmetry
operator, on account of the point symmetry of triangulenes. Intermolecular
hybridization among the zero modes of adjacent triangulenes, found for
instance in density functional theory calculations, can only be established by
third-neighbor hopping [13]. This stems from the fact that the zero modes are
hosted in the majority sublattice and the intermolecular binding sites belong
to the minority sublattice. Thus, both $t_{1}$ and $t_{2}$ are linked to the
binding sites where zero modes have zero probability amplitude and the leading
term that produces intermolecular hybridization is $t_{3}$.

It was found [13] that a good agreement between spin-unpolarized DFT and
tight-binding energy bands of $[n_{a},n_{b}]$-triangulene 2D crystals is
obtained if we take $t_{3}=0.1t_{1}$. These calculations display a wealth of
weakly dispersive half-filled narrow bands located around the Fermi energy,
well separated from higher/lower energy bands. Centrosymmetric
$[n,n]$-triangulene crystals have $2(n-1)$ narrow bands. For $n\geq 3$ they
always feature pairs of electron–hole conjugate low-energy flat-bands,
corresponding to states localized in supramolecular hexagonal rings. These
flat-bands are very different from $E=0$ bands, as they are not sublattice
polarized and feature intermolecular hybridization. Non-centrosymmetric
$[n_{a},n_{b}]$-triangulene crystals feature $|n_{a}-n_{b}|$ narrow bands at
$E=0$ due to the bipartite character of the lattice. More specifically, one is
able to emulate a graphene-like spectrum for $n_{a}=n_{b}=2$, spin-1 Dirac
electrons for $n_{a}=2,n_{b}=3$, $p_{x,y}$-orbital honeycomb physics for
$n_{a}=n_{b}=3$, and a gapped system with flat valence and conduction bands
for $n_{a}=n_{b}=4$. [13].

###  II.2 Reduced model

The total number of energy bands in the single-particle energy spectrum of
$[n_{a},n_{b}]$-triangulene ribbons is determined by the number of atoms
$N_{r}=\left[(n_{a}+2)^{2}+(n_{b}+2)^{2}-6\right]N_{rt}$ with $N_{rt}$ the
number of triangulene dimers in the ribbon unit cell. For the 2D case, the
$n_{a}+n_{b}-2$ weakly dispersive bands around $E=0$ are embedded inside a
large gap isolated from higher energy bands. Therefore, low-energy physics at
half-filling is expected to be governed by electrons occupying zero modes and
these include all $n_{a/b}-1$ zero modes from each $[n_{a/b}]$-triangulene
since intermolecular coupling is dominated by the small third-neighbor term.
The effective 2D Bloch Hamiltonian in the basis of the $n_{a}+n_{b}-2$ zero
modes as in [13] reads

$H_{\text{eff}}(\vec{k})=\begin{pmatrix}0_{n_{a}-1}&\tau(\vec{k})\\\
\tau^{\dagger}(\vec{k})&0_{n_{b}-1}\end{pmatrix}.$ (2)

This bipartite Hamiltonian describes sublattice polarized modes, $n_{a}-1$
($n_{b}-1$) of those populate sites belonging to the $A$ ($B$) sublattice [see
Fig. 1(b)]. The elements of the $(n_{a}-1)\times(n_{b}-1)$ off-diagonal matrix
$\tau(\vec{k})$ are given by

$\tau_{ab}(\vec{k})=\sum_{m=0,1,2}\tau^{ab}_{m},\qquad\tau^{ab}_{m}=e^{i\vec{k}.\vec{R}_{m}}t_{m}^{ab},$
(3)

determined using the primitive vectors
$\vec{R}_{1,2}=(n_{a}+n_{b}+1)\vec{a}_{1,2}$ that define the hexagonal lattice
of the triangulene dimer in terms of the graphene primitive vectors
$\vec{a}_{1,2}$. Intracell ($m=0$) matrix elements are determined by the
overlap between the components of the zero mode wavefunction at the binding
sites. Meanwhile, the collection of states in the same sublattice reflect a
rotational degree of freedom associated with $C_{3}$ symmetry. Inevitably,
intercell ($m=1,2$) contributions can be expressed in terms of the previous
ones using the phases ($\omega_{a/b}$) of the zero modes upon a $R_{2\pi/3}$
rotation operation as

$t_{1,2}^{ab}=e^{\pm i(\omega_{b}-\omega_{a})}t_{0}^{ab}.$ (4)

We thus see that intermolecular hoppings, represented in the basis of $C_{3}$
eigenstates, can acquire a Peierls-like phase factor [13].

It was found that the $n_{a}+n_{b}-2$ in-gap energy bands of the complete
tight-binding model of Eq. 1 are in excellent agreement with the bands
computed with the reduced model of Eq. 2. As we discuss now, the same
statement applies also for the case of triangulene ribbons.

###  II.3 Ribbon structures

In order to study the edge states we compute the electronic structure of
ribbons along the two main crystallographic directions of the honeycomb
lattice, that generate edges with zigzag and armchair geometry [35].

The Hamiltonian of the $[n_{a},n_{b}]$-triangulene ribbon with $N_{rt}$ dimers
in the ribbon unit cell can be expressed as the following banded block-
Toeplitz matrix

$H_{r}(k)=\begin{pmatrix}\tau_{r,0}&\tau_{r,1}(k)&&&&\\\
\tau_{r,1}^{\dagger}(k)&\tau_{r,0}&\tau_{r,1}^{*}(k)&&&\\\
&\tau_{r,1}^{\top}(k)&\tau_{r,0}&\tau_{r,1}(k)&\\\
&&\tau_{r,1}^{\dagger}(k)&\tau_{r,0}&\ddots\\\ &&&\ddots&\ddots\\\
\end{pmatrix}.$ (5)

where $\tau$ are matrices. Within the full tight-binding approach, the
dimension of each $\tau_{r,\sigma}$ entry is $N_{t_{a}}\times N_{t_{b}}$ with
$N_{t_{a/b}}$ the number of atoms in each $[n_{a/b}]$-triangulene, making the
Hamiltonian a total of
$(N_{t_{a}}+N_{t_{b}})N_{rt}\times(N_{t_{a}}+N_{t_{b}})N_{rt}$ matrix. In
contrast, within the reduced tight-binding approach, the dimension of the
$\tau_{r,\sigma}$ matrices reduces to $(n_{a}-1)\times(n_{b}-1)$ if we
consider the effective low-energy model on the basis of the zero modes. Block-
diagonal $\sigma=0$ entries account for couplings within the ribbon unit cell
and are block anti-diagonal following Eq. 2 while off-diagonal $\sigma=1$
matrices connect adjacent unit cells. Their corresponding elements for a
zigzag ribbon can be explicitly written in terms of the established entries of
the 2D effective Bloch Hamiltonian as

$\displaystyle\tau_{r,0}^{ab}$ $\displaystyle=\tau_{0}^{ab},$ (6)
$\displaystyle\tau_{r,1}^{ab}$
$\displaystyle=\tau_{1}^{ab}(k_{x}=0,k_{y}=0)+\tau_{2}^{ab}(k_{x},k_{y}=0),$

where the momentum $k_{x}$ has been renormalized to the width of each ribbon
unit cell.

In both edge configurations, the atoms on the triangular islands that populate
the boundaries belong to both sublattices. Since dimers connect two
triangulenes with opposite majority sublattices, sublattice imbalance between
edges is entirely suppressed in armchair ribbons as well as in zigzag ribbons
provided that $n_{a}=n_{b}$. The top (bottom) zigzag edge is composed of
$[n_{a}]$-triangulenes ($[n_{b}]$-triangulenes) with extra $(n_{a}-1)$
[$(n_{b}-1)$] atoms in the majority sublattice, while both armchair edges are
made of triangulene dimers with extra $|n_{a}-n_{b}|$ atoms in the majority
sublattice.

###  II.4 Evanescent Bloch waves for the calculation of edge states

We now review a method [41, 42, 43, 44, 45, 46] to obtain the energy
dispersion and the wave function of edge states, based on the generalization
of the Bloch states to the case of imaginary wave vectors. Conventional Bloch
states are eigenstates of the discrete translation operator that leaves the
crystal invariant. In a 1D crystal, relevant for our case, we define the
translation operator as $\mathcal{T}=\sum_{j\in\mathbf{Z}}\ket{j}\bra{j+1}$,
where $j$ identifies the unit cells of the crystal. We label the eigenstates
and eigelvalues as $\mathcal{T}\ket{z}=z\ket{z}$ with

$\ket{z}_{\gamma}=\sum_{j\in\mathbb{Z}}z^{j}\ket{j},z\in\mathbb{C}.$ (7)

Conventionally, the eigenvalues $z$ are phases $e^{ik}$, where $k$ is chosen
as to match the periodic boundary conditions. In this case, repeated
application of the translation operator results in a phase modulation of the
wave function, without affecting the amplitude. Analytical extension to the
complex plane permits one to define real eigenvalues, $z$, that result in
evanescent/exploding solutions, depending on whether $z<1$ or $z>1$. These
non-unitary representations of the states can only describe edge states.

Using the standard ansatz in $\langle j|\psi_{z}(j)\rangle=z^{j}|u(z)\rangle$,
we need to solve

$h(z)\ket{u(z)}=\epsilon(z)\ket{u(z)},$ (8)

where $h(z)$ is the usual square matrix of dimension $n$, the number of sites
of the crystal cell. The conventional Brillouin zone (BZ) is the unit circle
with $|z|=1$ in the complex plane.

For $|z|\neq 1$, $h(z)$ is non-Hermitian and its eigenstates have an
exponential behavior. The full set of solutions of the bulk equation is
obtained by requiring that the corresponding eigenspectrum is real. The
diagonalization of Eq. 8 yields two roots $\\{z_{n}\\}$ for a given energy
$\epsilon$. Since we are looking for solutions of the open system, the
boundary equation is introduced at this stage to determine which $\\{z_{n}\\}$
belong to the energy spectrum of $H$.

For a given $k_{x}$, the Hamiltonian of the ribbons defines a finite-size one-
dimensional chain. We can therefore, implement this approach for every $k_{x}$
and find thereby the analytic form of the edge modes of triangulene ribbons
assuming the generalized solutions only in the open boundary direction
($\gamma=x$) while using the Bloch theorem in the periodic direction to be
valid in the bulk region, thus obtaining the energy dispersion
$\epsilon(z_{n},k)$. The allowed $\\{z_{n}\\}$ solutions at the edge are
obtained by addressing the boundary equation, whose definition may not be
straightforward in models with multiple orbitals per unit cell. In the
following sections, we make this boundary equation self-evident by means of a
basis rotation via bonding and anti-bonding linear combination of the Wannier
states associated with the zero modes of each $[n]$-triangulene.

###  II.5 Topological invariant

The notion of the Zak phase for one-dimensional models is well defined as the
integral of the Berry connection across the Brillouin zone [47]. In two-
dimensional systems, this integration should be taken on a cut of a 2D
Brillouin zone in a direction transverse to the ribbon orientation.
Translation invariance along the direction $k_{\parallel}$ parallel to the
boundary guarantees a continuous parameter in space in which the state
$\ket{u(\vec{k})}$ can travel on a closed path. The parallel component is real
and still a good quantum number, while the normal component is complex
corresponding to the exponential decay factor [48]. Therefore, we can define
the Zak phase $\mathcal{Z}(k_{\parallel})$ of an isolated band as the
integration of the Berry connection across the 2D Brillouin zone along the
perpendicular direction $k_{\bot}$ in

$\mathcal{Z}(k_{\parallel})=i\oint\braket{u(\vec{k})}{\partial_{k_{\bot}}u(\vec{k})}dk_{\bot}.$
(9)

For finite but sufficiently long periodic systems ($N_{rt}\gg 1$), the
numerical computation of Eq. 9 can be carried out using a discrete set of
$k_{\bot}$-points in the first BZ as follows

$\mathcal{Z}(k_{\parallel})=-\operatorname{Im}\operatorname{ln}\prod_{j}^{N_{k_{\bot}}-1}\braket{u(k_{j})}{u(k_{j+1})}.$
(10)

Chiral and inversion symmetric models have quantized Zak phases that can only
take a trivial or non-trivial value ($0\,\wedge\,\pi\mod 2\pi$). According to
the bulk-boundary correspondence, there are low-energy excitations at the
interface between topologically distinct insulators with a trivial and a non-
trivial phase. This means that, when the phase is $\pi$-quantized, a ribbon
with open boundary conditions along $k_{\bot}$ (interface with vacuum) and
invariant by translation along $k_{\parallel}$ supports interface modes.

In the multiband formalism, the topological property of a bandgap can be
determined by the summation of the $\mathcal{Z}_{n}(k_{\parallel})$ phases of
the set of bands below this gap [49, 50]. Specifically for insulating
inversion symmetric models and as $k_{\bot}$ is integrated along the path
perpendicular to the edge, the total phase
$\sum_{n\in\text{occ}}\mathcal{Z}_{n}(k_{\parallel})/\pi$ of the occupied
bands is equivalent to the winding of the trajectory of
$\det\tau(k_{\parallel})$ and predicts the number of pairs of flat edge states
in the midgap [50, 51]. Here, we use this formalism to predict the topological
phase of all bandgaps for each triagulene ribbon, including the non-central
ones, where pairs of finite-energy dispersive edge states are found. The case
of degenerate Hamiltonians can be treated using the non-Abelian formalism
which measures the geometric phase associated with the matrix of a subspace of
the Hilbert space [52, 53]. However, this approach is not able to predict edge
states between degenerate bands. Herein, the presence of singular points of
degeneracy for each $k_{\parallel}$ was treated by making use of a smooth
gauge choice where phases are roughly uniformly distributed throughout the
loop. In the case where there is a finite number of degeneracies in the band
structure for each $k_{\parallel}$, this gauge allows for these singular
points to be excluded from the numerical calculation of
$\mathcal{Z}_{n}(k_{\parallel})$ without significantly affecting the results.
Additional considerations on this approach can be found in appendix B.

##  III $[2,2]$-triangulene ribbons

Figure 2: Full model for $[2,2]$-triangulene ribbons in the zigzag and
armchair configurations. The full tight-binding diagonalization reveals a
total of $N_{r}=26N_{rt}$ bands, as shown in the spectrum (a) for a
$[2,2]$-triangulene zigzag ribbon with $N_{rt}=20$ dimers, $2N_{rt}$ of those
isolated around zero-energy. A closeup near this point signals an identical
spectrum to the reported graphene ribbons, both in the (b) zigzag and (c)
armchair configurations.

The results for $[2,2]$-triangulene ribbons in both the zigzag and armchair
configurations with $N_{rt}=20$ triangulene dimers, with 26 sites per dimer,
in the unit cell are shown in Fig. 2. The full tight-binding models yield
$N_{r}=26N_{rt}$ bands, wherein $(n_{a}+n_{b}-2)N_{rt}=2N_{rt}$ of those are
located around the Fermi energy, well separated from the remainder $24N_{rt}$
higher/lower energy bands [see Fig. 2(b)].

The configuration for triangulene ribbons is established by replacing the unit
cell dimer of the triangular Bravais lattice with a triangulene dimer [see
Fig. 2(a)]. In fact, every $[2]$-triangulene hosts exactly one isolated zero
mode and intermolecular coupling between these states is solely assured by
$t_{3}$. Thus, one expects the low-energy bands for $[2,2]$-triangulene
ribbons to be isomorphic to the tight-binding model for a honeycomb ribbon
with nearest neighbor approximation and one orbital per site, provided that we
keep the analogy in the choice of open boundary conditions [54]. Nevertheless,
in constrast with the models of graphene with nearest-neighbor approximation,
the narrower low-energy bandwidth of triangulene ribbons is given by the
smaller amplitude of the third-neighbor hopping term. Therefore, full
concurrence between the low-energy bands of $[2,2]$-triangulene ribbons and
the $p_{z}$-orbitals of graphene ribbons can be obtained for an effective
hopping term of $\tilde{t}_{g}=\frac{t_{3}}{3}$ with energy bandwidth of
$6\tilde{t}_{g}$.

In the zigzag configuration, the boundary condition forces the wave functions
to vanish in a single sublattice at each edge such that sublattice symmetry is
broken and the two Dirac cones project onto inequivalent points $k_{x}=\pm
2\pi/3$ on the edge. This gives rise to two zero-energy edge states when
$|k_{x}|\in[2\pi/3,\pi]$. These states are localized at the edges, analog to
the extensively studied zigzag edge states in graphene [35, 54]. The single
band formalism was employed here to calculate the Zak phase of the band below
the midgap. We find a topological non-trivial phase $\mathcal{Z}=\pi$ for the
interval $|k|\in[2\pi/3,\pi]$ and a trivial phase otherwise, effectively
predicting the pair of zero-energy edge states. In the armchair configuration,
we also obtain a set of low-energy bands isomorphic to that of armchair
graphene ribbons [see Fig. 2(c,d)].

The analogy between the zigzag edge states of the $[2,2]$-triangulene crystals
and those of graphene is not perfect when we analyse the wave functions. In
the case of graphene, the amplitude of the edge mode associated to the zigzag
edge with $A$ atoms vanishes in the $B$ sublattice. In contrast, in the
triangulene case, the wave function of the edge states of the $A$-rich side
has a small weight in the $B$ rich triangulenes. This departure is only
captured in the complete model, but not in the reduced model. Therefore, it
has to be originated by a small influence of the higher/lower energy bands.

##  IV $[2,3]$-triangulene ribbons

In this section, we present the energy dispersion of $[2,3]$-triangulene
ribbons using the complete tight-binding model and compare with the reduced
effective model that only includes the zero modes. We also analyse the
localization and topological properties of edge modes found specifically in
zigzag ribbons, while considerations on armchair edge states are provided in
appendix A.

###  IV.1 Complete model

In Fig. 3 we show the low-energy bands corresponding to zigzag and armchair
ribbons, calculated with the complete tight-binding model. These low-energy
bands are identical to those obtained for the same geometries using the
minimal model, confirming that they are linear combinations of the zero modes
of the triangulenes.

We find three types of states for the $[2,3]$ ribbons. First, zero-energy
states associated to the sublattice imbalance of the $[2,3]$ dimer.
Specifically, we find as many $E=0$ states as $[2,3]$ dimers in the unit cell.
Therefore, only two of these states are localized at the edges. Second, an
electron-hole symmetric pair of dispersive edge states [see purple lines in
Fig. 3(a,b)], different from the graphene case in several counts:

  1. 1.

They are found both for zizag and armchair terminations in contrast with
graphene, where armchair boudaries do not have edge states;

  2. 2.

For the zigzag direction, these states are only localized in the edge
terminated with the $[2]$-triangulene [see Fig. 3a(i)];

  3. 3.

For the armchair configuration, edge states [colored bands in Fig. 3(b)] are
found in both sides of the ribbon and there are two pairs of electron-hole
symmetric edge states specifically at the zone boundary;

  4. 4.

The localization of the all edge states is maximal at the zone boundary, which
leads to a maximal electron-hole splitting, as a result of the inter-edge
nature of the states;

  5. 5.

The dispersive edge bands never reach $E=0$. In the next subsection we derive
an analytical formula for the dispersion of the zigzag edges states by
applying the evanescent Bloch wave method to the reduced tight-binding model;

  6. 6.

They cannot be predicted by the Zak phase, that is no longer quantized in the
absence of inversion symmetry, a characteristic of the $[2,3]$ non-
centrosymmetric triangulenes.

The third type of state that we found is the group of dispersive bulk modes
associated to the quantum confined subbands of the Dirac cone. Importantly,
because of the confinement induced by the finite width of the ribbons, there
is a gap at the $\Gamma$-point, in contrast with the 2D crystal. In the
neighborhood of $\Gamma$ we have found that the energy of the confined states
can be described by

$\epsilon(k_{x},n)=\pm\hbar
v_{F}\sqrt{k_{x}^{2}+\left(\frac{\pi\left(n-\frac{1}{2}\right)}{W\left(N_{rt}+1\right)}\right)^{2}},$
(11)

with $W$ the width of the ribbon and the Fermi velocity $\hbar
v_{F}^{[2,3]}=6\sqrt{\frac{3}{11}}t_{3}a_{cc}$ already established in the
literature [13]. In the next subsection we derive this formula, making use of
the reduced tight-binding model.

###  IV.2 Reduced model

Figure 3: Modelling $[2,3]$-triangulene ribbons with $N_{rt}=20$ dimers in the
unit cell. The numerical results of the energy dispersion for the (a) zigzag
and (b) armchair configurations show the presence of midgap edge states
highlighted in purple. In the former, the states localize in (i) the top edge
and (ii) decay through the bulk away from the zone boundary. The zigzag ribbon
unit cell effectively maps into a (c) diamond chain in the $C_{3}$ symmetric
basis and a (d) SSH model with dangling sites in the bonding/anti-bonding
combination basis, whose unit cells at the high-symmetry points are depicted
in (e).

We now shed light on some of the results of Fig. 3(a) using the reduced tight-
binding model. In the case of $[2,3]$-triangulenes, the minimal model
considers $n_{a}+n_{b}-2=3$ zero modes per triangulene dimer, one in the
$[2]$-triangulene, denoted as $\ket{a}$, two in the $[3]$-triangulene, denoted
as $\\{\ket{b_{1}},\ket{b_{2}}\\}$. Within the reduced tight-binding model,
the Hamiltonian matrix of a $[2,3]$ zigzag ribbon made of $N_{r}$ dimers, has
dimension $3*N_{r}$. The entries of this matrix are given by the intra- and
inter-cell contributions (see Eq. 3)

$\tau_{r,0}=\begin{pmatrix}0&\tau_{r,0}^{ab_{1}}&\tau_{r,0}^{ab_{2}}\\\
\tau_{r,0}^{ab_{1}\dagger}&0&0\\\
\tau_{r,0}^{ab_{2}\dagger}&0&0\end{pmatrix},\;\;\tau_{r,1}=\begin{pmatrix}0&0&0\\\
\tau_{r,1}^{ab_{1}}&0&0\\\ \tau_{r,1}^{ab_{2}}&0&0\end{pmatrix},$ (12)

whose entries follow

$\tau_{r,0}^{ab_{1,2}}=t_{[2,3]},\quad\tau_{r,1}^{ab_{1,2}}=\left(e^{i\theta}+e^{\mp
i\theta}e^{-ik_{x}}\right)t_{[2,3]},$ (13)

with $t_{[2,3]}=\frac{2t_{3}}{\sqrt{6\times 11}}$ and $\theta=\frac{2\pi}{3}$.
The structure of the $\tau_{r,\sigma}$ matrices permits one to map the low
energy states of the $[2,3]$ triangulene zigzag ribbon unit into a finite-size
diamond chain, where every site in the diamond lattice represents a zero mode
state. Within this effective model, hoppings have additional $e^{i\theta}$
Peierls-like factors, amounting to a magnetic flux within the unit cell of
$\phi_{uc}=\frac{\phi_{0}}{2\pi}2\theta=\frac{2}{3}\phi_{0}$ per plaquette,
normalized to $\phi_{0}=h/e$ [see Fig. 3(c)]. Naturally, the effective
magnetic flux within a plaquette is $k$-dependent due to the $e^{\pm ik_{x}}$
phase factors in $\tau_{r,1}$ matrices. We now perform a basis rotation using
bonding/anti-bonding linear combination between the two zero modes localized
in the B-type $[3]$-triangulene. This breaks $C_{3}$ symmetry between the
eigenfunctions. In the new basis set $\\{\ket{a},\ket{b_{+}},\ket{b_{-}}\\}$,
the $\tau_{r,\sigma}^{ab^{\pm}}$ matrices now satisfy

$\displaystyle\tau_{r,0}^{ab^{\pm}}$ $\displaystyle=(1\pm
1)t_{[2,3]},\quad\tau_{r,1}^{ab^{+}}=\left(-1-e^{-ik_{x}}\right)t_{[2,3]},$
(14) $\displaystyle\tau_{r,1}^{ab^{-}}$
$\displaystyle=\left(-1+e^{-ik_{x}}\right)\sqrt{3}it_{[2,3]}.$

As one of the intra-cell hopping terms vanishes, the antisymmetric combination
of states $\ket{b_{-}}$ become coupled only to adjacent triangulene dimers
[see Fig. 3(d)]. These states get fully decoupled from the chain at the
$\Gamma$-point where $\tau_{r,1}^{ab^{-}}=0$ and contribute with $N_{rt}$
zero-energy bands. At the $k_{x}=0$ point, the system behaves as an open
linear chain that couples the $2N_{rt}$ bonding states $\ket{a(b)_{+}}$
through an uniform hopping term
$\tau_{r,1}^{ab^{+}}=\tau_{r,0}^{ab^{+}}=2t_{[2,3]}$ [see Fig. 3(e)]. Note
that we have dropped the relative phase between $\tau_{r,0}^{ab^{+}}$ and
$\tau_{r,1}^{ab^{+}}$ since the closed loops disappear for $k_{x}=0$. Thus,
the energy of the dispersive subbands is given by

$\epsilon(k_{x}=0,k_{n})=4t_{[2,3]}\cos(k_{n}),$ (15)

where $k_{n}=\frac{\pi}{2N_{rt}+1}n$ and $n=1,\dots,2N_{rt}$. It can be seen
that Eq. 15 is consistent with Eq. 11 evaluated at $k_{x}=0$.

Analytical expressions can also be obtained at the zone boundary
$(k_{x}=\pi)$. Specifically, the three highly degenerate points with energies
$\epsilon(k_{x}=\pi)=\\{0,\pm 4t_{[2,3]}\\}$ are due to the trimerization
depicted in Fig. 3(e), on account of the vanishing $\tau_{r,1}^{ab^{+}}$
entries. As a result, their eigenfunctions have zero amplitude in the
$[2]$-triangulene of each trimer. Additionally, at the open boundary filled
with $[2]$-triangulenes, a dimer is formed leading to the chiral pair with
energies $\epsilon_{\text{edge}}(k_{x}=\pi)=\pm 2\sqrt{3}t_{[2,3]}$. In the
opposite boundary, a single decoupled site located in the $[3]$-triangulene
contributes with a zero-energy edge state.

The analytical form of the finite-energy edge state can be recovered by
assuming a generalized solution in the open $y$-direction, as introduced in
section II.4. We look for suitable solutions of the form
$\ket{\epsilon_{\text{edge}}}=e^{ik_{x}}z_{y}^{j}\ket{u(z)}$ in the basis set
of the zero modes $\\{\ket{a},\ket{b_{1}},\ket{b_{2}}\\}$ for the following
bulk Hamiltonian $h(k_{x},z_{y})$ of a $[2,3]$-triangulene crystal

$h(k_{x},z_{y})=t_{[2,3]}\begin{pmatrix}0&\tau_{ab_{1}}&\tau_{ab_{2}}\\\
\tau_{ab_{1}}^{\prime}&0&0\\\ \tau_{ab_{2}}^{\prime}&0&0\end{pmatrix},$ (16)

with entries

$\displaystyle\tau_{ab_{1,2}}$
$\displaystyle=1+2z_{y}\cos(R_{1_{x}}k_{x}\pm\theta),$ (17)
$\displaystyle\tau_{ab_{1,2}}^{\prime}$
$\displaystyle=1+2z_{y}^{-1}\cos(R_{1_{x}}k_{x}\pm\theta),$

and $R_{1_{x}}=(n_{a}+n_{b}+1)a/2$. The allowed $\\{z_{y_{n}}\\}$ solutions
are obtained by requiring that the eigensolutions of the Hamiltonian in Eq. 16
are also admissible in the open system. Zigzag open boundary conditions
require breaking in half the dimer unit cells at the edges. The corresponding
boundary equations can be easily identified in the rotated mapping version of
Fig. 3(d), where the number of intracell hopping terms is minimized. Given
that the chiral pair only localizes in the $[2]$-triangulene-terminated zigzag
edge, the relevant boundary equation imposes that the wave function vanishes
at $b_{+}$ sites outside the ribbon, that is
$\braket{b_{+}}{\epsilon_{\text{edge}}}=0$. This leads to the only admissible
solution $z_{y}=\cos(R_{1_{x}}k_{x})$, resulting in the following analytical
form of the energy dispersion for this pair of edge states

$\displaystyle\epsilon_{\text{edge}}=\pm\sqrt{6}t_{[2,3]}\sqrt{1-\cos^{2}(R_{1_{x}}k_{x})}=$
$\displaystyle=\pm\sqrt{6}t_{[2,3]}\sqrt{1-z_{y}(k_{x})^{2}},$ (18)

that shows how the splitting of the edge states is maximal when the
localization is maximal, $z_{y}=0$. This is at odds with conventional edge
states in zizgzag graphene, on account of the fact that the two edge states of
the $[2,3]$ zigzag ribbon live on the same edge. Therefore, their
hybridization is reduced when their localization decreases, which happens as
$k_{x}$ moves from the zone boundary towards the zone center. Eq. 18 is
consistent with the numerical calculations except around the $\Gamma$-point,
where the edge state becomes more extended and the effect of the other
boundary condition becomes relevant.

##  V $[3,3]$-triangulene ribbons

We present the results for the energy dispersion $[3,3]$-triangulene ribbons
using the tight-binding complete model and the effective reduced model. We
focus our analysis on the localization and topological properties of edge
modes found specifically in zigzag ribbons, while considerations on armchair
edge states are provided in appendix A.

###  V.1 Complete model

We now consider the edge states of $[3,3]$-ribbons. This parent 2D crystal has
four sets of low-energy bands that arise from the four zero modes of the unit
cell. Two of them are isomorphic to graphene bands, the other two are finite-
energy flat-bands, associated to ring states [21, 14, 13], degenerate with the
dispersive bands at the $\Gamma$-point. Therefore, we should expect edge
states associated to the graphene-like bands, as well as edge states
associated to the flat-bands. In the ribbon geometry, we find a total of
$4N_{rt}$ low-energy bands, where $N_{rt}$ is the number of $[3]$-triangulene
dimers in the unit cell, organized in four groups of low-energy bands plotted
in Fig. 4(a,b).

  * -

The first group of flat bands with finite energy, that belong to ring states,
i.e, to confined states in the smallest closed loops made of six
$[3]$-triangulenes in a hexagonal plaquette. The number of states of this type
is $2(N_{rt}-1)$ for zigzag ribbons and $2(N_{rt}-2)$ for armchair ribbons.

  * -

The second group has $2(N_{rt}-2)$ electron-hole symmetric dispersive bands.
Two Dirac cones at $k_{x}=\pm 2\pi/3$ in zigzag ribbons and valley admixing in
armchair ribbons are predicted in these and all centrosymetric triangulenes
where 3-fold symmetry is preserved. At the $\Gamma$-point, the bands of the
second group become becomes degenerate with the ring states, very much like in
the 2D crystal.

  * -

A third group of two (four) edge-state bands for the zigzag (armchair)
geometry, purple highlighted in Fig. 4(a) [4(b)] in the gap of the second
group of dispersive bands. For the zigzag configuration, the pair of bands
have zero-energy and are localized at the $[3]$-triangulene in both boundaries
when $0<|k|<2\pi/3$ [see Fig. 4a(i)]. This interval is complementary to the
one found in the edge states of $[2,2]$-triangulene ribbons and emulates the
behavior of edge states in bearded zigzag graphene [37]. In armchair ribbons,
this pair of dispersive states is also associated to edge states, maximally
localized at the zone boundary in both sides of the ribbon, populating the
outermost dimers (see appendix A).

  * -

The fourth group, highlighted in green in Fig. 4(a) [Fig. 4(b)], consists of
two degenerate pairs of states associated to the band touching points between
the first and second group of bands. These are present in both boundaries and
are compactly localized in one (two) edge dimers for $k_{x}=0$ ($k_{x}=\pi$).
In the vicinity, these states acquire a $k$-dependent $c_{y}$ decaying factor
before loosing its localization properties at the $k_{x}=\pi$ ($k_{x}=0$)
point [see Fig. 4a(ii)].

###  V.2 Reduced model

Individual $[3,3]$-triangulene dimers host $4$ zero modes that form the basis
set $\\{\ket{a_{1}},\ket{a_{2}},\ket{b_{1}},\ket{b_{2}}\\}$ for the internal
space of the effective zigzag ribbon Hamiltonian. This set can be identified
by their sublattice occupation and relative phase $\omega_{1,2}=\pm 2\pi/3$
between eigenstates of the rotation operator $R_{2\pi/3}$ within the same
sublattice. The finite elements of the effective zigzag ribbon Hamiltonian are
given by the intra and intercell matrices following

$\displaystyle\tau_{r,0}$
$\displaystyle=\begin{pmatrix}0&0&\tau_{r,0}^{a_{1}b_{1}}&\tau_{r,0}^{a_{1}b_{2}}\\\
0&0&\tau_{r,0}^{a_{2}b_{1}}&\tau_{r,0}^{a_{2}b_{2}}\\\
\tau_{r,0}^{a_{1}b_{1}\dagger}&\tau_{r,0}^{a_{2}b_{1}\dagger}&0&0\\\
\tau_{r,0}^{a_{1}b_{2}\dagger}&\tau_{r,0}^{a_{2}b_{2}\dagger}&0&0\end{pmatrix},$
(19) $\displaystyle\tau_{r,1}$ $\displaystyle=\begin{pmatrix}0&0&0&0\\\
0&0&0&0\\\ \tau_{r,1}^{a_{1}b_{1}}&\tau_{r,1}^{a_{2}b_{1}}&0&0\\\
\tau_{r,1}^{a_{1}b_{2}}&\tau_{r,1}^{a_{2}b_{2}}&0&0\end{pmatrix},$

whose entries are given by

$\displaystyle\tau_{r,0}^{a_{1,2}b_{1,2}}$
$\displaystyle=t_{[3,3]},\quad\tau_{r,1}^{a_{1}b_{1}}=\tau_{r,1}^{a_{2}b_{2}}=\left(1+e^{-ik_{x}}\right)t_{[3,3]},$
(20) $\displaystyle\tau_{r,1}^{a_{1,2}b_{2,1}}$ $\displaystyle=\left(e^{\pm
i\theta}+e^{-ik_{x}}e^{\mp i\theta}\right)t_{[3,3]},$

where $t_{[3,3]}=\frac{2t_{3}}{11}$. Very much like in the $[2,3]$ case, we
can map the minimal model of the zigzag ribbon unit cell into a real space
lattice, in this case a Creutz ladder, with $e^{\pm i\theta}$ Peierls factors
in the intercell hopping terms amounting to an effective magnetic flux in
alternating plaquettes [see Fig. 4(c)]. We perform a basis rotation using
bonding/anti-bonding linear combination between the two pairs of zero modes
localized in the same sublattice. This will break $C_{3}$ symmetry between the
eigenfunctions in the new basis set
$\\{\ket{a_{+}},\ket{a_{-}},\ket{b_{+}},\ket{b_{-}}\\}$ and the finite
$\tau_{r,\sigma}^{ab}$ matrices now follow

Figure 4: Model for $[3,3]$-triangulene ribbons with $N_{rt}=20$ dimers in the
unit cell. The numerical results of the energy dispersion for the (a) zigzag
and (b) armchair configurations show the presence of two types of midgap edge
states. In the former, states highlighted in purple [green] localize in a (i)
single triangulene [(ii) triangulene dimer] at the top and bottom boundaries.
The zigzag ribbon unit cell effectively maps into a (c) Creutz ladder in the
$C_{3}$ symmetric basis. The number of intracell hopping terms is reduced if
one considers (d) the bonding/anti-bonding combination basis, whose unit cells
at the high-symmetry points are depicted in (e).

$\displaystyle\tau_{r,0}^{a^{+}b^{+}}$
$\displaystyle=2t_{[3,3]};\quad\tau_{r,1}^{a^{\pm}b^{\mp}}=\mp\left(1-e^{-ik_{x}}\right)\sqrt{3}it_{[3]}/2;$
(21) $\displaystyle\tau_{r,1}^{a^{-}b^{-}}$
$\displaystyle=3\tau_{r,1}^{a^{+}b^{+}}=\left(1+e^{-ik_{x}}\right)t_{[3,3]}/2.$

In this new basis, the number of intracell hopping terms is reduced to a
single $\tau_{r,0}^{a^{+}b^{+}}$ and the ribbon unit cell effectively maps
into a SSH chain composed of sites belonging to bonding states
$\ket{a(b)_{+}}$ [see Fig. 4(d)]. Additionally, sites of anti-bonding states
$\ket{a(b)_{-}}$ are connected to the main chain through closed loops, whose
respective Peierls phase factors sum up to an effective zero transverse
magnetic flux.

At this stage, the presence of a flat-band and the emergence of finite-energy
edge states at a band touching point can be easily understood. When $k_{x}=0$,
anti-bonding states become fully dimerized and contribute with $N_{rt}$
energies $\epsilon(k_{x}=0)=\pm 3t_{[3,3]}$ to the flat-band. Meanwhile, the
bonding states participating in the open SSH chain with hopping terms
$\tau_{r,0}^{a^{+}b^{+}}$ and $\tau_{r,1}^{a^{+}b^{+}}$ follow the generic
form of the energy dispersion

$\epsilon(k_{x}=0,k_{n})=\\\
=\pm\sqrt{(\tau_{r,0}^{a^{+}b^{+}})^{2}+(\tau_{r,1}^{a^{+}b^{+}})^{2}+2\tau_{r,0}^{a^{+}b^{+}}\tau_{r,1}^{a^{+}b^{+}}\cos{(k_{n})}}\\\
=t_{[3,3]}\sqrt{5+4\cos(k_{n})}.$ (22)

There are $N_{rt}-1$ complete dimer unit cells in the chain so that
$k_{n}=\frac{\pi n}{N_{rt}}$ with $n=1,\dots,N_{rt}-1$. In this limit,
$\tau_{r,0}^{a^{+}b^{+}}>\tau_{r,1}^{a^{+}b^{+}}$ and the open chain is in the
topological regime, contributing to the missing chiral pair that completes the
set of eigenstates for $k_{x}=0$. This pair localizes in the edge sites mapped
as the bonding combination state of the two $C_{3}$ symmetric zero modes at
the top and bottom edge $[3]$-triangulenes. When $k_{x}=\pm 2\pi/3$, the
hopping terms $|\tau_{r,1}^{a^{-}b^{-}}|=|\tau_{r,1}^{a^{\pm}b^{\mp}}|$ and
the Peierls factors in the closed loops sum up to $\pi$, resulting in a
topological transition where the chiral edge state pair pick up a finite
opposite dispersion. At the zone boundary, $\tau_{r,1}^{a^{+}b^{+}}$ vanishes
and all bulk states are localized in two triangulene dimers contributing with
energies $\epsilon(k_{x}=\pi)=\\{\pm 3t_{[3,3]},\pm t_{[3,3]}\\}$ to the flat
and dispersive bands, respectively. There are $N_{rt}-1$ of such 4-site
clusters which effectively removes from each band one state of the accounted
at $k_{x}=0$ [see Fig. 4(e)]. The missing four states stem from the additional
dimerization at both boundaries responsible for the twice degenerate finite-
energy edge states with
$\epsilon_{\text{edge}}(k_{x}=\pi)=\pm\sqrt{3}t_{[3,3]}$.

We recover the analytical form of the finite-energy edge states assuming
solutions of the form
$\ket{\epsilon_{\text{edge}}}=e^{ik_{x}}z_{y}^{j}\ket{u(z)}$ in the rotated
basis set $\\{\ket{a_{+}},\ket{a_{-}},\ket{b_{+}},\ket{b_{-}}\\}$ for the
following bulk Hamiltonian $h(k_{x},z_{y})$ of a $[3,3]$-triangulene crystal

$h(k_{x},z_{y})=t_{[3,3]}\begin{pmatrix}0&0&\tau_{a_{+}b_{+}}&\tau_{a_{+}b_{-}}\\\
0&0&\tau_{a_{-}b_{+}}&\tau_{a_{-}b_{-}}\\\
\tau_{a_{+}b_{+}}^{\prime}&\tau_{a_{-}b_{+}}^{\prime}&0&0\\\
\tau_{a_{+}b_{-}}^{\prime}&\tau_{a_{-}b_{-}}^{\prime}&0&0\end{pmatrix},$ (23)

with entries

$\displaystyle\tau_{a_{-}b_{-}}$
$\displaystyle=2+\cos(R_{1_{x}}k_{x})z_{y}=3\tau_{a_{+}b_{+}}-2,$ (24)
$\displaystyle\tau_{a_{\pm}b_{\mp}}$
$\displaystyle=\pm\sqrt{3}\sin(R_{1_{x}}k_{x})z_{y},$
$\displaystyle\tau_{a_{-}b_{-}}^{\prime}$
$\displaystyle=2+\cos(R_{1_{x}}k_{x})z_{y}^{-1}=3\tau_{a_{+}b_{+}}^{\prime}-2,$
$\displaystyle\tau_{a_{\pm}b_{\mp}}^{\prime}$
$\displaystyle=\pm\sqrt{3}\sin(R_{1_{x}}k_{x})z_{y}^{-1}.$

The allowed $\\{z_{y_{n}}\\}$ solutions are obtained by requiring that the
eigensolutions of the Hamiltonian in Eq. 23 also satisfy the open boundary
equations. In the rotated mapping version of Fig. 4(d), one can easily
identify that the open conditions impose that the wave function vanishes at
$b_{+}$ sites for the 0th dimer and at $a_{+}$ sites for the $(N_{rt}+1)$
dimer, that is $\braket{a/b_{+}}{\epsilon_{\text{edge}}}=0$. This leads to two
admissible solutions $z_{y}=\\{\cos(R_{1_{x}}k_{x}),\sec(R_{1_{x}}k_{x})\\}$
that recover the decaying behavior from the top and bottom edge. The
analytical form of the energy dispersion for the finite-energy edge states in
$[3,3]$-triangulene ribbons follows

$\epsilon_{\text{edge}}=\pm\sqrt{3}t_{[3,3]}\sqrt{2-\cos(2R_{1_{x}}k_{x})}.$
(25)

All analytical results provided in the former discussion are in agreement with
the numerical calculations of zigzag ribbons considering the reduced model but
not with the complete model. In the latter, we notice an overall broadening of
the narrow bands around zero-energy and a small dispersion of the flat-bands
in the vicinity of $\Gamma$, implying a small influence of higher/lower energy
bands.

###  V.3 Zak phase analysis

For $[3,3]$-triangulene zigzag ribbons, the loop around the BZ in the
direction perpendicular to the edge has at most four singularity points of
degeneracy: two between the dispersive bands at $\\{K,K^{\prime}\\}$ and two
between the dispersive and flat-bands at the $\Gamma$-point. These degeneracy
points were disregarded in the calculation of the Zak phase for each band (see
methods in appendix B). The pair of zero-energy states in the midgap is
predicted by the sum of the phases for the two lowest bands
$\mathcal{Z}(k_{x})=\mathcal{Z}_{1}+\mathcal{Z}_{2}$. Results show a
topological non-trivial phase for $k_{x}\in[-2\pi/3,2\pi/3]$ and a trivial
phase otherwise.

The topology of the lowest and highest bandgaps were also evaluated. In the
former, the total phase is given by the Zak phase of the lowest band
$\mathcal{Z}=\mathcal{Z}_{1}$ and in the latter one should account for the
three bands below the gap $\mathcal{Z}=\sum_{j=1}^{3}\mathcal{Z}_{j}$. Both
phases are found to be always $\pi$-quantized, thus predicting the two pairs
of dispersive edge states.

##  VI $[4,4]$-triangulene ribbons

In this section we introduce the results for the energy dispersion of
$[4,4]$-triangulene ribbons using the tight-binding complete model and the
effective reduced model on the basis of the zero modes. Unlike the three cases
discussed so far, $[4,4]$-triangulene crystals have a gap at $E=0$. We focus
our analysis on the localization of the edge modes found in both
configurations and on the topological properties of the bandgaps in the band
structure of zigzag ribbons, whose characterization can be equally reproduced
in armchair ribbons.

Figure 5: Model for $[4,4]$-triangulene ribbons with $N_{rt}=20$ dimers in the
unit cell. The numerical results of the energy dispersion for the (a) zigzag
and (b) armchair configurations show the presence of two types of midgap edge
states, the ones highlighted in purple [green] have dispersionless
[dispersive] zero [finite] edge states and localize in (i) a single zigzag
edge triangulene or (ii) armchair edge dimer with small weights in adjacent
triangulenes. The ribbon unit cell effectively maps into a (c) 2-orbital
kagomé ribbon in the corner localized basis. Green highlight shows edge
localization for $\upsilon=0$.

###  VI.1 Complete model

The energy bands of $[4,4]$-triangulene ribbons in Fig. 5(a,b) feature two
copies of the Kagome ribbon bulk bands separated by an energy gap at half-
filling. The group of flat-bands, whose states localize in rings composed of
six $[4]$-triangulenes, is degenerate at the $\Gamma$-point with a group of
dispersive graphene-like bands. At the Dirac points between the latter groups,
twice degenerate dispersive edge states emerge and are maximally localized at
$k=\pi$ in the outermost edge dimers for zigzag boundaries, surviving for
$|k|\in[2\pi/3,\pi]$. Armchair finite-energy edge states have a larger decay
through the bulk even at the zone boundary.

Contrary to the $[3,3]$ case, the two identified band touching points between
the $2(N_{rt}-1)$ [$2(N_{rt}-2)$] flat-band groups and the $2(N_{rt}-1)$
dispersive band groups at the $\Gamma$-point in the zigzag (armchair)
configuration do not generate finite-energy dispersive edge states. The
missing $2$ ($4$) states are found in the gap between the electron-hole
symmetric bands, with no dispersion, both for zigzag and armchair ribbons. In
the zigzag case, these edge states have $E=0$, whereas in the armchair case
this non-dispersive band features a very small width-independent energy
splitting, likely due to intraedge hybdirization. As can be seen in Fig.
5(i,ii), their $k$-independent wave function is compactly localized in the
majority sublattice at the top and bottom zigzag edge $[4]$-triangulene and in
the armchair edge dimer with small weight in adjacent $[4]$-triangulenes.

###  VI.2 Reduced model

The degeneracy of the zero modes in $[n]$-triangulenes for $n>2$ allow one to
choose a suitable representation for its set. In the particular case of
$[4]$-triangulenes accommodating $n_{a/b}-1=3$ zero modes, these can be
represented as the eigenstates of the 3-fold symmetry operator of Fig.1(b) or
rearranged into three orbitals related by a $2\pi/3$ rotation [14]. We use the
analogue of a tight-binding ring of three sites where each site represents a
state that is essentially populating one corner of the $[4]$-triangulene. The
new corner localized basis $\ket{c_{j}}$ can be constructed from the $C_{3}$
symmetric zero modes $\ket{\omega}$ with $\omega=\\{0,\pm 2\pi/3\\}$ from the
linear combination $\ket{\omega}=\frac{1}{\sqrt{3}}\sum_{j=1}^{3}e^{i\omega
j}\ket{c_{j}}$. The finite elements of the effective zigzag ribbon Hamiltonian
in this basis are given by the intra and intercell contributions

$\tau_{r,0}=\begin{pmatrix}0^{[n_{a}-1]}&\tau_{r,0}^{c}\\\
\tau_{r,0}^{c}{}^{\dagger}&0^{[n_{b}-1]}\end{pmatrix},\quad\tau_{r,1}=\begin{pmatrix}0^{[n_{a}-1]}&0\\\
\tau_{r,1}^{c}{}^{\dagger}&0^{[n_{b}-1]}\end{pmatrix},$ (26)

whose entries follow

$\displaystyle\tau_{r,0}^{c}$
$\displaystyle=t_{[4,4]}\begin{pmatrix}1&\upsilon&\upsilon\\\ \upsilon&0&0\\\
\upsilon&0&0\end{pmatrix},$ (27) $\displaystyle\tau_{r,1}^{c}$
$\displaystyle=t_{[4,4]}\begin{pmatrix}0&\upsilon&\upsilon e^{-ik_{x}}\\\
\upsilon&1&\upsilon(1+e^{-ik_{x}})\\\ \upsilon
e^{-ik_{x}}&\upsilon(1+e^{-ik_{x}})&e^{-ik_{x}}\end{pmatrix},$

where $t_{[4,4]}=0.35t_{3}$ and $\upsilon t_{[4,4]}=0.033t_{3}$ describe first
and second neighboring overlap between adjacent corner states, respectively.
The non-vanishing second neighbor term is a byproduct of the non-zero weights
of the corner localized states at sites close to the other two corners. Next
order terms are smaller than $\upsilon t_{[4,4]}$ and can be discarded.

The zigzag unit cell of the reduced model is depicted in Fig. 5(c). For
$\upsilon=0$, every corner localized state only hybridizes with a single
corner state on the adjacent triangulene. In this dimerized limit, a
dispersionless valence bond solid is formed from the corresponding bonding and
antibonding states. These states are centered at the midpoint between
neighboring corners and are weakly coupled by a non-zero $\upsilon t_{[4,4]}$,
effectively tailoring kagome ribbon lattices. Indeed, the low-energy spectrum
gives a pair of kagome subbands from these bonding and antibonding states.

As expected, the effective hopping terms do not connect corner states within
each [4]-triangulene nor belonging to the same sublattice. Therefore, each
triangle [see Fig. 5(c)] belongs to a different sublattice, preserving the
chiral symmetry of the model. Zigzag edge states obey this symmetry as they
remain zero-energy states with increasing $\upsilon$. In the armchair
configuration and when $\upsilon=0$, two edge states with support in both
sublattices emerge on the same edge. A non-zero $\upsilon$ term couples these
states, leading to energy splitting of order $\upsilon^{2}$. Higher order
terms make the edge state band dispersive with small bandwidth at least of the
order of $\upsilon^{4}$. The limit of $\upsilon\gg 1$ does not lead to
dimerization and therefore there is no gap reopening. For each $k_{x}$ in the
ribbon, $\upsilon=0$ implies a dimerized chain similar to the $t_{1}=0$ limit
of the SSH chain.

###  VI.3 Zak phase analysis

We ought to evaluate the topology of both the insulating midgap and the two
bandgaps between the Dirac points of the $[4,4]$ zigzag ribbon band structure.
We identify six singularity points of degeneracy: four between the dispersive
bands at $\\{K,K^{\prime}\\}$ and two between the dispersive and flat-bands at
the $\Gamma$-point. The overlap terms at these points were dropped from the
numerical calculation of the Zak phase for the corresponding bands.

The pair of zero-energy edge states in the midgap is predicted by the sum of
the phases for the three lowest bands
$\mathcal{Z}=\sum_{j=1}^{3}\mathcal{Z}_{j}$. Results show a topological non-
trivial phase for all $k_{x}$.

The topology of the lowest and highest bandgaps were also evaluated, having
into account the presence of degeneracies at the Dirac cones. For the lowest
bandgap, the total phase is given by the Zak phase of the lowest band
$\mathcal{Z}=\mathcal{Z}_{1}$ and for the highest bandgap one should account
for the five bands below the gap $\mathcal{Z}=\sum_{j=1}^{5}\mathcal{Z}_{j}$.
Both phases are found to be $\pi$ for $|k_{x}|\in[2\pi/3,\pi]$ and trivial
when $k_{x}\in[-2\pi/3,2\pi/3]$, thus effectively predicting the two pairs of
dispersive edge states.

##  VII Conclusions

We have carried out an extensive study on the properties of edge states on a
class of models that generalize the conventional honeycomb tight-binding
Hamiltonian to the case where each site has several degenerate modes,
eigenstates of the $C_{3}$ symmetry operator. This model is relevant to
describe many systems, including triangulene crystals, the main motivation for
our work. We have discussed in detail four classes of triangulene crystals by
examining the band structure of one-dimensional triangulene ribbons. The
construction of simpler one-dimensional effective models with $k$-dependent
hopping parameters, such as the diamond chain for the $[2,3]$-triangulene and
the creutz ladder for the $[3,3]$-triangulene, allowed one to understand the
regimes where these edge states emerge. We have found exotic edge phenomena,
absent in the conventional honeycomb tight-binding model of graphene, deriving
from the non-trivial features of the two-dimensional bands of triangulenes:

  * -

A pair of edge states living in the same $[2]$-triangulene edge of the $[2,3]$
zigzag ribbon: a consequence of the zero mode mismatch between triangulenes
results in edge dimerization for the top edge and an additional zero-energy
edge state in the bottom edge at the zone boundary.

  * -

Edge states occurring at flat-parabolic touching-points and bearded edge
states when degenerate with a bulk state from the Dirac point of $[3,3]$
ribbons, resulting from edge dimerization on both boundaries and the non-
trivial regime of an SSH chain, respectively.

  * -

Zero-energy edge states in the bandgap at half-filling of the $[4,4]$ bands,
that arise as a consequence of unpaired zero-modes in both edges.

  * -

Edge states surviving in armchair ribbons whenever the dimensions of the
triangulene dimer $n_{a}+n_{b}>4$.

We have assessed the topological origin of zigzag edge states using the
single-band Zak phase formalism. Although the bulk-boundary correspondence for
finite-energy dispersive edge states remains an open question, we argue that
the presence of all edge states in $[3,3]$ and $[4,4]$ zigzag ribbons was
successfully linked to non-trivial Zak phase values, including the dispersive
finite-energy edge states appearing between degenerate bands. The same
approach can be applied to armchair edge states for the right choice of unit
cell and $k_{\parallel}$-direction.

Our results reveal the prominence of peculiar properties in $C_{3}$ symmetric
multi-orbital honeycomb lattices. Full concurrence between the complete and
reduced tight-binding versions correlates the appearance of edge states to the
low-energy bands associated with the zero modes of single triangulenes. We
argue that the agreement between the complete and the mininmal model stresses
the crucial role played by the interplay between the $C_{3}$ symmetry of the
high-fold fermions in each site of the honeycomb lattice and the $C_{3}$
symmetry of the lattice. Therefore, we expect that our analys will be relevant
for other multi-orbital models with $C_{n}$-symmetry on other crystalline
arrangements.

###### Acknowledgements.

L.M. acknowledges the financial support from FCT (grant No. PTDC/FIS-
MAC/2085/2020). L.M. and J.F.-R. acknowledge support from the European Union
(Grant FUNLAYERS-101079184). We acknowledge fruitful discussions with João
Henriques and António Costa. J.F.-R. acknowledges financial support from
Generalitat Valenciana (Prometeo2021/017 and MFA/2022/045), from SNF Sinergia
(Grant Pimag, CRSII5-205987) and MICIN-Spain (Grants PID2019-109539GB-C41 and
PRTR-C1y.I1). L.M. and R.G.D. developed their work within the scope of
Portuguese Institute for Nanostructures, Nanomodelling and Nanofabrication
(i3N) Projects No. UIDB/50025/2020, No. UIDP/50025/2020, and No.
LA/P/0037/2020, financed by national funds through the Fundação para a Ciência
e Tecnologia (FCT) and the Ministério da Educação e Ciência (MEC) of Portugal.

##  Appendix A Edge states in armchair terminated ribbons

Figure 6: Profile of the wavefunction for the dispersive edge states found in
(a) $[2,3]$ armchair ribbons at (i) $k_{y}=\pi$ and (ii) at $k_{y}=0$; (b)
$[3,3]$ armchair ribbons at (i), (iii) $k_{y}=\pi$ and at (ii), (iv)
$k_{y}=0$. We show the probability distribution of the lowest energy edge
state in the unit cell of the armchair ribbon periodic in the $y$-direction
with $N_{rt}=21$ and $N_{rt}=41$ dimers for a(i).

The observed valley mixing on armchair graphene largely influences the
$\pi$-electron states near the Fermi level and accounts for the absence of
edge states [55]. Herein, the isospin of the fermions within each individual
triangulene for $n_{a/b}>2$ recovers edge localization, despite featuring no
intra-edge sublattice imbalance (edges terminate with A-B triangulene dimers).
In the following appendix, we briefly discuss the identified midgap edge
states in $[2,3]$ and $[3,3]$ triangulene armchair ribbons.

In Fig. 6(a) we represent the lower energy state of the electron-hole pair of
dispersive edges states found in the gap between the dispersive and the flat
subbands of the [2,3] armchair ribbon. Since inter-edge sublattice imbalance
is absent in the armchair configuration, the edge state can be found in both
sides of the ribbon and accounts for the double electron-hole edge state pair
at the zone boundary. As one can see in Fig. 6a(i) for an armchair ribbon with
$N_{rt}=21$ and a larger ribbon with $N_{rt}=41$ dimers, this state is mostly
localized in the outermost $[2]$-triangulene and two adjacent
$[3]$-triangulenes. Similar to its zigzag counterpart, the localization length
of this state increases as one approaches the $\Gamma$-point [see Fig.
6a(ii)].

The probability distribution of the edge states found in the three gaps of
$[3,3]$ armchair ribbons is illustrated in Fig. 6(b). We find a total of two
doubly degenerate electron-hole pairs of dispersive edge states. The first
pair is localized in the midgap of the bulk dispersive subbands while the
second pair can be found in the gap between the dispersive and flat subbands.
These modes are present for all momenta, are maximally localized at the zone
boundary in both sides of the ribbon [see Fig. 6b(i,iii)] and the exponential
decay through the bulk increases closer to the $\Gamma$-point [see Fig.
6b(ii,vi)]. The dispersive edge states in the central gap populate the
outermost dimers while the lowest energy pair localizes in the second to the
last row of dimers, with a gradual decrease towards the bulk.

##  Appendix B Numerical calculation of the total Zak phase for [3,3] and
[4,4] ribbons

Figure 7: Numerical calculation of the total Zak phase for the identified
bandgaps on the band structure of (a) $[3,3]$ and (b) $[4,4]$ triangulene
ribbons. Inset tables display the $\mathcal{Z}_{n}$ for each band numbered in
increasing energy order for the relevant intervals of momentum $k_{x}$. The
phase of the highest band is not shown here as is irrelevant for the
calculation of $\mathcal{Z}$.

We have used the single band formalism to predict the topological phases of
all bandgaps of $[3,3]$ and $[4,4]$ triangulene ribbons. This often poses some
challenges in the presence of degenerate bands if no considerations are taken
on the gauge choice. We have treated the presence of singular points of
degeneracy for each $k_{\parallel}$ as follows: we make use of a smooth gauge
choice where phases are roughly uniformly distributed throughout the loop.
This distribution becomes less uniform close to topological transition points.
The gauge transformation is possible by finding, for each band, the non-zero
component of $\ket{u(k)}$ that has the smallest change throughout the
$k$-scanning in the Wilson loop and choosing the gauge where this component is
real. The advantage of using a smooth gauge is that it allows for the
calculation of the Zak phase when there is finite number of band crossings.
These points can be dropped out numerically, but one must take into account
the change in the ordering of the bands at these crossing points using for
example a largest overlap criterium to correctly order the bands. For both the
[3,3] and [4,4] ribbons, the bands considered for the estimation of the total
Zak phase have at most three singularity points of degeneracy: two between the
dispersive bands at $k_{\bot}=2\pi/3$ and one between dispersive and a flat-
band at $k_{\bot}=0$. Disregarding the overlaps elements in these points will
not drastically affect the calculation of the total phase.

We also have to note that the value of Zak phase of each band may have a
difference of $\pm\pi$ for different choices of the bulk unit cell. We make a
choice of unit cell dimer such that it allows the full reconstruction of the
lattice including its specific zigzag edge.

The total Zak phases for the three identified bandgaps of the $[3,3]$ and
$[4,4]$ ribbons are depicted in Fig. 7(a,b), respectively. These were
determined by the sum of the Zak phases for each band below the corresponding
bandgap, whose values can be found in the inset tables in Fig. 7. Specifically
for the lowest bandgap, the total Zak phase is equivalent to $\mathcal{Z}_{1}$
of the lowest energy band (which is the only band found below this gap). We
can see that the non central bandgaps of [3,3] ribbons are always $\pi$ while
the midgap takes a non-trivial topological phase when
$k_{x}\in[-2\pi/3,2\pi/3]$, agreeing with the localization of the states found
in these gaps. On the contrary, only the midgap of $[4,4]$ ribbons is
$\mathcal{Z}=\pi$ independent of $k_{x}$. The remaining gaps have non-trivial
topology in the complementary interval $|k_{x}|\in[2\pi/3,\pi]$, which
effectively predicts the pair of dispersive edge states.

Singular points of degeneracy are present when $k_{x}=\\{0,\pm 2\pi/3\\}$. In
[3,3] ribbons, $\mathcal{Z}_{n}(0)$ possesses two band crossings at $k_{y}=0$
between bands $1,2$ and $3,4$ and $\mathcal{Z}_{n}(\pm 2\pi/3)$ with $n=2,3$
has two degenerate points at $k_{y}=\pm\pi$. In the case of [4,4] ribbons, we
find two degenerate points between bands $2,3$ and $4,5$ that will affect the
overlap elements for $k_{y}=0$ and we also find four band crossings for
$k_{x}=\pm 2\pi/3$ between the two highest and lowest bands at $k_{y}=\pm\pi$.
We have disregarded the overlap elements of all these degeneracy points which
causes the total phase $\mathcal{Z}(0,\pm 2\pi/3)$ to have a minimal deviation
from the $0\,\wedge\,\pi\mod 2\pi$ quantization expected for inversion
symmetric models. This error gets mitigated if we consider a smaller $\Delta
k_{y}$.

## References

  * Novoselov  _et al._ [2005] K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang, M. I. Katsnelson, I. V. Grigorieva, S. V. Dubonos, and A. A. Firsov, Nature 438, 197 (2005).
  * Yarmohammadi and Zareyan [2016] M. Yarmohammadi and M. Zareyan, Chinese Physics B 25, 068105 (2016).
  * Katsnelson  _et al._ [2006] M. I. Katsnelson, K. S. Novoselov, and A. K. Geim, Nature Physics 2, 620 (2006).
  * Beenakker [2008] C. Beenakker, Reviews of Modern Physics 80, 1337 (2008).
  * Katsnelson [2006] M. Katsnelson, The European Physical Journal B 51, 157 (2006).
  * Nair  _et al._ [2008] R. R. Nair, P. Blake, A. N. Grigorenko, K. S. Novoselov, T. J. Booth, T. Stauber, N. M. Peres, and A. K. Geim, Science 320, 1308 (2008).
  * Brey and Fertig [2006] L. Brey and H. A. Fertig, Physical Review B 73, 235411 (2006).
  * Xu  _et al._ [2015] S.-Y. Xu, I. Belopolski, N. Alidoust, M. Neupane, G. Bian, C. Zhang, R. Sankar, G. Chang, Z. Yuan, C.-C. Lee, S.-M. Huang, H. Zheng, J. Ma, D. S. Sanchez, B. Wang, A. Bansil, F. Chou, P. P. Shibayev, H. Lin, S. Jia, and M. Z. Hasan, Science 349, 613 (2015).
  * Rhim  _et al._ [2020] J.-W. Rhim, K. Kim, and B.-J. Yang, Nature 584, 59 (2020).
  * Oh  _et al._ [2022] C.-g. Oh, D. Cho, S. Y. Park, and J.-W. Rhim, Communications Physics 5, 320 (2022).
  * Bradlyn  _et al._ [2016] B. Bradlyn, J. Cano, Z. Wang, M. Vergniory, C. Felser, R. J. Cava, and B. A. Bernevig, Science 353, aaf5037 (2016).
  * Hasan  _et al._ [2021] M. Z. Hasan, G. Chang, I. Belopolski, G. Bian, S.-Y. Xu, and J.-X. Yin, Nature Reviews Materials 6, 784 (2021).
  * Ortiz  _et al._ [2022] R. Ortiz, G. Catarina, and J. Fernández-Rossier, 2D Materials 10, 015015 (2022).
  * Tamaki  _et al._ [2020] G. Tamaki, T. Kawakami, and M. Koshino, Physical Review B 101, 205311 (2020).
  * Fernández-Rossier and Palacios [2007] J. Fernández-Rossier and J. J. Palacios, Physical Review Letters 99, 177204 (2007).
  * Ortiz  _et al._ [2019] R. Ortiz, R. Á. Boto, N. García-Martínez, J. C. Sancho-García, M. Melle-Franco, and J. Fernández-Rossier, Nano Letters 19, 5991 (2019).
  * Tran  _et al._ [2017] V.-T. Tran, J. Saint-Martin, P. Dollfus, and S. Volz, AIP Advances 7, 075212 (2017).
  * Mizoguchi  _et al._ [2021] T. Mizoguchi, Y. Kuno, and Y. Hatsugai, Physical Review B 104, 035161 (2021).
  * Weeks and Franz [2010] C. Weeks and M. Franz, Physical Review B 82, 085310 (2010).
  * Tsai  _et al._ [2015] W.-F. Tsai, C. Fang, H. Yao, and J. Hu, New Journal of Physics 17, 055016 (2015).
  * Wu  _et al._ [2007] C. Wu, D. Bergman, L. Balents, and S. Das Sarma, Physical Review Letters 99, 070401 (2007).
  * Congjun and Sarma [2008] W. Congjun and S. D. Sarma, Physical Review B 77, 235107 (2008).
  * Fang  _et al._ [2024] T. Fang, T. Zhang, T. Hu, and Z. Wang, Nano Letters 24, 3059 (2024).
  * Delgado  _et al._ [2023] A. Delgado, C. Dusold, J. Jiang, A. Cronin, S. G. Louie, and F. R. Fischer, arXiv: 2301.06171 (2023).
  * Mishra  _et al._ [2020] S. Mishra, D. Beyer, K. Eimre, R. Ortiz, J. Fernández-Rossier, R. Berger, O. Gröning, C. A. Pignedoli, R. Fasel, X. Feng, and P. Ruffieux, Angewandte Chemie International Edition 59, 12041 (2020).
  * Jacob  _et al._ [2021] D. Jacob, R. Ortiz, and J. Fernández-Rossier, Physical Review B 104, 075404 (2021).
  * Henriques and Fernández-Rossier [2023] J. C. G. Henriques and J. Fernández-Rossier, Physical Review B 108, 155423 (2023).
  * Mishra  _et al._ [2021] S. Mishra, G. Catarina, F. Wu, R. Ortiz, D. Jacob, K. Eimre, J. Ma, C. A. Pignedoli, X. Feng, P. Ruffieux, J. Fernández-Rossier, and R. Fasel, Nature 598, 287 (2021).
  * Catarina and Fernández-Rossier [2022] G. Catarina and J. Fernández-Rossier, Physical Review B 105, L081116 (2022).
  * Martínez-Carracedo  _et al._ [2023] G. Martínez-Carracedo, L. Oroszlány, A. García-Fuente, L. Szunyogh, and J. Ferrer, Physical Review B 107, 035432 (2023).
  * Saleem  _et al._ [2024] Y. Saleem, T. Steenbock, E. R. J. Alhadi, W. Pasek, G. Bester, and P. Potasz, Nano Letters 24, 7417 (2024).
  * Zhou and Liu [2020] Y. Zhou and F. Liu, Nano Letters 21, 230 (2020).
  * Sethi  _et al._ [2021] G. Sethi, Y. Zhou, L. Zhu, L. Yang, and F. Liu, Physical Review Letters 126, 196403 (2021).
  * Catarina  _et al._ [2023] G. Catarina, J. C. G. Henriques, A. Molina-Sánchez, A. T. Costa, and J. Fernández-Rossier, Physical Review Research 5, 043226 (2023).
  * Nakada  _et al._ [1996] K. Nakada, M. Fujita, G. Dresselhaus, and M. S. Dresselhaus, Physical Review B 54, 17954 (1996).
  * Sierda  _et al._ [2023] E. Sierda, X. Huang, D. I. Badrtdinov, B. Kiraly, E. J. Knol, G. C. Groenenboom, M. I. Katsnelson, M. Rösner, D. Wegner, and A. A. Khajetoorians, Science 380, 1048 (2023).
  * Plotnik  _et al._ [2013] Y. Plotnik, M. C. Rechtsman, D. Song, M. Heinrich, J. M. Zeuner, S. Nolte, Y. Lumer, N. Malkova, J. Xu, A. Szameit, Z. Chen, and M. Segev, Nature Materials 13, 57 (2013).
  * Jacqmin  _et al._ [2014] T. Jacqmin, I. Carusotto, I. Sagnes, M. Abbarchi, D. Solnyshkov, G. Malpuech, E. Galopin, A. Lemaître, J. Bloch, and A. Amo, Physical Review Letters 112, 116402 (2014).
  * Gomes  _et al._ [2012] K. K. Gomes, W. Mar, W. Ko, F. Guinea, and H. C. Manoharan, Nature 483, 306 (2012).
  * Chen  _et al._ [2023] A. Chen, H. Brand, T. Helbig, T. Hofmann, S. Imhof, A. Fritzsche, T. Kießling, A. Stegmaier, L. K. Upreti, T. Neupert, T. Bzdušek, M. Greiter, R. Thomale, and I. Boettcher, Nature Communications 14, 622 (2023).
  * Bellec  _et al._ [2014] M. Bellec, U. Kuhl, G. Montambaux, and F. Mortessagne, New Journal of Physics 16, 113023 (2014).
  * Kohmoto and Hasegawa [2007] M. Kohmoto and Y. Hasegawa, Physical Review B 76, 205402 (2007).
  * Laude  _et al._ [2009] V. Laude, Y. Achaoui, S. Benchabane, and A. Khelif, Physical Review B 80, 092301 (2009).
  * Alase  _et al._ [2017] A. Alase, E. Cobanera, G. Ortiz, and L. Viola, Physical Review B 96, 195133 (2017).
  * Cobanera  _et al._ [2018] E. Cobanera, A. Alase, G. Ortiz, and L. Viola, Physical Review B 98, 245423 (2018).
  * Kunst  _et al._ [2019] F. K. Kunst, G. van Miert, and E. J. Bergholtz, Physical Review B 99, 085427 (2019).
  * Zak [1989] J. Zak, Physical Review Letters 62, 2747 (1989).
  * Delplace  _et al._ [2011] P. Delplace, D. Ullmo, and G. Montambaux, Physical Review B 84, 195452 (2011).
  * Xiao  _et al._ [2015] M. Xiao, G. Ma, Z. Yang, P. Sheng, Z. Q. Zhang, and C. T. Chan, Nature Physics 11, 240 (2015).
  * Wang  _et al._ [2016] L. Wang, R.-Y. Zhang, M. Xiao, D. Han, C. T. Chan, and W. Wen, New Journal of Physics 18, 103029 (2016).
  * Milićević  _et al._ [2017] M. Milićević, T. Ozawa, G. Montambaux, I. Carusotto, E. Galopin, A. Lemaître, L. Le Gratiet, I. Sagnes, J. Bloch, and A. Amo, Physical Review Letters 118, 107403 (2017).
  * Sarandy and Lidar [2006] M. S. Sarandy and D. A. Lidar, Physical Review A 73, 062101 (2006).
  * Sugawa  _et al._ [2021] S. Sugawa, F. Salces-Carcoba, Y. Yue, A. Putra, and I. B. Spielman, npj Quantum Information 7, 144 (2021).
  * Lado  _et al._ [2015] J. Lado, N. García-Martínez, and J. Fernández-Rossier, Synthetic Metals 210, 56 (2015).
  * Wakabayashi  _et al._ [2009] K. Wakabayashi, Y. Takane, M. Yamamoto, and M. Sigrist, Carbon 47, 124 (2009).