# Extended unitarity and absence of skin effect in periodically driven systems

Aditi Chakrabarty<EMAIL_ADDRESS>Department of Physics and
Astronomy, National Institute of Technology, Rourkela, Odisha-769008, India
Sanjoy Datta<EMAIL_ADDRESS>Department of Physics and Astronomy,
National Institute of Technology, Rourkela, Odisha-769008, India

###### Abstract

One of the most striking features of non-Hermitian quasiperiodic systems with
arbitrarily small asymmetry in the hopping amplitudes and open boundaries is
the accumulation of all the bulk eigenstates at one of the edges of the
system, termed in literature as the skin effect, below a critical strength of
the potential. In this Letter, we uncover that a time-periodic drive in such
systems can eliminate the SE up to a finite strength of this asymmetry.
Remarkably, the critical value for the onset of SE is independent of the
driving frequency and approaches to the static behavior in the thermodynamic
limit. We find that the absence of SE is intricately linked to the emergence
of extended unitarity in the delocalized phase, providing dynamical stability
to the system. Interestingly, under periodic boundary condition, our non-
Hermitian system can be mapped to a Hermitian analogue in the large driving
frequency limit that leads to the extended unitarity irrespective of the
hopping asymmetry and the strength of the quasiperiodic potential, in stark
contrast to the static limit. Additionally, we numerically verify that this
behavior persists Based on our findings, we propose a possible experimental
realization of our driven system, which could be used as a switch to control
the light funneling mechanism.

Introduction- In recent years, non-Hermitian Hamiltonians have garnered
widespread attention due to their ability to accurately mimic experimental
procedures involving interactions with the environment. They also offer
possibilities for exotic phases of quantum matter that are typically absent in
their Hermitian counterparts. Unlike Hermitian systems, the reality of
eigenenergies and stable unitary dynamics is guaranteed only in specific
classes of non-Hermitian Hamiltonians that possess either
$\mathcal{PT}$-symmetry [1, 2] or pseudo-Hermiticity [3, 4].
Among the non-Hermitian systems, the paradigmatic Hatano-Nelson (HN)
Hamiltonian with asymmetric hopping amplitudes [5, 6] stands out as another
class of non-Hermitian systems of particular interest. This interest stems
from a unique phenomenon known as the skin effect (SE), which refers to the
exponential localization of all bulk eigenstates at the edges of a lattice
with an open boundary [7, 8, 9]. Interestingly, the spectral behavior in these
systems is drastically sensitive to the choice of boundary conditions [10].
Moreover, recent studies have illustrated the existence of SE accompanying a
delocalization-localization (DL) phase transition in such non-Hermitian
quasicrystals [11].
In parallel, the exploration of a diversified range of periodically driven
systems has been triggered due to the presence of rich and intriguing features
typically absent in their temporally static counterparts [12, 13, 14, 15, 16].
The underlying quantum mechanics and the general understanding of such systems
influenced by some external time-periodic drive (generally known as Floquet
systems) have gained interest in recent years due to their applications in
ultrafast spintronics [17, 18], quantum optics [19], ultra-cold atomic systems
[20, 21], and trapped ions [22]. On the other hand, the interplay between
periodic driving and non-Hermiticity has led to several novel findings in
recent years [23, 24, 25, 26, 27, 28, 29, 30]. It has been demonstrated that a
time-periodic drive can induce SE in systems with an on-site loss [31].
Additionally, Floquet engineering has been utilized to control the direction
of edge modes [32]. Furthermore, it has been pointed out that a time-periodic
drive can stabilize the dynamics of a two-level non-Hermitian Rabi model [33],
where the Floquet quasienergies turn out to be real, leading to extended
unitarity (EU) at the end of each complete driving period.
In this Letter, we unveil for the first time that a time-periodic drive can
lead to the EU condition in HN quasicrystals, with a concurrent disappearance
of the SE up to a finite strength of the non-reciprocity in the hopping
amplitudes. In the limit of high frequency in the drive, we illustrate
analytically that such a time-periodic modulation reduces the non-Hermitian
system to a Hermitian equivalent counterpart under the periodic boundary
condition (PBC), giving rise to completely real Floquet quasienergies
irrespective of the strength of the potential and the asymmetry in the hopping
amplitude, at each stroboscopic time period. This is in stark contrast to
static HN quasicrystals, where the energy spectrum undergoes a complex to real
transition at a critical value of the potential, although both these systems
manifest a DL transition. Surprisingly, however, under an open boundary
condition (OBC), the EU condition persists only up to a critical strength of
the non-reciprocity in the delocalized regime. Interestingly, we demonstrate
that the emergence of EU destroys the SE. Moreover, contrary to the static
limit, where the energies are real in the SE phase, we find that in time-
periodic systems, it can exist when the Floquet quasienergy spectrum is
complex. These counter-intuitive results on the relationship between the SE
and the eigenenergies completely alter our understanding of conventional HN
systems.
In practice, from the perspective of experimentalists, it is noteworthy that
the dynamics in such systems has found wide realizations in photonics and in
quantum systems by enhancing the extent of optical sensing [34, 35].
Additionally, the dynamics in a Hamiltonian with asymmetric hopping has
recently been studied experimentally in photonic lattices to focus an incident
light at a desired location, irrespective of the point of excitation, a
phenomenon termed as light funneling [36]. In this Letter, we propose an
experimental setup that can exploit our findings to tune the light funneling
effect.
Model and methods- We consider a time-dependent version of the single-particle
quasiperiodic non-Hermitian Hamiltonian [37, 11] defined as,

$\displaystyle\mathcal{H}(t)=\displaystyle\sum_{n=1}^{N-1}(Je^{h\text{cos}(\omega
t)}c^{\dagger}_{n+1}c_{n}+Je^{-h\text{cos}(\omega t)}c^{\dagger}_{n}c_{n+1})$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}+\sum_{n=1}^{N}V\text{cos}(2\pi\alpha
n)c^{\dagger}_{n}c_{n},~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}$ (1)

where $c^{\dagger}_{n}$ and $c_{n}$ denote the fermionic creation and
annihilation operators, $N$ represents the number of sites in the lattice,
where $n$ is the site index. The lattice size $L=Na$, where $a$ is the lattice
period of translation (considered to be 1 in arb. units). In the prototypical
static HN Hamiltonian, the asymmetric hopping of the fermions towards the left
and the right is incorporated in the tight-binding notation using an imaginary
magnetic vector potential by the terms $e^{-h}$ and $e^{h}$ respectively. In
this Letter, we consider continuous temporally cosine modulated non-reciprocal
hopping amplitudes as indicated in the first term of the above Hamiltonian.
The second term characterizes the on-site quasiperiodic potential of the
Aubry-André-Harper (AAH) type, where $\alpha$ defines the incommensurability,
set as $(\sqrt{5}-1)/2$ throughout this work. It is evident that after a
stroboscopic period $T$, $\mathcal{H}(t)=\mathcal{H}(t+T)$.
For the time-periodic Hamiltonian as given in Eq. 1, the Floquet theory has
been instrumental in determining the states after a time $T$. According to the
Floquet theory, the Floquet propagator for one complete period and an initial
time $t_{0}=0$ is defined as,

$\displaystyle
U(T,0)=\mathcal{T}e^{-(i/\hbar)\int_{0}^{T}\mathcal{H}(t)dt}=e^{-i\mathcal{H}_{F}T/\hbar},$
(2)

where $\mathcal{T}$ takes care of the time-ordering of the Hamiltonians at
different instants of time. In the above equation, $\mathcal{H}_{F}$ is the
Floquet Hamiltonian, whose eigenstates can be obtained by the exact
diagonalization of $U(T,0)$. The reduced Planck’s constant ($\hbar$) is
considered to be of unit magnitude throughout this work. In general, $U(T,0)$
is non-unitary when the elemental Hamiltonian is non-Hermitian, and can be
constructed using the biorthogonal formalism [38].

---
Figure 1: Projected $\left\langle\text{BIPR}\right\rangle$ as a function of
$h$ and $V$ for the driven HN system at different driving frequencies: (a)
$\omega=\pi/2$, (b) $\omega=\pi$, (c) $\omega=2\pi$, and (d) $\omega=4\pi$.
The blue region in all the phase diagrams indicate the delocalized phase. The
red solid line depicts the critical point for DL transition in the static HN
Hamiltonian. The black dotted lines manifests the critical value under the
drive as obtained analytically from the expression in Eq. 10. The DL
transition is demonstrated for a lattice with 144 sites and under the PBC. The
white data markers superposed in all the phase diagrams separate the regions
with and without the SE, obtained under OBC. The spectral behavior in the
different regimes under both PBC and OBC are also indicated.

The eigenvectors and eigenvalues of $U(T,0)$ on exact diagonalization yield
the eigenspectrum of $\mathcal{H}_{F}$ given as,

$\displaystyle
U(T,0)=\displaystyle\sum_{n}E_{n}\ket{\psi_{nR}}\bra{\psi_{nL}},\text{and}~{}E_{n}=e^{-i\epsilon_{n}^{F}T},$
(3)

where $E_{n}$, $\ket{\psi_{nR}}$ and $\ket{\psi_{nL}}$ are the eigenenergies
and the right and left eigenvectors respectively. The Floquet quasienergies
$\epsilon_{n}^{F}$ satisfy the relation
$\mathcal{H}_{F}{\psi_{n}}=\epsilon_{n}^{F}{\psi_{n}}$ (defined modulo
$\hslash\omega$). The solutions to the time dependent Hamiltonians to
ascertain the physical properties can mostly be attained numerically using the
Floquet eigenstates and quasienergies as described above.
To identify the delocalized and localized phases in the system, the concept of
Inverse Participation Ratio (IPR) [39, 40] is widely used in the literature.
The concept of IPR has been extended in non-Hermitian systems, where a new
measure of bidirectional-IPR(BIPR) for an eigenstate labelled ‘$j$’ has been
introduced recently [41] and is defined as,

$\displaystyle
BIPR_{j}=\frac{\displaystyle\sum_{n=1}^{N}|\psi_{nL}^{j}\psi_{nR}^{j}|^{2}}{\Big{(}\displaystyle\sum_{n=1}^{N}|{\psi_{nL}^{j}\psi_{nR}^{j}}|\Big{)}^{2}}~{}~{}~{}~{}~{},$
(4)

where the sum of the weights of the wavevectors is over all the lattice sites
indicated by $n$. $\left\langle BIPR\right\rangle$ signifies the average of
BIPR over all the eigenstates. The $\left\langle BIPR\right\rangle$ is
$\mathcal{O}(L^{-1})$ when the states are completely delocalized, and
$\mathcal{O}(1)$ in the localized regime. For all the subsequent details and
the findings on our periodically driven system, we have considered $J=1$ (in
arb. units), $L=144$, and the Trotter time step as $\Delta t=0.001$, unless
specifically stated.

---
Figure 2: The Floquet quasienergy spectrum under the PBC (in red) and OBC (in
blue) for the driven system at a low value of the imaginary vector potential,
i.e., $h=0.1$ and $V=0.5$ in (a) $\omega=\pi/2$ and (c) $\omega=4\pi$. (b,d)
depict lattice-site resolved $|\psi|^{2}$ for all the eigenstates
corresponding to the two parameters in (a) and (c) respectively, indicating
the absence of SE under the OBC. $L=144$ in all the cases.
---
Figure 3: The Floquet quasienergy spectrum under the PBC (in red) and OBC (in
blue) for the driven system at $h=0.2$ and $V=0.5$ in (a) $\omega=\pi/2$ and
(c) $\omega=4\pi$. (b,d) demonstrates the lattice-site resolved $|\psi|^{2}$
for different eigenstates for the two parameters corresponding to (a) and (c),
demonstrating the presence of SE under the OBC. The lattice consists of $144$
sites in all the cases.

Results and discussions- The aim of this work is achieved in two steps. We
first explore the behavior of DL transition of the time-dependent HN system
under PBC in the parameter space of $h$ and $V$ for different driving
frequencies. In the subsequent analysis, we assess the sensitivity of the
energy spectrum to the boundary conditions alongwith the investigation of the
existence of SE.
Figs. 1(a-d) illustrate the phase diagram of the time-periodic version of the
HN Hamiltonian as discussed in Eq. 1 under the PBC. In sharp contrast to the
static system, interestingly, we observe from Figs. 1(a-d) that the DL
transition (indicated by a transition from blue to a different color) does not
follow the critical value of the undriven system (indicated by the solid red
line) determined from the condition $V_{c}=\text{max}[Je^{h},Je^{-h}]$, as
demonstrated in Ref. [37]. From the phase diagrams, it is evident that the
transition from delocalized to localized states strongly depends on the
driving frequency.
In the following, we demonstrate that the critical point of the DL transition
under the PBC can be obtained analytically in the regime of a large frequency
of the drive. The analytical expression of the phase boundary is derived by
retrieving an effective Floquet Hamiltonian upon time-averaging the original
Hamiltonian as discussed in earlier works [42, 43, 44]. The effective
Hamiltonian, which can only be obtained when the frequency of the drive is
large enough, is written as follows,

