# Topology in a Su–Schrieffer–Heeger plasmonic crystal

D. A. Miranda  0000-0002-9784-3340 Centro de Física (CF-UM-UP) and
Departamento de Física, Universidade do Minho, P-4710-057 Braga, Portugal  Y.
V. Bludov  0000-0001-9648-1459 Centro de Física (CF-UM-UP) and Departamento de
Física, Universidade do Minho, P-4710-057 Braga, Portugal  N. Asger Mortensen
0000-0001-7936-6264 POLIMA—Center for Polariton-driven Light–Matter
Interactions, University of Southern Denmark, Campusvej 55, DK-5230 Odense M,
Denmark  Danish Institute for Advanced Study, University of Southern Denmark,
Campusvej 55, DK-5230 Odense M, Denmark  Nuno M. R. Peres  0000-0002-7928-8005
Centro de Física (CF-UM-UP) and Departamento de Física, Universidade do Minho,
P-4710-057 Braga, Portugal  POLIMA—Center for Polariton-driven Light–Matter
Interactions, University of Southern Denmark, Campusvej 55, DK-5230 Odense M,
Denmark  International Iberian Nanotechnology Laboratory (INL), Av Mestre José
Veiga, 4715-330 Braga, Portugal

###### Abstract

In this paper we study the topology of the bands of a polaritonic crystal
composed of graphene and of a metallic grating. We first derive a
Kronig–Penney type of equation for the polaritonic bands as function of the
Bloch wavevector and discuss the propagation of the surface plasmon polaritons
on the polaritonic crystal using a transfer-matrix approach. Secondly, we
reformulate the problem as a tight-binding model that resembles the
Su–Schrieffer-Heeger (SSH) Hamiltonian, one difference being that the hopping
amplitudes are, in this case, energy dependent. In possession of the tight-
binding equations it is a simple task to determine the topology of the bands.
Similarly to the SSH model, we show that there is a tunable parameter that
induces topological phase transitions from trivial to non-trivial. In our
case, it is the distance $d$ between the graphene sheet and the metallic
grating. We note that $d$ is a parameter that can be easily tuned
experimentally simply by controlling the thickness of the spacer between the
grating and the graphene sheet. It is then experimentally feasible to engineer
devices with the required topological properties.

##  I Introduction

It is well known [1, 2, 3] that continuous metallic interfaces and metallic
nanoparticles sustain collective charge excitations, dubbed plasmon
polaritons. In the case of nanoparticles, this type of excitations are
localized at their external boundary [4, 5]. For metallic films, the charge
excitations are also localized at the metal vacuum interface but, contrary to
metallic nanoparticles, they are free to travel along the interface. This type
of excitation, half matter, half light, is dubbed surface plasmon polariton.
In noble metals, such as silver and gold, plasmon polaritons are sustained in
the visible and near-IR spectral range. For lower frequencies the plasmon
polaritons become severely damped due to Ohmic losses [6]. Often, devices
based on plasmon polaritons explore the highly confined electric fields for
applications such as sensing, optoelectronic applications and probing [7, 8,
9].

The severe aforementioned Ohmic losses call for other materials that can
sustain undamped plasmon polaritons. It has been shown that graphene fills
this gap by exhibiting both propagating and localized surface plasmons in the
mid-infrared (IR) spectrum, that can be tuned by gating or doping [10, 11, 12,
13]. When encapsulated in hexagonal boron nitride, surface plasmon polaritons
in graphene show higher propagation lengths and stronger electromagnetic field
confinement [14, 15, 16]. This property made graphene suitable for sensing and
optoelectronic applications in the mid-IR spectral range [17, 18, 19].

Figure 1: Representation of a plasmonic crystal periodic in the $x$ direction.
A graphene sheet is separated from a metal grating by a dielectric of
thickness $d$ and dielectric constant $\epsilon_{2}$. Bellow the graphene
sheet, another semi-infinite dielectric, with dielectric constant
$\epsilon_{1}$, acts as substrate. The metal grating can be divided in two
parts: a metallic one with length $l_{g}$, named gated region, followed by a
stripe where the metal is absent, and with length $l_{u}$, named ungated
region. The unit cell of the plasmonic crystal has length $l=l_{g}+l_{u}$.

It has been shown that plasmons in noble metal based structures can, in
certain circumstances, sustain topological plasmons, as is the case of
structured silver films where topological plasmonic vortices have been
observed [20]. Topology in structured materials is actually a quite general
feature of matter, exhibited in systems as diverse as photonic and phononic
systems [21, 22, 23, 24, 25, 26, 27], and cold atom lattices [28, 29], and
being rather ubiquitous in condensed matter systems [30, 31, 32, 33], where
the crystalline structure plays an essential role in determining the presence
or absence of topological energy bands.

Graphene itself is not a topological material (although it was shown that a
hexagonal array of metallic nanoparticles can be [34]). Therefore, the
question arises if structured graphene can exhibit topological properties. In
particular, if plasmon polaritons in this material can present a polaritonic
band structure with topological bands. A preliminary and fully numerical study
suggests that this is indeed the case [35]. Graphene can be structured in
several different ways. Two possibilities are patterning a periodic array of
graphene ribbons [36, 37, 38, 39] or coupling a continuous graphene sheet with
a metallic grating [40, 41, 42, 43]. The advantage of the second approach is
twofold. First, no additional scattering channels are introduced in the system
as opposed to the patterned case, which leads to additional scattering of the
surface plasmons at the edges of the ribbons. Secondly, the micro-fabrication
process is much facilitated. We, therefore, focus our attention on the second
system.

It is well known [41, 42, 44] that graphene surface plasmons behave
considerably different when in the vicinity of a metal substrate, in
comparison to usual graphene plasmons when graphene is deposited on a
dielectric substrate alone. In the latter case the plasmon are dubbed graphene
plasmons (or ungated plasmons), whereas in the first case they are dubbed
screened graphene plasmons (occasionally also refereed to as acoustic
plasmons, gated plasmons, or even image polaritons [45]). This distinction is
made because, in the presence of a metal, the electric field is screened,
leading to a softening of the plasmon dispersion, that becomes linear with the
wave vector [44], as opposed to the usual graphene plasmons which disperse
with the square root of the wave vector [44]. Thus, graphene in the vicinity
of a metallic grating is a hybrid system with alternating regions of gated and
ungated plasmons. We call this system a plasmonic crystal. The unit cell is
then composed of a region of gated plasmons (acoustic graphene plasmons with
$\omega\propto q$) followed by a region of ungated ones (graphene plasmons
with $\omega\propto\sqrt{q}$). The fact that the dispersion of the two types
of plasmons is different, as noted before, has the consequence that both types
of plasmons have different wave vector for the plasmon polaritons of a give
frequency. From the point of view of a theoretical description, the envisioned
system is much like a Kronig–Penney model [46] for electrons moving in
periodic arrangement of electrostatic potentials. Indeed, when we describe the
system in terms of electrostatic potentials and electric current densities,
the transcendental equation leading to the polaritonic bands is formally
identical to that of the Kronig–Penney model. However, the Kronig–Penney model
formalism is not simple to analyze due to its infinite number of bands. Also,
from the point of view of topology, the original formulation of the model
lacks the form of an electronic tight-binding Hamiltonian, for which well know
tools exist for describing the system’s topology. One such electronic model is
the one-dimensional Su–Schrieffer–Heeger (SSH) model [47, 33]. Electrons
described by the SSH model are known to show topological behavior when certain
conditions are met for the hopping parameters. The question now arises if it
is possible to describe our system, a periodic array of gated and ungated
graphene plasmons, with the SSH model formalism. Our analysis below shows that
this is indeed the case. With this mapping in hands, the description of
topology becomes much simpler than in the original Kronig–Penney formulation.

This paper is organized has follows. In Sec. II we introduce the polaritonic
crystal, as realized in Ref. 41 alongside the mathematical formalism.
Thereafter, we derive the plasmonic energy dispersion and the crystal
propagation properties. In Sec. III the SSH tight-binding representation of
the plasmonic crystal is introduced. There, we investigate how the distance
between the metal grating and graphene affects the energy-dependent tight-
binding coefficients and the topology associated with the plasmonic bands.
Additionally, we derive an analytical expression for the particular distance
that results in the energy gaps closing, as well as an analytical
approximation for the bands. Finally, we offer our conclusions and detail
important calculations in the Appendixes.

##  II Plasmons in a graphene-based SSH-like grating

In this system, experimentally studied in Ref. 41, we have a graphene sheet
sandwiched between two dielectrics, characterized by dielectric constants
$\epsilon_{1}$ and $\epsilon_{2}$. On the upper-most dielectric, at a distance
$d$ from the graphene sheet, there is a metal grating. The grating is a one-
dimensional (1D) periodic structure composed of two alternating portions: a
metallic section of length $l_{g}$ (denominated gated region) and an empty
portion of length $l_{u}$(named ungated region), such that the unit cell of
length $l=l_{u}+l_{g}$ can be defined. The presence of a conductor
significantly alters the properties of plasmons on graphene. For instance, the
wavevectors of the gated ($q_{g})$ and ungated $(q_{u})$ regions differ, being
given by [44]

$q_{g}(\omega,d)=\omega\sqrt{\frac{\hbar\epsilon_{2}}{4E_{F}c\alpha d}},$ (1a)
$q_{u}(\omega)=\frac{\hbar(\epsilon_{1}+\epsilon_{2})\omega^{2}}{4E_{F}c\alpha},$
(1b)

respectively, where $E_{F}$ is the Fermi energy of graphene (considered to be
equal in both regions), $\alpha$ is the fine-structure constant, $c$ is the
speed of light in vacuum, and $\omega$ is the plasmon frequency. We note that
Eqs. (1a) and (1b) are analytical approximations to the exact dispersion,
which can be known only numerically. To obtain them, it is necessary to assume
that (i) the frequency regime is such that the conductivity of graphene is
dominated by its Drude conductivity
$\sigma(\omega)=\frac{4iE_{F}}{\hbar(\omega+i\gamma)}\alpha\epsilon_{0}c$,
(ii) absorption is zero ($\gamma=0)$, and (iii) we can approximate the
electrostatic limit, $q_{u/g}\gg\frac{\sqrt{\epsilon_{j}}\omega}{c}$, where
$j=1,2$ is the index labeling the dielectrics. Also, we consider values of $d$
such that $q_{g}d\ll 1$, which has been considered in previous experiments
[48]. The calculations leading to Eqs. (1a) and (1b) are detailed in Appendix
A.

