# Tunable exciton polaritons in band-gap engineered hexagonal boron nitride

Pedro Ninhos  0000-0002-1143-7457 POLIMA—Center for Polariton-driven
Light–Matter Interactions, University of Southern Denmark, Campusvej 55,
DK-5230 Odense M, Denmark  Christos Tserkezis  0000-0002-2075-9036
POLIMA—Center for Polariton-driven Light–Matter Interactions, University of
Southern Denmark, Campusvej 55, DK-5230 Odense M, Denmark  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
POLIMA—Center for Polariton-driven Light–Matter Interactions, University of
Southern Denmark, Campusvej 55, DK-5230 Odense M, Denmark  Centro de Física
(CF-UM-UP) and Departamento de Física, Universidade do Minho, P-4710-057
Braga, Portugal  International Iberian Nanotechnology Laboratory (INL), Av
Mestre José Veiga, 4715-330 Braga, Portugal

###### Abstract

We show that hexagonal boron nitride (hBN), a two-dimensional insulator, when
subjected to an external superlattice potential forms a new paradigm for
electrostatically tunable excitons in the near- and mid-ultraviolet (UV). The
imposed potential has three consequences: (i) it renormalizes the effective
mass tensor, leading to anisotropic effective masses; (ii) it renormalizes the
band gap, eventually reducing it; (iii) it reduces the exciton binding
energies. All these consequences depend on a single dimensionless parameter,
which includes the product of strength of the external potential with its
period. In addition to the excitonic energy levels, we compute the optical
conductivity along two orthogonal directions, and from it the absorption
spectrum. The results for the latter show that our system is able to mimic a
grid polarizer. These characteristics make one-dimensional hBN superlattices a
viable and unexplored platform for fine-tuned polaritonics in the UV to
visible spectral range.

##  I Introduction

The properties of two-dimensional (2D) materials are characterized by a
richness of physical phenomena, encompassing unusual elastic [1], electronic
[2], and optical properties [3]. An exemplary case of enhanced light–matter
interactions in 2D materials is provided by graphene [4], whose transmission
in the visible range of the spectrum is characterized by the value of the fine
structure constant alone, a unique case in materials. Concerning the
electronic properties in 2D materials, we can find semi-metals [5], conductors
[6], semi-conductors and insulators [7, 8], strange metals [9], metallic
ferromagnets [10, 11], and superconductors [12, 13, 14, 15]. Therefore, 2D
materials have become a condensed-matter physicist’s playground due to this
wealth of physical properties. Some systems, such as hexagonal boron nitride
(hBN), have an interesting story: its first application was as an
encapsulating medium for graphene and other 2D materials, due to its ultra-
flat surface and low density of defects [16]. Such use allows to increase the
mobility of charge carriers [17] and produce well defined (large quality
factor $Q$) electromagnetic (EM) resonances in polaritonic devices. Later, it
became clear that hBN is interesting in its own right when light–matter
interaction is considered [18, 19]. The interaction with infrared radiation
due to polar phonons in hBN gives rise to large-$Q$ polaritonic resonances
used for example in gas sensing [20]. At the same time, the interaction of
ultraviolet (UV) radiation with the electrons in hBN has opened the field of
UV photonics [21].

Figure 1: (a) A hBN monolayer over a substrate, subjected to an external
periodic potential with period $L$ along the $x$-direction, under a time-
periodic drive of a linearly polarized optical field. (b) Brillouin zone (BZ)
of the original monolayer hBN, and (c) the BZ of the superlattice resulting
from applying the 1D potential of (a).

2D materials can be combined by stacking several layers, forming van der Waals
heterostructures; materials consisting of few layers have proven to be a
fertile theoretical and experimental playground for the study of
unconventional physics [22]. On one hand, heterostructures are promising for
tailoring optoelectronic properties of 2D insulators (such as hBN) [18]. On
the other hand, homostructures reveal quite interesting features, distinct
from those of the initial isolated constituents. From the example of AB-
stacking graphene bilayer, it was realized that a band gap can be opened and
controlled electrostatically, using a voltage bias [23, 24]. In twisted
bilayer graphene, it was observed that a superlattice potential has a
tremendous impact on the electronic and optical properties of 2D materials
[25, 26, 27]. Armed with this information, it is possible to envision a
combination of both electrostatic and superlattice potentials such that the
electronic and optical properties of 2D materials can be externally
controlled. One can manipulate the electronic and optical properties of a
material by band-gap engineering [28]. For example, it is possible to change
the gap of bulk hBN by femtosecond laser pulses [29]. For monolayer hBN, the
band gap can be reduced by cladding the hBN between two graphene layers [30],
or by functionalization [31].

In this context, it is pertinent to ask if one can manipulate the physical
properties of hBN, the 2D material considered is this work, by forming a
superlattice, resulting from combining a hBN lattice with an electrostatic
periodic potential [Fig. 1(a)]. As a pristine 2D lattice, hBN is a wide band
gap insulator in which boron (B) and nitrogen (N) atoms are arranged in a
honeycomb lattice, bounded by strong covalent bonds. Unlike graphene, hBN
presents a large band gap at the K and K’-points of the Brillouin zone (BZ)
[see Fig. 1(b)], and the low-energy dispersion is quadratic and not linear, as
it is in graphene. In a recent work [32], a direct measurement of the density
of states of monolayer hBN by scanning tunneling microscopy revealed a band
gap of $6.8\pm 0.2$ eV, while optical spectroscopy revealed an exciton binding
energy of $0.7\pm 0.2$ eV. With hBN being a large-band-gap insulator,
excitonic effects play an important role in its optical properties in the deep
UV range [33]. This was shown when an experimental measurement of
photoluminescence of monolayer hBN on quartz [19] reported an optical
resonance around 6.1 eV, that can only be explained by taking excitonic
effects into account.

When studying the sub-gap optical properties of hBN, we cannot ignore the fact
that the exciton can excite an EM mode confined at the surface of the 2D
material, called exciton polariton. The intriguing properties of polaritons
allow for the detection of small changes in the refractive index of the
surrounding material, rendering the 2D material viable for the fabrication of
optoelectronic sensors. Given the energy scales of the band gap and exciton
binding energies of hBN, one can think of it as a good candidate for the
detection of bio-molecules that present resonances in the UV, such as cyclic
$\beta$-helical peptides, relevant for the biophysics of the human body [34].

In this work, we propose means of manipulating the optoelectronic properties
of hBN by applying a one-dimensional (1D) periodic potential on a monolayer,
as represented in Fig. 1(a). As we will see, the resultant superlattice has
distinct electronic properties from the pristine hBN monolayer. The layer is
deposited on a quartz substrate, and is illuminated from above with
monochromatic light, linearly polarized along either the $x$ or the $y$-axis.
This driving field will allow the formation of excitons that, by coupling with
light, will give rise to exciton polaritons propagating along the sample.

This paper is organized as follows: In Sec. II we derive the renormalized low-
energy Hamiltonian, from which we can extract the emergent electronic
structure of hBN with a 1D periodic potential on top. In Sec. III we introduce
the formalism adopted for the calculation of the excitonic optical response of
the system under study. In Sec. IV we provide results for the new excitonic
physics, and the consequences it brings to the sub-gap optical response,
namely the optical conductivity and optical absorption. We finish the section
with a study of the spatial confinement of transverse magnetic (TM) exciton
polaritons. In Sec. V we summarize our main results and offer our conclusions.

##  II Renormalized Hamiltonian

The electronic properties of independent electrons in 2D materials are
described in terms of tight-binding models. In general, these models use a
fairly large number of bands, which prevents analytical calculations of the
materials’ properties. Fortunately, the physics of excitons is often
determined by effective models with only two bands, whose Hamiltonian can be
written using Pauli matrices. These are the cases of graphene and hBN. The
difference between the two systems lies in the presence of a band gap in the
latter. With all this in mind, an effective, low-energy, independent-electron
Hamiltonian can be written as [35]

$H_{0}=-i\hbar
v_{F}\left(\sigma_{x}\partial_{x}+\sigma_{y}\partial_{y}\right)+M\sigma_{z}\,,$
(1)

where $v_{F}$ is the Fermi velocity, $E_{g}=2M$ is the band gap obtained
within density-functional theory (DFT) [33], and the $\sigma$’s are the Pauli
matrices. The low-energy dispersion relation reads

$E_{s}(\mathbf{q})=s\frac{1}{2}\sqrt{E_{g}^{2}+4\hbar^{2}v_{F}^{2}(q_{x}^{2}+q_{y}^{2})}\approx
s\left(\frac{E_{g}}{2}+\frac{\hbar^{2}|\mathbf{q}|^{2}}{2m^{*}}\right)\,,$ (2)

where the wavevector $\mathbf{q}$ is measured from the $\mathbf{K}$ point in
the BZ (see Fig. 1), while $m^{*}=E_{g}/2v_{F}^{2}$ is the common effective
mass of electrons in the conduction band ($s=+1$) and holes in the valence
band ($s=-1$), respectively. The approximation of parabolic bands is valid
when the band gap is large and when looking at electronic states with small
momentum (measured from the point around which we are expanding the
Hamiltonian).

Since we want to describe the effect of an external potential on the
electronic properties of electrons in hBN (including its excitonic
properties), we introduce an external electrostatic potential $V(x)$, so that
the total Hamiltonian reads

$H=H_{0}+V(x).$ (3)