$\displaystyle\mathcal{H}_{F}=\frac{1}{T}\int_{0}^{T}\mathcal{H}(t)dt.$ (5)

To evaluate $\mathcal{H}_{F}$ as given in Eq. 5, we first consider the hopping
towards the right, which is given by,

$\displaystyle\frac{1}{T}\int_{0}^{T}\displaystyle\sum_{n}(Je^{h\text{cos}(\omega
t)}c^{\dagger}_{n+1}c_{n})dt.$ (6)

In the above equation $\int_{0}^{T}e^{h\text{cos}(\omega t)}dt$ can be
identified as the generating function of the modified Bessel equation, i.e.,
$\int e^{\frac{z}{2}(x+\frac{1}{x})}dx$ with $z=h$ and $x=\text{exp}(i\omega
t)$. Integration over the entire time period yields the effective Floquet
Hamiltonian for this part, and is given as,

$\displaystyle\displaystyle\sum_{n}\Big{(}J\mathcal{I}_{0}(h)\Big{)}c^{\dagger}_{n+1}c_{n},$
(7)

where $\mathcal{I}_{0}$ is the zeroth order modified Bessel function. It is
interesting to note that the above expression is independent of $\omega$.
Similarly, for the hopping in the opposite direction (towards the left), the
other part of the effective Hamiltonian can be written as,

$\displaystyle\displaystyle\sum_{n}\Big{(}J\mathcal{I}_{0}(-h)\Big{)}c^{\dagger}_{n}c_{n+1}.$
(8)

It is easy to see that $\mathcal{I}_{0}(h)=\mathcal{I}_{0}(-h)$. The effective
Floquet Hamiltonian can then be explicitly written as,

$\displaystyle\mathcal{H}_{F}=\displaystyle\sum_{n}\Big{(}J\mathcal{I}_{0}(h)\Big{)}(c^{\dagger}_{n+1}c_{n}+c^{\dagger}_{n}c_{n+1})+$
$\displaystyle\sum_{n}V\text{cos}(2\pi\alpha
n)c^{\dagger}_{n}c_{n}.~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}$ (9)

Remarkably, from Eq. 9, it becomes quite evident that the effective Floquet
Hamiltonian for the HN systems with a drive in the magnetic vector potential
reduces to a Hermitian AAH Hamiltonian with a rescaled hopping amplitude given
as $J\mathcal{I}_{0}(h)$ in the regime of high frequency. Thus, we expect the
DL transition to occur at a critical value of the quasiperiodic potential
determined by the self-duality of the AAH Hamiltonian in the Hermitian limit,
expressed as,

$\displaystyle
V_{c}=2J^{\prime},~{}~{}\text{where}~{}~{}J^{\prime}=J\mathcal{I}_{0}(h).$
(10)

In Figs. 1(a-d), we have indicated the analytically obtained phase boundary in
Eq. 10 by a black dotted line. It is clear that as we approach to a higher
value of the driving frequency ($\omega=4\pi$), the numerically determined
critical value for the DL transition agrees excellently with the analytical
result in the entire parameter space. However, for lower driving frequencies,
the numerical critical points follow the analytical expression only upto a
certain strength of the asymmetry. In addition, the effect of the periodic
drive in such systems is to shift the DL transition to a lower value of
$V_{c}$, as compared to its static counterpart (demonstrated by the red line
in Figs. 1(a-d)). Moreover, our analytical result suggests that the Floquet
quasienergies in the entire parameter space of $h$ and $V$ should be real for
high driving frequencies under PBC, due to the exact mapping of the original
non-Hermitian system to the Hermitian AAH Hamiltonian. We have numerically
verified this assertion. Suprisingly, however, the reality of the
quasienergies persists for all frequencies of the drive. This is in complete
contrast to the static HN Hamiltonian under the PBC, where the DL transition
is concurrent with a spectral transition from complex to real. In our case,
the existence of the real quasienergies irrespective of the strength of the
potential (refer to Fig. 6 of the supplemental material for more details) and
non-Hermiticity, which is the hallmark for EU, suggests that our system
becomes dynamically stable after a stroboscopic driving period.
In non-Hermitian systems a change in the boundary condition from PBC to OBC
drastically alters the behaviour of electronic states and eigenenergies,
especially in the delocalized phase. In the static HN counterparts, it has
been well demonstrated that under the OBC, in the delocalized regime, the
eigenenergies lie on the real axis [9], alongwith the appearance of SE.
Naturally, one of the most important question that arises is whether such a
correspondence between the quasienergies and SE exist under the time-periodic
drive. We address this question in the subsequent discussion.

Figure 4: System-size behavior of $h_{c}$ for the onset of SE as a function of
inverse lattice size in different regimes of the Floquet drive, i.e., at
$\omega=\pi/2$ (green), $\omega=4\pi$ (blue) and $\omega=8\pi$ (red). The data
are obtained for $L$= 55, 89, 144 and 233 under the OBC.
---

In analogy with the static system, one can naively anticipate that the SE
would persist in the entire delocalized regime (represented in blue in Figs
1(a-d)) retaining the reality of the Floquet quasienergies. However, we
demonstrate that this expectation is in complete contradiction when the system
is driven. From Figs. 2 (a,c) it is evident that for weak asymmetry in the
hopping, the periodic drive leads to real Floquet quasienergies under both the
periodic and open boundary conditions. However, remarkably, we find that under
the OBC, the time-periodic drive demolishes the SE regardless of the driving
frequency (Figs. 2 (b,d)), contrary to the undriven systems where the SE
exists for an arbitrarily weak strength of the imaginary magnetic vector
potential.
To further understand the dependence of the SE on $h$, we consider a slightly
greater amplitude of the asymmetricity in the Hamiltonian. In this case, the
quasienergy spectrum under the PBC still exhibits real eigenvalues, similar to
a Hermitian system, as demonstrated in Figs. 3 (a,c), manifesting EU as
indicated previously. This is in accordance with the finding in Eq. 9.
However, with an increase in the value in $h$, the spectrum under OBC lies on
the complex plane and does not satisfy the EU condition. Surprisingly, when
the spectrum becomes complex under OBC, the system exhibits SE irrespective of
the driving frequency, as illustrated in Figs. 3 (b,d). This is in stark
contrast to the static systems. From these findings, it is clear that there
exists a critical value $h_{c}$ at which the SE appears. We have numerically
determined $h_{c}$ in the entire parameter space of $h$ and $V$. Our findings
are summarized in Figs. 1(a-d) for different frequencies of the drive. It is
remarkable that $h_{c}$ is independent of the driving frequency and the
quasiperiodic potential. The numerical value of $h_{c}$ is estimated to be
$h\simeq 0.20$ with an error of $\pm 0.01$. The preceding discussions clearly
indicate that the emergence of EU and disappearance of SE are closely linked
in such driven systems. The determination of $h_{c}$ for the different driving
frequencies has been elucidated in Fig. 7 of the supplemental material.
It is important to note that the critical strength of the asymmetricity for
the onset of SE in such driven systems depends strongly on the size of the
lattice. In order to get a comprehensive idea on this system size dependence,
in Fig. 4 we present $h_{c}$ as a function of $1/L$ for different driving
frequencies. It is evident that $h_{c}$ scales inversely with the system size.
In addition, $h_{c}$ becomes independent of the driving frequency when the
size of the lattice is large enough. Furthermore, it can be easily observed
that $h_{c}$ approaches to the static limit of HN Hamiltonian for
$L\to\infty$.

Figure 5:  Schematic of light funneling using two oppositely directed static
HN chains. The interface is formed at $n=0$. The chains on the left and right
side of the interface have greater unidirectionality towards the right and
left directions respectively. The output light traverses through the funnel-
like structure where it is collected.
---

Switch for light funneling- In this discussion, we propose a possible
experimental implementation of the system considered in this work from the
point view of its application in controlling the light funneling effect.
Similar to the recent realizations of the non-reciprocal lattices in photonic
systems as demonstrated in Refs. [45, 36], we consider a set-up with two
optical fibers with anisotropic hopping in either directions, seperated by an
interface as illustrated in Fig. 5. The driven HN system can be replicated by
changing the anisotropy in the hopping with the help of a frequency dependent
beam-splitter [46], and selectively using the desired beam of light. Such an
unidirectionality in the two fibers causes the incident light that is impinged
on the lattice to be pushed towards the interface due to SE, resulting in a
funneling effect of the output light. Since the value of $h$ can be tuned with
the anisotropic beam splitter and our results suggest the existence of $h_{c}$
for the onset of SE, it is easy to understand that the SE can be turned on or
off by controlling the beam splitter, effectively acting as a switch to tune
the light-funneling effect.
Conclusions- In conclusions, this Letter demonstrates that a time-periodic
drive can fundamentally alter the behaviour of the paradigmatic quasiperiodic
HN systems. By introducing a drive, we unfold that under PBC, the Floquet
Hamiltonian in the large frequency limit becomes equivalent to a Hermitian AAH
counterpart, giving rise to the extended unitarity and real Floquet
quasienergies in the entire paramter space of the Hamiltonian, in stark
contrast to the static limit. We find that the driven system undergoes a DL
transition similar to the static case, albeit at a lower strength of the
potential for a given asymmtetry in the hopping, determined by the self-
duality condition of the effective Hermitian AAH Hamiltonian. Under OBC,
however, we find that the EU condition survives only up to a critical value of
the asymmetry in the hopping, $h_{c}$, along with the disappearance of SE
within that asymmetry. In complete contrast to the static limit, the SE
appears above this critical value with complex Floquet quasienergies.
Remarkably, for a given system size, $h_{c}$ is found to be independent of the
driving frequency and quasiperiodic potential. Furthermore, we find that
$h_{c}$ scales inversely with the system size approaching to the static
behaviour in the thermodynamic limit. Our work deepens the understanding of
the DL transition under the PBC, Floquet quasienergies and its connection to
the SE under the OBC in HN systems. Finally, we provide an experimental set-up
that can exploit the findings of this Letter to control the light-funelling
mechanism.
A.C. acknowledges the financial support received from the Council of
Scientific $\&$ Industrial Research (CSIR)-HRDG, India, via File No.-
09/983(0047)/2020-EMR-I. The data computation in this work were carried out in
the cluster procured from SERB (DST), India (Grant No. EMR/2015/001227) and
using the High Performance Computing (HPC) facilities of the National
Institute of Technology Rourkela.