Our goal in this section is to use the transfer-matrix formalism to understand
the propagation of graphene surface plasmons and to obtain the energy
dispersion, associating the Bloch momentum $q$, arising from the periodic
arrangement of the metal grating, with the energy $\hbar\omega$.

###  II.1 Fields, boundary conditions and transfer matrices

We begin by defining the expressions of the electrostatic potential describing
propagating plasmons in the ungated and gated regions

$\phi_{u}(x)=\phi_{u}^{(+,n)}e^{iq_{u}(x+l_{u}-nl)}+\phi_{u}^{(-,n)}e^{-iq_{u}(x+l_{u}-nl)},$
(2a)
$\phi_{g}(x)=\phi_{g}^{(+,n)}e^{iq_{g}(x-nl)}+\phi_{g}^{(-,n)}e^{-iq_{g}(x-nl)},$
(2b)

where the $+(-)$ superscript correspond to the wave component traveling to the
right (left). We can think of the gated regions as rectangular barriers of
length $l_{g}$ that begin at $x=nl$ and end at $x=nl+l_{g}$, and of ungated
regions, that begin at $nl-l_{u}$ and end at $nl$, where $n\in\mathbb{Z}$ (for
details see Fig. 2). We can then make use of appropriate boundary conditions
to relate the electrostatic potential between gated and ungated regions. At
the boundaries between regions we impose: (i) the continuity of the
electrostatic potential $\phi(x)$ and (ii) the continuity of the current
density $\sigma\frac{\partial\phi(x)}{\partial x}$, where $\sigma$ is the
graphene conductivity (Drude conductivity) defined at the beginning of the
section. At $x=-l_{u}+nl$ we have:

$\phi_{g}^{(+,n-1)}e^{iq_{g}l_{g}}+\phi_{g}^{(-,n-1)}e^{-iq_{g}l_{g}}=\phi_{u}^{(+,n)}+\phi_{u}^{(-,n)},$
(3a)
$\displaystyle\sigma_{g}q_{g}\left(\phi_{g}^{(+,n-1)}e^{iq_{g}l_{g}}-\phi_{g}^{(-,n-1)}e^{-iq_{g}l_{g}}\right)=$
$\displaystyle\sigma_{u}q_{u}\left(\phi_{u}^{(+,n)}-\phi_{u}^{(-,n)}\right),$
(3b)

while at $x=nl$:

$\phi_{u}^{(+,n)}e^{iq_{u}l_{u}}+\phi_{u}^{(-,n)}e^{-iq_{u}l_{u}}=\phi_{g}^{(+,n)}+\phi_{g}^{(-,n)},$
(4a)
$\displaystyle\sigma_{u}q_{u}\left(\phi_{u}^{(+,n)}e^{iq_{u}l_{u}}-\phi_{u}^{(-,n)}e^{-iq_{u}l_{u}}\right)=$
$\displaystyle\sigma_{g}q_{g}\left(\phi_{g}^{(+,n)}-\phi_{g}^{(-,n)}\right).$
(4b)

Figure 2: Scheme of an unit cell beginning at $x=nl-l_{u}$ and finishing at
$x=nl+l_{g}$ from a front view. The gated region is comprised between
$nl<x<nl+l_{g}$ and the ungated region is comprised between $nl-l_{u}<x<nl$.
The transfer matrix relating the wave amplitudes to the left of the gated
region with the wave amplitudes to the right of the gated region is ${\bf
M}_{g}$. The transfer matrix relating the wave amplitudes to the left of the
gated region with the wave amplitudes to the right of the ungated region is
${\bf M}_{u}$. The unit cell transfer matrix is then defined as ${\bf M}={\bf
M}_{u}{\bf M}_{g}$.

We then define the gated and ungated region transfer matrices, according to
Fig. 2:

$\displaystyle\left(\begin{array}[]{c}\phi_{g}^{(+,n)}\\\
\phi_{g}^{(-,n)}\end{array}\right)={\bf
M}_{u}\left(\begin{array}[]{c}\phi_{u}^{(+,n)}\\\
\phi_{u}^{(-,n)}\end{array}\right)$ (5e)
$\displaystyle\left(\begin{array}[]{c}\phi_{u}^{(+,n)}\\\
\phi_{u}^{(-,n)}\end{array}\right)={\bf
M}_{g}\left(\begin{array}[]{c}\phi_{g}^{(+,n-1)}\\\
\phi_{g}^{(-,n-1)}\end{array}\right)$ (5j)

The transfer matrix ${\bf M}_{g}$ can be obtained from Eqs. (3a) and (3b) and
will have the form

${\bf
M}_{g}=\tfrac{1}{2}\left(\begin{array}[]{cc}\left(1+\frac{\sigma_{g}q_{g}}{\sigma_{u}q_{u}}\right)e^{iq_{g}l_{g}}&\left(1-\frac{\sigma_{g}q_{g}}{\sigma_{u}q_{u}}\right)e^{-iq_{g}l_{g}}\\\
\\\ \vspace{-7mm}\hfil\\\
\left(1-\frac{\sigma_{g}q_{g}}{\sigma_{u}q_{u}}\right)e^{iq_{g}l_{g}}&\left(1+\frac{\sigma_{g}q_{g}}{\sigma_{u}q_{u}}\right)e^{-iq_{g}l_{g}}\end{array}\right).$
(6a) In a similar manner, the transfer matrix ${\bf M}_{u}$ can be obtained
from Eqs. (4a) and (4b) and represented as ${\bf
M}_{u}=\tfrac{1}{2}\left(\begin{array}[]{cc}\left(1+\frac{\sigma_{u}q_{u}}{\sigma_{g}q_{g}}\right)e^{iq_{u}l_{u}}&\left(1-\frac{\sigma_{u}q_{u}}{\sigma_{g}q_{g}}\right)e^{-iq_{u}l_{u}}\\\
\\\ \vspace{-7mm}\hfil\\\
\left(1-\frac{\sigma_{u}q_{u}}{\sigma_{g}q_{g}}\right)e^{iq_{u}l_{u}}&\left(1+\frac{\sigma_{u}q_{u}}{\sigma_{g}q_{g}}\right)e^{-iq_{u}l_{u}}\end{array}\right).$
(6b)

###  II.2 Dispersion and energy bands

Coefficients of the electrostatic potential across one period of the
structure, at $x=(n-1)l$ and $x=nl$, can be related by combining Eqs. (6a) and
(6b). Consequently, we must then have:

$\left(\begin{array}[]{c}\phi_{g}^{(+,n)}\\\
\phi_{g}^{(-,n)}\end{array}\right)={\bf
M}\left(\begin{array}[]{c}\phi_{g}^{(+,n-1)}\\\
\phi_{g}^{(-,n-1)}\end{array}\right).$ (7)

The matrix ${\bf M}={\bf M}_{u}{\bf M}_{g}$ is the unit cell transfer matrix,
which is unitary, i.e., ${\rm det}\left({\bf M}\right)=1$. The Bloch momentum
$q$ is introduced by invoking Bloch’s theorem:

$\left(\begin{array}[]{c}\phi_{g}^{(+,n)}\\\
\phi_{g}^{(-,n)}\end{array}\right)=e^{iqnl}\left(\begin{array}[]{c}\phi_{g}^{(+,0)}\\\
\phi_{g}^{(-,0)}\end{array}\right).$ (8)

Non-trivial solutions of the system of Eqs. (7) and (8) for the coefficients
$\phi_{g}^{(+,0)},\phi_{g}^{(-,0)}$ are obtained if the determinant vanishes,
i.e.,

${\rm det}\left({\bf M}-e^{iql}{\bf I}\right)=0,$ (9)

where ${\bf I}$ is the identity matrix. After taking into account the property
of unitary matrices, Eq. (9) can be represented as [49]

$2\cos(ql)={\rm Tr}({\bf M}).$ (10)

This condition yields a Kronig–Penney-like dispersion relation for the
graphene plasmons which, after taking the explicit form of matrices ${\bf
M}_{u}$ and ${\bf M}_{g}$ into account, can be written in the form

${\rm cos}(ql)={\rm cos}(q_{u}l_{u}){\rm cos}(q_{g}l_{g})-Z{\rm
sin}(q_{u}l_{u}){\rm sin}(q_{g}l_{g}),$ (11)

where
$Z=\left[\frac{\sigma_{u}q_{u}}{2\sigma_{g}q_{g}}+\frac{\sigma_{g}q_{g}}{2\sigma_{u}q_{u}}\right]$.
The numerical solution of the dispersion relation (11) is shown in Fig. 3. The
periodicity of the structure results in a spectrum that consists in energy
bands (where propagation of waves is allowed) and forbidden gaps (where the
propagation is exponentially suppressed). Notice that in the low-frequency
region the dispersion is linear (as shown in inset of Fig. 3). In this regime,
we can expand Eq. (11) into a series with respect to the frequency $\omega$,
obtaining the following approximation

$\left|q\right|\approx\omega\frac{l_{g}}{l}\sqrt{\frac{\hbar\epsilon_{2}}{4E_{F}c\alpha
d}\left(1+\frac{\sigma_{g}l_{u}}{\sigma_{u}l_{g}}\right)}.$ (12)

Figure 3: First three energy bands as a function of the phase $ql$. The
parameters used in this numerical calculation were $l=500$ nm, $l_{g}=240$ nm,
$l_{u}=260$ nm, $d=100$ nm, and $E_{F}=0.45$ eV. The dielectrics encapsulating
the graphene sheet were considered air ($\epsilon_{1}=1)$ and h-BN
($\epsilon_{2}=3.5)$, such that the metal grating is on top of h-BN. The inset
in the figure shows the dispersion relation closer to $ql=0$ where it is
possible to see that, for very small $q$, the dispersion is linear.

###  II.3 Propagation of graphene plasmons