The potential varies only along the $x$-direction and has period $L$, that is,
$V(L+x)=V(x)$. This external potential then creates a superlattice, leading to
a mini-BZ of size $[-\pi/L,\pi/L[$, as in Fig. 1(c), and the appearance of a
set of mini-bands. This problem is similar to graphene under an external
potential [36], but with an extra term that confers a mass to the electrons
thus opening an energy gap at $K$ of the original BZ. By introducing the
external potential $V(x)$, we add an extra periodicity in the system in
addition to the original lattice periodicity (already encoded in the original
low-energy Hamiltonian). As noted above, we form a 1D superlattice whose
reciprocal lattice vectors are given by $\mathbf{G}_{\ell}=\ell G_{0}\hat{x}$,
with $G_{0}=2\pi/L$ and $\ell$ an integer.

It would be convenient if we could remove the potential from the Hamiltonian
(3), while maintaining the Dirac-like form of $H_{0}$. Fortunately, this is
possible via a unitary transformation written in terms of $V(x)$ and given by
[36]

$U_{\alpha}=\frac{1}{\sqrt{2}}\left(\begin{array}[]{cc}e^{-i\alpha(x)/2}&-e^{i\alpha(x)/2}\\\
e^{-i\alpha(x)/2}&e^{i\alpha(x)/2}\end{array}\right),$ (4)

where $\alpha(x)=(2/\hbar
v_{F})\int_{0}^{x}\mathrm{d}x^{\prime}V(x^{\prime})$. Focusing on electronic
states whose wavevector $\mathbf{k}=\mathbf{q}+\ell\mathbf{G}_{0}/2$, with
$|\mathbf{q}|\ll G_{0}$, we obtain an analytical expression for the
renormalized low-energy Hamiltonian, that is, the initial low-energy
Hamiltonian but with renormalized parameters. For that, we use the following
states as the basis functions

$\left(\begin{smallmatrix}1\\\ \\\
0\end{smallmatrix}\right)e^{i(\mathbf{q}+\ell\mathbf{G}/2)\cdot\mathbf{r}}\quad\text{
and }\quad\left(\begin{smallmatrix}0\\\ \\\
1\end{smallmatrix}\right)e^{i(\mathbf{q}-\ell\mathbf{G}/2)\cdot\mathbf{r}},$
(5)

and expand the exponential in a Fourier series as

$e^{i\alpha(x)}=\sum_{\ell=-\infty}^{\infty}J_{\ell}\left(\beta\right)e^{i\ell
G_{0}x},\quad\beta\equiv\frac{V_{0}L}{\pi\hbar v_{F}},$ (6)

where $J_{\ell}$ is the Bessel function of the first kind of order $\ell$.
After applying the unitary transformation, computing the matrix elements for
the previous basis states, and changing to a basis made of a symmetric and
anti-symmetric combination of the previous states, we obtain the renormalized
low-energy Hamiltonian

$\bar{H}=\frac{\ell}{2}G_{0}\hbar v_{F}+\hbar
v_{F}[\sigma_{x}q_{x}+\sigma_{y}J_{\ell}(\beta)q_{y}]-\sigma_{z}MJ_{\ell}(\beta).$
(7)

Here, the first term simply shifts the energies. Since $\ell$ is an integer,
we have now multiple bands, whose emergence can also be predicted by the
theory of nearly-free electrons. The comparison between the two treatments is
done in App. A. In Eq. (7), we notice that for a quasiparticle propagation
along the $y$-direction the group velocity is renormalized as
$v_{F}\rightarrow J_{\ell}(\beta)v_{F}$. Since $|J_{\ell}(\beta)|<1$, the
slope of the dispersion along the $y$-direction is smaller than along the
$x$-direction. This gives rise to an anisotropic Dirac cone. Also, the band
gap $E_{g}=2M$ changes as $E_{g}\rightarrow E_{g}|J_{\ell}(\beta)|$.
Therefore, it is possible to reduce the pristine gap of hBN by choosing an
appropriate value for $\beta$. With this mechanism we can control both the gap
and the excitonic states, which now become the excitons of a synthetic
anisotropic material, that resembles naturally occurring anisotropic
materials, such as phosphorene.

Next, we explore the consequences of the introduction of the 1D potential in
the electronic properties of hBN. From the renormalized low-energy Hamiltonian
(7), we infer the low-energy dispersion (up to a global shift in the energies)

$E_{s}\left(\mathbf{q}\right)=s\frac{1}{2}\sqrt{E_{g}^{2}J_{\ell}^{2}(\beta)+4\hbar^{2}v_{F}^{2}\left[q_{x}^{2}+q_{y}^{2}J_{\ell}^{2}(\beta)\right]},$
(8)

which is now anisotropic, while we recover the isotropic Eq. (2) for
$\beta\rightarrow 0$ and $\ell=0$. This deformed paraboloid has two distinct
effective masses $m_{x}=m^{*}|J_{\ell}(\beta)|$ and
$m_{y}=m^{*}/|J_{\ell}(\beta)|$ along the $x$ and $y$-directions,
respectively. The eigenproblem
$\bar{H}|s,\mathbf{q}\rangle=E_{s}(\mathbf{q})|s,\mathbf{q}\rangle$ has an
analytical solution for the spinors

$\displaystyle|+,\mathbf{q}\rangle$
$\displaystyle=\left(\begin{array}[]{c}e^{-i\theta_{q}\left(\beta\right)}\sin\left(\frac{\xi_{q}}{2}\right)\\\
\cos\left(\frac{\xi_{q}}{2}\right)\end{array}\right),\,$ (9c)
$\displaystyle|-,\mathbf{q}\rangle$
$\displaystyle=\left(\begin{array}[]{c}\cos\left(\frac{\xi_{q}}{2}\right)\\\
-e^{i\theta_{q}\left(\beta\right)}\sin\left(\frac{\xi_{q}}{2}\right)\end{array}\right),$
(9f)

where $\xi_{q}=\arctan\left(\frac{2\hbar
v_{F}q\left(\beta\right)}{J_{\ell}(\beta)E_{g}}\right)$, with
$q\left(\beta\right)=\sqrt{q_{x}^{2}+J_{\ell}^{2}(\beta)q_{y}^{2}}$, and
$\theta_{q}\left(\beta\right)=\arctan\left[J_{\ell}(\beta)\tan\theta\right]$,
with $\tan(\theta)=q_{y}/q_{x}$. These spinors will be useful later on to
compute the optical conductivity.

##  III Excitonic contributions to sub-gap optical response

The presence of the external potential $V(x)$ influences considerably the
electronic properties of hBN, by rendering the low-energy dispersion
anisotropic. This brings consequences to the excitonic physics, and
subsequently to the optical response. There are several approaches to the
calculation of the excitonic properties of a semiconductor or a band
insulator. Excitons with a radius large enough compared to the lattice spacing
are denoted Wannier excitons. The general approach to determine the exciton
energy levels and wavefunctions is the solution of the Bethe–Salpeter equation
[37] in momentum space. This microscopic approach is powerful but often
computationally demanding. It would then be desirable to have in reach a
simpler approach that still captures the essence of the exciton physics.
Starting from a microscopic approach [38], if some —relatively flexible—
conditions are met, the Fourier transform of the Bethe–Salpeter equation
transforms into an effective Schrödinger equation with an electrostatic
attractive potential called the Wannier equation. The method of the Wannier
equation has been benchmarked against the Bethe–Salpeter equation with good
results, as the latter equation accounts for fine details only [38].

The Wannier equation can be solved, among others, with the variational method,
which allows for a good approximation of the energies and approximate
analytical expressions for the excitonic wavefunctions. This approach has
previously been applied to black phosphorus [39] and TMDs [40, 41], with
significant accuracy.

###  III.1 Variational method

The variational method [42] consists in using a trial wavefunction, that
depends on a set of parameters, to compute the expectation value of the
Hamiltonian $H$ and minimize it with respect to said parameters. Symbolically,
we define a state $|\psi_{\alpha}(\mathbf{r})\rangle$, with $\alpha$
representing the set of unknown parameters. Then we compute
$E_{\alpha}=\langle\psi_{\alpha}(\mathbf{r})|H|\psi_{\alpha}(\mathbf{r})\rangle$
and solve $\partial E/\partial\alpha=0$ for $\alpha$. Finally, we return to
$E_{\alpha}$ and substitute the optimal values for $\alpha$.

The Hamiltonian $H$ pertains to the system comprised of an electron–hole pair,

$H=\frac{p_{x}^{2}}{2\mu_{x}}+\frac{p_{y}^{2}}{2\mu_{y}}+V_{\text{RK}}(r),$
(10)

where the first two terms correspond to the kinetic contribution, with
$\mu_{x}$ and $\mu_{y}$ being the electron–hole reduced masses along the $x$
and $y$-directions, respectively, and the last term corresponding to the
electrostatic attraction between the electron and the hole in the 2D material,
here given by the electrostatic Rytova–Keldysh (RK) potential [37]

$V_{\text{RK}}(r)=-\frac{e^{2}}{4\pi\varepsilon_{0}r_{0}}\frac{\pi}{2}\left[\mathbf{H}_{0}\left(\frac{\kappa
r}{r_{0}}\right)-Y_{0}\left(\frac{\kappa r}{r_{0}}\right)\right].$ (11)

Here, $r_{0}$ is a screening length parameter, $\kappa$ is the average between
the relative permittivity of the substrate $\varepsilon_{2}$ and the cladding
material surrounding the 2D material $\varepsilon_{1}$, and $\varepsilon_{0}$
is the vacuum permittivity. Finally, $\mathbf{H}_{0}$ is the Struve function
of order zero and $Y_{0}$ is the Bessel function of the second kind of order
zero. The exciton dynamics is determined by the Wannier equation
$E_{b}\psi=H\psi$, which is a Schrödinger-like equation, where $H$ is the
Hamiltonian given in Eq. (10), $\psi$ is the exciton wavefunction and
$E_{b}=E-\Delta$ is the exciton binding energy, while $E$ is the energy level
of the exciton and $\Delta$ the band gap. To obtain the energy levels of the
exciton we have to compute the binding energies and the band gap. For the
former, we apply the variational method to the first four excitonic states
$1s$, $2x$, $2y$ and $2s$, and use trial wavefunctions that resemble the
solutions of the 2D hydrogen atom, but taking into account the effect of the
anisotropy in the mass tensor. The form of each wavefunction is given in App.
C.1, and the details of the calculations for finding the variational
parameters are given in [39]. The band gap is the sum of the renormalized DFT
gap and the electron–electron exchange correction, both contributions
renormalizing with the potential (see App. B for details).

###  III.2 Optical conductivity

Subjecting the excitonic material to an external and monochromatic EM field
with frequency $\omega$, the optical conductivity, within the linear response
and in the dipole approximation, is given by [43]

$\sigma(\omega)=\frac{e^{2}}{8\pi^{3}i\hbar}\sum_{\nu}\left[\frac{E_{\nu}|\Omega_{\nu}|^{2}}{E_{\nu}-(\hbar\omega+i\Gamma)}+\frac{E_{\nu}|\Omega_{\nu}|^{2}}{E_{\nu}+\hbar\omega+i\Gamma}\right]$
(12)

where $E_{\nu}$ is the energy of the excitonic state $\nu$. The broadening
parameter $\Gamma\ll\omega$ is the sum of the non-radiative decay due to
phonons, the non-radiative decay due to impurities and the radiate decay, but
its calculation is outside the scope of this work. Therefore, we made an
estimation based on previous works [44]. Furthermore, $\Omega_{\nu}$ is the
Berry curvature of the excitonic state $\nu$, being computed by
$\Omega_{\nu}=\sum_{\mathbf{q}}\Psi_{\nu}(\mathbf{q})\langle-,\mathbf{q}|\mathbf{r}\cdot\hat{\mathbf{e}}|+,\mathbf{q}\rangle$,
where $\Psi_{\nu}(\mathbf{k})$ is the wavefunction of the excitonic state
$\nu$ in momentum space, and the braket is the optical matrix element of the
dipole operator between the spinors given in Eq. (9). The unit vector
$\hat{\mathbf{e}}$ is dictated by the polarization of the wave, and herein we
consider light polarized along either the $x$ and the $y$-directions. As the
system is anisotropic, the conductivity will be different depending on the
polarization of the incident EM field. As explained in App. C, we have two
distinct conductivities $\sigma_{xx}$ and $\sigma_{yy}$ for light polarized
along the $x$ and $y$-directions, respectively. Knowing the optical
conductivity, we compute the sub-gap optical absorption through

$\mathcal{A}=1-\mathcal{R}-\real\left\\{\sqrt{\frac{\varepsilon_{2}}{\varepsilon_{1}}}\right\\}\mathcal{T},$
(13)

where the reflection and transmission coefficients $\mathcal{R}$ and
$\mathcal{T}$ follow from Fresnel’s equations for normal incidence on a
conductive material cladded between two dielectrics [45]. We note that the
relative permittivities that enter in Eq. (13) are the dynamic dielectric
functions of the respective materials. In the electron–electron or
electron–hole interaction, the relevant quantities are the static
permittivities of the materials. In the end, the dependence is only on the
average between the static permittivities, that we denoted by $\kappa$ in Eq.
(11).

##  IV Results and discussion

We have seen that the 1D external potential influences the electronic
properties of hBN. In particular, it renormalizes the gap and the effective
masses, depending on the strength and period of the potential. In fact, by
inspection of the renormalized low-energy Hamiltonian in Eq. (7), we can claim
that the renormalized parameters depend only on the product $V_{0}L$ between
the strength and period of the potential. Hence, for a complete analysis, it
only suffices to study the system varying $\beta$. The surrounding medium is
air ($\varepsilon_{1}=1$), while for the substrate we choose silicon dioxide
(SiO2). The static dielectric function of SiO2 is $\varepsilon_{2}(\omega\to
0)=3.9$, and in the UV range it is essentially constant, with
$\varepsilon_{2}(\omega)\approx 2.16$.

###  IV.1 Excitonic binding energies

As we have seen in Sec. II, the effective masses renormalize as
$m_{x}=m^{*}|J_{0}(\beta)|$ and $m_{y}=m^{*}/|J_{0}(\beta)|$. This means that
as $J_{0}(\beta)\to 0$, the dispersion becomes more and more anisotropic. To
study the effects of this anisotropy in the electronic properties of the
system at hand, we compute the band gap and the binding energies of the
excitonic states $1s$, $2x$, $2y$ and $2s$ for varying $\beta$, from the case
$\beta=0$ without potential, to $\beta=2.3$ sufficiently close to the zero of
the Bessel function $J_{0}(\beta)$. We stress that the band gap includes the
electron–electron exchange correction, that also renormalizes with $\beta$
(see App. B). The results are displayed in Fig. 2. It can be seen in Fig. 2(a)
that by manipulating the value of $\beta$, we can tune the band gap. When
$\beta$ reaches the zero of $J_{0}(\beta)$, the gap closes. Regarding the
binding energies, Fig. 2(b) shows that the anisotropy lifts the degeneracy
between the states $2x$ and $2y$. Also, the binding energies of the states
$1s$, $2x$ and $2s$ decrease as we increase $\beta$, but that of the $2y$
state is non-monotonic. Nonetheless, the state $2y$ is always the first
excited state. Initially, the second excited state is $2x$, but for
$\beta\approx 1.5$, the $2s$ takes its place.

Figure 2: (a) Band gap and (b) binding energies of the excitonic states $1s$
(black continuous line), $2x$ (red dashed line), $2y$ (red continuous line)
and $2s$ (black dashed line) as a function of $\beta$. In both panels, the top
horizontal axis displays the value of the effective mass anisotropy defined as
$\Delta\mu=(m_{y}-m_{x})/m^{*}$. For hBN, the parameters are [33]:
$E_{g}=3.92$ eV, $\hbar v_{F}=5.06$ eV$\cdot$Å, and $r_{0}=10$ Å. For the
surrounding medium we have air with $\varepsilon_{1}=1$, while for the quartz
substrate $\varepsilon_{2}=3.9$. Figure 3: Sub-gap optical conductivity of hBN
under a 1D periodic potential, as depicted in Fig. 1. Panels (a) and (b) show
the real and imaginary parts of the conductivity, respectively, for light
polarized along the $x$-axis. Panels (c) and (d) show the real and imaginary
parts of the conductivity, respectively, for light polarized along the
$y$-axis. The in-set in panels (a) and (b) show a zoom-in of the real and
imaginary parts of the conductivity, respectively, for the cases $\beta=0$ and
$\beta=0.53$. Results for different $\beta$ have been rescaled as noted in
each panel, for clarity. We used $\Gamma=10$ meV, and for the quartz substrate
$\varepsilon_{2}=2.16$.

###  IV.2 Optical conductivity and Absorption

A recent experimental work [19] showed that single-layer hBN displays optical
activity in the UV range. The measured photoluminescence exhibits two sharp
peaks at a frequency around $6.1$ eV. These peaks are excitonic in nature, and
their theoretical prediction in [44] agrees very well with the experimental
measurement. Here we apply the same theory as in [44] to compute the optical
absorption of monolayer hBN subjected to the external periodic potential. But
first, we computed the optical conductivity using Eq. (12) for selected values
of $\beta$, to obtain the results displayed in Figs. 3a-3d. Since the optical
conductivity is complex, we plot separately the real part in Figs. 3a and 3c,
and the imaginary part in Figs. 3b and 3d. In the calculation we focus on the
excitonic states $1s$ and $2s$, while the states $2x$ and $2y$ do not
contribute to the conductivity. We also considered both cases of light
polarized along the $x$ and $y$ axes. We observe a tendency of a redshift in
the peaks as we increase $\beta$, and consequently the anisotropy of the
system. We note that the energy level $E_{\nu}$ of the excitonic state $\nu$
is given by $E_{\nu}=E_{b,\nu}+\Delta$, where $E_{b,\nu}$ is the respective
binding energy and $\Delta$ the band gap, quantities computed in the previous
section. For light polarized along the $x$-axis, we observe in Figs. 3a and 3c
that the peaks of the conductivity coming from the $1s$ and $2s$ states are
enhanced as we increase $\beta$. On the other hand, the peaks in the
conductivity for light polarized along the $y$-direction are attenuated.

Figure 4: Sub-gap optical absorption of hBN under a 1D periodic potential.
Panels (a) and (b) show the absorption for light polarized along the $x$ and
$y$-directions, respectively. Each color corresponds to a different value of
$\beta$, where $\beta=0$ corresponds to the case with no potential. We used
the same values for $\beta$ of Fig. 3. For each realization, two peaks are
observed. The one occurring at lower energy corresponds to the $1s$ state,
while the other corresponds to the $2s$ state. Panel (c) shows the values of
the absorption peaks due to the excitonic states $1s$ (black lines) and $2s$
(red lines) as a function of $\beta$ for light polarized along the $x$
(continuous lines) and $y$ (dashed lines) directions, respectively. Note that
for each value of $\beta$ the peaks occur at different energies, as one can
see in panels (a) and (b).

One can notice also in Fig. 3b that for the case $\beta=2.3$, the excitonic
states $1s$ and $2s$ become so close, that the imaginary part of the
conductivity $\imaginary\\{\sigma(\omega)\\}$ has only one zero. We expect
that this influences the character of the polaritons, since the sign of
$\imaginary\\{\sigma(\omega)\\}$ dictates which kind of polaritons, transverse
magnetic (TM) or transverse electric (TE) the 2D material supports, as we will
explain in Sec. IV.3.

We can now compute the optical absorption at normal incidence using Eq. (13).
Figs. 4a and 4b display the optical absorption for the same values of $\beta$
as in Fig. 3, for light polarized along the $x$ and $y$-directions,
respectively. Interestingly, the behaviour of the absorption is not as simple
as that of the conductivity. For instance, let us look at Fig. 4a. The peak
due to the $1s$ state seems to attenuate from right to left, as we increase
$\beta$. In contrast, the peak due to the $2s$ state seems to intensify. To
study this behaviour more carefully, in Fig. 4c we plot the value of the
absorption peaks due to the $1s$ (black lines) and $2s$ (red lines) states
separately, for polarization along the $x$ (continuous lines) and $y$
directions (dashed lines), as a function of $\beta$, which is allowed to vary
continuously. Indeed, we see that for polarization along the $x$ axis the $1s$
peak decreases monotonically as we approach the zero of $J_{0}(\beta)$, while
the $2s$ peak tends to increase, until it reaches a maximum right before the
zero of the Bessel function. Interestingly, we can have a $2s$ peak more
intense than a $1s$ peak, a feature that can be observed in Figs. 4a (notice
the green and red curves) and 4c. This phenomenon has already been observed
experimentally [46] and explained theoretically [47]. For light polarized
along the $y$-direction, one can observe that the absorption for light
polarized along the $y$-direction is not as intense as for light polarized
along the $x$-direction. Furthermore, from Fig. 4c, we can note that it is now
the $2s$ peak that decreases monotonically, while the $1s$ peak is maximum for
$\beta\approx 0.9$.

###  IV.3 Exciton polaritons

Polaritons are light–matter hybrids, arising from the strong-coupling of
elementary excitations of a 2D material with light. Polaritonic modes are
characterized by an evanescent decay away from the material and a strong
enhancement of the EM field in its vicinity. There are many kinds of
polaritons, like phonon polaritons, plasmon polaritons, magnon polaritons and
exciton polaritons [48]. Here we are concerned only with exciton polaritons.
The dispersion relation of the polariton depends on the polarization of the EM
field. For TE modes, the electric field is in the plane of the 2D material,
and perpendicular to the wavevector of the polariton, while the magnetic field
lies on the $xz$-plane. The dispersion relation for TE modes is

$\kappa_{1}+\kappa_{2}-i\omega\mu_{0}\sigma_{yy}(\omega)=0,$ (14)

where $\varepsilon_{j}$, $j=1,2$, is the relative electric permittivity of the
two media cladding the hBN layer in the frequency region herein considered,
$\kappa_{j}=\sqrt{q^{2}-\varepsilon_{j}\omega^{2}/c^{2}}$, $\omega$ is the
frequency of the polariton, $q$ is the in-plane wavenumber, and $\mu_{0}$ is
the vacuum magnetic permeability.

For TM modes, it is the magnetic field that lies on the plane of the 2D
material, while the electric field is in the $xz$-plane. Their dispersion
relation is

$\frac{\varepsilon_{1}}{\kappa_{1}}+\frac{\varepsilon_{2}}{\kappa_{2}}+i\frac{\sigma_{xx}(\omega)}{\varepsilon_{0}\omega}=0.$
(15)

As a rule of thumb, TE modes are supported when the imaginary part of the
conductivity is negative, whereas TM modes are supported when the imaginary
part of the conductivity is positive. The solutions of the previous two
equations, that can only be obtained numerically, gives the dispersion
relation $\omega(q)$, where the wavenumber has now an imaginary part, since
the conductivity is complex. It is important to highlight that in Eqs. (14)
and (15) two different conductivities appear. In isotropic systems, such as
the case of monolayer hBN, the conductivities $\sigma_{xx}$ and $\sigma_{yy}$
coincide. However, they are in general different, as we have seen in the
previous section.

Figure 5: The dispersion relation $\omega(q)$ for exciton polaritons. We solve
numerically Eq. (15) for different values of $\beta$. The solution is the
frequency as a function of the wavenumber, that has a real part (solid blue
line) and an imaginary part (dashed blue line). For comparison, we show the
dispersion relation for free propagating light in vacuum (solid black line)
and in SiO2 (dashed black line). The case $\beta=0$ represents monolayer hBN
without potential on top.

Regarding the TE modes, we observe that they are either very poorly confined
and present a similar dispersion to free light as expected [45, 49], or they
do not exist at all [that is, Eq. (14) does not present any solution]. Hence,
we omit the analysis for TE modes, as they are not relevant. We focus then on
TM modes solely. To obtain the dispersion relation $\omega(q)$ for TM
polarized polaritons, we solve numerically Eq. (15) in the regions where
$\imaginary\\{\sigma(\omega)\\}>0$, for different values of $\beta$. In
general, the imaginary part of the conductivity admits more than one zero, and
in our case, we have observed in Sec. IV.2 that it has either two or one. The
main region of frequencies of interest is that between the first and second
zeros, in the case $\imaginary\\{\sigma(\omega)\\}$ does have two zeros. After
the third zero, we are very near the bottom of the conduction band, where the
sub-gap conductivity we computed herein is not accurate. To compute the
conductivity above the gap, one needs to also consider unbound states [50].

The dispersion relation of TM modes is presented in Figs. 5a-c, for three
different values of $\beta$. The case $\beta=0$ (Fig. 5a) represents hBN with
no potential on top, and we also choose an intermediate case $\beta=0.5$ (Fig.
5b), and an extreme case $\beta=2.3$ (Fig. 5c). The case without potential and
the intermediate case share similar features, namely the dispersive behaviour
very distinct from free light, with the real part of the wavenumber
$\real\\{q\\}$ reaching a high maximum value at a specific frequency, and the
high degree of confinement given the high values that $\imaginary\\{q\\}$ can
reach. However, the $\beta=2.3$ case of Fig. 5c is qualitatively different.
This is due to the facts that in this case the imaginary part of the
conductivity has only one zero and the conductivity peaks at higher values. As
a consequence, the dispersion behaves like free light in SiO2 for low
frequencies and then behaves like the previous two cases, with $\real\\{q\\}$
achieving a maximum around $36$ $\mu$m-1. Then, the dispersion bends until
$\real\\{q\\}$ reaches a minimum of 14 $\mu$m-1, and afterwards it increases
back again [51]. The degree of confinement is also lower in this case, as can
be seen by the different wavenumber scales (nm-1 in Figs. 5a and 5b, $\mu$m-1
in Fig. 5c). To gain intuition on the confinement of the polariton, we show in
Fig. 7 a spatial map of the EM field of a TM mode in the vicinity of the
interface between air and the substrate.

##  V Conclusion and Outlook

With the advent of the experimental realization of 2D materials, the way of
manipulating the electronic properties of 2D materials by applying an
additional periodic modulation has gained significant attention.
Notwithstanding the current popularity of this topic, the idea of forming
superlattices to manipulate the properties of materials dates back to the end
of the last century, with studies on plasmonic excitations in externally
modulated 2D electron-gases [52, 53]. The periodic modulation can result from
combining two monolayers to form a bilayer structure and applying a twist to
form a moiré superlattice [54, 55], or even by combining bulk systems, like
for example two hBN slabs, bringing their surfaces together and applying a
twist [56].

In our case, the external potential can be imposed by a periodic patterned
back/top gate. The tuning of the potential has a high degree of flexibility,
since our predictions depend only on the product of the strength of the
potential with its period, but not on both independently. In a typical
semiconductor, the exciton size is of the order of 5 nm. For the period, we
need to choose a high enough value compared to the exciton radius, so of the
order of 40 nm. For this length scale and to have $\beta\sim 1$, we need a
potential amplitude of the order of 40 meV, which is feasible in experiments.
However, our theory does not make predictions when $\beta$ matches exactly the
zero of the Bessel function, since the renormalized effective mass along $x$
vanishes and that along $y$ diverges, rendering the Wannier equation
description inapplicable. Also, the simplicity of the potential allowed for an
analytical expression for the renormalized low-energy Hamiltonian, but it
comes with the disadvantage that a cosine potential is more difficult to
implement experimentally. We remark that our theory also applies to a square
periodic potential, more common in experiments. Hence, our results widely
apply to a variety of periodic potential.

Our work presents a way of changing the frequency of the excitonic resonances
of hBN, as an alternative to the usual one of changing the substrate. The
redshift in the resonances as we increase the strength (or period) of the
potential makes it possible to make hBN optically non-trivial in between the
near- and the deep-UV. Furthermore, the system response is different depending
on the polarization of the incident field, which is the essential
characteristic of a grid wire polarizer.

The study of excitonic and optical properties in twisted bilayer systems is a
natural extension of this work. A key feature in these systems is the
emergence of flat electronic bands, enhancing interactions and promoting many-
body correlated states, which offer a plethora of exotic phenomena in physics
[13, 15]. There are already studies of excitons in these systems [57, 58] that
predict that the optical spectrum of the bilayer system is very different from
that of the isolated monolayer. Furthermore, plasmon polaritons have been
studied in twisted bilayer graphene as well [59], where the appearance of
electronic flat bands leads to plasmon polaritons with almost constant
dispersion. However, polaritonics in the UV range within the realm of twisted
bilayer systems remains largely unexplored. With the experimental and
technological tools available nowadays, one can expect solid developments in
the field of 2D moiré structures in general, and in particular in the study of
excitons [60, 61].

## Acknowledgements

N. A. M. is a VILLUM Investigator supported by VILLUM FONDEN (Grant No.
16498). 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 projects PTDC/FIS-
MAC/2045/2021 and EXPL/FIS-MAC/0953/ 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 Nearly-free electron theory

As an independent check of the previous results, we now set up a quasi-free
electron calculation, where the electrons described by $H_{0}$ feel the effect
of the external potential $V(x)$ as a perturbation. The description of the
effect of an external potential on free electrons is a well known problem in
condensed-matter physics [62]. Here we apply the same reasoning to electrons
in hBN. Both the wave functions and the energy bands can be obtained from [63]

$\left[E-E_{s}\left(\mathbf{q}\right)\right]c(s,\mathbf{q})=\sum_{s^{\prime},\mathbf{G}}U_{\mathbf{G}}\langle
s,\mathbf{q}|s^{\prime},\mathbf{q}-\mathbf{G}\rangle
c(s^{\prime},\mathbf{q}-\mathbf{G}),$ (16)

where $E_{s}\left(\mathbf{q}\right)$ is the unperturbed dispersion relation
(2), while the braket uses the spinors of the unperturbed low-energy
Hamiltonian, and $U_{\mathbf{G}}$ are the Fourier components of the potential.
In our case, $U_{\mathbf{G}}=V_{0}/2$ for $\mathbf{G}=\pm G_{0}\hat{x}$, and
$U_{\mathbf{G}}=0$ otherwise. The constants $c(s,\mathbf{q})$ are the
amplitudes of superposition of unperturbed electronic states allowing to
retrieve the electronic wave function. In principle, Eq. (16) is exact.
However, often the summation on the right-hand side in truncated for including
a finite number of reciprocal lattice vectors. Here, we include only three
reciprocal lattice vectors $G=0,\pm 2\pi/L$. This has two consequences: (i)
only the two lowest energy bands (the top valence band and the bottom
conduction bands) are well described and (ii) the value of $\beta$ should be
smaller than 1. In Fig. 6 we compare the results from the exact approach
described in the previous section and quasi-free-electron description
introduced here. The agreement is good.

Figure 6: Renormalized original low-energy dispersion, obtained from the exact
Hamiltonian Eq. (7), and from solving numerically Eq. (16), in blue and dashed
yellow lines, respectively. The particular example with $M=50$ meV is for
$\beta=0.8$. The right side of the plot represents the dispersion as a
function of $q_{x}$ for $q_{y}=0$, while the left one is for dependence on
$q_{y}$ with $q_{x}=0$. The horizontal line is a to guide the eye, so that it
is easier to notice that the dispersion is anisotropic. For this plot we chose
the actual value for $\hbar v_{F}$ for hBN, but a smaller value for the gap,
for illustration reasons.

These results show that, at least for $\beta<1$, we can very well approach
this problem from a nearly-free-electron perspective, and the unitary
transformation is not needed. One of the advantages of using this approach is
that we work directly with the eigenfunctions of $H_{0}$. However, we observed
that this method breaks down if one chooses a combination of $V_{0}$ and $L$
such that $\beta$ assumes values comparable or higher than 1. We stress that
the purpose of this section was to conciliate the exact result from Sec. II
and the nearly-free electron theory, and in what conditions the predictions of
the two different methods agree. From here onward, we work solely with the
exact result.

##  Appendix B Electron-electron exchange correction

The Bethe–Salpeter equation is an eigenproblem for the wavefunction
$\phi_{\nu}(\mathbf{k})$ of the exciton in momentum space and its energy
$E_{\nu}$, and it writes [44]

$\displaystyle
E_{\nu}\phi_{\nu}(\mathbf{k})=\phi_{\nu}(\mathbf{k})(E_{\mathbf{k},+}-E_{\mathbf{k},-})$
$\displaystyle-\frac{1}{A}\sum_{\mathbf{p}}\phi_{\nu}(\mathbf{p})V_{\mathbf{k}-\mathbf{p}}F_{+--+}(\mathbf{k},\mathbf{p},\mathbf{k}-\mathbf{p})$
(17)
$\displaystyle+\frac{1}{A}\phi_{\nu}(\mathbf{k})\sum_{\mathbf{p}}V_{\mathbf{p}}\left[F_{----}(\mathbf{p}+\mathbf{k},\mathbf{k},\mathbf{p})-F_{-+-+}(\mathbf{p}+\mathbf{k},\mathbf{k},\mathbf{p})\right],$

where the form factor

$\displaystyle
F_{\lambda_{1}\lambda_{2}\lambda_{3}\lambda_{4}}\left(\mathbf{p},\mathbf{k},\mathbf{q}\right)=\langle\lambda_{2},\mathbf{p}-\mathbf{q}|\lambda_{3},\mathbf{p}\rangle\langle\lambda_{1},\mathbf{k}+\mathbf{q}|\lambda_{4},\mathbf{k}\rangle$
(18)

is defined in terms of products of the spinors (9). Now we elucidate the
physical meaning of the different terms in Eq. (17). The first term on the
right-hand side accounts for the gain in energy by exciting an electron in the
valence band to the conduction band, without accounting for electron–hole
interactions. The second term accounts for the interaction between the
electron in the conduction band and the hole in the valence band. The third
term corresponds to the electron–electron exchange corrections to the bands,
and it is the one in which we are interested now. To compute the exchange
correction to the band gap, we compute the sum over momentum in the third term
of the Bethe–Salpeter equation (17) for $\mathbf{k}=0$. Denoting the exchange
correction to the gap by $\Delta_{\text{ex}}$, we have

$\Delta_{\text{ex}}=\frac{1}{A}\sum_{\mathbf{p}}V_{\mathbf{p}}\left[F_{----}(\mathbf{p},0,\mathbf{p})-F_{-+-+}(\mathbf{p},0,\mathbf{p})\right].$
(19)

The form factors evaluate to
$F_{----}(\mathbf{p},0,\mathbf{p})=|\langle-,0|-,\mathbf{p}\rangle|^{2}=\cos^{2}(\xi_{p}/2)$
and
$F_{-+-+}(\mathbf{p},0,\mathbf{p})=|\langle+,0|-,\mathbf{p}\rangle|^{2}=\sin^{2}(\xi_{p}/2)$,
so that

$\displaystyle\Delta_{\text{ex}}$
$\displaystyle=\frac{1}{A}\sum_{\mathbf{p}}V_{\mathbf{p}}\cos(\xi_{p})=\frac{1}{A}\sum_{\mathbf{p}}V_{\mathbf{p}}\frac{1}{\sqrt{1+\left(\frac{2\hbar
v_{F}p(\beta)}{J_{0}(\beta)E_{g}}\right)^{2}}}.$ (20)

To proceed, we need to convert the sum into an integral, and the expression
for $V_{\mathbf{p}}$, which is simply the Fourier transform of the
Rytova–Keldysh potential. It reads

$V_{\mathbf{p}}=\frac{e^{2}}{2\varepsilon_{0}p(r_{0}p+\kappa)},$ (21)

and turning the sum into an integral and introducing the auxiliary function
$\gamma(\theta)=\sqrt{\cos^{2}(\theta)+J_{0}^{2}(\beta)\sin^{2}(\theta)}$ we
get

$\displaystyle\Delta_{\text{ex}}=\frac{e^{2}}{8\pi^{2}\varepsilon_{0}\kappa}\int_{0}^{2\pi}\int_{0}^{\infty}\frac{\mathrm{d}p\mathrm{d}\theta}{\left(1+\frac{r_{0}p}{\kappa}\right)\sqrt{1+\left(\frac{2\hbar
v_{F}p\gamma(\theta)}{k_{0}}\right)^{2}}},$ (22)

where $k_{0}=|J_{0}(\beta)|E_{g}/2\hbar v_{F}$ is the characteristic
wavenumber that depends on $\beta$. The radial integral can be performed
analytically by making the change of variables
$\tan(u)=p\gamma(\theta)/k_{0}$, leading to the expression

$\Delta_{\text{ex}}=\frac{e^{2}}{8\pi^{2}\varepsilon_{0}}\frac{k_{0}^{2}r_{0}}{\kappa^{2}}\int_{0}^{2\pi}\mathrm{d}\theta\,\frac{\ln\left[\tfrac{1+\alpha(\theta)+\sqrt{1+\alpha^{2}(\theta)}}{1+\alpha(\theta)-\sqrt{1+\alpha^{2}(\theta)}}\right]}{\alpha(\theta)\sqrt{1+\alpha^{2}(\theta)}}$
(23)

with $\alpha(\theta)\equiv r_{0}k_{0}\gamma(\theta)/\kappa$. We are not aware
of any analytical solution, but the angular integral can be evaluated
numerically.

##  Appendix C Calculation of the Optical conductivity

We give details on the calculation of the Berry curvature

$\Omega_{\nu}=\sum_{\mathbf{q}}\Psi_{\nu}(\mathbf{q})\langle-,\mathbf{q}|\mathbf{r}\cdot\hat{\mathbf{e}}|+,\mathbf{q}\rangle.$
(24)

The main elements that we need are the wavefunction of the excitonic state in
momentum space, and the optical matrix element, which is the same for all
states.

###  C.1 Fourier transform of excitonic wavefunctions

For the variational method, we used the following trial wavefunctions:

$\Psi_{1s}\left(x,y\right)=C_{1s}e^{-\sqrt{a_{1s}x^{2}+b_{1s}y^{2}}},$ (25a)
$\Psi_{2x}\left(x,y\right)=C_{2x}xe^{-\sqrt{a_{2x}x^{2}+b_{2x}y^{2}}},$ (25b)
$\Psi_{2y}\left(x,y\right)=C_{2y}ye^{-\sqrt{a_{2y}x^{2}+b_{2y}y^{2}}},$ (25c)
$\Psi_{2s}\left(x,y\right)=C_{2s}\text{{\hbox{\left(1-d\sqrt{a_{2s}x^{2}+b_{2s}y^{2}}\right)}}}e^{-\sqrt{a_{2s}x^{2}+b_{2s}y^{2}}}.$
(25d)

For each one of the states $1s$, $2x$, $2y$, and $2s$, we have two variational
parameters denoted by $a$ and $b$, which account for the anisotropy of the
Hamiltonian. In the isotropic case, $a=b$ for each state and then the
variational wavefunctions coincide with those of the 2D hydrogen atom. We
stress that the major difference between 2D Hydrogen and pristine hBN is the
form of the electrostatic potential. The prefactors denoted by $C$ are
normalization constants, and all states must be orthogonal to each other. By
symmetry, the $s$ states and the $p$ states are orthogonal, and the states $p$
are also orthogonal to each other. Orthogonality between the $1s$ and $2s$
states can be achieved by imposing
$\int\mathrm{d}\mathbf{r}\,\Psi_{2s}^{*}\Psi_{1s}=0$, obtaining in this way an
expression for $d$, that depends on $a_{2s}$ and $b_{2s}$, which in turn
depend on the previously calculated $a_{1s}$ and $b_{1s}$. The details of how
the variational parameters for each state are calculated can be consulted in
[39].

The momentum-resolved wavefunction of the state $\nu$ is obtained by Fourier
transforming the respective wavefunction defined in real space. For the ground
state $1s$ the 2D Fourier transform
$\Psi\left(\mathbf{k}\right)=\int\mathrm{d}\mathbf{r}\,e^{-i\mathbf{k\cdot\mathbf{r}}}\,\Psi\left(\mathbf{r}\right)$
gives

$\displaystyle\Psi_{1s}\left(\mathbf{k}\right)=\frac{2\sqrt{2\pi}}{\left(a_{1s}b_{1s}\right)^{1/4}}\left(1+\frac{k_{x}^{2}}{a_{1s}}+\frac{k_{y}^{2}}{b_{1s}}\right)^{-3/2},$
(26a) where we performed the change of variables
$x=\rho/\sqrt{a_{1s}}\cos(\theta)$, $y=\rho/\sqrt{b_{1s}}\sin(\theta)$, and
computed the radial and angular integrals analytically. We note that even
though the expression is analytical, it has a non-trivial angular dependence.
The excited states considered in this work have also analytical expressions
for the Fourier transforms $\Psi_{2x}(\mathbf{k})=-C_{2x}\frac{3i\pi
k_{x}/a_{2x}}{\left(1+\frac{k_{x}^{2}}{a_{2x}}+\frac{k_{y}^{2}}{b_{2x}}\right)^{5/2}},$
(26b) $\Psi_{2y}(\mathbf{k})=-C_{2y}\frac{3i\pi
k_{y}/b_{2y}}{\left(1+\frac{k_{x}^{2}}{a_{2y}}+\frac{k_{y}^{2}}{b_{2y}}\right)^{5/2}},$
(26c)
$\Psi_{2s}(\mathbf{k})=C_{2s}2\pi\frac{1-2d+(1+d)\left(\frac{k_{x}^{2}}{a_{2s}}+\frac{k_{y}^{2}}{b_{2s}}\right)}{\left(1+\frac{k_{x}^{2}}{a_{2s}}+\frac{k_{y}^{2}}{b_{2s}}\right)^{5/2}}\,,$
(26d)

where calculation details are similar to the case of the state $1s$.

###  C.2 Optical matrix element

The optical matrix element can be obtained by using the relation

$\langle-,\mathbf{q}|\mathbf{r}|+,\mathbf{q}\rangle=\frac{\langle-,\mathbf{q}|[H,\mathbf{r}]|+,\mathbf{q}\rangle}{E_{-,\mathbf{q}}-E_{+,\mathbf{q}}},$
(27)

and the canonical commutation relations to evaluate the commutator, to yield

$\displaystyle\langle-,\mathbf{q}|\mathbf{r}|+,\mathbf{q}\rangle=$
$\displaystyle\frac{i\hbar
v_{F}}{E_{-,\mathbf{q}}-E_{+,\mathbf{q}}}\Big{[}\cos^{2}\left(\tfrac{\xi_{q}}{2}\right)\left[\hat{x}-iJ_{\ell}\left(\beta\right)\hat{y}\right]$
(28)
$\displaystyle-e^{-2i\theta_{q}\left(\beta\right)}\sin^{2}\left(\tfrac{\xi_{q}}{2}\right)\left[\hat{x}+iJ_{\ell}\left(\beta\right)\hat{y}\right]\Big{]}\,.$

In the isotropic case, the previous expression dictates that only the $s$ and
$d$ states contribute to the optical conductivity. The term with cosine
contributes only for $s$ states, while the term with sine contributes only for
$d$ states. However, in our case, the imaginary exponential has a non-trivial
dependence on the angle $\theta$. Explicitly,
$\exp\left[-2i\theta_{k}(\beta)\right]=\exp\left\\{-2i\arctan\left[J_{\ell}(\beta)\arctan(\theta)\right]\right\\}$.
In this way, the term with the sine gives a contribution to the Berry
curvature for $s$ states in the anisotropic case. So we write the Berry
curvature $\Omega_{\nu}=\Omega^{\cos}_{\nu}+\Omega^{\sin}_{\nu}$ as a sum of
two contributions. We also note that in the optical matrix element, the terms
along $y$-direction have an extra factor of $J_{\ell}\left(\beta\right)$,
resulting in different conductivities depending on the polarization of the
incident field.

To compute the Berry curvature (24), we convert the sum over momentum into an
integral, for which the following auxiliary function will be useful

$\sigma\left(\theta\right)=\sqrt{\frac{\cos^{2}(\theta)}{a_{1s}}+\frac{\sin^{2}(\theta)}{b_{1s}}}\,.$
(29)

With these functions, we can write the quantities

$\displaystyle E_{-,\mathbf{q}}-E_{+,\mathbf{q}}$
$\displaystyle=\sqrt{E_{g}^{2}J_{\ell}^{2}(\beta)+4\hbar^{2}v_{F}^{2}q^{2}\gamma^{2}\left(\theta\right)},$
(30) $\displaystyle\cos^{2}\left(\tfrac{\xi_{q}}{2}\right)$
$\displaystyle=\frac{1}{2}+\frac{1}{2\sqrt{1+\left(\frac{2\hbar
v_{F}q\gamma\left(\theta\right)}{J_{\ell}\left(\beta\right)E_{g}}\right)^{2}}},$
$\displaystyle\sin^{2}\left(\tfrac{\xi_{q}}{2}\right)$
$\displaystyle=\frac{1}{2}-\frac{1}{2\sqrt{1+\left(\frac{2\hbar
v_{F}q\gamma\left(\theta\right)}{J_{\ell}\left(\beta\right)E_{g}}\right)^{2}}},$

that appear in the optical matrix element, and for the ground state
wavefunction

$\Psi_{1s}\left(\mathbf{k}\right)=\frac{2\sqrt{2\pi}}{\left(a_{1s}b_{1s}\right)^{1/4}}\frac{1}{\left[1+k^{2}\sigma^{2}(\theta)\right]^{3/2}}.$
(31)

###  C.3 Optical conductivity

Focusing on the cosine term and light polarized along the $x$-direction, we
can write

$\displaystyle\Omega^{\cos}_{1s}=\frac{-2\sqrt{2\pi}i\hbar
v_{F}}{(a_{1s}b_{1s})^{-1/4}}\int_{0}^{2\pi}\mathrm{d}\theta\left[\frac{\hbar
v_{F}\gamma(\theta)-E_{g}|J_{\ell}(\beta)|\sigma(\theta)}{E_{g}^{2}J_{\ell}^{2}(\beta)\sigma^{3}(\theta)-4\hbar^{2}v_{F}^{2}\sigma(\theta)\gamma^{2}(\theta)}\frac{E_{g}\hbar
v_{F}|J_{\ell}(\beta)|\gamma(\theta)\arcsec\left(\frac{E_{g}|J_{\ell}(\beta)|\sigma(\theta)}{2\hbar
v_{F}\gamma(\theta)}\right)}{\left(E_{g}^{2}J_{\ell}^{2}(\beta)\sigma^{2}(\theta)-4\hbar^{2}v_{F}^{2}\gamma^{2}(\theta)\right)^{3/2}}\right]\,,$
(32)

where the radial integral was performed analytically, whereas the angular
integral can only be worked out numerically. The term that goes with the sine
leads to

$\displaystyle\Omega^{\sin}_{1s}$ $\displaystyle=\frac{-2\sqrt{2\pi}i\hbar
v_{F}}{(a_{1s}b_{1s})^{-1/4}}\int_{0}^{2\pi}\mathrm{d}\theta\left[\frac{-\hbar
v_{F}\gamma(\theta)}{E_{g}^{2}J_{\ell}^{2}(\beta)\sigma^{3}(\theta)-4\hbar^{2}v_{F}^{2}\sigma(\theta)\gamma^{2}(\theta)}+\frac{E_{g}\hbar
v_{F}|J_{\ell}(\beta)|\gamma(\theta)\arcsec\left(\frac{E_{g}|J_{\ell}(\beta)|\sigma(\theta)}{2\hbar
v_{F}\gamma(\theta)}\right)}{\left(E_{g}^{2}J_{\ell}^{2}(\beta)\sigma^{2}(\theta)-4\hbar^{2}v_{F}^{2}\gamma^{2}(\theta)\right)^{3/2}}\right]e^{-2i\theta_{k}\left(\beta\right)}\,,$
(33)

that needs to be evaluated numerically as well. The Berry curvature is then
simply $\Omega^{x}_{1s}=\Omega^{\cos}_{1s}+\Omega^{\sin}_{1s}$, where we put
an $x$ in the superscript to stress that the calculation was for light
polarized along the $x$-direction. For light polarized along the
$y$-direction, we can identify from the optical matrix element (28) an extra
factor of $J_{m}(\beta)$, and also a relative sign between the cosine and sine
terms, so that the Berry curvature is now
$\Omega^{y}_{1s}=J_{m}(\beta)(\Omega^{\cos}_{1s}-\Omega^{\sin}_{1s})$. For
symmetry reasons, neither the state $2x$ nor $2y$ contribute to the optical
conductivity. However, the excited state $2s$ contributes to the conductivity.
Its contribution is calculated in the same way as for the $1s$ state, however,
even though the radial integral can be computed analytically, its expression
is quite involved, and is not relevant for the discussion. At the end, the
angular integral has to be computed numerically anyway. Finally, the optical
conductivity is computed through

$\sigma_{jj}(\omega)=\frac{e^{2}}{8\pi^{3}i\hbar}\sum_{\nu=1s,2s}\frac{E_{\nu}|\Omega^{j}_{\nu}|^{2}}{E_{\nu}-(\hbar\omega+i\Gamma)}+\frac{E_{\nu}|\Omega^{j}_{\nu}|^{2}}{E_{\nu}+\hbar\omega+i\Gamma}$
(34)

where $j=x,y$ for light polarized along the $x$ and $y$-directions,
respectively.

##  Appendix D Spatial map of the EM field

In the figure below we show a spatial map of the EM field of a TM polariton in
the vicinity of the hBN monolayer. As one can see, the intensity of both
magnetic and electric fields decay away from the 2D material located at the
plane $z=0$. We can also infer from Fig. 7 that the field in the substrate is
more confined than in air. This degree of spatial confinement is in the order
of the nanometers ($\sim 0.5$ nm).

Figure 7: Spatial map of the EM field of a TM mode. We make an intensity plot
representing the intensity of the single component magnetic field, and a
stream plot representing the off-plane components of electric field and
respective intensities. The top half represents the field in air, while the
bottom half represents the field in SiO2. For the frequency wavenumber pair,
we picked the point in the dispersion of Fig. 5(b) of highest (real)
wavenumber. This plot does not take into account the imaginary part of the
wavenumber. Parameters used: $q=4.45$ nm-1, $\hbar\omega=5.68$ eV,
$\varepsilon_{1}=1$ and $\varepsilon_{2}=2.16$.

## References

  * Jiang  _et al._ [2020] H. Jiang, L. Zheng, Z. Liu, and X. Wang, Two-dimensional materials: From mechanical properties to flexible mechanical sensors, InfoMat 2, 1077–1094 (2020).
  * Dutta  _et al._ [2023] T. Dutta, N. Yadav, Y. Wu, G. J. Cheng, X. Liang, S. Ramakrishna, A. Sbai, R. Gupta, A. Mondal, Z. Hongyu, and A. Yadav, Electronic properties of 2D materials and their junctions, Nano Mater. Sci. (2023), DOI: 10.1016/j.nanoms.2023.05.003.
  * Jia [2019] B. Jia, 2D optical materials and the implications for photonics, APL Photonics 4, 080401 (2019).
  * Gan  _et al._ [2012] X. Gan, K. F. Mak, Y. Gao, Y. You, F. Hatami, J. Hone, T. F. Heinz, and D. Englund, Strong enhancement of light–matter interaction in graphene coupled to a photonic crystal nanocavity, Nano Lett. 12, 5626–5631 (2012).
  * Geim and Novoselov [2007] A. K. Geim and K. S. Novoselov, The rise of graphene, Nat. Mater. 6, 183–191 (2007).
  * Miller [2023] J. L. Miller, A family of two-dimensional conductors comes into bloom, Phys. Today 76(6), 12 (2023).
  * Song  _et al._ [2013] X. Song, J. Hu, and H. Zeng, Two-dimensional semiconductors: recent progress and future perspectives, J. Mater. Chem. C 1, 2952–2969 (2013).
  * Huang  _et al._ [2022] X. Huang, C. Liu, and P. Zhou, 2D semiconductors for specific electronic applications: from device to system, npj 2D Mater. Appl. 6, 51 (2022).
  * Cao  _et al._ [2020] Y. Cao, D. Chowdhury, D. Rodan-Legrain, O. Rubies-Bigorda, K. Watanabe, T. Taniguchi, T. Senthil, and P. Jarillo-Herrero, Strange metal in magic-angle graphene with near Planckian dissipation, Phys. Rev. Lett. 124, 076801 (2020).
  * Jiang  _et al._ [2018] Z. Jiang, P. Wang, X. Jiang, and J. Zhao, MBene (MnB): a new type of 2D metallic ferromagnet with high Curie temperature, Nanoscale Horiz. 3, 335–341 (2018).
  * Yu  _et al._ [2022] S. Yu, X. Li, J. Tang, Y. Wang, F. Xu, and X. Wang, Recent advances in two-dimensional ferromagnetism: strain-, doping-, structural- and electric field-engineering toward spintronic applications, Sci. Technol. Adv. Mater. 23, 140–160 (2022).
  * Qiu  _et al._ [2021] D. Qiu, C. Gong, S. Wang, M. Zhang, C. Yang, X. Wang, and J. Xiong, Recent advances in 2D superconductors, Adv. Mater. 33, 2006124 (2021).
  * Cao  _et al._ [2018a] Y. Cao, V. Fatemi, S. Fang, K. Watanabe, T. Taniguchi, E. Kaxiras, and P. Jarillo-Herrero, Unconventional superconductivity in magic-angle graphene superlattices, Nature 556, 43–50 (2018a).
  * Park  _et al._ [2021] J. M. Park, Y. Cao, K. Watanabe, T. Taniguchi, and P. Jarillo-Herrero, Tunable strongly coupled superconductivity in magic-angle twisted trilayer graphene, Nature 590, 249–255 (2021).
  * Cao  _et al._ [2018b] Y. Cao, V. Fatemi, A. Demir, S. Fang, S. L. Tomarken, J. Y. Luo, J. D. Sanchez-Yamagishi, K. Watanabe, T. Taniguchi, E. Kaxiras, R. C. Ashoori, and P. Jarillo-Herrero, Correlated insulator behaviour at half-filling in magic-angle graphene superlattices, Nature 556, 80–84 (2018b).
  * Abidi  _et al._ [2019] I. H. Abidi, N. Mendelson, T. T. Tran, A. Tyagi, M. Zhuang, L.-T. Weng, B. Özyilmaz, I. Aharonovich, M. Toth, and Z. Luo, Selective defect formation in hexagonal boron nitride, Adv. Opt. Mater. 7, 1900397 (2019).
  * Banszerus  _et al._ [2016] L. Banszerus, M. Schmitz, S. Engels, M. Goldsche, K. Watanabe, T. Taniguchi, B. Beschoten, and C. Stampfer, Ballistic transport exceeding 28 $\mu\text{m}$ in CVD grown graphene, Nano Lett. 16, 1387–1391 (2016).
  * Rizzo  _et al._ [2023] D. J. Rizzo, J. Zhang, B. S. Jessen, F. L. Ruta, M. Cothrine, J. Yan, D. G. Mandrus, S. E. Nagler, T. Taniguchi, K. Watanabe, M. M. Fogler, A. N. Pasupathy, A. J. Millis, A. Rubio, J. C. Hone, C. R. Dean, and D. N. Basov, Polaritonic probe of an emergent 2D dipole interface, Nano Lett. 23, 8426–8435 (2023).
  * Elias  _et al._ [2019] C. Elias, P. Valvin, T. Pelini, A. Summerfield, C. J. Mellor, T. S. Cheng, L. Eaves, C. T. Foxon, P. H. Beton, S. V. Novikov, B. Gil, and G. Cassabois, Direct band-gap crossover in epitaxial monolayer boron nitride, Nat. Commun. 10, 2639 (2019).
  * Yang  _et al._ [2023] W. Yang, L. K. Ang, W. Zhang, J. Han, and Y. Xu, High sensitivity gas sensor based on surface exciton polariton enhanced photonic spin Hall effect, Opt. Express 31, 27041–27053 (2023).
  * Li  _et al._ [2022] S. Li, A. Pershin, G. Thiering, P. Udvarhelyi, and A. Gali, Ultraviolet quantum emitters in hexagonal boron nitride from carbon clusters, J. Phys. Chem. Lett. 13, 3150–3157 (2022).
  * Geim and Grigorieva [2013] A. K. Geim and I. V. Grigorieva, Van der Waals heterostructures, Nature 499, 419 (2013).
  * Henriques  _et al._ [2022a] J. C. G. Henriques, B. Amorim, R. M. Ribeiro, and N. M. R. Peres, Excitonic response of AA’ and AB stacked hBN bilayers, Phys. Rev. B 105, 115421 (2022a).
  * Castro  _et al._ [2007] E. V. Castro, K. S. Novoselov, S. V. Morozov, N. M. R. Peres, J. M. B. L. dos Santos, J. Nilsson, F. Guinea, A. K. Geim, and A. H. C. Neto, Biased bilayer graphene: Semiconductor with a gap tunable by the electric field effect, Phys. Rev. Lett. 99, 216802 (2007).
  * Bistritzer and MacDonald [2011] R. Bistritzer and A. H. MacDonald, Moiré bands in twisted double-layer graphene, Proc. Nat. Acad. Sci. U. S. A. 108, 12233–12237 (2011).
  * Lopes dos Santos  _et al._ [2007] J. M. B. Lopes dos Santos, N. M. R. Peres, and A. H. Castro Neto, Graphene bilayer with a twist: Electronic structure, Phys. Rev. Lett. 99, 256802 (2007).
  * Hennighausen and Kar [2021] Z. Hennighausen and S. Kar, Twistronics: a turning point in 2D quantum materials, Electron. Struct. 3, 014004 (2021).
  * Chaves  _et al._ [2020] A. Chaves, J. G. Azadani, H. Alsalman, D. R. da Costa, R. Frisenda, A. J. Chaves, S. H. Song, Y. D. Kim, D. He, J. Zhou, A. Castellanos-Gomez, F. M. Peeters, Z. Liu, C. L. Hinkle, S.-H. Oh, P. D. Ye, S. J. Koester, Y. H. Lee, P. Avouris, X. Wang, and T. Low, Bandgap engineering of two-dimensional semiconductor materials, npj 2D Mater. Appl. 4, 29 (2020).
  * [29] L. Xue, S. Liu, Y. Hang, A. M. Summers, D. J. Wilson, X. Wang, P. Chen, T. G. Folland, J. A. Hachtel, H. Shi, S. Hosseini-Zavareh, S. R. Das, S. Lei, Z. Zhang, C. M. Sorensen, W. Guo, J. D. Caldwell, J. H. Edgar, C. I. Blaga, and C. A. Trallero-Herrero, Unraveling ultrafast photoionization in hexagonal boron nitride, arXiv.2101.10429 .
  * Thygesen [2017] K. S. Thygesen, Calculating excitons, plasmons, and quasiparticles in 2D materials and van der Waals heterostructures, 2D Materials 4, 022004 (2017).
  * Bhattacharya  _et al._ [2012] A. Bhattacharya, S. Bhattacharya, and G. P. Das, Band gap engineering by functionalization of BN sheet, Phys. Rev. B 85, 035415 (2012).
  * Román  _et al._ [2021] R. J. P. Román, F. J. R. C. Costa, A. Zobelli, C. Elias, P. Valvin, G. Cassabois, B. Gil, A. Summerfield, T. S. Cheng, C. J. Mellor, P. H. Beton, S. V. Novikov, and L. F. Zagonel, Band gap measurements of monolayer h-BN and insights into carbon-related point defects, 2D Materials 8, 044001 (2021).
  * Ferreira  _et al._ [2019] F. Ferreira, A. J. Chaves, N. M. R. Peres, and R. M. Ribeiro, Excitons in hexagonal boron nitride single-layer: a new platform for polaritonics in the ultraviolet, J. Opt. Soc. Am. B 36, 674–683 (2019).
  * Fears  _et al._ [2013] K. P. Fears, D. Y. Petrovykh, S. J. Photiadis, and T. D. Clark, Circular dichroism analysis of cyclic $\beta$-helical peptides adsorbed on planar fused quartz, Langmuir 29, 10095–10101 (2013).
  * Castro Neto  _et al._ [2009] A. H. Castro Neto, F. Guinea, N. M. R. Peres, K. S. Novoselov, and A. K. Geim, The electronic properties of graphene, Rev. Mod. Phys. 81, 109 (2009).
  * Park  _et al._ [2008a] C.-H. Park, L. Yang, Y.-W. Son, M. L. Cohen, and S. G. Louie, New generation of massless Dirac fermions in graphene under external periodic potentials, Phys. Rev. Lett. 101, 126804 (2008a).
  * Quintela  _et al._ [2022] M. F. C. M. Quintela, J. C. G. Henriques, L. G. M. Tenório, and N. M. R. Peres, Theoretical methods for excitonic physics in 2D materials, phys. stat. solid. (b) 259, 2200097 (2022).
  * Have  _et al._ [2019] J. Have, G. Catarina, T. G. Pedersen, and N. M. R. Peres, Monolayer transition metal dichalcogenides in strong magnetic fields: Validating the Wannier model using a microscopic calculation, Phys. Rev. B 99, 035416 (2019).
  * Gomes  _et al._ [2021] J. N. S. Gomes, C. Trallero-Giner, and M. I. Vasilevskiy, Variational calculation of the lowest exciton states in phosphorene and transition metal dichalcogenides, J. Phys.: Cond. Matter 34, 045702 (2021).
  * Martins Quintela and Peres [2020] M. F. C. Martins Quintela and N. M. R. Peres, A colloquium on the variational method applied to excitons in 2D materials, Eur. Phys. J. B 93, 222 (2020).
  * Henriques  _et al._ [2021] J. C. G. Henriques, N. A. Mortensen, and N. M. R. Peres, Analytical description of the $1s$ exciton linewidth temperature dependence in transition metal dichalcogenides, Phys. Rev. B 103, 235402 (2021).
  * Sakurai and Napolitano [2020] J. J. Sakurai and J. Napolitano,  _Modern Quantum Mechanics_, 3rd ed. (Cambridge University Press, 2020).
  * Pedersen [2015] T. G. Pedersen, Intraband effects in excitonic second-harmonic generation, Phys. Rev. B 92, 235432 (2015).
  * Henriques  _et al._ [2019] J. C. G. Henriques, G. B. Ventura, C. D. M. Fernandes, and N. M. R. Peres, Optical absorption of single-layer hexagonal boron nitride in the ultraviolet, J. Phys.: Cond. Matter 32, 025304 (2019).
  * BLudov  _et al._ [2013] Y. V. BLudov, A. Ferreira, N. M. R. Peres, and M. I. Vasilevskiy, A primer on surface plasmon-polaritons in graphene, Int. J. Mod. Phys. B 27, 1341001 (2013).
  * Ju  _et al._ [2017] L. Ju, L. Wang, T. Cao, T. Taniguchi, K. Watanabe, S. G. Louie, F. Rana, J. Park, J. Hone, F. Wang, and P. L. McEuen, Tunable excitons in bilayer graphene, Science 358, 907–910 (2017).
  * Henriques  _et al._ [2022b] J. C. G. Henriques, I. Epstein, and N. M. R. Peres, Absorption and optical selection rules of tunable excitons in biased bilayer graphene, Phys. Rev. B 105, 045411 (2022b).
  * Basov  _et al._ [2021] D. N. Basov, A. Asenjo-Garcia, P. J. Schuck, X. Zhu, and A. Rubio, Polariton panorama, Nanophotonics 10, 549 (2021).
  * Gonçalves and Peres [2016] P. A. D. Gonçalves and N. M. R. Peres,  _An Introduction to Graphene Plasmonics_ (World Scientific, Singapore, 2016).
  * Leppenen  _et al._ [2021] N. V. Leppenen, L. E. Golub, and E. L. Ivchenko, Sommerfeld enhancement factor in two-dimensional Dirac materials, Phys. Rev. B 103, 235311 (2021).
  * Wolff  _et al._ [2018] C. Wolff, K. Busch, and N. A. Mortensen, Modal expansions in periodic photonic systems with material loss and dispersion, Phys. Rev. B 97, 104203 (2018).
  * van Zyl and Zaremba [1999] B. P. van Zyl and E. Zaremba, Thomas–Fermi–Dirac–von Weizsäcker hydrodynamics in laterally modulated electronic systems, Phys. Rev. B 59, 2079 (1999).
  * Fessatidis and Cui [1991] V. Fessatidis and H. L. Cui, Hydrodynamic theory of magnetoplasmons in a modulated two-dimensional electron gas, Phys. Rev. B 43, 11725 (1991).
  * Chu  _et al._ [2020] Y. Chu, L. Liu, Y. Yuan, C. Shen, R. Yang, D. Shi, W. Yang, and G. Zhang, A review of experimental advances in twisted graphene moiré superlattice, Chin. Physics B 29, 128104 (2020).
  * Yan  _et al._ [2019] W. Yan, L. Meng, Z. Meng, Y. Weng, L. Kang, and X.-a. Li, Probing angle-dependent interlayer coupling in twisted bilayer WS2, J. Phys. Chem. C 123, 30684–30688 (2019).
  * Lee  _et al._ [2021] H. Y. Lee, M. M. Al Ezzi, N. Raghuvanshi, J. Y. Chung, K. Watanabe, T. Taniguchi, S. Garaj, S. Adam, and S. Gradečak, Tunable optical properties of thin films controlled by the interface twist angle, Nano Lett. 21, 2832–2839 (2021).
  * Brem and Malic [2023] S. Brem and E. Malic, Bosonic delocalization of dipolar moiré excitons, Nano Lett. 23, 4627–4633 (2023).
  * Roman-Taboada  _et al._ [2023] P. Roman-Taboada, E. Obregon-Castillo, A. R. Botello-Mendez, and C. Noguez, Excitons in twisted ${\mathrm{aa}}^{’}$ hexagonal boron nitride bilayers, Phys. Rev. B 108, 075109 (2023).
  * Stauber and Kohler [2016] T. Stauber and H. Kohler, Quasi-flat plasmonic bands in twisted bilayer graphene, Nano Lett. 16, 6844–6849 (2016).
  * Wu  _et al._ [2022] B. Wu, H. Zheng, S. Li, J. Ding, J. He, Y. Zeng, K. Chen, Z. Liu, S. Chen, A. Pan, and Y. Liu, Evidence for moiré intralayer excitons in twisted WSe2/WSe2 homobilayer superlattices, Light Sci. Appl. 11, 166 (2022).
  * Zhao  _et al._ [2023] S. Zhao, Z. Li, X. Huang, A. Rupp, J. Göser, I. A. Vovk, S. Y. Kruchinin, K. Watanabe, T. Taniguchi, I. Bilgin, A. S. Baimuratov, and A. Högele, Excitons in mesoscopically reconstructed moiré heterostructures, Nat. Nanotechnol. 18, 572–579 (2023).
  * Ashcroft and Mermin [1976] N. W. Ashcroft and N. D. Mermin,  _Solid State Physics_ (Harcourt, Orlando, 1976).
  * Park  _et al._ [2008b] C.-H. Park, L. Yang, Y.-W. Son, M. L. Cohen, and S. G. Louie, Anisotropic behaviours of massless Dirac fermions in graphene under periodic potentials, Nat. Phys. 4, 213–217 (2008b).


Undergraduate Texts in Mathematics. Springer New York, 2008.