## References

  * Bender and Boettcher [1998] C. M. Bender and S. Boettcher, Real spectra in non-hermitian hamiltonians having pt symmetry, Phys. Rev. Lett. 80, 5243 (1998).
  * Bender  _et al._ [2002] C. M. Bender, M. V. Berry, and A. Mandilara, Generalized pt symmetry and real spectra, Journal of Physics A: Mathematical and General 35, L467 (2002).
  * Mostafazadeh [2002] A. Mostafazadeh, Pseudo-hermiticity versus symmetry: The necessary condition for the reality of the spectrum of a non-hermitian hamiltonian, J. Math. Phys. 43, 205 (2002).
  * Kawabata  _et al._ [2019] K. Kawabata, K. Shiozaki, M. Ueda, and M. Sato, Symmetry and topology in non-hermitian physics, Phys. Rev. X 9, 041015 (2019).
  * Hatano and Nelson [1996] N. Hatano and D. R. Nelson, Localization transitions in non-hermitian quantum mechanics, Phys. Rev. Lett. 77, 570 (1996).
  * Hatano and Nelson [1998] N. Hatano and D. R. Nelson, Non-hermitian delocalization and eigenfunctions, Phys. Rev. B 58, 8384 (1998).
  * Lee [2016] T. E. Lee, Anomalous edge state in a non-hermitian lattice, Phys. Rev. Lett. 116, 133903 (2016).
  * Lee  _et al._ [2019] C. H. Lee, L. Li, and J. Gong, Hybrid higher-order skin-topological modes in nonreciprocal systems, Phys. Rev. Lett. 123, 016805 (2019).
  * Okuma  _et al._ [2020] N. Okuma, K. Kawabata, K. Shiozaki, and M. Sato, Topological origin of non-hermitian skin effects, Phys. Rev. Lett. 124, 086801 (2020).
  * Xiong [2018] Y. Xiong, Why does bulk boundary correspondence fail in some non-hermitian topological models, Journal of Physics Communications 2, 035043 (2018).
  * Longhi [2021] S. Longhi, Phase transitions in a non-hermitian aubry-andré-harper model, Phys. Rev. B 103, 054203 (2021).
  * Lignier  _et al._ [2007] H. Lignier, C. Sias, D. Ciampini, Y. Singh, A. Zenesini, O. Morsch, and E. Arimondo, Dynamical control of matter-wave tunneling in periodic potentials, Phys. Rev. Lett. 99, 220403 (2007).
  * Ivanov  _et al._ [2008] V. V. Ivanov, A. Alberti, M. Schioppo, G. Ferrari, M. Artoni, M. L. Chiofalo, and G. M. Tino, Coherent delocalization of atomic wave packets in driven lattice potentials, Phys. Rev. Lett. 100, 043602 (2008).
  * Poli  _et al._ [2011] N. Poli, F.-Y. Wang, M. G. Tarallo, A. Alberti, M. Prevedelli, and G. M. Tino, Precision measurement of gravity with cold atoms in an optical lattice and comparison with a classical gravimeter, Phys. Rev. Lett. 106, 038501 (2011).
  * Chen  _et al._ [2011] Y.-A. Chen, S. Nascimbène, M. Aidelsburger, M. Atala, S. Trotzky, and I. Bloch, Controlling correlated tunneling and superexchange interactions with ac-driven optical lattices, Phys. Rev. Lett. 107, 210405 (2011).
  * Bitter and Milner [2016] M. Bitter and V. Milner, Experimental observation of dynamical localization in laser-kicked molecular rotors, Phys. Rev. Lett. 117, 144104 (2016).
  * Oka and Kitamura [2019] T. Oka and S. Kitamura, Floquet engineering of quantum materials, Annual Review of Condensed Matter Physics 10, 387 (2019).
  * Losada  _et al._ [2019] J. M. Losada, A. Brataas, and A. Qaiumzadeh, Ultrafast control of spin interactions in honeycomb antiferromagnetic insulators, Phys. Rev. B 100, 060410 (2019).
  * Eckhardt  _et al._ [2022] C. J. Eckhardt, G. Passetti, M. Othman, C. Karrasch, F. Cavaliere, M. A. Sentef, and D. M. Kennes, Quantum floquet engineering with an exactly solvable tight-binding chain in a cavity, Communications Physics 5, 122 (2022).
  * Meinert  _et al._ [2016] F. Meinert, M. J. Mark, K. Lauber, A. J. Daley, and H.-C. Nägerl, Floquet engineering of correlated tunneling in the bose-hubbard model with ultracold atoms, Phys. Rev. Lett. 116, 205301 (2016).
  * Sandholzer  _et al._ [2022] K. Sandholzer, A.-S. Walter, J. Minguzzi, Z. Zhu, K. Viebahn, and T. Esslinger, Floquet engineering of individual band gaps in an optical lattice using a two-tone drive, Phys. Rev. Res. 4, 013056 (2022).
  * He  _et al._ [2021] R. He, M.-Z. Ai, J.-M. Cui, Y.-F. Huang, Y.-J. Han, C.-F. Li, G.-C. Guo, G. Sierra, and C. E. Creffield, Riemann zeros from floquet engineering a trapped-ion qubit, npj Quantum Information 7, 109 (2021).
  * Della Valle and Longhi [2013] G. Della Valle and S. Longhi, Spectral and transport properties of time-periodic $\mathcal{PT}$-symmetric tight-binding lattices, Phys. Rev. A 87, 022119 (2013).
  * Liu  _et al._ [2022] C.-H. Liu, H. Hu, and S. Chen, Symmetry and topological classification of floquet non-hermitian systems, Phys. Rev. B 105, 214305 (2022).
  * Zhou [2021] L. Zhou, Floquet engineering of topological localization transitions and mobility edges in one-dimensional non-hermitian quasicrystals, Phys. Rev. Res. 3, 033184 (2021).
  * Huang  _et al._ [2016] Y. Huang, Z.-q. Yin, and W. L. Yang, Realizing a topological transition in a non-hermitian quantum walk with circuit qed, Phys. Rev. A 94, 022302 (2016).
  * Koutserimpas and Fleury [2018] T. T. Koutserimpas and R. Fleury, Nonreciprocal gain in non-hermitian time-floquet systems, Phys. Rev. Lett. 120, 087401 (2018).
  * Höckendorf  _et al._ [2020] B. Höckendorf, A. Alvermann, and H. Fehske, Topological origin of quantized transport in non-hermitian floquet chains, Phys. Rev. Res. 2, 023235 (2020).
  * Banerjee and Narayan [2020] A. Banerjee and A. Narayan, Controlling exceptional points with light, Phys. Rev. B 102, 205423 (2020).
  * Banerjee and Sengupta [2023] T. Banerjee and K. Sengupta, Emergent conservation in the floquet dynamics of integrable non-hermitian models, Phys. Rev. B 107, 155117 (2023).
  * Ke  _et al._ [2023] S. Ke, W. Wen, D. Zhao, and Y. Wang, Floquet engineering of the non-hermitian skin effect in photonic waveguide arrays, Phys. Rev. A 107, 053508 (2023).
  * Li  _et al._ [2023] Y. Li, C. Lu, S. Zhang, and Y.-C. Liu, Loss-induced floquet non-hermitian skin effect, Phys. Rev. B 108, L220301 (2023).
  * Gong and Wang [2015] J. Gong and Q.-H. Wang, Stabilizing non-hermitian systems by periodic driving, Phys. Rev. A 91, 042135 (2015).
  * McDonald and Clerk [2020] A. McDonald and A. A. Clerk, Exponentially-enhanced quantum sensing with non-hermitian lattice dynamics, Nature Communications 11, 5382 (2020).
  * De Carlo  _et al._ [2022] M. De Carlo, F. De Leonardis, S. R.A., L. Colatorti, and V. Passaro, Non-hermitian sensing in photonics and electronics: A review., Sensors 22, 3977 (2022).
  * Weidemann  _et al._ [2020] S. Weidemann, M. Kremer, T. Helbig, T. Hofmann, A. Stegmaier, M. Greiter, R. Thomale, and A. Szameit, Topological funneling of light, Science 368, 311 (2020).
  * Jiang  _et al._ [2019] H. Jiang, L.-J. Lang, C. Yang, S.-L. Zhu, and S. Chen, Interplay of non-hermitian skin effects and anderson localization in nonreciprocal quasiperiodic lattices, Phys. Rev. B 100, 054301 (2019).
  * Brody [2013] D. C. Brody, Biorthogonal quantum mechanics, Journal of Physics A: Mathematical and Theoretical 47, 035305 (2013).
  * Evers and Mirlin [2000] F. Evers and A. D. Mirlin, Fluctuations of the inverse participation ratio at the anderson transition, Phys. Rev. Lett. 84, 3690 (2000).
  * Wessel and Milat [2005] S. Wessel and I. Milat, Quantum fluctuations and excitations in antiferromagnetic quasicrystals, Phys. Rev. B 71, 104427 (2005).
  * Wang  _et al._ [2019] P. Wang, L. Jin, and Z. Song, Non-hermitian phase transition and eigenstate localization induced by asymmetric coupling, Phys. Rev. A 99, 062112 (2019).
  * Goldman and Dalibard [2014] N. Goldman and J. Dalibard, Periodically driven quantum systems: Effective hamiltonians and engineered gauge fields, Phys. Rev. X 4, 031027 (2014).
  * Marin Bukov and Polkovnikov [2015] L. D. Marin Bukov and A. Polkovnikov, Universal high-frequency behavior of periodically driven systems: from dynamical stabilization to floquet engineering, Advances in Physics 64, 139 (2015).
  * Eckardt and Anisimovas [2015] A. Eckardt and E. Anisimovas, High-frequency approximation for periodically driven quantum systems from a floquet-space perspective, New Journal of Physics 17, 093039 (2015).
  * Longhi [2020] S. Longhi, Stochastic non-hermitian skin effect, Opt. Lett. 45, 5250 (2020).
  * Makarov [2021] D. N. Makarov, Theory of a frequency-dependent beam splitter in the form of coupled waveguides, Scientific Reports 11, 5014 (2021).

Supplemental material for “Extended unitarity and absence of skin effect in
periodically driven systems”

S.1. Existence of extended unitarity in the localized regime

Figure 6: The Floquet quasienergy spectrum in the complex plane under the PBC
(in red) and OBC (in blue) at $h=0.1$ and $V=2.5$ (localized regime) at: (a)
$\omega=\pi/2$, (b) $\omega=4\pi$. $L=144$ as considered in the main text.
---

Figure 7: $|\psi|^{2}$ for all the eigenstates projected over the lattice
sites at (a,e) $\omega=\pi/2$, (b,f) $\omega=\pi$ (c,g) $\omega=2\pi$, (d,h)
$\omega=4\pi$. The figures in the upper panels correspond to the value of $h$
just before the onset of SE, whereas the ones in the lower panel correspond to
the value of $h$ at which the SE appears (determined within an error of $\pm
0.01$). We have considered the OBC in a lattice with 144 sites.

In Figs. 1(a-d) of the main text, we have illustrated the phase diagrams of
the driven HN Hamiltonian considered in our work. It is well known that the
delocalization-localization phase transition in the static HN systems are
accompanied by a complex to real transition in the energy spectrum under the
periodic boundary condition (PBC). Under the open boundary condition (OBC),
such systems exhibit the skin effect (SE) as discussed in the main text.
However, as we have demonstrated that in the driven systems and under the PBC,
the Floquet quasienergies are always real in the delocalized regime, satisying
the extended unitarity (EU) condition as shown in Figs. 2(a,c) and Figs.
3(a,c) in the main text for $h=0.1$ and $h=0.2$ respectively. This is in
contrary to the behavior of static HN systems. Moreover, under the OBC, the
Floquet quasienergy spectra changes from real (in the presence(absence) of
EU(SE)) to complex (in the absence(presence) of EU(SE)). To assess the
behavior of the EU in the localized regime of the phase diagram, we have
presented the quasienergy spectrum under both the PBC and OBC in Figs. 6(a-b).
It is evident that for any frequency in the drive, the Floquet quasienergies
satisfy the EU condition in the localized regime of the phase diagram under
both the PBC and OBC. We have verified that this EU condition persists
irrespective of the strength of the quasiperiodic potential and asymmetric
hopping amplitude.

S.2. Demonstration of the SE in the driven HN systems

To find out the critical value of $h_{c}$ for the onset of SE, we resort to
the conventional method of determining whether all the eigenstates localize at
one end of the lattice using the wave-function probabilities as demonstrated
in Figs. 7(a-h). In an unidirectional system, all the eigenstates under the
OBC pile up at one of the edges towards which there is a greater
directionality of the fermionic hopping. Therefore, in our system, since the
amplitude of hopping towards the right dominates over the left hopping
amplitude, we say that our system possesses SE when all the states are
localized towards the right end. It is clear from the upper panel that Figs.
7(a-d), a few eigenstates are loocalized in the bulk of the system. Therefore
we do not consider them as skin modes. We have clearly illustrated the value
of $h_{c}$ at $0.20\pm 0.01$ for different driving frequencies in a lattice
with 144 sites as discussed in the main text.


# Emergent intersubband-plasmon-polaritons of Dirac electrons under one-
dimensional superlattices

Minwoo Jung<EMAIL_ADDRESS>Department of Physics, Cornell University,
Ithaca, New York, 14853, USA  Gennady Shvets<EMAIL_ADDRESS>School of
Applied and Engineering Physics, Cornell University, Ithaca, New York 14853,
USA

###### Abstract

Artifical superlattice (SL) potentials have been employed extensively for band
structure engineering of two-dimensional (2D) Dirac electron gas in graphene
Park  _et al._ (2008, 2009); Brey and Fertig (2009); Forsythe  _et al._
(2018); Li  _et al._ (2021). While such engineered electronic band structures
can modify optical or plasmonic properties of graphene, an emergent
polaritonic behavior beyond weak perturbative effects (e.g. anisotropic Drude
weights Brey  _et al._ (2020)) has not been reported. Here, we show that an
extreme modulation of one-dimensional (1D) SL potentials in monolayer graphene
deforms the underlying Dirac band dispersion and introduces ladder-like energy
levels near the Fermi surface, which result in emergent intersubband
polaritonic responses in optical conductivity. In our proposed system, hBN-
encapsulated graphene is placed on top of a 1D periodic metagate. In addition,
a backgate placed beneath the metagate is used as the second gate, further
modulating carrier density on regions in graphene that are not directly
screened by the metagate. With a strong carrier density modulation, graphene
electrons experience an array of deep potential wells, and at large enough
momenta perpendicular to the modulation direction, Dirac electrons are
waveguided via total internal reflections. These waveguided modes appear as
flat subbands with nearly equispaced energy levels. As a result, there arise
hybrid intersubband-polaritons with ultra-strong coupling in plasmonic
dispersions. Our study opens up an avenue for exploring emergent polaritons in
two-dimensional materials with gate-tunable electronic band structures.

††preprint: APS/123-QED