With both ${\bf M}_{u}$ and ${\bf M}_{g}$, we have all necessary tools to
describe the propagation of the graphene plasmon polaritons through the
system. The transmission coefficient $T_{N}$ for a system with $N$ identical
unit cells is given by [49]

$T_{N}=\frac{1}{1+|{\bf M}_{12}^{N}|^{2}}$ (13)

where ${\bf M}^{N}$ is the transfer matrix of the system with $N$ unit cells,
which is the product of $N$ matrices ${\bf M}$. The element ${\bf M}_{12}^{N}$
is given by Chebyshev’s identity [49] as:

${\bf M}_{12}^{N}={\bf M}_{12}\frac{{\rm sin}(Nql)}{{\rm sin}(ql)}.$ (14)

The quantity $ql$ is a phase related to the Bloch momentum and obeys Eq. (10).
Thus, ${\bf M}_{12}^{N}$ is easily obtainable after multiplying ${\bf
M}_{u},{\bf M}_{g}$ and the computation of $T_{N}$ yields:

$\displaystyle T_{N}$ $\displaystyle=\frac{1}{1+\frac{{\rm sin^{2}}(Nql)}{{\rm
sin^{2}}(ql)}|{\bf M}_{12}|^{2}}$ (15) $\displaystyle=\frac{1}{1+\frac{{\rm
sin^{2}}(Nql)}{{\rm sin^{2}}(ql)}\left|\frac{i{\rm
sin}(q_{u}l_{g})e^{iq_{g}l_{u}}\left(\sigma_{g}q_{g}+\sigma_{u}q_{u}\right)\left(\sigma_{g}q_{g}-\sigma_{u}q_{u}\right)}{2\sigma_{g}q_{g}\sigma_{u}q_{u}}\right|^{2}}.$

In Fig. 4 we represent the transmission coefficient $T_{N}$ as a function of
the energy $\hbar\omega$ for $N=100$ and $N=25$ unit cells. As it can be seen
from Fig. 4, for the frequency ranges corresponding to the gaps in the
spectrum (see Fig. 3), transmittance turns to be zero, i.e., $T_{N}=0$. Energy
regions with zero transmittance correspond to the condition $2<{\rm Tr}({\bf
M})$, as it follows from Eq. (10). For these frequencies the periodic
structure serves as a plasmonic Bragg mirror. At the same time, for the
frequency ranges corresponding to allowed bands in the spectrum, the nonzero
transmittance exhibits periodical oscillations. Moreover, for frequencies
corresponding to Bloch wavevector $ql=m\pi/N$, where $1\leq m\leq N-1$, the
system becomes fully transparent and $T_{N}=1$.

Figure 4: System’s transmission coefficient $T_{N}$ as a function of the
energy $\hbar\omega$ for $100$ unit cells (blue) and $25$ unit cells (red).
The parameters used in this figure were $l=500$ nm, $l_{g}=240$ nm,
$l_{u}=260$ nm, $d=100$ nm and $E_{F}=0.45$ eV. The dielectrics encapsulating
the graphene sheet were considered air ($\epsilon_{1}=1)$ and h-BN
($\epsilon_{2}=3.5)$.

##  III Plasmonic SSH tight-binding model

Figure 5: Scheme of the SSH representation of the plasmonic crystal. The two
sub-lattices are depicted with different colors: the purple one has lattice
sites at $x=nl+lg$ and the green one at $x=nl$, where $n\in\mathbb{Z}$. Since
these lattice sites are defined at boundaries between gated and ungated
regions, the sub-indexes $u,g$ are irrelevant.

The goal in this section is to represent the electrostatic potential
describing the plasmons in the tight-binding formalism. This representation is
rather useful for the discussion of the topological description of the bands.
The symmetry of the problem (absence of inversion symmetry and two segments
composing one unit cell) invites us to make an analogy with a system that
possesses sub-lattice symmetry, like the SSH model. To do this, we need to
define the two sub-lattices: the first one is defined by the points at
$x=nl+lg$, depicted in purple in Fig. 5 and is associated with the end of
gated regions. The other sub-lattice is defined by points at $x=nl$, depicted
in green, and is associated with the end of ungated regions. To derive the
tight-binding representation we then need to find expressions for the fields
$\phi_{u}(x)$, $\phi_{g}(x)$, $\sigma_{u}\frac{\partial\phi_{u}(x)}{\partial
x}$, and $\sigma_{g}\frac{\partial\phi_{g}(x)}{\partial x}$ at each lattice
site as a function of the previous and the next sites using the transfer-
matrix formalism. For example, we want to obtain
$C\phi(nl+l_{g})=A\phi(nl)+B\phi(nl+l)$ and
$C\phi(nl)=A\phi(nl+l_{g})+B\phi(nl-l_{u})$, where $A$, $B$, and $C$ are the
tight-binding coefficients to be determined.

###  III.1 Energy-dependent hopping amplitudes

The detailed derivation of the tight-binding equations can be found in
Appendix B, where we show how to write the following equations relating
lattice sites from two different sub-lattices:

$\phi(nl+l_{g})\frac{\sigma_{g}q_{g}}{{\rm sin}(q_{g}l_{g})}+\phi(nl-
lu)\frac{\sigma_{u}q_{u}}{{\rm
sin}(q_{u}l_{u})}=\phi(nl)\left[\frac{\sigma_{u}q_{u}{\rm
cos}(q_{u}l_{u})}{{\rm sin}(q_{u}l_{u})}+\frac{\sigma_{g}q_{g}{\rm
cos}(q_{g}l_{g})}{{\rm sin}(q_{g}l_{g})}\right]$ (16a)
$\phi(nl)\frac{\sigma_{g}q_{g}}{{\rm
sin}(q_{g}l_{g})}+\phi(nl+l)\frac{\sigma_{u}q_{u}}{{\rm
sin}(q_{u}l_{u})}=\phi(nl+l_{g})\left[\frac{\sigma_{u}q_{u}{\rm
cos}(q_{u}l_{u})}{{\rm sin}(q_{u}l_{u})}+\frac{\sigma_{g}q_{g}{\rm
cos}(q_{g}l_{g})}{{\rm sin}(q_{g}l_{g})}\right]$ (16b)

These equations relate the electrostatic potential in one sub-lattice with the
electrostatic potential at nearest-neighbor lattice points, contained in the
other sub-lattice. The functions multiplying the electrostatic potentials play
the role of hopping amplitudes in electronic tight-binding equations. To make
Eqs. (16a) and (16b) look more similar to the SSH tight-binding Hamiltonian,
we introduce

$\displaystyle A$ $\displaystyle=$ $\displaystyle\frac{\sigma_{g}q_{g}}{{\rm
sin}(q_{g}l_{g})}$ (17a) $\displaystyle B$ $\displaystyle=$
$\displaystyle\frac{\sigma_{u}q_{u}}{{\rm sin}(q_{u}l_{u})}$ (17b)
$\displaystyle C$ $\displaystyle=$
$\displaystyle\left[\frac{\sigma_{u}q_{u}{\rm cos}(q_{u}l_{u})}{{\rm
sin}(q_{u}l_{u})}+\frac{\sigma_{g}q_{g}{\rm cos}(q_{g}l_{g})}{{\rm
sin}(q_{g}l_{g})}\right]$ (17c) $\displaystyle\psi_{n}$ $\displaystyle=$
$\displaystyle\phi(nl)$ (17d) $\displaystyle\xi_{n}$ $\displaystyle=$
$\displaystyle\phi(nl+l_{g})$ (17e)

which conveniently allows us to rewrite Eqs. (16a) and (16b) as:

$\displaystyle C\xi_{n}=A\psi_{n}+B\psi_{n+1}$ (18a) $\displaystyle
C\psi_{n}=A\xi_{n}+B\xi_{n-1}$ (18b)

These coefficients are energy dependent through $q_{u}$, $q_{g}$,
$\sigma_{u}$, and $\sigma_{g}$, and in this regard they differ from the usual
SSH model where $A$, $B$, and C are constants. In order to vary these
coefficients, it is then necessary to vary the gate tension applied to the
metallic rods. Bloch’s theorem allows us to write $\psi_{n}=\psi_{0}e^{iqnl}$
and $\xi_{n}=\xi_{0}e^{-iqnl}$. In matrix notation, Eq. (18) then becomes:

$\left(\begin{array}[]{cc}0&A+Be^{-iql}\\\
A+Be^{iql}&0\end{array}\right)\left(\begin{array}[]{c}\psi_{0}\\\
\xi_{0}\end{array}\right)=C\left(\begin{array}[]{c}\psi_{0}\\\
\xi_{0}\end{array}\right)$ (19)

The matrix with the coefficients $A$ and $B$ can be written as a linear
combination of the Pauli matrices $\sigma_{x}$ and $\sigma_{y}$, as usual in
these type of problems [47]. That is, we can write

$\begin{array}[]{c}\left[A(\omega,d)+B(\omega){\rm
cos}(ql)\right]\sigma_{x}+B(\omega){\rm sin}(ql)\sigma_{y}\\\ \equiv
d_{x}(\omega,d,q)\sigma_{x}+d_{y}(\omega,q)\sigma_{y}.\end{array}$ (20)

The coefficient $C$, playing the role of the eigenvalue on the SSH model, has
the form

$C(\omega,q,d)=\pm\sqrt{A^{2}(\omega,d)+B^{2}(\omega)+2A(\omega,d)B(\omega){\rm
cos}(ql)},$ (21)

where the explicit dependencies of the coefficients on $\omega$, $q$, and $d$
are kept as arguments of $A$, $B$, and $C$. The two possible solutions of
$C(\omega,q,d)$ correspond to two separate bands when plotted against $ql$.
This allows us to make a parallel with the original SSH [47], in which
$C(\omega,q,d)$ plays the role of the SSH model’s energy and Eq. (20) plays
the role of the Hamiltonian. We make use of this parallelism to study the
topological features of the plasmonic SSH.

###  III.2 Topological characterization in momentum space