Polaritons in 2D semiconductor materials (e.g. plasmon polaritons in graphene
Woessner  _et al._ (2015); Alonso-Gonzalez  _et al._ (2017) or exciton-trion
polaritons in transitional metal dichalcogenides Sidler  _et al._ (2017);
Koksal  _et al._ (2021)) promise an ideal platform for novel opto-electronic
devices, supported by various ways to control the carrier densities in these
2D materials. Photo-carrier injection via high-intensity pulses is useful to
probe ultrafast transient responses of such polaritons Ni  _et al._ (2016);
Tan  _et al._ (2020); Sternbach  _et al._ (2020), and capacitive field-effect
gating is exploited for active control over steady-state responses. While a
uniform gate can tune overall polariton dispersions due to carrier density
dependence of Drude weight in graphene Woessner  _et al._ (2015); Alonso-
Gonzalez  _et al._ (2017) and exciton/trion oscillator strengths in
transitional metal dichalcogenides Rana  _et al._ (2021), more exotic electro-
optic controls can be achieved by using a metagate–a gating structure with
spatially varying capacitance. Especially in graphene plasmonics, spatial
modulation of carrier densities via metagate-tuning has been considered for
various applications such as local phase modulation Woessner  _et al._ (2017),
topological phase switching Fan  _et al._ (2019), and Bloch polariton steering
Xiong  _et al._ (2021).

Figure 1: SL engineering with double(metagate/backgate)-gating. a Schematic of
the device design; the inset on the exposed graphene plane is an example of a
SL potential ($h_{t}=5$nm, $h_{b}=10$nm, $h_{g}=10$nm, $h_{0}=150$nm,
$L=300$nm, $S=80$nm, $V_{M}=1$V and $V_{B}=-9$V). b Electrostatic simulation
showing the electric field along the direction perpendicular to graphene
sheet, $E_{z}$, and the electric field lines in $x-z$ plane; the role of the
backgate is to further modulate the carrier density locally on the region that
is not screened by the metagate. c The doped carrier density $n(x)$ (top), the
corresponding Fermi level $E_{F}(x)$ (middle), and the SL potential $U_{E}(x)$
(bottom) for different backgate voltages ($V_{M}=1$V is for all three cases);
dotted red: $V_{B}=12$V, dashed blue: $-6$V, and solid green: $-9$V.

Owing to quantum capacitance effect Fang  _et al._ (2007), spatially modulated
carrier densities $n(\mathbf{r})$ under electrochemical equilibrium gives rise
to SL electric potentials $U_{E}(\mathbf{r})$ for Dirac electrons, given as

$\hbar v_{F}\sqrt{\pi n(\mathbf{r})}+U_{E}(\mathbf{r})=\mu_{0},$ (1)

where $v_{F}$ is the Fermi velocity of Dirac electrons and $\mu_{0}$ is the
electrochemical potential. Under a periodic SL potential, the conical Dirac
dispersion deforms into minibands or subbands Park  _et al._ (2008, 2009);
Brey and Fertig (2009); Forsythe  _et al._ (2018); Li  _et al._ (2021). These
SL-induced subbands has been probed through emergent electronic transport
properties Forsythe  _et al._ (2018); Li  _et al._ (2021), but an emergent
optical or polaritonic phenomenon directly stemming from this band structure
deformation hasn’t been reported. Even though the above-mentioned works
Woessner  _et al._ (2017); Fan  _et al._ (2019); Xiong  _et al._ (2021)
utilized spatial modulation of $n(\mathbf{r})$, their phenomena are fully
captured at the level of a simple optical conductivity model with locally
varying Drude weights computed from the conical Dirac dispersion before SL-
induced deformation. We note that a recent study Brey  _et al._ (2020)
investigated plasmonic responses in graphene under a strict form of SL
potentials $U_{E}(x)\propto\cos{(G_{0}x)}$, and reported anisotropic plasmonic
dispersions and increased plasmon damping due to subband transitions. While
such marginal changes in graphene plasmons are certianly direct consequences
of SL engineering, these effects don’t necessarily introduce a new type of
polaritonic excitation.

In this study, we show that a 1D SL potential in a 2D Dirac electron gas leads
to coherent intersubband resonances in the optical conductivity. Along with
plasmon-polaritons, these newly emergent intersubband-polaritons form hybrid
intersubband-plasmon-polaritons (HIPPs). Each SL unit cell carries a square
potential well deep enough to host multiple tightly-confined bound states with
negligible coupling to bound states in adjacent unit cells, thereby forming
flat subbands in the deformed band structure. Notably, the massless nature of
Dirac electrons plays a pivotal role in the ultrastrong resonant enhancement
of oscillator strengths of the intersubband transition (ISBT) at certain
quantized frequencies, as the subbands are given at ladder-like equispaced
energy levels unlike the quadratic energy levels of massive particles in a
square potential well. Also, in our proposed system, the HIPP dispersions can
be systemically controlled via the combination of a metagate and a backgate,
where the metagate controls the underlying plasmon-polariton dispersions and
the backgate controls the intersubband resonance frequencies and the Rabi-
splitting strengths. Our results demonstrates that SL engineering in 2D
materials, a versatile technique to study exotic electronic transport
properties, can be used for the search of novel polaritonic materials as well.

Figure 1a depicts the geometry of the considered device platform. The
combination of a patterned metagate and an additional normal gate has been
used in several other works Li  _et al._ (2021); Fan  _et al._ (2019) to
attain better tunability of graphene carrier densities. We put the flat gate
as a backgate beneath the metagate, so that the backgate is used to control
the carrier densities only at a particular region on graphene that is not
screened by the metagate. In this way, if the duty cycle of air gaps in the
metagate grating is not so large ($S/L<50\%$, roughly speaking), the baseline
Drude weight ($\propto E_{F}=\hbar v_{F}\sqrt{\pi n}$) of graphene optical
conductivity is mostly determined by the metagate voltage $V_{M}$, and the SL
modulation depth is fine-tuned with the backgate voltage $V_{B}$. The
electrostatic simulation results shown in Fig. 1b-c are carried out by an
iterative solver in COMSOL Multiphysics to match the nonlinear boundary
condition given in Eq. (1) Jung  _et al._ (2018), and we consider the density-
dependent renormalization of the Fermi velocity:
$v_{F}(n)=[0.85+0.035\ln(n_{0}/n)]\times 10^{6}$m$/$s, where
$n_{0}=10^{15}$cm-2 Yu  _et al._ (2013). Figure 1b clearly shows that the
electric field lines from the backgate is mostly screened by the metagate, but
a few of them would penetrate through the air gap and reach to graphene. As a
result, the doped carrier density in the region above the metallic grating is
flat and controlled by $V_{M}$, whereas the region above the air gaps can even
be almost depleted when a large enough negative voltage $V_{B}$ is applied to
the backgate (Fig. 1c). Since we consider graphene to be grounded $\mu_{0}=0$,
the SL potential is given by simply negating the Fermi level
$U_{E}(x)=-E_{F}(x)$. The resulting shape of $U_{E}(x)$ can be viewed as a
periodic array of wide square potential wells with narrower separating walls,
where the height of separating walls is directly controlled by $V_{B}$.

Figure 2: Ladder-like flat subbands of Dirac electrons via TIREDD. a Left:
Subbands with $k_{y}=0$ appear as simply folded conical dispersions due to
Klein tunneling. Filled bands ($E<\mu_{0}=0$) are colored in blue, and
unoccupied bands are in green. The red dots are Dirac crossing points. Right:
With a sufficiently large value of $k_{y}$ and a strongly modulated SL
potential, several subbands around the Fermi surface appear as flat bands
along $k_{x}$. Not only that, the level spacing of these flat subbands is
almost uniform due to the nature of Dirac electrons that have linear
dispersion. b The energy diagram on top shows the energy levels of several
example bound states–(i) a mode right below the Fermi surface, (ii) the lowest
energy bound state formed by conduction band TIREDD, and (iii) the highest
energy bound state formed by valence band TIREDD–relative to the SL potential
landscape ($V_{B}=-9$V). The following graphs are the wavelength amplitudes
$\psi^{\dagger}\psi$ of these three example bound states. Each potential well
in each unit cell confines the bound states so tightly that the tunneling into
the adjacent unit cells is vanishing, which explains the appearance of the
flat subbands. c Schematics for the mechanism of TIREDD for the waveguiding of
conduction band electrons (top) and of valence band electrons (bottom). Each
circle represent the iso-energy contour of a Dirac cone in the momentum space,
and its radius $k_{r}(x)=|E-U_{E}(x)|/\hbar v_{F}$ is determined by the
difference of the eigenenergy level (dashed lines) and the SL potential (green
lines) at a local position. Electron waves can propagate along $x$-axis if
$k_{y}$ is less than $k_{r}$, but it evanescently decays along $x$-axis if
$k_{y}$ is greater than $k_{r}$. Thus, when an electron travels from a domain
with $k_{y}<k_{r}$ to another domain with $k_{y}>k_{r}$, there occurs TIREDD.
These TIREDD-based bound states can be formed both in potential wells (top;
for conduction band electron) and in potential barriers (bottom; for valence
band electron).

The electronic subband structures for Dirac electrons in 1D SL potentials
periodic along the $x$-axis are calculated by solving Dirac Hamiltonian
equation $\left(v_{F}\nu
p_{x}\sigma_{x}+v_{F}p_{y}\sigma_{y}+U_{E}\mathcal{I}\right)\ket{\mathbf{k};j}=E_{\mathbf{k};j}\ket{\mathbf{k};j}$
with the Bloch ansatz
$\braket{\mathbf{r}}{\mathbf{k};j}=e^{i\mathbf{k}\cdot\mathbf{r}}\psi_{\mathbf{k};j}(x)$

$\begin{bmatrix}U_{E}(x)&\hbar v_{F}\nu(k_{x}-i\partial_{x})-i\hbar
v_{F}k_{y}\\\ \hbar v_{F}\nu(k_{x}-i\partial_{x})+i\hbar
v_{F}k_{y}&U_{E}(x)\end{bmatrix}\psi_{\mathbf{k};j}(x)=E_{\mathbf{k};j}\psi_{\mathbf{k};j}(x),$
(2)

where $\nu$ is the valley index ($+1$ for $\mathbf{K}$-valley and $-1$ for
$\mathbf{K}^{\prime}$-valley), $j\in\mathbb{Z}$ is the subband index, and the
Bloch wavefunction satisfies the periodicity condition
$\psi_{\mathbf{k};j}(x)=\psi_{\mathbf{k};j}(x+L)$. Here,
$\mathbf{k}=k_{x}\hat{\mathbf{x}}+k_{y}\hat{\mathbf{y}}$ is the electron
wavevector, where $k_{x}$ is the Bloch wavenumber parallel to the periodic
direction and $k_{y}$ is the wavenumber perpendicular to the SL modulation.
The same band structure is repeated for each spin subspace. With vanishing
$k_{y}=0$, the conical linear dispersion is exactly preserved even in the
presence of strongly modulated $U_{E}(x)$ (Fig. 2a left) due to Klein
tunneling Katsnelson  _et al._ (2006); Young and Kim (2009). At sufficiently
large $k_{y}$, in contrast, there occur several flat subbands around and below
the Fermi surface (Fig. 2a right), as the bound states are formed inside the
potential wells via the total internal reflection of electrons with Dirac
dispersion (TIREDD). As clearly shown in Fig. 2b, these bound states are
tightly confined within each unit cell and their wavefunction amplitudes
$\psi^{\dagger}\psi$ decay rapidly to zero across the separating potential
walls so that Dirac electrons are waveguided along $y$-axis with negligible
tunneling across the adjacent wells. The negligible coupling between the bound
states in the adjacent unit cells manifests as flat subbands in the band
diagrams.

The TIREDD condition is met when a Dirac electron at a given energy $E$
travels with a nonzero $k_{y}$ from a domain, where the propagation along
$x$-axis is allowed, to another domain, where the propagation is not allowed.
At a local position along the SL potential, the iso-energy contour of the
conial Dirac dispersion is given as a circle in the momentum space with a
radius of $k_{r}(x)=|E-U_{E}(x)|/\hbar v_{F}$. Therefore, the local wavenumber
along $x$-axis $\sqrt{k_{r}(x)^{2}-k_{y}^{2}}$ becomes imaginary if $k_{y}$
exceeds $k_{r}(x)$, and the electron wavefunction decays exponentially along
$x$-axis as $\psi\propto e^{-\alpha x}$, where the decaying factor is given as
$\alpha=\sqrt{k_{y}^{2}-k_{r}^{2}}$. It is evident that this evanescent decay
along $x$-axis is not possible with $k_{y}=0$, which, along with the linear
nature of Dirac dispersion, leads to Klein tunneling as discussed earlier.
Notably, as illustrated in Fig. 2c, TIREDD occurs not only when the electron
travels from a lower potential domain to a higher potential domain–a more
familiar picture–but it also happens in the opposite case when the electron
travels from a higher potential domain to a lower potential domain. The latter
case corresponds to TIREDD of an antiparticle in the high-energy physics
language, or to TIREDD of the valence band electron in graphene. As a result,
bound states out of valence band electrons are formed around the potential
barriers, see the state (iii) in Fig. 2b. Each flat subband in the band
diagram in Fig. 2a corresponds to a TIREDD-based bound state, and the flatness
depends on whether the width of decaying barriers, the regions satisfying
$k_{y}>k_{r}(x)$, is sufficiently thicker than the inverse of decaying factor
$\alpha$. We note that the mean free path of electrons in hBN-encapsulated
graphene can be as high as $1\mu$m already at the room temperature and near
$10\mu$m at cryogenic temperatures $T<100$K Wang  _et al._ (2013). Thus, we
expect that the TIREDD condition can be realized in the realistic experiments,
since the mean free path in high-quality samples is much longer than the width
of the SL potential wells.

Figure 3: The massless dispersion of Dirac electrons makes the ISBT energy to
be nearly uniform over a broad region Fermi surface. The electronic subband
structure is shown along $k_{y}$ at a fixed $k_{x}=0$ ($V_{B}=-9V$ case). Each
vertical black bar is given as a guide to eyes for denoting a vertical
transition ($\Delta j=1$) from an occupied state below the Fermi surface to a
state above the Fermi surface, and all bars have have the same length. In this
case, the ISBT energy appears as very uniform roughly within
$6\pi/L<|k_{y}|<15\pi/L$, which is almost 45% of the whole area of the Fermi
surface $|k_{y}|<20\pi/L$.

Another notable feature of these TIREDD-based flat subbands is the ladder-like
energy level spacing around the Fermi surface. With a rough approximation (for
more detailed analysis, see the Supplementary Information), the bound state
energy levels of Dirac electrons in a square potential well (for now, let’s
consider the conduction band TIREDD only) are given as $E_{j}\sim\hbar
v_{F}\sqrt{(\pi j/W)^{2}+k_{y}^{2}}+U_{0}$, where $W$ is the width of
potential well, $U_{0}$ is the potential inside the well, and $j$ is the band
index. Even though this expression is not completely linear in $j$, it quickly
approaches to the asymptotic linear relation $E_{j+1}-E_{j}\sim\pi\hbar
v_{F}/W$ when $(E_{j}-U_{0})$ is only twice greater than $\hbar v_{F}k_{y}$.
Thus, the ISBT energy is maintained as nearly uniform over a substantial
portion of the Fermi surface (see Fig. 3), which resonantly enhances the
oscillator strength of the ISBT at certain quantized frequencies given at
integer multiples of $\pi v_{F}/W$. In Fig. 3, six or seven bands are
altogether contributing to the ISBT nearly at the same frequency.

Figure 4: Graphene optical conductivity and ISBT resonances. a Real part of
the conductivity
Re$[\sigma_{xx}(\mathbf{q}=q\hat{\mathbf{x}},\mathbf{q}^{\prime}=q^{\prime}\hat{\mathbf{x}};\omega)]$
calculated for $q=q^{\prime}$; in order to visualize the features at higher
frequencies better, we plotted $\omega\times\sigma_{xx}$. With strongly
modulated SL potentials ($V_{B}=-6$V, $-9$V), there appear several Lorenzian
peaks that correspond to the ISBT. $\Delta j$ refers to the difference between
subband indices of two bands involved in the ISBT. As a comparison, the
conductivity for uniformly doped ($E_{F}=0.15$eV) graphene is also provided. b
Both real (left) and imaginary (right) part of the conductivity calculated at
$q=q^{\prime}=0$ (top) and $q=q^{\prime}=\frac{\pi}{2L}$ (bottom); solid thin
black: no modulation ($E_{F}=0.15$eV), dotted red: $V_{B}=12$V, dashed blue:
$-6$V, and solid thick green: $-9$V.

Ultrastrong ISBT optical responses manifest as resonant features in graphene
optical conductivity, which is calculated by the Kubo formula under random
phase approximation (RPA) Brey  _et al._ (2020); Ramezanali  _et al._ (2009),
given as

$\frac{\sigma_{xx}(\mathbf{q},\mathbf{q}^{\prime};\omega)}{\sigma_{0}}=i\pi
g_{s}g_{\nu}\hbar^{2}v_{F}^{2}\int\frac{d^{2}\mathbf{k}}{(2\pi)^{2}}\int\frac{d^{2}\mathbf{k}^{\prime}}{(2\pi)^{2}}\sum_{j,j^{\prime}}\frac{f(E_{\mathbf{k};j})-f(E_{\mathbf{k}^{\prime};j^{\prime}})}{E_{\mathbf{k};j}-E_{\mathbf{k}^{\prime};j^{\prime}}}\frac{\bra{\mathbf{k};j}e^{-i\mathbf{q}\cdot\mathbf{r}}\sigma_{x}\ket{\mathbf{k}^{\prime};j^{\prime}}\bra{\mathbf{k}^{\prime};j^{\prime}}e^{i\mathbf{q}^{\prime}\cdot\mathbf{r}}\sigma_{x}\ket{\mathbf{k};j}}{\hbar(\omega+i\gamma)+E_{\mathbf{k};j}-E_{\mathbf{k}^{\prime};j^{\prime}}},$
(3)

where $\sigma_{0}=2e^{2}/h$ is the conductance quantum, $g_{s}=2$ and
$g_{\nu}=2$ are the spin and valley degeneracy,
$f(E)=1/\left[\exp(-E/k_{B}T)+1\right]$ is the Fermi-Dirac distribution
($\mu_{0}=0$), and each 2D momentum integration is done as $\int
d^{2}\mathbf{k}=\int_{-\pi/L}^{\pi/L}dk_{x}\int_{-\infty}^{\infty}dk_{y}$. To
avoid any confusion for notations, we use $q$ to denote the momentum of the
optical field, whereas we have used $k$ to denote the momentum of Dirac
electrons. Since our system is periodic in $x$-axis, the conductivity is
vanishing except when $\mathbf{q}-\mathbf{q}^{\prime}$ is an integer multiple
of $G_{0}\hat{\mathbf{x}}$ ($G_{0}=\frac{2\pi}{L}$); therefore, the surface
current density response is given as
$K_{x}(\mathbf{q})=\sum_{\mathbf{q}^{\prime}\in\\{\mathbf{q}+mG_{0}\hat{\mathbf{x}}|m\in\mathbb{Z}\\}}\sigma_{xx}(\mathbf{q},\mathbf{q}^{\prime})E_{x}(\mathbf{q}^{\prime})$.
In this work, we only consider the plasmonic excitations along $x-$axis
($\mathbf{q}=q\hat{\mathbf{x}}$,
$\mathbf{q}^{\prime}=q^{\prime}\hat{\mathbf{x}}$) with transverse-magnetic
polarization ($B_{x}=E_{y}=0$) Luo  _et al._ (2013). For the results plotted
in Fig. 4, we assumed a cryogenic temperature of $T=60$K and a plasmonic
scattering rate (Drude loss) of $\gamma=2\pi\times 0.2$THz, which are well
within the experimentally attainable ranges Ni  _et al._ (2018). We calculated
the conductivity $\sigma_{xx}$ for each of three SL potentials given in Fig.
1c and also for an unmodulated graphene with uniform doping of $E_{F}=0.15$eV
as a comparison. Even in the no-modulation case, a resonant behavior is found
along a linear line $v_{F}q=\omega$. This corresponds to the intraband
transitions that occur when the phase velocity of the optical field matches
with the Fermi velocity of Dirac electrons. This velocity-matching effect
becomes one of the key distinguishing factors between the simple Drude
conductivity model and the nonlocal RPA conductivity model in uniformly doped
graphene Lundeberg  _et al._ (2017). When the SL modulation depth is moderate
($V_{B}=12$V case), only the first ($\Delta j=1$; $\Delta j$ refers to the
difference between two subband indices) ISBT resonance appears faintly, while
the conductivity at higher frequency is almost the same to the no-modulation
case. With much stronger modulation depths, however, the first ISBT resonance
peak becomes the most prominent feature and the higher-order ($\Delta
j=2,3,...$) ISBT resonance peaks also become visible.

Figure 5: HIPP dispersion with ultra-strong coupling and far-field detection
of HIPPs. a Density of states or
$-\text{Im}\left[\text{Tr}\left([\epsilon(\mathbf{q},\omega)]^{-1}\right)\right]$
for visualizing the HIPP dispersion. The white dashed lines in $V_{B}=-9$V
panel are denoting the same frequencies of the ISBT resonances ($\Delta
j=1,2$, and 3) shown in Fig. 4a. b Reflection spectra for the normal incidence
of light polarized along $x$-axis; solid thin black: no modulation
($E_{F}=0.15$eV), dotted red: $V_{B}=12$V, dashed blue: $-6$V, and solid thick
green: $-9$V.

Figure 5a shows the resulting HIPP dispersions featuring the ultra-strong
coupling between the ISBT and the underlying plasmon-polaritons. If the system
is spatially homogeneous, the polariton eigenmodes in graphene appear as the
zeros of the scalar dynamical dielectric function
$\epsilon(q,\omega)=1-q^{2}\sigma(q,\omega)/i\omega C(q,\omega)$ Ramezanali
_et al._ (2009); Torre  _et al._ (2017). Here, $C(q,\omega)$ is the dynamical
capacitance of the system Xiong  _et al._ (2019), which connects the dynamic
carrier density oscillation $\delta n(q,\omega)$ and the dynamic electric
potential field on graphene $\delta U_{E}(q,\omega)$: $\delta
U_{E}=\frac{e^{2}}{C}\delta n$. In essence, the dynamic capacitance $C$
encodes the information about the dielectric environment around graphene, in
contrast to the conductivity $\sigma$ that encodes the dielectric property of
graphene itself. With a periodic modulation along $x$-axis like in our system,
the dynamical dielectric function is given as a matrix form

$[\epsilon(\mathbf{q},\omega)]_{m,m^{\prime}}=\delta_{m,m^{\prime}}-\sum_{l}C^{-1}(\mathbf{q}_{m},\mathbf{q}_{l};\omega)\frac{\mathbf{q}_{l}\cdot\mathbf{q}_{m^{\prime}}}{i\omega}\sigma_{xx}(\mathbf{q}_{l},\mathbf{q}_{m^{\prime}};\omega),$
(4)

where $m,m^{\prime},l\in\mathbb{Z}$ are integer indices,
$\mathbf{q}_{m}=\mathbf{q}+mG_{0}\hat{x}$ is the harmonic overtone of the
polariton Bloch wavevector $\mathbf{q}$, and $C^{-1}$ is the inverse dynamic
capacitance that governs a linear relation $\delta
U_{E}(\mathbf{q})=e^{2}\sum_{\mathbf{q}^{\prime}}C^{-1}(\mathbf{q},\mathbf{q}^{\prime})\delta
n(\mathbf{q}^{\prime})$ (for more detailed explanation, see the Supplementary
Information). Then, the polariton Bloch eigenmodes appear as the zeros of the
determinant of the dynamical dielectric function matrix Torre  _et al._
(2017); Jung  _et al._ (2018); Xiong  _et al._ (2021, 2019)—i.e. the matrix
$[\epsilon(\mathbf{q},\omega)]$ becomes non-invertible. Therefore, in Fig. 5a,
we plotted the density of states approximated as
$DOS(\mathbf{q},\omega)=-\text{Im}\left[\text{Tr}\left([\epsilon(\mathbf{q},\omega)]^{-1}\right)\right]$
to visualize the HIPP dispersion for polartions propagating along $x$-axis
($\mathbf{q}=q\hat{\mathbf{x}}$).

With a moderate depth of the SL modulation ($V_{B}=12$V case), the HIPP
dispersion is similar to the plasmon-polariton dispersion with no modulation,
and the ISBT feature is very subtle. As the SL modulation gets deeper, there
emerge several HIPP branches resulting from the hybridization between the
underlying plasmon-polartions and the ISBT resonances. With an extreme
modulation ($V_{B}=-9$V case), we observe a huge Rabi-splitting ($\sim 2$THz)
between the lowest branch and the second lowest, which is even comparable to
the ISBT frequency itself ($\sim 2$THz). In such an ultra-strong coupling
regime, a recent study reported that the electronic band structure of the
material could be modified in return due to the vacuum fluctuation of the
strongly-interacting polaritonic modes Kurman and Kaminer (2020). We believe
that our system would exhibit a similar behavior, but we didn’t consider such
additional corrections in this work. Figure 5b illustrates that this HIPP
phenomenon can be detected in the far-field reflection as well. Each peak in
the reflection spectra corresponds to the $q=0$ mode along the second lowest
HIPP branch. In the reflection calculation, we assumed that the backgate
substrate is silicon doped with a carrier density of $10^{15}$cm-3. The
diverging behavior at $\omega\rightarrow 0$ for no-modulation and $V_{B}=12$V
cases is due to the Drude response of the silicon backgate.