In this section, we investigate the topology of our plasmonic SSH model,
focusing our attention on the behavior of the first plasmonic energy gap,
alongside the coefficient $C(\omega,q,d)$. We note that when there is a gap in
the energy bands, there is also a gap in the plot $C(\omega,q,d)$ versus $ql$.
As we vary the distance from the metal $d$ we note that, as the first band gap
closes, the gap in $C(\omega,q,d)$ also closes at $ql=\pm\pi$. This is exactly
what happens in the electronic SSH model, but with the energy $E$ instead of
$C(\omega,q,d)$ [47]. We then expect that, if $d_{c}$ is the distance from the
metal that closes the energy gap and the gap in $C(\omega,q,d)$, values of $d$
greater and smaller than $d_{c}$ will correspond to different topological
phases. It is then necessary to find an analytical expression for $d_{c}$, as
we show next for the first energy band gap. Since both bands of
$C(\omega,q,d)$ touch at $ql=\pm\pi$ when the first energy gap closes, we can
rewrite Eq. (21) as

$\begin{array}[]{c}C(\omega,\pm\frac{\pi}{l},d)=\pm\sqrt{A^{2}(\omega,d)+B^{2}(\omega)-2A(\omega,d)B(\omega)}=0\\\
\\\
C(\omega,\pm\frac{\pi}{l},d)=\pm\left[A(\omega,d)-B(\omega)\right]=0\end{array}$
(22)

Then it is clear that, for this condition to hold, $A=B$. By taking the
explicit form of the coefficients as in Eq. (17) into consideration, the
conditions $A=B$ and $C=0$ can be rewritten as:

$C=A{\rm cos}(q_{g}l_{g})+A{\rm cos}(q_{u}l_{u})=0$ (23a)
$A=B=\frac{\sigma_{g}q_{g}}{{\rm sin}(q_{g}l_{g})}=\frac{\sigma_{u}q_{u}}{{\rm
sin}(q_{u}l_{u})}$ (23b)

From Eq. (23a), we see that the cosines are equal but with opposite signs

${\rm cos}(q_{g}l_{g})=-{\rm cos}(q_{u}l_{u}).$ (24)

The sine terms ${\rm sin}(q_{g}l_{g})$ and ${\rm sin}(q_{u}l_{u})$ are then
either identical or have the same modulus but opposite sign as well. From Eq.
(23b), we see that the sines must have the same sign, since $\sigma_{g}q_{g}$
and $\sigma_{u}q_{u}$ are strictly positive quantities and, to fulfill
condition (23b), they should satisfy

$\sigma_{g}q_{g}=\sigma_{u}q_{u}.$ (25)

Angles that have the same sine and opposite cosine sum to an odd multiple of
$\pi$:

$q_{g}l_{g}+q_{u}l_{u}=m\pi,\quad m=1,3,5\ldots$ (26)

Now, from Eqs. (25) and (26) we can derive an analytical expression for the
distance $d_{c}$ from the metal that makes the two bands of $C$ touch. This is
straightforward if the explicit expressions for $q_{g}$ and $q_{u}$ are
considered, given by Eqs. (1a) and (1b). Thus, for equal conductivities
$\sigma_{g}=\sigma_{u}$, condition (25) transforms into $q_{u}=q_{g}$, which
is fulfilled for a certain frequency $\omega_{*}$ (such that $\hbar\omega_{*}$
is the energy at the point where the two bands touch) [see Fig. 6(a)]. By
combining Eqs. (1a) and (1b), an expression for frequency $\omega_{*}$ can be
determined

$\omega_{*}\left(d\right)=\sqrt{\frac{4E_{F}c\alpha\epsilon_{2}}{\hbar
d(\epsilon_{1}+\epsilon_{2})^{2}}}.$ (27)

This frequency $\omega_{*}$ varies with thickness $d$ as $\sim 1/\sqrt{d}$, as
it is shown in Fig. 6(b) by the blue line (left axis). At the same time,
condition (26) is fulfilled for $m=1$ and for critical thickness $d_{c}$ as
expected [crossing of orange line in Fig. 6(b) with horizontal dashed line].
Substituting Eqs. (1a) and (1b) for frequency $\omega_{*}$ into condition
(26), we obtain an equation for the critical thickness in the form

$\omega_{*}\left(d_{c}\right)\sqrt{\frac{\hbar\epsilon_{2}}{4E_{F}c\alpha
d_{c}}}l_{g}+\frac{\hbar(\epsilon_{1}+\epsilon_{2})}{4E_{F}c\alpha}\left[\omega_{*}\left(d_{c}\right)\right]^{2}l_{u}=m\pi,$
(28)

from which we obtain

$d_{c}^{(m)}=\frac{\epsilon_{2}l}{m\pi(\epsilon_{1}+\epsilon_{2})}.$ (29)

Figure 6: (a) Dispersion relations of gated [Eq. (1a), green line] or ungated
[Eq. (1b), red line] plasmons for $d=100$ nm. Dashed horizontal line depicts
frequency $\omega_{*}$, for which $q_{u}=q_{g}$; (b) Dependence of frequency
$\omega_{*}$ (left vertical axis) and phase over period
$q_{u}l_{u}+g_{g}l_{g}$ (right axis) upon distance $d$. Dashed horizontal line
correspond to phase $q_{u}l_{u}+g_{g}l_{g}=\pi$, while vertical dashed line
depicts critical thickness $d_{c}$. In all panels other parameters are:
$E_{F}=0.45$ eV, $\epsilon_{1}=1$, $\epsilon_{2}=3.5$, $l_{u}=260$ nm, and
$d=100$ nm .

This expression for $d_{c}$ is valid for all odd band gaps. By substituting
the critical distance given by Eq. (29) into Eq. (27), we find the energies
for which the odd band gaps close as only a function of $m$ and the parameters
of the system. We observe that if $m=1$, we have the value of
$\hbar\omega_{*}$ at which the first band gap closes and, if $m=3$, we have
the value of $\hbar\omega_{*}$ at which the third band gap closes, and so on.
At the end, we see that $m$ plays the role of the band gap index. We would
like to point out that this formalism also applies to those band gaps that
close at $ql=0$, like the second and fourth ones. The only difference is that,
for this case, the cosines are equal in sign and modulus, while the sines have
equal modulus but opposite sign. We would arrive at the conclusion that
$q_{g}l_{g},q_{u}l_{u}$ sum to an even multiple of $\pi$, obeying
$q_{g}l_{g}+q_{u}l_{u}=p\pi$ with $p=2,4,6\ldots$. Thus, Eq. (29) applies to
all band gaps if we redefine $m=1,2,3,4\ldots$.

With the value of $d_{c}$ at hand, the remaining task is to identify which
topological phase corresponds to $d>d_{c}$ and which corresponds to $d<d_{c}$
for the first energy band. This information can be obtained with the functions
$d_{x}$ and $d_{y}$ as defined in Eq. (20). From the electronic SSH model [47]
we know that, as $ql$ varies from $-\pi$ to $\pi$, these functions describe a
curve in the $d_{x}-d_{y}$ plane. This curve contains all information we need
to state on which topological phase the system is in: if the curve winds
around the origin, the system is said to be in the topological phase; if the
curve does not wind around the origin, the system is said to be in the trivial
phase; if the curve passes exactly through the origin, the system lies between
topological phases and the energy band gap of the SSH model is closed. So, we
expect that after choosing a value of $d$ and plotting $d_{x}(\omega,d,q)$ and
$d_{y}(\omega,q)$ for the corresponding energies of the first energy band, as
$ql$ varies from $-\pi$ to $\pi$, we can find in which topological phase the
system is in.

The results of this discussion are shown in Fig. 7. The first row corresponds
to $d<d_{c}$, the second row to $d=d_{c}$ and the third row to $d>d_{c}$. In
the first column, the first two energy bands are depicted. In the second
column, the coefficient $C(\omega,q,d)$ is plotted against $ql$. In the third
column, we plot the parametric curve in the $d_{x}-d_{y}$ plane. Initially,
looking at only the energy bands and the coefficient $C(\omega,q,d)$, it is
not possible to distinguish between topological phases. However, the curve in
the third column for $d<d_{c}$ does not wind around the origin, while the
curve for $d>d_{c}$ does. This indicates that the first band has trivial
topology for $d<d_{c}$ and non-trivial topology for $d>d_{c}$. In the second
row, we use $d=d_{c}$ to show that indeed the energy gap and the
$C(\omega,q,d)$ gap close simultaneously. Additionally, the corresponding
$d_{x}-d_{y}$ curve passes exactly at the origin.

Figure 7: Topological characterization of the first energy band of the
plasmonic SSH crystal. The parameters for this representations were $l=500$
nm, $l_{g}=240$ nm, $l_{u}=260$ nm, $E_{F}=0.45$ eV, $\epsilon_{1}=1$ and
$\epsilon_{2}=3.5$. For this choice of parameters, $d_{c}\approx 123.79$ nm,
for which the gap closes at $\hbar\omega\approx 0.060$ eV. The distance to the
metal was chosen to be $d=100$ nm for the first row and $d=150$ nm for the
second row. ${\bf a)},{\bf d)},{\bf g)}$ Energy bands as a function of $ql$;
${\bf b)},{\bf e)},{\bf h)}$ Coefficient $C(\omega,q,d)$ as a function of
$ql$; ${\bf c)},{\bf f)},{\bf i)}$ Parametric plot in the plane $d_{x}-d_{y}$
as $ql$ varies from $-\pi$ to $\pi$. We emphasize that for each $ql$, the
corresponding value of $\hbar\omega$ in the first energy band has to be
computed into $d_{x}(\omega,d,q),d_{y}(\omega,q)$. For $d=100$ nm, the
topology of the first band is trivial. For $d=d_{c}$, we find ourselves at the
boundary between topological phases. For $d>d_{c}$, the topology of the first
band is non-trivial.

###  III.3 Defect states