The emergent HIPPs found in our proposed system have several unique features
compared to the usual intersubband-polaritons or intersubband-plasmon-
polaritons studied in other platforms. First, the quantum well structure is
given along the direction of the polariton propagation. Accordingly, the ISBT
of our system occurs through in-plane electric fields along $x$-axis, which
allows the far-field detection even with normal incidence of light. In
contrast, conventional ISBT structures are based on vertical engineering of
quantum wells Kurman and Kaminer (2020); Kyriienko and Shelykh (2012); Zanotto
_et al._ (2010), and the optical coupling requires out-of-plane electric
fields. Second, as discussed earlier, the linear dispersion of Dirac electrons
allows multiple (6$\sim$7) bands with equi-spaced energy levels to resonantly
build up the ISBT strength. This resonant enhancement from multiple ladder-
like bands below the Fermi surface is still possible with a quadratic
dispersion of massive particles, but it will require a harmonic potential
instead of a square potential. However, in realistic material platforms, it
would remain as an extremely challenging task to engineer the quadratic shape
of the harmonic potential precisely enough to maintain the uniform energy
level spacing upto the 6th or 7th energy levels. Third, both the ISBT quanutm
well structures and the plasmon-polaritons that couple to the ISBT are hosted
simultaneously by monolayer graphene. This also contributes to the ultra-
strong coupling, since the plasmon-polariton field strengths are by nature
maximum at the plane of graphene. Similarly, in conventional vertically
confined quantum well structures, when the ground state subband is populated,
the 2D electron gas is naturally formed, hosting plasmon-polaritons confined
around it Kyriienko and Shelykh (2012). But, both the ISBT ground state
wavefunction and the plasmon-polariton fields have finite widths along
$z$-axis, unlike our system where the subband states are confined at an
atomically thin layer.

We emphasize again that the HIPPs shown in Fig. 5 operate in the ultra-
strongly coupled regime, featuring a giant Rabi-splitting that becomes
comparable to the ISBT frequency. In this regime, several quantum
electrodynamic phenomena can arise, such as material bandgap renormalization
Kurman and Kaminer (2020) or anti-resonant coupling that breaks the rotating
wave approximation and the Kubo conductivity formula Anappara  _et al._
(2009). Therefore, more precise determination of the HIPP dispersion would
require a full quantum description of the ultra-strong coupling physics, which
we leave as a future work. Lastly, even apart from the HIPP physics, the
TIREDD-induced ladder-like energy bands themselves can be useful for high
harmonic generations. Nonlinear optical responses can be resonantly enhanced
by engineering the equi-spaced energy level of subbands in quantum well
structures Lee  _et al._ (2014). As mentioned eariler, the linear dispersion
of Dirac electrons naturally ensures the equi-spaced energy levels of the
bound states in the square potential well. Thus, the 1D SL potential in
graphene can be also used as a novel material platform for nonlinear optics.

In conclusion, our study suggests that the SL engineering in 2D materials can
lead to the discovery of novel polariton phenomena emerging from the deformed
electronic band structures. The modified subband structure of Dirac electrons
under a 1D SL adds a completely new dimension to the polartion composition,
leading to the formation of the HIPPs. This emergent HIPP is easily tunable by
the double-gating scheme, provides a way to detect the SL-induced band
structure changes with a far-field optical measurement, and becomes suitable
for the study of quantum and nonlinear optics based on ultra-strong light-
matter interaction. Introducing a 2D SL Forsythe  _et al._ (2018) or
patterning other 2D materials beyond graphene kun Shi  _et al._ (2019) could
lead to more opportunities to study emergent polaritons with other novel
formation mechanisms.

###### Acknowledgements.

This work was supported by the Army Research Office (ARO) under a Grant No.
W911NF-16-1-0319, and by the National Science Foundation (NSF) under the
Grants No. DMR-1741788 and DMR-1719875. M.J. was also supported in part by
Cornell Fellowship and in part by the Kwanjeong Fellowship from Kwanjeong
Educational Foundation.

## References

  * Park  _et al._ (2008) C.-H. Park, L. Yang, Y.-W. Son, M. L. Cohen, and S. G. Louie, Nat. Phys. 4, 213–217 (2008).
  * Park  _et al._ (2009) C.-H. Park, Y.-W. Son, L. Yang, M. L. Cohen, and S. G. Louie, Phys. Rev. Lett. 103, 046808 (2009).
  * Brey and Fertig (2009) L. Brey and H. A. Fertig, Phys. Rev. Lett. 103, 046809 (2009).
  * Forsythe  _et al._ (2018) C. Forsythe, X. Zhou, K. Watanabe, T. Taniguchi, A. Pasupathy, P. Moon, M. Koshino, P. Kim, and C. R. Dean, Nat. Nanotechnol. 13, 566–571 (2018).
  * Li  _et al._ (2021) Y. Li, S. Dietrich, C. Forsythe, T. Taniguchi, K. Watanabe, P. Moon, and C. R. Dean, Nat. Nanotechnol. 16, 525–530 (2021).
  * Brey  _et al._ (2020) L. Brey, T. Stauber, L. Martín-Moreno, and G. Gómez-Santos, Phys. Rev. Lett. 124, 257401 (2020).
  * Woessner  _et al._ (2015) A. Woessner, M. B. Lundeberg, Y. Gao, A. Principi, P. Alonso-Gonzalez, M. Carrega, K. Watanabe, T. Taniguchi, G. Vignale, M. Polini, J. Hone, R. Hillenbrand, and F. H. L. Koppens, Nature Mater. 14, 421 (2015).
  * Alonso-Gonzalez  _et al._ (2017) P. Alonso-Gonzalez, A. Y. Nikitin, Y. Gao, A. Woessner, M. B. Lundeberg, A. Principi, N. Forcellini, W. Yan, S. Velez, A. J. Huber, K. Watanabe, T. Taniguchi, F. Casanova, L. E. Hueso, M. Polini, J. Hone, F. H. L. Koppens, and R. Hillenbrand, Nat. Nanotech. 12, 31 (2017).
  * Sidler  _et al._ (2017) M. Sidler, P. Back, O. Cotlet, A. Srivastava, T. Fink, M. Kroner, E. Demler, and A. Imamoglu, Nat. Phys. 13, 255 (2017).
  * Koksal  _et al._ (2021) O. Koksal, M. Jung, C. Manolatou, A. N. Vamivakas, G. Shvets, and F. Rana, Phys. Rev. Res. 3, 033064 (2021).
  * Ni  _et al._ (2016) G. X. Ni, L. Wang, M. D. Goldflam, M. Wagner, Z. Fei, A. S. McLeod, M. K. Liu, F. Keilmann, B. Özyilmaz, A. H. C. Neto, J. Hone, M. M. Fogler, and D. N. Basov, Nat. Photon. 10, 244 (2016).
  * Tan  _et al._ (2020) L. B. Tan, O. Cotlet, A. Bergschneider, R. Schmidt, P. Back, Y. Shimazaki, M. Kroner, and A. İmamoğlu, Phys. Rev. X 10, 021011 (2020).
  * Sternbach  _et al._ (2020) A. J. Sternbach, S. H. Chae, S. Latini, A. A. Rikhter, Y. Shao, B. Li, D. Rhodes, B. Kim, P. J. Schuck, X. Xu, X.-Y. Zhu, R. D. Averitt, J. Hone, M. M. Fogler, A. Rubio, and D. N. Basov, Science 371, 617 (2020).
  * Rana  _et al._ (2021) F. Rana, O. Koksal, M. Jung, G. Shvets, A. N. Vamivakas, and C. Manolatou, Phys. Rev. Lett. 126, 127402 (2021).
  * Woessner  _et al._ (2017) A. Woessner, Y. Gao, I. Torre, M. B. Lundeberg, C. Tan, K. Watanabe, T. Taniguchi, R. Hillenbrand, J. Hone, M. Polini, and F. H. L. Koppens, Nat. Photonics 11, 421 (2017).
  * Fan  _et al._ (2019) Z. Fan, S. Dutta-Gupta, R. Gladstone, S. Trendafilov, M. Bosch, M. Jung, G. R. S. Iyer, A. J. Giles, M. Shcherbakov, B. Feigelson, J. D. Caldwell, M. Allen, J. Allen, and G. Shvets, Nanophotonics 8, 1417 (2019).
  * Xiong  _et al._ (2021) L. Xiong, Y. Li, M. Jung, C. Forsythe, S. Zhang, A. S. McLeod, Y. Dong, S. Liu, F. L. Ruta, C. Li, K. Watanabe, T. Taniguchi, M. M. Fogler, J. H. Edgar, G. Shvets, and C. R. D. adn DN Basov, Sci. Adv. 7, eabe8087 (2021).
  * Fang  _et al._ (2007) T. Fang, A. Konar, H. Xing, and D. Jena, Appl. Phys. Lett. 91, 092109 (2007).
  * Jung  _et al._ (2018) M. Jung, Z. Fan, and G. Shvets, Phys. Rev. Lett. 121, 086807 (2018).
  * Yu  _et al._ (2013) G. L. Yu, R. Jalil, B. Belle, A. S. Mayorov, P. Blake, F. Schedin, S. V. Morozov, L. A. Ponomarenko, F. Chiappini, S. Wiedmann, U. Zeitler, M. I. Katsnelson, A. K. Geim, K. S. Novoselov, , and D. C. Elias, Proc. Natl. Acad. Sci. U.S.A. 110, 3282 (2013).
  * Katsnelson  _et al._ (2006) M. I. Katsnelson, K. S. Novoselov, and A. K. Geim, Nat. Phys. 2, 620 (2006).
  * Young and Kim (2009) A. F. Young and P. Kim, Nat. Phys. 5, 222 (2009).
  * Wang  _et al._ (2013) L. Wang, I. Meric, P. Huang, Q. Gao, Y. Gao, H. Tran, T. Taniguchi, K. Watanabe, L. Campos, D. Muller, J. Guo, P. Kim, J. Hone, K. Shepard, and C. R. Dean, Science 342, 614 (2013).
  * Ramezanali  _et al._ (2009) M. R. Ramezanali, M. M. Vazifeh, R. Asgari, M. Polini, and A. H. MacDonald, J. Phys. A: Math. Theor. 42, 214015 (2009).
  * Luo  _et al._ (2013) X. Luo, T. Qiu, W. Lu, and Z. Ni, Mater. Sci. Eng. R Rep. 74, 351 (2013).
  * Ni  _et al._ (2018) G. X. Ni, A. S. 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 577, 530 (2018).
  * Lundeberg  _et al._ (2017) M. B. Lundeberg, Y. Gao, R. Asgari, C. Tan, B. V. D. adn Marta Autore adn Pablo Alonso-González, A. Woessner, K. Watanabe, T. Taniguchi, R. Hillenbrand, J. Hone, M. Polini, and F. Koppens, Science 357, 187 (2017).
  * Torre  _et al._ (2017) I. Torre, M. I. Katsnelson, A. Diaspro, V. Pellegrini, and M. Polini, Phys. Rev. B 96, 035433 (2017).
  * Xiong  _et al._ (2019) L. Xiong, C. Forsythe, M. Jung, A. S. McLeod, S. S. Sunku, Y. M. Shao, G. X. Ni, A. J. Sternbach, S. Liu, J. H. Edgar, E. J. Mele, M. M. Fogler, G. Shvets, C. R. Dean, and D. N. Basov, Nat. Commun. 10, 4780 (2019).
  * Kurman and Kaminer (2020) Y. Kurman and I. Kaminer, Nat. Phys. 16, 868 (2020).
  * Kyriienko and Shelykh (2012) O. Kyriienko and I. A. Shelykh, J. or Nanophotonics 6, 061804 (2012).
  * Zanotto  _et al._ (2010) S. Zanotto, G. Biasiol, R. Degl’Innocenti, L. Sorba1, and A. Tredicucci, Appl. Phys. Lett. 97, 231123 (2010).
  * Anappara  _et al._ (2009) A. A. Anappara, S. D. Liberato, A. Tredicucci, C. Ciuti, G. Biasiol, L. Sorba, and F. Beltram1, Phys. Rev. B 79, 201303(R) (2009).
  * Lee  _et al._ (2014) J. Lee, M. Tymchenko, C. Argyropoulos, P.-Y. Chen, F. Lu, F. Demmerle, G. Boehm, M.-C. Amann, A. Alù, and M. A. Belkin, Nature 511, 65 (2014).
  * kun Shi  _et al._ (2019) L. kun Shi, J. Ma, and J. C. W. Song, 2D Mater. 7, 015028 (2019).

Supplementary Materials: Emergent intersubband-plasmon-polaritons of Dirac
electrons under one-dimensional superlattice

##  I Bound states of Dirac electrons in square potential wells

Here, we provide an analytical solution for bound states based on total
internal reflection of electrons with Dirac dispersion (TIREDD) in a square
potential well or a square potential barrier. Suppose a square potential well
($U_{0}<U_{1}$) or a square potential barrier ($U_{0}>U_{1}$) given as

$U_{E}(x)=U_{0}\text{ }(0<x<W),\text{ }U_{1}\text{ }(\text{elsewhere}).$ (S1)

Then, we can set an ansatz for a bound state with eigenenergy $E$ as

$\psi(\mathbf{r})=e^{ik_{y}y}\times\begin{cases}L\begin{bmatrix}1\\\
\frac{-i\alpha+ik_{y}}{(E-U_{1})/\hbar v_{F}}\end{bmatrix}e^{\alpha
x}&(x<0)\\\ A\begin{bmatrix}1\\\ \frac{q+ik_{y}}{(E-U_{0})/\hbar
v_{F}}\end{bmatrix}e^{iqx}+B\begin{bmatrix}1\\\
\frac{-q+ik_{y}}{(E-U_{0})/\hbar v_{F}}\end{bmatrix}e^{-iqx}&(0<x<W)\\\
R\begin{bmatrix}1\\\ \frac{i\alpha+ik_{y}}{(E-U_{1})/\hbar
v_{F}}\end{bmatrix}e^{-\alpha x}&(x>W)\end{cases},$ (S2)