One of the advantages of characterization of structures by their topological
properties is the ability to predict the existence of surface and defect
states. In this subsection, we consider periodical structures with defects,
which give rise to the presence of localized states. In Fig. 8 we show the
spectrum of defect states in the structure, obtained from the periodic
plasmonic crystal (considered in previous sections) by increasing twice the
length of one of the gated regions. The presence of a single defect in the
structure allows for the existence of localized defect states of two types:
even and odd, whose spectra are shown in Fig. 8(a) by blue and red lines,
respectively. Both defect states appear inside the first gap of the energy
spectrum of the periodical structure [whose edges are depicted by black lines
in Fig. 8(a)]. Even defect states exist for $d>d_{c}$, for which the first
band is characterized by non-trivial topology, while odd defect states exist
for $d<d_{c}$, for which the first band has trivial topology. Also, the
spatial distribution of an even defect state is characterized by a reflective
symmetry of the scalar potential $\phi(x)=\phi(-x)$ and reflective anti-
symmetry of the current density $j_{x}(x)=-j_{x}(-x)$ with respect to the
symmetry plane $x=0$ [blue lines in Figs. 8(b) and 8(c)]. On the contrary, the
spatial distribution of an odd defect state [red lines in Figs. 8(b) and 8(c)]
is characterized by a reflective anti-symmetry of the scalar potential
$\phi(x)=-\phi(-x)$ and reflective symmetry of the current density
$j_{x}(x)=j_{x}(-x)$. The existence of defect states is also possible in a
different version of structure, where one of the ungated regions has twice its
usual length instead [see inset in Fig. 9]. As it can be seen from Fig. 9, the
situation is reversed: even defect states (depicted by blue line) exist in the
gap for $d<d_{c}$, while odd defect states (depicted by red line) exist when
$d>d_{c}$. Thus, we have the opposite situation if we compare it to the case
with a defect in the central gated region [compare with Fig. 8(a)], where even
and odd states exist for $d>d_{c}$ and $d<d_{c}$, respectively.

Figure 8: Defect states in the structure, where one (central) of the gated
regions has double thickness, compared to that of the periodic structure. (a)
Frequencies of first energy band gap edges (black lines), even (blue line) and
odd (red line) defect state versus thickness $d$; (b,c) Dependence of scalar
potential [panel (b)] or current density [panel (c)] upon coordinate $x$ of
even (blue lines) and odd (red line) defect states. In panels (b) and (c),
even and odd defect states are calculated for $d=150$ nm and $d=100$ nm,
respectively. Also in panels (b) and (c), ungated (gated) regions are shaded
by blue (grey) colors. As in all other panels, the parameters of the structure
are $E_{F}=0.45$ eV, $\epsilon_{1}=1$, $\epsilon_{2}=3.5$. Figure 9: Even
(blue line) and odd (red line) defect states in the structure, where one
(central) of the ungated regions has double thickness, compared to that in
periodic structure with parameters $E_{F}=0.45$ eV, $\epsilon_{1}=1$,
$\epsilon_{2}=3.5$. The schematics of considered structure is shown in the
inset, where ungated (gated) regions are shaded by blue (grey) colors.
Frequencies of first energy band gap edges are depicted by black lines.

###  III.4 Analytical approximation for Dirac bands

Another interesting aspect of the plasmonic SSH model is the fact that, when
the energy gaps close, we obtain Dirac cones. It is possible to obtain an
analytical approximations for the bands by inspecting the relationship between
coefficients $A$ and $B$. As stated in the previous section, for odd band gaps
we have

$\begin{array}[]{c}A=B=\frac{\sigma_{g}q_{g}}{{\rm
sin}(q_{g}l_{g})}=\frac{\sigma_{u}q_{u}}{{\rm sin}(q_{u}l_{u})}\\\ \\\ {\rm
sin}(q_{g}l_{g})={\rm sin}(q_{u}l_{u})\end{array}$ (30)

meaning that $\sigma_{g}q_{g}=\sigma_{u}q_{u}$. This is also true for even
band gaps, since in these cases the conditions are

$\begin{array}[]{c}A=-B=\frac{\sigma_{g}q_{g}}{{\rm
sin}(q_{g}l_{g})}=-\frac{\sigma_{u}q_{u}}{{\rm sin}(q_{u}l_{u})}\\\ \\\ {\rm
sin}(q_{g}l_{g})=-{\rm sin}(q_{u}l_{u})\end{array}$ (31)

Thus, the relationship $\sigma_{g}q_{g}=\sigma_{u}q_{u}$ is always true, given
that $q_{g},q_{u}$ are evaluated at the energy where the bands touch and at
$d=d_{c}$. In this scenario, we can alter the form of the Kronig–Penney-like
dispersion found for the graphene plasmons in Eq. (11) by substituting
$\sigma_{g}q_{g}=\sigma_{u}q_{u}$

$\begin{array}[]{c}\begin{array}[]{c}{\rm cos}(ql)={\rm cos}(q_{u}l_{u}){\rm
cos}(q_{g}l_{g})-{\rm sin}(q_{u}l_{u}){\rm sin}(q_{g}l_{g})\end{array}\\\
={\rm cos}(q_{u}l_{u}+q_{g}l_{g})\end{array}$ (32)

The general solution to this equation has the form

$\pm
ql+2n\pi=\frac{\hbar(\epsilon_{1}+\epsilon_{2})\omega^{2}}{4E_{F}c\alpha}l_{u}+\omega\sqrt{\frac{\hbar\epsilon_{2}}{4E_{F}c\alpha
d_{c}}}l_{g},$ (33)

where $n\in\mathbb{Z}$ and the explicit expressions for $q_{u}$ and $q_{g}$
were used. Solving the quadratic equation for $\hbar\omega$ and preserving
only the solutions corresponding to positive energies we arrive at

$\hbar\omega=\frac{\hbar^{2}}{(\epsilon_{1}+\epsilon_{2})l_{u}E_{F}c\alpha}\sqrt{\frac{\epsilon_{2}}{d_{c}}}l_{g}\left(\sqrt{1+\frac{4d_{c}}{\epsilon_{2}}\frac{(\epsilon_{1}+\epsilon_{2})l_{u}(\pm
ql+2n\pi)}{l_{g}^{2}}}-1\right).$ (34)

Again, we emphasize that
$\frac{\sigma_{u}q_{u}}{2\sigma_{g}q_{g}}+\frac{\sigma_{g}q_{g}}{2\sigma_{u}q_{u}}$
= 1 only at the point where the bands touch and thus Eq. (34) is only an
approximation. At the vicinity of the first gap (which is closed in this case)
we use an expression for the critical thickness (29) for $m=1$, as well as
define Block wavevector in the limits $0\leq ql\leq 2\pi$. In this situation,
two branches of frequencies can be obtained from Eq. (34) as

$\displaystyle\omega=\omega_{*}\left(d_{c}^{(1)}\right)\frac{\sqrt{l_{g}^{2}+\frac{ql}{\pi}4ll_{u}}-l_{g}}{2l_{u}},$
(35)
$\displaystyle\omega=\omega_{*}\left(d_{c}^{(1)}\right)\frac{\sqrt{l_{g}^{2}+\left(2-\frac{ql}{\pi}\right)4ll_{u}}-l_{g}}{2l_{u}}.$
(36)

At the vicinity of the second gap, the expression for the critical thickness
(29) should be used with $m=2$, and the wavevector is defined inside the
limits $-\pi\leq ql\leq\pi$. In this case, Eq. (34) for the two branches of
frequencies can be rewritten as

$\displaystyle\omega=\omega_{*}\left(d_{c}^{(2)}\right)\frac{\sqrt{l_{g}^{2}+\left(1\pm\frac{ql}{2\pi}\right)4ll_{u}}-l_{g}}{2l_{u}}.$
(37)

Close to the Dirac points, the square root in the expression above can be
Taylor expanded such that the relationship between $\omega$ and $ql$ is
linear. For the first Dirac point, Eqs. (35) and (36) in the vicinity of
$ql\simeq\pi$ will be represented as

$\omega=\omega_{*}\left(d_{c}^{(2)}\right)\left[1\pm\left(\frac{ql}{\pi}-1\right)\frac{l}{l+l_{u}}\right].$
(38)

For the second Dirac point, Eq. (37) in the vicinity of $ql\simeq 0$ can be
expanded as

$\omega=\omega_{*}\left(d_{c}^{(2)}\right)\left(1\pm\frac{ql}{2\pi}\frac{l}{l+l_{u}}\right).$
(39)

In Fig. 10 we superpose the numerical calculations and the analytical
approximation for the first and second energy band gaps. As expected, the
result in Eq. (34) is better suited for $ql$ around the Dirac points.

Figure 10: Energy bands and analytical approximation for $\mathbf{a)}$ first
band gap closed, represented from $ql=[\pi,2\pi]$ and $\mathbf{b)}$ second
band gap closed, represented from $ql=[-\pi,\pi]$. The analytical
approximation is very good for $ql$ close to the Dirac point and diverges
significantly from the numerical calculations close to the band edges. This is
due to the fact that, in the approximation, the function
$\frac{\sigma_{u}q_{u}}{2\sigma_{g}q_{g}}+\frac{\sigma_{g}q_{g}}{2\sigma_{u}q_{u}}$
was evaluated only at the energy where the Dirac point occurs. The parameters
used to produce these figures were $l=500$ nm, $l_{g}=240$ nm, $l_{u}=260$ nm,
$E_{F}=0.45$ eV, $\epsilon_{1}=1$ and $\epsilon_{2}=3.5$. For this choice, the
first gap closes at $\hbar\omega\approx 0.060$ eV for $d_{c}\approx 123.79$ nm
. The second gap closes at $\hbar\omega\approx 0.085$ eV for $d_{c}\approx
61.89$ nm.

##  IV Conclusions

In this work we investigated a graphene-based polaritonic crystal. We showed
that the plasmon energies obey a Kronig–Penney dispersion and calculated their
propagation on the crystal with the transfer-matrix formalism. Upon
realization that the crystal possesses sub-lattice symmetry, we managed to
derive tight-binding equations for the electrostatic potential $\phi(x)$
describing the plasmons. At the end, these equations resemble the SSH model
tight-binding, one difference being that the hopping amplitudes are energy
dependent. This representation is convenient because it allows us to bring the
traditional formalism used to describe topological insulators in one dimension
to plasmonics. By investigating the energy dispersion and the tight-binding
coefficients we demonstrate that, by tuning the distance $d$ between graphene
and the metal grating, it is possible to associate trivial and non-trivial
topology to the energy bands of the polaritonic crystal. Additionally, we
derived an analytical expression for $d_{c}$, the value of $d$ that closes the
energy band gaps. This expression can be used alongside numerical methods to
easily identify the boundary between topological phases and to simplify the
calculations involving the energy-dependent hopping amplitudes. Lastly, an
analytical approximation for the energy bands was obtained by slightly
modifying the Kronig–Penney-like dispersion. We show that this approximation
describes the energy bands well if the energy gaps are closed, specially
around Dirac points. We believe the findings of our work have the potential to
facilitate future studies on topological plasmonics and provide a reliable
framework to investigate similar polaritonic crystals analytically.