where $q=\sqrt{\left(\frac{E-U_{0}}{\hbar v_{F}}\right)^{2}-k_{y}^{2}}$ is the
momentum along $x$-axis within the well/barrier, and
$\alpha=\sqrt{k_{y}^{2}-\left(\frac{E-U_{1}}{\hbar v_{F}}\right)^{2}}$ is the
decaying factor outside the well/barrier. Both $q$ and $\alpha$ are real;
therefore, this ansatz is possible only when $|E-U_{1}|<\hbar
v_{F}k_{y}<|E-U_{0}|$.

Figure S1: a. Plotting the inverse tangent of the right hand side of Eq. (S3);
$U_{0}=-0.2$eV, $U_{1}=-0.02$eV (or $U_{0}=-0.02$eV, $U_{1}=-0.2$eV; both
yield the same result), and $W=200$nm ($v_{F}=1.1\times 10^{6}$m/s). b.
Dispersion of bound states in $k_{y}$ for a potential well (left:
$U_{0}=-0.2$eV, $U_{1}=-0.02$eV) and for a potential barrier (left:
$U_{0}=-0.02$eV, $U_{1}=-0.2$eV). The dased lines refer to $\hbar
v_{F}k_{y}=|E-U_{0}|$ and $\hbar v_{F}k_{y}=|E-U_{1}|$.

By imposing the continuity of $\psi$ and the continuity of probability current
$\mathbf{J}=\psi^{\dagger}\sigma_{x}\psi\hat{\mathbf{x}}+\psi^{\dagger}\sigma_{y}\psi\hat{\mathbf{y}}$,
we obtain the following condition for the eigenenergy $E$:

$\tan(qW)=\frac{q\alpha}{\frac{E-U_{0}}{\hbar v_{F}}\frac{E-U_{1}}{\hbar
v_{F}}-k_{y}^{2}}.$ (S3)

Figure S1a shows the inverse tangent of the right hand side of the above
equation. Since the right hand side is vanisihng in most region, the bound
state energy condition simply reduces to $\tan(qW)\sim 0$. Therefore, we get
$E_{j}\sim\hbar v_{F}\sqrt{(\pi j/W)^{2}+k_{y}^{2}}+U_{0}$, which we discussed
in the main text. Figure S1b shows that a potential barrier also can host
bound states via TIREDD of the valence band electrons, as discussed in the
main text.

##  II Numerical method for the calculation of HIPP dispersion and normal
reflection spectrum

In this section, we elaborate on the numerical method used for the calculation
of HIPP dispersion and normal reflection spectrum shown in the main text. As
discussed in the main text, we only consider the transverse magnetic modes
that can be described with $E_{x}$, $E_{z}$ and $B_{y}$. Also, we consider the
dispersion in the momentum along $x$-axis ($\mathbf{q}=q\hat{\mathbf{x}}$).

Recall that the metagate periodicity is $L$, and the width of air gaps in the
metagate is $S$. For convenience, let’s define several variables to describe
the plane wave solutions in each of the layers—the air above all layers ($A$),
hBN layers ($B$), the oxide layer ($O$), and the substrate ($S$):
$q_{m}=q+\frac{2\pi m}{L}$, $q_{\omega}=\frac{\omega}{c}$,
$\kappa_{m}^{A}=-i\sqrt{q_{\omega}^{2}-q_{m}^{2}}$,
$\kappa_{m}^{B}=-i\sqrt{\epsilon_{xy}^{B}q_{\omega}^{2}-\frac{\epsilon_{xy}^{B}}{\epsilon_{z}^{B}}q_{m}^{2}}$,
$\kappa_{m}^{O}=-i\sqrt{\epsilon^{O}q_{\omega}^{2}-q_{m}^{2}}$, and
$\kappa_{m}^{S}=-i\sqrt{\epsilon^{S}q_{\omega}^{2}-q_{m}^{2}}$, where
$m\in\mathbb{Z}$ is an integer index. To describe the modes ($E_{x}$, $E_{z}$
and $B_{y}$) in the air gaps of the metagate ($M$):
$\eta_{\mu}=\frac{\pi\mu}{S}$, $\phi_{\mu}(x)=\sqrt{2-\delta_{\mu
0}}\cos\left[\eta_{\mu}(x-\frac{L-S}{2})\right]$, and
$\kappa_{\mu}^{M}=-i\sqrt{q_{\omega}^{2}-\eta_{\mu}^{2}}$, where $\mu$ is a
non-negative integer index.

Then, we can set an ansatz for the HIPP mode as below.

(i) In the air above the top hBN layer ($z>h_{t}$):

$\begin{split}cB_{y}&=\sum_{m}e^{iq_{m}x}\left(R_{m}e^{-\kappa_{m}^{A}(z-h_{t})}+I_{m}e^{\kappa_{m}^{A}(z-h_{t})}\right)\\\
k_{\omega}E_{x}&=-i\sum_{m}e^{iq_{m}x}\kappa_{m}^{A}\left(-R_{m}e^{-\kappa_{m}^{A}(z-h_{t})}+I_{m}e^{\kappa_{m}^{A}(z-h_{t})}\right)\\\
k_{\omega}E_{z}&=-\sum_{m}e^{iq_{m}x}q_{m}\left(R_{m}e^{-\kappa_{m}^{A}(z-h_{t})}+I_{m}e^{\kappa_{m}^{A}(z-h_{t})}\right).\end{split}$
(S4)

Here, $R_{0}$ is the reflection coefficient, in the presence of the normal
($q=0$) incident lght $I_{m}=\delta_{m0}$. An eigenmode exists even with
vanighing external drive term $I_{m}=0$.

(ii) In the top hBN layer ($0<z<h_{t}$):

$\begin{split}cB_{y}&=\sum_{m}e^{iq_{m}x}\left(A_{m}\cosh(\kappa_{m}^{B}z)+B_{m}\sinh(\kappa_{m}^{B}z)\right)\\\
k_{\omega}E_{x}&=\frac{-i}{\epsilon_{xy}^{B}}\sum_{m}e^{iq_{m}x}\kappa_{m}^{B}\left(A_{m}\sinh(\kappa_{m}^{B}z)+B_{m}\cosh(\kappa_{m}^{B}z)\right)\\\
k_{\omega}E_{z}&=\frac{-1}{\epsilon_{z}^{B}}\sum_{m}e^{iq_{m}x}q_{m}\left(A_{m}\cosh(\kappa_{m}^{B}z)+B_{m}\sinh(\kappa_{m}^{B}z)\right)\end{split}$
(S5)

(iii) In the bottom hBN layer ($-h_{b}<z<0$):

$\begin{split}cB_{y}&=\sum_{m}e^{iq_{m}x}\left(C_{m}\cosh(\kappa_{m}^{B}z)+B_{m}\sinh(\kappa_{m}^{B}z)\right)\\\
k_{\omega}E_{x}&=\frac{-i}{\epsilon_{xy}^{B}}\sum_{m}e^{iq_{m}x}\kappa_{m}^{B}\left(C_{m}\sinh(\kappa_{m}^{B}z)+B_{m}\cosh(\kappa_{m}^{B}z)\right)\\\
k_{\omega}E_{z}&=\frac{-1}{\epsilon_{z}^{B}}\sum_{m}e^{iq_{m}x}q_{m}\left(C_{m}\cosh(\kappa_{m}^{B}z)+B_{m}\sinh(\kappa_{m}^{B}z)\right)\end{split}$
(S6)

(iv) In the air gap of the metagate ($-h_{b}-h_{m}<z<-h_{b}$):

$\begin{split}cB_{y}&=\sum_{\mu}\phi_{\mu}(x)\left(D_{\mu}\cosh(\kappa_{\mu}^{M}(z+h_{b}))+E_{\mu}\sinh(\kappa_{\mu}^{M}(z+h_{b}))\right)\\\
k_{\omega}E_{x}&=-i\sum_{\mu}\phi_{\mu}(x)\kappa_{\mu}^{M}\left(D_{\mu}\sinh(\kappa_{\mu}^{M}(z+h_{b}))+E_{\mu}\cosh(\kappa_{\mu}^{M}(z+h_{b}))\right)\\\
k_{\omega}E_{z}&=i\sum_{\mu}\partial_{x}\phi_{\mu}(x)\left(D_{\mu}\cosh(\kappa_{\mu}^{M}(z+h_{b}))+E_{\mu}\sinh(\kappa_{\mu}^{M}(z+h_{b}))\right)\end{split}$
(S7)

(v) In the oxide layer ($-h_{b}-h_{m}-h_{o}<z<-h_{b}-h_{m}$):

$\begin{split}cB_{y}&=\sum_{m}e^{iq_{m}x}\left(F_{m}\cosh(\kappa_{m}^{O}(z+h_{b}+h_{m}))+G_{m}\sinh(\kappa_{m}^{O}(z+h_{b}+h_{m}))\right)\\\
k_{\omega}E_{x}&=\frac{-i}{\epsilon^{O}}\sum_{m}e^{iq_{m}x}\kappa_{m}^{O}\left(F_{m}\sinh(\kappa_{m}^{O}(z+h_{b}+h_{m}))+G_{m}\cosh(\kappa_{m}^{O}(z+h_{b}+h_{m}))\right)\\\
k_{\omega}E_{z}&=\frac{-1}{\epsilon^{O}}\sum_{m}e^{iq_{m}x}q_{m}\left(F_{m}\cosh(\kappa_{m}^{O}(z+h_{b}+h_{m}))+G_{m}\sinh(\kappa_{m}^{O}(z+h_{b}+h_{m}))\right)\end{split}$
(S8)

(vi) In the substrate ($z<-h_{b}-h_{m}-h_{o}$):

$\begin{split}cB_{y}&=\sum_{m}e^{iq_{m}x}H_{m}e^{\kappa_{m}^{S}(z+h_{b}+h_{m}+h_{o})}\\\
k_{\omega}E_{x}&=\frac{-i}{\epsilon^{S}}\sum_{m}e^{iq_{m}x}\kappa_{m}^{S}H_{m}e^{\kappa_{m}^{S}(z+h_{b}+h_{m}+h_{o})}\\\
k_{\omega}E_{z}&=\frac{-1}{\epsilon^{S}}\sum_{m}e^{iq_{m}x}q_{m}H_{m}e^{\kappa_{m}^{S}(z+h_{b}+h_{m}+h_{o})}\end{split}$
(S9)

At $z=0$, we can define the dynamic electric potential field on graphene
$\delta U_{E}$ from $-\partial_{x}(\delta U_{E}/(-e))=E_{x}(z=0)$, and the
dynamic carrier density oscillation $\delta n$ from $-e\delta
n=D_{z}(z=0^{+})-D_{z}(z=0^{-})$

$\begin{split}k_{\omega}\delta
U_{E}&=\frac{-e}{\epsilon_{xy}^{B}}\sum_{m}\frac{\kappa_{m}^{B}}{q_{m}}B_{m}e^{iq_{m}x}\\\
k_{\omega}\delta
n&=\frac{1}{e}\sum_{m}q_{m}(A_{m}-C_{m})e^{iq_{m}x}\end{split}$ (S10)

Now, we match the boundary conditions (continuity of $E_{x}$ and continuity of
$B_{y}$ with no free current density). At $z=0$, due to the current density at
graphene, we get $B_{y}(z=0^{+})-B_{y}(z=0^{-})=-\mu_{0}\sigma*E_{x}(z=0)$.
Again, for convenience, let’s define several vector/matrix notations:

$\\{A_{m}\\},\\{B_{m}\\},...,\\{I_{m}\\},\\{R_{m}\\}\rightarrow A,B,...,I,R$,

$[\mathcal{C}_{t}]_{mm^{\prime}}=\delta_{mm^{\prime}}\cosh(\kappa_{m}^{B}h_{t})$,
$[\mathcal{S}_{t}]_{mm^{\prime}}=\delta_{mm^{\prime}}\sinh(\kappa_{m}^{B}h_{t})$,

$[\mathcal{C}_{b}]_{mm^{\prime}}=\delta_{mm^{\prime}}\cosh(\kappa_{m}^{B}h_{b})$,
$[\mathcal{S}_{b}]_{mm^{\prime}}=\delta_{mm^{\prime}}\sinh(\kappa_{m}^{B}h_{b})$,

$[\mathcal{C}_{m}]_{\mu\mu^{\prime}}=\delta_{\mu\mu^{\prime}}\cosh(\kappa_{\mu}^{M}h_{m})$,
$[\mathcal{S}_{m}]_{\mu\mu^{\prime}}=\delta_{\mu\mu^{\prime}}\sinh(\kappa_{\mu}^{M}h_{m})$,

$[\mathcal{C}_{o}]_{mm^{\prime}}=\delta_{mm^{\prime}}\cosh(\kappa_{m}^{O}h_{o})$,
$[\mathcal{S}_{o}]_{mm^{\prime}}=\delta_{mm^{\prime}}\sinh(\kappa_{m}^{O}h_{o})$,

$[\mathcal{K}^{A}]_{mm^{\prime}}=\delta_{mm^{\prime}}\kappa_{m}^{A}$,
$[\mathcal{K}^{B}]_{mm^{\prime}}=\delta_{mm^{\prime}}\kappa_{m}^{B}/\epsilon_{xy}^{B}$,
$[\mathcal{K}^{M}]_{\mu\mu^{\prime}}=\delta_{\mu\mu^{\prime}}\kappa_{\mu}^{M}$,
$[\mathcal{K}^{O}]_{mm^{\prime}}=\delta_{mm^{\prime}}\kappa_{m}^{O}/\epsilon^{O}$,
$[\mathcal{K}^{S}]_{mm^{\prime}}=\delta_{mm^{\prime}}\kappa_{m}^{S}/\epsilon^{S}$,
$[\mathcal{Q}]_{mm^{\prime}}=\delta_{mm^{\prime}}q_{m}$,

$[\Sigma]_{mm^{\prime}}=-i\sigma(q_{n},q_{m},\omega)$, and $[\mathcal{T}]_{\mu
m}=\frac{1}{S}\int_{S}dx\psi_{\mu}^{*}(x)e^{iq_{m}x}$. Here, $\mathcal{T}$ is
the basis transformation matrix (from plane waves to the eigenmodes in the air
gap). From the boundary conditions, we get:

(i) Between the air and the top hBN:

$\begin{split}\mathcal{C}_{t}A+\mathcal{S}_{t}B&=R+I\\\
\mathcal{K}^{B}(\mathcal{S}_{t}A+\mathcal{C}_{t}B)&=\mathcal{K}^{A}(-R+I)\end{split}$
(S11)

(ii) Between the top and the bottom hBN:

$C-A=\frac{1}{\omega}\Sigma\mathcal{K}^{B}B$ (S12)

(iii) Between the bottom BN and the air gap of the metagate:

$\begin{split}D&=\mathcal{T}(\mathcal{C}_{b}C-\mathcal{S}_{b}B)\\\
\mathcal{K}^{M}E&=\mathcal{T}\mathcal{K}^{B}(\mathcal{C}_{b}B-\mathcal{S}_{b}C)\\\
\frac{S}{L}\mathcal{T}^{\dagger}\mathcal{K}^{M}E&=\mathcal{K}^{B}(\mathcal{C}_{b}B-\mathcal{S}_{b}C)\end{split}$
(S13)

(iv) Between the air gap of the metagate and the oxide:

$\begin{split}\mathcal{C}_{m}D-\mathcal{S}_{m}E&=\mathcal{T}F\\\
\mathcal{K}^{M}(\mathcal{C}_{m}E-\mathcal{S}_{m}D)&=\mathcal{T}\mathcal{K}^{O}G\\\
\frac{S}{L}\mathcal{T}^{\dagger}\mathcal{K}^{M}(\mathcal{C}_{m}E-\mathcal{S}_{m}D)&=\mathcal{K}^{O}G\end{split}$
(S14)

(v) Between the oxide and the substrate:

$\begin{split}\mathcal{C}_{o}F-\mathcal{S}_{o}G&=H\\\
\mathcal{K}^{O}(\mathcal{C}_{o}G-\mathcal{S}_{o}F)&=\mathcal{K}^{S}H\end{split}$
(S15)

From Eq. (S15), we can eliminate $H$ to obtain

$G=(\mathcal{K}^{O}\mathcal{C}_{o}+\mathcal{K}^{S}\mathcal{S}_{o})^{-1}(\mathcal{K}^{O}\mathcal{S}_{o}+\mathcal{K}^{S}\mathcal{C}_{o})F.$
(S16)

Then, from Eq. (S14) and Eq. (S16), we can eliminate $F$ and $G$ to obtain:

$E=(\mathcal{V}\mathcal{K}^{M}\mathcal{C}_{m}+\mathcal{S}_{m})^{-1}(\mathcal{V}\mathcal{K}^{M}\mathcal{S}_{m}+\mathcal{C}_{m})D,$
(S17)

where
$\mathcal{V}=\frac{S}{L}\mathcal{T}(\mathcal{K}^{O}\mathcal{S}_{o}+\mathcal{K}^{S}\mathcal{C}_{o})^{-1}(\mathcal{K}^{O}\mathcal{C}_{o}+\mathcal{K}^{S}\mathcal{S}_{o}){K^{O}}^{-1}\mathcal{T}^{\dagger}$.
Further, from Eq. (S13) and Eq. (S17), we can eliminate $D$ and $E$ to obtain:

$C=(\mathcal{W}\mathcal{C}_{b}+\mathcal{K}^{B}\mathcal{S}_{b})^{-1}(\mathcal{W}\mathcal{S}_{b}+\mathcal{K}^{B}\mathcal{C}_{b})B=\left[{\mathcal{K}^{B}}^{-1}\mathcal{S}_{b}{\mathcal{C}_{b}}^{-1}+(\mathcal{C}_{b}\mathcal{W}\mathcal{C}_{b}+\mathcal{K}^{B}\mathcal{C}_{b}\mathcal{S}_{b})^{-1}\right]\mathcal{K}^{B}B,$
(S18)

where
$\mathcal{W}=\frac{S}{L}\mathcal{T}^{\dagger}\mathcal{K}^{M}(\mathcal{V}\mathcal{K}^{M}\mathcal{C}_{m}+\mathcal{S}_{m})^{-1}(\mathcal{V}\mathcal{K}^{M}\mathcal{S}_{m}+\mathcal{C}_{m})\mathcal{T}$
or, equivalently,
$\mathcal{W}=\frac{S}{L}\mathcal{T}^{\dagger}\left[\mathcal{K}^{M}\mathcal{S}_{m}{\mathcal{C}_{m}}^{-1}+(\mathcal{C}_{b}\mathcal{V}\mathcal{C}_{b}+{\mathcal{K}^{M}}^{-1}\mathcal{C}_{m}\mathcal{S}_{m})^{-1}\right]\mathcal{T}$
(the latter form makes it clear that $\mathcal{W}$ is a Hermitian matrix).
From Eq. (S11), we can eliminate eliminate $R$ to obtain:

$A=2\mathcal{K}^{A}(\mathcal{K}^{A}\mathcal{C}_{t}+\mathcal{K}^{B}\mathcal{S}_{t})^{-1}I-\left[{\mathcal{K}^{B}}^{-1}\mathcal{S}_{t}{\mathcal{C}_{t}}^{-1}+(\mathcal{C}_{t}\mathcal{K}^{A}\mathcal{C}_{t}+\mathcal{K}^{B}\mathcal{C}_{t}\mathcal{S}_{t})^{-1}\right]\mathcal{K}^{B}B,$
(S19)

or we can eliminate $A$ to obtain:

$R=(\mathcal{K}^{A}\mathcal{C}_{t}+\mathcal{K}^{B}\mathcal{S}_{t})^{-1}\left[(\mathcal{K}^{A}\mathcal{C}_{t}-\mathcal{K}^{B}\mathcal{S}_{t})I-\mathcal{K}^{B}B\right].$
(S20)

If there is no external drive input $I=0$, Eq. (S18), and Eq. (S19) altogether
gives:

$\begin{split}C-A&=\mathcal{Y}\mathcal{K}^{B}B\text{, where}\\\
\mathcal{Y}=\left[{\mathcal{K}^{B}}^{-1}\mathcal{S}_{t}{\mathcal{C}_{t}}^{-1}+(\mathcal{C}_{t}\mathcal{K}^{A}\mathcal{C}_{t}+\mathcal{K}^{B}\mathcal{C}_{t}\mathcal{S}_{t})^{-1}\right]&+\left[{\mathcal{K}^{B}}^{-1}\mathcal{S}_{b}{\mathcal{C}_{b}}^{-1}+(\mathcal{C}_{b}\mathcal{W}\mathcal{C}_{b}+\mathcal{K}^{B}\mathcal{C}_{b}\mathcal{S}_{b})^{-1}\right].\end{split}$
(S21)

Then, by examining Eq. (S10) and recalling the definition of the dynamic
capacitance $\delta
U_{E}(\mathbf{q})=e^{2}\sum_{\mathbf{q}^{\prime}}C^{-1}(\mathbf{q},\mathbf{q}^{\prime})\delta
n(\mathbf{q}^{\prime})$, we arrive at

$C^{-1}(q_{m},q_{m^{\prime}})=\frac{[\mathcal{Y}^{-1}]_{mm^{\prime}}}{q_{m}q_{m^{\prime}}}=[\mathcal{Q}^{-1}\mathcal{Y}^{-1}\mathcal{Q}^{-1}]_{mm^{\prime}}.$
(S22)

The inverse of dynamic capacitance is equal to the dynamic coulomb interaction
Jung  _et al._ (2018); Torre  _et al._ (2017); Xiong  _et al._ (2019). It is
easily checked that, in a suspended graphene in vacuum
($h_{t},h_{b}\rightarrow\infty$,
$\epsilon_{xy}^{B}=\epsilon_{z}^{B}=\epsilon_{0}$), the above expression
reduces to the expression obtained by taking the 2D Fourier transform of
Coulomb potential
$e^{2}C^{-1}(q_{m},q_{m^{\prime}})\rightarrow\delta_{mm^{\prime}}\frac{e^{2}}{2\epsilon_{0}|q_{m}|}$.
Thus, the dynamical dielectric function given in the main text is reduced as a
compact matrix form:

$\epsilon(\mathbf{q},\omega)=\mathbf{1}-\frac{1}{\omega}\mathcal{Q}^{-1}\mathcal{Y}^{-1}\Sigma\mathcal{Q}=\mathcal{Q}^{-1}\left[\mathbf{1}-\frac{1}{\omega}\mathcal{Y}^{-1}\Sigma\right]\mathcal{Q}.$
(S23)

Therefore, the density of state, which is approximated as the imaginary part
of the inverse of the dynamical dielectric function, is given as

$DOS(\mathbf{q},\omega)=-\text{Im}\left[\text{Tr}\left([\epsilon(\mathbf{q},\omega)]^{-1}\right)\right]=-\text{Im}\left[\text{Tr}\left(\left[\mathbf{1}-\frac{1}{\omega}\mathcal{Y}^{-1}\Sigma\right]^{-1}\right)\right].$
(S24)

This is the quantity plotted in the main text for the figure containing the
HIPP dispersions.

Now, we can combine Eq. (S12), Eq. (S18), and Eq. (S19) to obtain:

$R=(\mathcal{K}^{A}\mathcal{C}_{t}+\mathcal{K}^{B}\mathcal{S}_{t})^{-1}\left[(\mathcal{K}^{A}\mathcal{C}_{t}-\mathcal{K}^{B}\mathcal{S}_{t})-2\left(\mathcal{Y}-\frac{1}{\omega}\Sigma\right)^{-1}\mathcal{K}^{A}(\mathcal{K}^{A}\mathcal{C}_{t}+\mathcal{K}^{B}\mathcal{S}_{t})^{-1}\right]I.$
(S25)

As mentioned earlier, the reflection upon normal incidence is calculated as
$|R_{0}|^{2}$ with $I_{m}=\delta_{m0}$ and $q=0$.

Figure S2: a SL potential $U_{E}(x)$ (bottom) for different backgate voltages
($V_{M}=1$V is for all three cases); dotted red: $V_{B}=7$V, dashed blue:
$0$V, and solid green: $-2$V. ($h_{t}=5$nm, $h_{b}=10$nm, $h_{g}=10$nm,
$h_{0}=150$nm, $L=300$nm, $S=150$nm). b Reflection spectra for the normal
incidence of light polarized along $x$-axis; solid thin black: no modulation
($E_{F}=0.15$eV), dotted red: $V_{B}=7$V, dashed blue: $0$V, and solid thick
green: $-2$V. c Density of states or
$-\text{Im}\left[\text{Tr}\left([\epsilon(\mathbf{q},\omega)]^{-1}\right)\right]$
for visualizing the HIPP dispersion.

##  III HIPP dispersions under a different superlattice design

In this section, we show that the appearance of the hybrid intersubband-
plasmon-polaritons (HIPPs) shown in the main text is not contingent upon a
specific set of parameter conditions, by providing the HIPP dispersions for a
different superlattice design. In the main text, the results were shown for a
system with a periodicity of $L=300$nm and the air gap width of $S=80$nm.
Here, in Fig. S2, we provided the same calculations for $S=150$nm. For a
similar degree of $U_{E}$ modulation, the intersubband transition (ISBT)
frequencies are slightly blue-shifted, compared to the results in the main
text, since the potential well width $L-S$ is now deceased. Other than such
small details, the HIPP phenomenon is qualitatively the same. Therefore, the
experimental verification of this HIPP emergence under 1D SL in graphene would
be universally possible for nearly any choice of parameter conditions.