## Acknowledgements

D. A. M. acknowledges the project PTDC/FIS-MAC/2045/2021 for a research grant.
The Center for Polariton-driven Light–Matter Interactions (POLIMA) is funded
by the Danish National Research Foundation (Project No. DNRF165). N. M. R. P.
acknowledges support by the Portuguese Foundation for Science and Technology
(FCT) in the framework of the Strategic Funding UIDB/04650/2020, COMPETE 2020,
PORTUGAL 2020, FEDER, and through project PTDC/FIS-MAC/2045/2021. N. M. R. P.
also acknowledges the Independent Research Fund Denmark (grant no.
2032-00045B) and the Danish National Research Foundation (Project No.
DNRF165).

##  Appendix A Approximations for gated and ungated momenta

In the non-retarded (electrostatic) limit electromagnetic fields in graphene
and surrounding media are defined by Poisson equation

$\epsilon_{0}{\rm div}\left[\epsilon\left(z\right){\rm
grad}\phi\left(x,z\right)\right]=-\rho_{2D}\delta\left(z\right),$ (40)

where $\phi\left(x,z\right)$ is the scalar potential,

$\epsilon\left(z\right)=\left\\{\begin{array}[]{c}\epsilon_{1},\qquad z<0\\\
\epsilon_{2},\qquad z>0\end{array}\right.$ (41)

is the dielectric permeability, and $\rho_{2D}$ is the two-dimensional charge
density in graphene. In the right-hand part of Eq. (40), the delta-function
describes the fact that charges in graphene are concentrated on a two-
dimensional plane, arranged at $z=0$. Beyond the graphene plane, the Poisson
equation can be solved separately in the spatial domains $z>0$ and $z<0$,
where the dielectric permeabilities are homogeneous, and the solution will
have the form

$\phi\left(x,z,t\right)=\left\\{\begin{array}[]{c}\phi\left(0,0,0\right)e^{iq_{u}x-i\omega
t}e^{q_{u}z},\qquad z<0\\\ \phi\left(0,0,0\right)e^{iq_{u}x-i\omega
t}e^{-q_{u}z},\qquad z>0\end{array}\right.$ (42)

Here the scalar potential is continuous at $z=0$, and chose of signs in the
exponential terms, containing $\pm q_{u}z$ guarantees that the surface wave is
evanescent in the $z$ direction, and decays with increasing distance from the
graphene plane at $z=0$. Also solution (42) implies that surface plasmons
propagate along the $x$-axis with the wavevector $q_{u}>0$ and is oscillating
in time with frequency $\omega$.

If the Poisson equation (40) is integrated along $z$-coordinate in the limits
$z\in\left[-0,+0\right]$, it transforms into

$\epsilon_{2}\left.\frac{\partial\phi}{\partial
z}\right|_{z=+0}-\epsilon_{1}\left.\frac{\partial\phi}{\partial
z}\right|_{z=-0}=-\frac{\rho_{2D}}{\epsilon_{0}}.$ (43)

Substituting Eq. (42) into Eq. (43), we obtain

$q_{u}\left(\epsilon_{2}+\epsilon_{1}\right)\phi\left(0,0,0\right)e^{iq_{u}x-i\omega
t}=\frac{\rho_{2D}}{\epsilon_{0}}.$ (44)

The correspondent two-dimensional charge density can be obtained from the
continuity equation

$\frac{\partial\rho_{2D}}{\partial t}\delta\left(z\right)+{\rm div}\,{\bf
j}=0,$ (45)

where three-dimensional current j in our particular case possesses only
$x$-component

$j_{x}=-\sigma\left(\omega\right)\delta\left(z\right)\frac{\partial\phi}{\partial
x}.$ (46)

where
$\sigma(\omega)=\frac{4iE_{F}}{\hbar(\omega+i\gamma)}\alpha\epsilon_{0}c$ is
the Drude 2D conductivity of graphene, the dominant conductivity term from THz
to mid-IR [44]. Here $E_{F}$ and $\gamma$ stand for the Fermi energy and
inverse scattering time, respectively. In this paper we consider the case,
where losses in graphene are absent ($\gamma\equiv 0$). Substituting current
density from Eq. (46) into the continuity equation (45) and integrating it in
time, we obtain the 2D charge density

$\rho_{2D}=\frac{q_{u}^{2}\sigma\left(\omega\right)}{i\omega}\phi\left(0,0,0\right)e^{iq_{u}x-i\omega
t}.$ (47)

This equation, being substituted into Eq. (43), gives the dispersion relation
for ungated plasmons

$q_{u}=\frac{\hbar(\epsilon_{1}+\epsilon_{2})\omega^{2}}{4cE_{F}\alpha}.$ (48)

If metal is present at a distance $d$, the scalar potential turns put to be
more complicated

$\displaystyle\phi\left(x,z,t\right)=\phi\left(0,0,0\right)e^{iq_{g}x-i\omega
t}\times$ $\displaystyle\left\\{\begin{array}[]{c}e^{q_{g}z},\qquad z<0\\\
\frac{{\rm sinh}(q_{g}[d-z])}{{\rm sinh}(q_{g}d)},\qquad 0<z<d\\\ 0,\qquad
z>d\end{array}\right.$ (52)

In this case the 2D charge density will be the same as that in Eq. (47)
[except $q_{u}$ should be substituted by $q_{g}$]. Substitution of the 2D
charge density and the particular solution of Poisson equation (52) into Eq.
(43) results in the dispersion relation for gated plasmons

$\epsilon_{1}+\frac{\epsilon_{2}}{{\rm
tanh}(q_{g}d)}=\frac{q_{a}\sigma(\omega)}{i\omega\epsilon_{0}}.$ (53)

For small distances between graphene and metal, $dq<<1$ , ${\rm
tanh}(q_{g}d)\simeq q_{g}d$ and $\epsilon_{1}\ll\epsilon_{2}/(q_{g}d)$. In
this case the dispersion relation for gated plasmons can be expressed as

$q_{g}=\omega\sqrt{\frac{\hbar\epsilon_{2}}{4E_{F}c\alpha d}}.$ (54)

##  Appendix B  Derivation of energy-dependent hopping amplitudes

###  B.1 Transfer matrix

We begin by determining the general transfer matrix relating, in one specific
region labeled by $j=u,g$, the functions at $x$ and $x+\Delta x$:

$\left(\begin{array}[]{c}\phi_{j}(x+\Delta x)\\\
\sigma_{j}\phi^{\prime}_{j}(x+\Delta
x)\end{array}\right)=\mathbf{W}_{j}(\Delta
x)\left(\begin{array}[]{c}\phi_{j}(x)\\\
\sigma_{j}\phi^{\prime}_{j}(x)\end{array}\right)$ (55)

where $\phi^{\prime}_{j}(x)=\frac{\partial\phi_{j}(x)}{\partial x}$ and
$\mathbf{W}_{j}(\Delta x)$ is the corresponding transfer matrix in the region
labeled by $j$. The transfer-matrix coefficients can be found by expanding the
expression above, using the definitions of $\phi(x)$ in Eq. (2b):

$\begin{array}[]{c}\phi_{j}^{+}e^{iq_{j}x}e^{iq_{j}\Delta
x}+\phi_{j}^{-}e^{-iq_{j}x}e^{-iq_{j}\Delta x}=\mathbf{W}_{j}(\Delta
x)_{11}\left(\phi_{j}^{+}e^{iq_{j}x}+\phi_{j}^{-}e^{-iq_{j}x}\right)+\mathbf{W}_{j}(\Delta
x)_{12}i\sigma_{j}q_{j}\left(\phi_{j}^{+}e^{iq_{j}x}-\phi_{j}^{-}e^{-iq_{j}x}\right)\\\
\\\ i\sigma_{j}q_{j}\left(\phi_{j}^{+}e^{iq_{j}x}e^{iq_{j}\Delta
x}-\phi_{j}^{-}e^{-iq_{j}x}e^{-iq_{j}\Delta x}\right)=\mathbf{W}_{j}(\Delta
x)_{21}\left(\phi_{j}^{+}e^{iq_{j}x}+\phi_{j}^{-}e^{-iq_{j}x}\right)+\mathbf{W}_{j}(\Delta
x)_{22}i\sigma_{j}q_{j}\left(\phi_{j}^{+}e^{iq_{j}x}-\phi_{j}^{-}e^{-iq_{j}x}\right)\end{array}$
(56)

Rearranging terms related to the same coefficients:

$\begin{array}[]{c}\phi_{j}^{+}e^{iq_{j}x}e^{iq_{j}\Delta
x}+\phi_{j}^{-}e^{-iq_{j}x}e^{-iq_{j}\Delta
x}=\phi_{j}^{+}e^{iq_{j}x}\left(\mathbf{W}_{j}(\Delta
x)_{11}+\mathbf{W}_{j}(\Delta
x)_{12}i\sigma_{j}q_{j}\right)+\phi_{j}^{-}e^{-iq_{j}x}\left(\mathbf{W}_{j}(\Delta
x)_{11}-\mathbf{W}_{j}(\Delta x)_{12}i\sigma_{j}q_{j}\right)\\\ \\\
\phi_{j}^{+}e^{iq_{j}x}e^{iq_{j}\Delta
x}-\phi_{j}^{-}e^{-iq_{j}x}e^{-iq_{j}\Delta
x}=\phi_{j}^{+}e^{iq_{j}x}\left(\frac{\mathbf{W}_{j}(\Delta
x)_{21}+\mathbf{W}_{j}(\Delta
x)_{22}i\sigma_{j}q_{j}}{i\sigma_{j}q_{j}}\right)-\phi_{j}^{-}e^{-iq_{j}x}\left(\frac{\mathbf{W}_{j}(\Delta
x)_{22}i\sigma_{j}q_{j}-\mathbf{W}_{j}(\Delta
x)_{21}}{i\sigma_{j}q_{j}}\right)\end{array}$ (57)

The exponentials of $\Delta x$ must then obey:

$\begin{array}[]{c}e^{iq_{j}\Delta x}=\mathbf{W}_{j}(\Delta
x)_{11}+\mathbf{W}_{j}(\Delta
x)_{12}i\sigma_{j}q_{j}=\frac{\mathbf{W}_{j}(\Delta
x)_{21}+\mathbf{W}_{j}(\Delta x)_{22}i\sigma_{j}q_{j}}{i\sigma_{j}q_{j}}\\\
\\\ e^{-iq_{j}\Delta x}=\mathbf{W}_{j}(\Delta x)_{11}-\mathbf{W}_{j}(\Delta
x)_{12}i\sigma_{j}q_{j}=\frac{\mathbf{W}_{j}(\Delta
x)_{22}i\sigma_{j}q_{j}-\mathbf{W}_{j}(\Delta
x)_{21}}{i\sigma_{j}q_{j}}\end{array}$ (58)

Expanding the exponentials into cosine and sine and solving the system of
equations, it follows that for arbitrary $\Delta x$ the transfer-matrix
coefficients are:

$\begin{array}[]{cc}\mathbf{W}_{j}(\Delta x)_{11}={\rm cos}(q_{j}\Delta
x),&\mathbf{W}_{j}(\Delta x)_{12}=\frac{{\rm sin}(q_{j}\Delta
x)}{\sigma_{j}q_{j}}\\\ \mathbf{W}_{j}(\Delta
x)_{21}{\rm=-\sigma_{j}q_{j}sin}(q_{j}\Delta x),&\mathbf{W}_{j}(\Delta
x)_{22}={\rm cos}(q_{j}\Delta x)\end{array}$ (59)

###  B.2 $\sigma_{g}\phi^{\prime}_{g}(x)$ as a function of
$\phi_{g}(nl),\phi_{g}(nl+l_{g})$:

For the gated region in the $n^{th}$ unit cell we define:

$\left(\begin{array}[]{c}\phi_{g}(nl+l_{g})\\\
\sigma_{g}\phi^{\prime}_{g}(nl+l_{g})\end{array}\right)=\mathbf{W}_{g}(l_{g})\left(\begin{array}[]{c}\phi_{g}(nl)\\\
\sigma_{g}\phi^{\prime}_{g}(nl)\end{array}\right)$ (60)
$\left(\begin{array}[]{c}\phi_{g}(x)\\\
\sigma_{g}\phi^{\prime}_{g}(x)\end{array}\right)=\mathbf{W}_{g}(x-nl)\left(\begin{array}[]{c}\phi_{g}(nl)\\\
\sigma_{g}\phi^{\prime}_{g}(nl)\end{array}\right)$ (61)

Using the first row of the transfer matrix, Eq. (60) leads to

$\phi_{g}(nl+l_{g})={\rm cos}(q_{g}l_{g})\phi_{g}(nl)+\frac{{\rm
sin}(q_{g}l_{g})}{\sigma_{g}q_{g}}\sigma_{g}\phi^{\prime}_{g}(nl)$ (62)

which can be rearranged as

$\sigma_{g}\phi^{\prime}_{g}(nl)=\frac{\sigma_{g}q_{g}\phi_{g}(nl+l_{g})}{{\rm
sin}(q_{g}l_{g})}-\frac{{\rm\sigma_{g}q_{g}cos}(q_{g}l_{g})\phi_{g}(nl)}{{\rm
sin}(q_{g}l_{g})}.$ (63)

Using the second row of Eq. (61) leads to:

$\sigma_{g}\phi^{\prime}_{g}(x)={\rm-\sigma_{g}q_{g}sin}(q_{g}\left[x-nl\right])\phi_{g}(nl)+{\rm
cos}(q_{g}\left[x-nl\right])\sigma_{g}\phi^{\prime}_{g}(nl)$ (64)

which can be substituted into Eq. (63). After some manipulation and using
trigonometric identities, this yields:

$\sigma_{g}\phi^{\prime}_{g}(x)=-\phi_{g}(nl)\frac{\sigma_{g}q_{g}{\rm
cos}(q_{g}\left[l_{g}+nl-x\right])}{{\rm
sin}(q_{g}l_{g})}+\phi_{g}(nl+l_{g})\frac{\sigma_{g}q_{g}{\rm
cos}(q_{g}\left[x-nl\right])}{{\rm sin}(q_{g}l_{g})}.$ (65)

###  B.3 $\sigma_{u}\phi^{\prime}_{u}(x)$ as a function of
$\phi_{u}(nl+l),\phi_{u}(nl+l_{g})$:

For the ungated region in the $n$th unit cell we define:

$\left(\begin{array}[]{c}\phi_{u}(nl+l)\\\
\sigma_{u}\phi^{\prime}_{u}(nl+l)\end{array}\right)={\bf
W}_{u}(l_{u})\left(\begin{array}[]{c}\phi_{u}(nl+l_{g})\\\
\sigma_{u}\phi^{\prime}_{u}(nl+l_{g})\end{array}\right),$ (66)
$\left(\begin{array}[]{c}\phi_{u}(x)\\\
\sigma_{u}\phi^{\prime}_{u}(x)\end{array}\right)={\bf W}_{u}(x-nl-
l_{g})\left(\begin{array}[]{c}\phi_{u}(nl+l_{g})\\\
\sigma_{u}\phi^{\prime}_{u}(nl+l_{g})\end{array}\right).$ (67)

This situation is exactly like the one for the gated region, but with the
changes

$\begin{array}[]{rcl}{\rm cos}(q_{g}l_{g}),{\rm
sin}(q_{g}l_{g})&\longrightarrow&{\rm cos}(q_{u}l_{u}),{\rm
sin}(q_{u}l_{u})\\\ \phi_{g}(nl)&\longrightarrow&\phi_{u}(nl+l_{g})\\\
\phi_{g}(nl+l_{g})&\longrightarrow&\phi_{u}(nl+l)\\\ {\bf
W}_{g}(x-nl)&\longrightarrow&{\bf W}_{u}(x-nl-l_{g})\end{array}$

With these substitutions, we find

$\sigma_{u}\phi^{\prime}_{u}(x)=-\phi_{u}(nl+l_{g})\frac{\sigma_{u}q_{u}{\rm
cos}(q_{u}\left[nl+l-x\right])}{{\rm
sin}(q_{u}l_{u})}+\phi_{u}(nl+l)\frac{\sigma_{u}q_{u}{\rm
cos}(q_{u}\left[x-nl-l_{g}\right])}{{\rm sin}(q_{u}l_{u})}.$ (68)

###  B.4 The tight-binding equation

With Eqs. (65) and (68) we managed to write, for a given region, a function
describing the plasmon in an arbitrary position $x$ in terms of functions at
the boundaries. Now, we equate both Eqs. (65) and (68) at $x=nl$ and
$x=nl+l_{g}$ to obtain the tight-binding equations (remember that the indexes
$u$ and $g$ are irrelevant on the wavefunctions now): At $x=nl$ (note that at
this point, $\sigma_{g}\phi^{\prime}_{g}(x)$ is evaluated in the $n$th unit
cell while $\sigma_{u}\phi^{\prime}_{u}(x)$ is evaluated at the $(n-1)$th unit
cell):

$\sigma_{g}\phi^{\prime}_{g}(nl)=\sigma_{u}\phi^{\prime}_{u}([n-1]l)$ (69)
$-\phi(nl)\frac{\sigma_{g}q_{g}{\rm cos}(q_{g}l_{g})}{{\rm
sin}(q_{g}l_{g})}+\phi(nl+l_{g})\frac{\sigma_{g}q_{g}}{{\rm
sin}(q_{g}l_{g})}=-\phi_{u}([n-1]l+l_{g})\frac{\sigma_{u}q_{u}}{{\rm
sin}(q_{u}l_{u})}+\phi_{u}([n-1]l+l)\frac{\sigma_{u}q_{u}{\rm
cos}(q_{u}\left[x-(n-1)l-l_{g}\right])}{{\rm sin}(q_{u}l_{u})}$ (70)
$\phi(nl+l_{g})\frac{\sigma_{g}q_{g}}{{\rm sin}(q_{g}l_{g})}+\phi(nl-
lu)\frac{\sigma_{u}q_{u}}{{\rm
sin}(q_{u}l_{u})}=\phi(nl)\left[\frac{\sigma_{u}q_{u}{\rm
cos}(q_{u}l_{u})}{{\rm sin}(q_{u}l_{u})}+\frac{\sigma_{g}q_{g}{\rm
cos}(q_{g}l_{g})}{{\rm sin}(q_{g}l_{g})}\right]$ (71)

At $x=nl+l_{g}$:

$\sigma_{g}\phi^{\prime}_{g}(nl+l_{g})=\sigma_{u}\phi^{\prime}_{u}(nl+l_{g})$
(72) $-\phi(nl)\frac{\sigma_{g}q_{g}}{{\rm
sin}(q_{g}l_{g})}+\phi(nl+l_{g})\frac{\sigma_{g}q_{g}{\rm
cos}(q_{g}l_{g})}{{\rm
sin}(q_{g}l_{g})}=-\phi(nl+l_{g})\frac{\sigma_{u}q_{u}{\rm
cos}(q_{u}l_{u})}{{\rm sin}(q_{u}l_{u})}+\phi(nl+l)\frac{\sigma_{u}q_{u}}{{\rm
sin}(q_{u}l_{u})}$ (73)

where we used the first representation of $\sigma_{u}\phi^{\prime}_{u}(x)$.

$\phi(nl)\frac{\sigma_{g}q_{g}}{{\rm
sin}(q_{g}l_{g})}+\phi(nl+l)\frac{\sigma_{u}q_{u}}{{\rm
sin}(q_{u}l_{u})}=\phi(nl+l_{g})\left[\frac{\sigma_{u}q_{u}{\rm
cos}(q_{u}l_{u})}{{\rm sin}(q_{u}l_{u})}+\frac{\sigma_{g}q_{g}{\rm
cos}(q_{g}l_{g})}{{\rm sin}(q_{g}l_{g})}\right]$ (74)

With Eqs. (71) and (74) we have the two desired tight-binding equations.

## References

  * Polman and Atwater [2005] A. Polman and H. A. Atwater, Materials Today 8, 56 (2005).
  * Perenboom, Wyder, and Meier [1981] J. A. A. J. Perenboom, P. Wyder, and F. Meier, Physics Reports 78, 173 (1981).
  * Bashevoy  _et al._ [2006] M. V. Bashevoy, F. Jonsson, A. V. Krasavin, N. I. Zheludev, Y. Chen, and M. I. Stockman, Nano Letters 6, 1113 (2006).
  * Mie [1908] G. Mie, Annalen der Physik 330, 377 (1908).
  * Fan, Zheng, and Singh [2014] X. Fan, W. Zheng, and D. J. Singh, Light: Science $\&$ Applications 3, e179 (2014).
  * Khurgin [2015] J. B. Khurgin, Nature Nanotechnology 10, 2 (2015).
  * Zeng  _et al._ [2014] S. Zeng, D. Baillargeat, H.-P. Ho, and K.-T. Yong, Chemical Society Reviews 43, 3426 (2014).
  * Yu  _et al._ [2019] H. Yu, Y. Peng, Y. Yang, and Z.-Y. Li, npj Computational Materials 5, 45 (2019).
  * Gadelha  _et al._ [2021] A. C. Gadelha, D. A. Ohlberg, C. Rabelo, E. G. Neto, T. L. Vasconcelos, J. L. Campos, J. S. Lemos, V. Ornelas, D. Miranda, R. Nadas, F. C. Santana, K. Watanabe, T. Taniguchi, B. van Troeye, M. Lamparski, V. Meunier, V.-H. Nguyen, D. Paszko, J.-C. Charlier, L. C. Campos, L. G. Cançado, G. Medeiros-Ribeiro, and A. Jorio, Nature 590, 405 (2021).
  * Jablan, Buljan, and Soljačić [2009] M. Jablan, H. Buljan, and M. Soljačić, Phys. Rev. B 80, 245435 (2009).
  * Grigorenko, Polini, and Novoselov [2012] A. N. Grigorenko, M. Polini, and K. S. Novoselov, Nature Photonics 6, 749 (2012).
  * Ni  _et al._ [2018] G. Ni, A. McLeod, Z. Sun, L. Wang, L. Xiong, K. W. Post, S. S. Sunku, B.-Y. Jiang, J. Hone, C. R. Dean, M. M. Fogler, and D. N. Basov, Nature 557, 530 (2018).
  * Fei  _et al._ [2012] Z. Fei, A. S. Rodin, G. O. Andreev, W. Bao, A. S. McLeod, M. Wagner, L. M. Zhang, Z. Zhao, M. Thiemens, G. Dominguez, M. M. Fogler, A. H. Castro Neto, C. N. Lau, F. Keilmann, and D. N. Basov, Nature 487, 82–85 (2012).
  * Woessner  _et al._ [2015] A. Woessner, M. B. Lundeberg, Y. Gao, A. Principi, P. Alonso-González, M. Carrega, K. Watanabe, T. Taniguchi, G. Vignale, M. Polini, J. Hone, R. Hillenbrand, and F. H. L. Koppens, Nature Materials 14, 421 (2015).
  * Brar  _et al._ [2013] V. W. Brar, M. S. Jang, M. Sherrott, J. J. Lopez, and H. A. Atwater, Nano Letters 13, 2541 (2013).
  * Koppens, Chang, and García de Abajo [2011] F. H. L. Koppens, D. E. Chang, and F. J. García de Abajo, Nano Letters 11, 3370 (2011).
  * Rodrigo  _et al._ [2015] D. Rodrigo, O. Limaj, D. Janner, D. Etezadi, F. J. García de Abajo, V. Pruneri, and H. Altug, Science 349, 165 (2015).
  * Vasić, Isić, and Gajić [2013] B. Vasić, G. Isić, and R. Gajić, Journal of Applied Physics 113, 013110 (2013).
  * Li  _et al._ [2014] Y. Li, H. Yan, D. B. Farmer, X. Meng, W. Zhu, R. M. Osgood, T. F. Heinz, and P. Avouris, Nano Letters 14, 1573 (2014).
  * Dai  _et al._ [2020] Y. Dai, Z. Zhou, A. Ghosh, R. S. K. Mong, A. Kubo, C.-B. Huang, and H. Petek, Nature 588, 616 (2020).
  * Gong  _et al._ [2021] Y. Gong, L. Guo, S. Wong, A. J. Bennett, and S. S. Oh, Scientific Reports 11, 1055 (2021).
  * Kim and Rho [2020] M. Kim and J. Rho, Nanophotonics 9, 3227 (2020).
  * Süsstrunk and Huber [2015] R. Süsstrunk and S. D. Huber, Science 349, 47 (2015).
  * Nash  _et al._ [2015] L. M. Nash, D. Kleckner, A. Read, V. Vitelli, A. M. Turner, and W. T. M. Irvine, Proceedings of the National Academy of Sciences 112, 14495 (2015).
  * Vila, Pal, and Ruzzene [2017] J. Vila, R. K. Pal, and M. Ruzzene, Physical Review B 96, 134307 (2017).
  * Chen  _et al._ [2018] H. Chen, H. Nassar, A. N. Norris, G. K. Hu, and G. L. Huang, Physical Review B 98, 094302 (2018).
  * Shi  _et al._ [2021] X. Shi, I. Kiorpelidis, R. Chaunsali, V. Achilleos, G. Theocharis, and J. Yang, Physical Review Research 3, 033012 (2021).
  * Zhang  _et al._ [2018] D.-W. Zhang, Y.-Q. Zhu, Y. X. Zhao, H. Yan, and S.-L. Zhu, Advances in Physics 67, 253 (2018).
  * Wintersperger  _et al._ [2020] K. Wintersperger, C. Braun, F. N. Ünal, A. Eckardt, M. D. Liberto, N. Goldman, I. Bloch, and M. Aidelsburger, Nature Physics 16, 1058 (2020).
  * Konig  _et al._ [2007] M. Konig, S. Wiedmann, C. Brune, A. Roth, H. Buhmann, L. W. Molenkamp, X.-L. Qi, and S.-C. Zhang, Science 318, 766 (2007).
  * Fu, Kane, and Mele [2007] L. Fu, C. L. Kane, and E. J. Mele, Physical Review Letters 98, 106803 (2007).
  * Chen  _et al._ [2009] Y. L. Chen, J. G. Analytis, J.-H. Chu, Z. K. Liu, S.-K. Mo, X. L. Qi, H. J. Zhang, D. H. Lu, X. Dai, Z. Fang, S. C. Zhang, I. R. Fisher, Z. Hussain, and Z.-X. Shen, Science 325, 178 (2009).
  * Su, Schrieffer, and Heeger [1979] W. P. Su, J. R. Schrieffer, and A. J. Heeger, Physical Review Letters 42, 1698 (1979).
  * Zhang  _et al._ [2023] L. Zhang, X.-M. Wang, X.-M. Qiu, Z. Wang, and J.-Y. Yan, Physical Review B 108, 085402 (2023).
  * Jin  _et al._ [2017] D. Jin, T. Christensen, M. Soljačić, N. X. Fang, L. Lu, and X. Zhang, Physical Review Letters 118, 245301 (2017).
  * Ju  _et al._ [2011] L. Ju, B. Geng, J. Horng, C. Girit, M. Martin, Z. Hao, H. A. Bechtel, X. Liang, A. Zettl, and F. Wang, Nature Nanotechnology 6, 630 (2011).
  * Nikitin  _et al._ [2012] A. Y. Nikitin, F. Guinea, F. J. Garcia-Vidal, and L. Martin-Moreno, Physical Review B 85, 081405(R) (2012).
  * Strait  _et al._ [2013] J. H. Strait, P. Nene, W.-M. Chan, C. Manolatou, S. Tiwari, F. Rana, J. W. Kevek, and P. L. McEuen, Physical Review B 87, 241410(R) (2013).
  * Zhao and Zhang [2015] B. Zhao and Z. M. Zhang, ACS Photonics 2, 1611 (2015).
  * Dias and Peres [2017] E. J. C. Dias and N. M. R. Peres, ACS Photonics 4, 3071 (2017).
  * Bylinkin  _et al._ [2019] A. Bylinkin, E. Titova, V. Mikheev, E. Zhukova, S. Zhukov, M. Belyanchikov, M. Kashchenko, A. Miakonkikh, and D. Svintsov, Physical Review Applied 11, 054017 (2019).
  * Rappoport  _et al._ [2021] T. G. Rappoport, Y. V. Bludov, F. H. L. Koppens, and N. M. R. Peres, ACS Photonics 8, 1817 (2021).
  * Guo and Argyropoulos [2023] T. Guo and C. Argyropoulos, Journal of Applied Physics 134, 050901 (2023).
  * Gonçalves and Peres [2016] P. A. D. Gonçalves and N. M. R. Peres,  _An Introduction To Graphene Plasmonics_ (World Scientific Publishing Company, 2016).
  * Menabde  _et al._ [2022] S. G. Menabde, J. T. Heiden, J. D. Cox, N. A. Mortensen, and M. S. Jang, Nanophotonics 11, 2433 (2022).
  * Kronig and Penney [1931] R. D. L. Kronig and W. G. Penney, Proceedings of the Royal Society of London A 130, 499 (1931).
  * Asbóth, Oroszlány, and Pályi [2016] J. Asbóth, L. Oroszlány, and A. Pályi,  _A Short Course on Topological Insulators: Band Structure and Edge States in One and Two Dimensions_, Lecture Notes in Physics (Springer, 2016).
  * Iranzo  _et al._ [2018] D. A. Iranzo, S. Nanot, E. J. C. Dias, I. Epstein, C. Peng, D. K. Efetov, M. B. Lundeberg, R. Parret, J. Osmond, J.-Y. Hong, J. Kong, D. R. Englund, N. M. R. Peres, and F. H. L. Koppens, Science 360, 291 (2018).
  * Markos and Soukoulis [2008] P. Markos and C. Soukoulis,  _Wave Propagation: From Electrons to Photonic Crystals and Left-Handed Materials_ (Princeton University Press, 2008).